Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc.qmd

---
title-block-banner: true
title: "Analysis of Ki67/PDGFRβ colocalization in defined ROIs of the ipsilateral hemisphere"
subtitle: "Data analysis notebook"
date: today
date-format: full
author: 
  - name: "Daniel Manrique-Castano"
    orcid: 0000-0002-1912-1764
    degrees:
      - PhD
    affiliation: 
      - name: Univerisity Laval 
        department: Psychiatry and Neuroscience
        group: Laboratory of neurovascular interactions 
note: "GitHub: https://daniel-manrique.github.io/"
keywords: 
  - Ki67
  - PDGFRβ
  - Brain injury
  - Cell proliferation
  - Bayesian modeling 
   
license: "CC BY"

format:
   pdf: 
    toc: true
    number-sections: true
    colorlinks: true
   html:
    code-fold: true
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: left
    number-sections: true
    theme: spacelab

knitr:
  opts_chunk: 
    warning: false
    message: false
    
csl: science.csl
bibliography: references.bib
---

# Preview

This notebook reports the analysis of Ki67/PDGFRβ colocalization in defined ROIs of the ipsilateral hemisphere following cerebral ischemia.

**Parent dataset:** Ki67, PDGFRβ, and CD31 stained ischemic hemispheres imaged at 20x using confocal microscopy. Samples are grouped at 0 (Sham), 3, and 7 days post-ischemia (DPI). The raw images and pre-processing scripts to generate the analyzed Z-projected images are available at the Zenodo repository (10.5281/zenodo.10553084) under the name `Confocal_20x_ROIs_Ki67-Pdgfrb-CD31.zip`.

**Working dataset**: The `Data_Raw/Confocal_20x_ROIs_Ki67-Pdgfrb-CD31/Image.csv`data frame contains the cell detection and colocalization analysis performed on CellProfiler [@stirling2021]. The CellProfiler pipeline is available at https://osf.io/79wzq/.

# Install and load required packages

Install and load all required packages. Please uncomment (delete #) the line code if installation is required. Load the installed libraries each time you start a new R session.

```{r}
#| label: Install_Packages
#| include: true
#| warning: false
#| message: false

#install.packages("devtools")
#library(devtools)

#install.packages(c("bayesplot", "bayestestR", "brms","dplyr", "easystats", "ggplot","gtsummary", "modelbased", "modelr", "modelsummary", "patchwork", "poorman","plyr", "tidybayes", "tidyverse"))


library(bayesplot)
library(bayestestR)
library(brms)
library(dplyr)
library(easystats)
library(ggplot2)
library(gtsummary)
library(modelbased)
library(modelr)
library(modelsummary)
library(patchwork)
library(poorman)
library(plyr)
library(tidybayes)
library(tidyverse)
```

# Visual themes

We create a visual theme to use in our plots.

```{r}
#| label: Plot_Theme
#| include: true
#| warning: false
#| message: false
  
Plot_theme <- theme_classic() +
  theme(
      plot.title = element_text(size=18, hjust = 0.5, face="bold"),
      plot.subtitle = element_text(size = 10, color = "black"),
      plot.caption = element_text(size = 12, color = "black"),
      axis.line = element_line(colour = "black", size = 1.5, linetype = "solid"),
      axis.ticks.length=unit(7,"pt"),
     
      axis.title.x = element_text(colour = "black", size = 16),
      axis.text.x = element_text(colour = "black", size = 16, angle = 0, hjust = 0.5),
      axis.ticks.x = element_line(colour = "black", size = 1),
      
      axis.title.y = element_text(colour = "black", size = 16),
      axis.text.y = element_text(colour = "black", size = 16),
      axis.ticks.y = element_line(colour = "black", size = 1),
      
      legend.position="right",
      legend.direction="vertical",
      legend.title = element_text(colour="black", face="bold", size=12),
      legend.text = element_text(colour="black", size=10),
      
      plot.margin = margin(t = 10,  # Top margin
                             r = 2,  # Right margin
                             b = 10,  # Bottom margin
                             l = 10) # Left margin
      ) 
```

# Load the data set

We load the dataset and handle it the subset the columns of interest.

```{r}
#| label: Ki67-Pdgfrb_Load 
#| include: true
#| warning: false
#| message: false
#| cache: true

# We load the dataset in case is not present in the R environment
Ki67_Cells <- read.csv(file = "Data_Raw/Confocal_20x_ROIs_Ki67-Pdgfrb-CD31/Image.csv", header = TRUE)

## We subset the relevant columns (cell number)
Ki67_Data <- subset(Ki67_Cells, select = c("FileName_CD31_Raw", "Count_All_Pdgfrb_Ki67_Colocalized","Count_Vascular_Pdgfrb_Ki61_Resized", "Count_Ki67_Filtered", "Count_Pdgfrb_Filled"))

## And extract metadata from the image name
Ki67_Data  <- cbind(Ki67_Data, do.call(rbind , strsplit(Ki67_Data$FileName_CD31_Raw, "[_\\.]"))[,1:4])

Ki67_Data <- subset(Ki67_Data, select = -c(FileName_CD31_Raw))

## We Rename the relevant columns 
colnames(Ki67_Data) <- c("Total_Ki67_Pdgfrb", "Vascular_Ki67_Pdgfrb", "Count_Ki67_Filtered","Count_Pdgfrb_Filled", "MouseID", "DPI", "Condition", "Region")

## We set the factors
Ki67_Data$DPI <- factor(Ki67_Data$DPI, levels = c("3D", "7D"))
Ki67_Data$Region <- factor(Ki67_Data$Region, levels = c("Peri", "Str", "Ctx"))


# Create an additional DPI variable (numeric)

DPI_mapping <- c("3D" = "3", "7D" = "7")
Ki67_Data$DPI_Cont <- as.numeric(DPI_mapping[as.character(Ki67_Data$DPI)])

write.csv(Ki67_Data, "Data_Processed/Confocal_20x_ROIs_Ki67-Pdgfrb-CD31/Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc.csv")

gt::gt(Ki67_Data[1:10,])
```

The dataset available at 10.5281/zenodo.10553084 includes images at 0D (sham animals) that were used to train the Pixel classification in ilastik [@berg2019]. There are not Ki67/PDGFRβ colocalized cells at 0 DPI. Also, includes images at 14D and 30D that were taken just for visualization and reference but will not be included in the analysis. 


# Exploratory data visualization

We perform the exploratory visualization for the dataset. We'll focus on the total colocalization of Ki67/PDGFRβ and the colocalization in vascular PDGFRβ cells.

```{r}
#| label: fig-Ki67_Coloc_Exploratory
#| include: true
#| warning: false
#| message: false
#| fig-cap: Exploratory data visualization 
#| fig-height: 5
#| fig-width: 12

set.seed(8807)

Ki67_Coloc_Sctr <- 
  ggplot(
    data  = Ki67_Data, 
    aes(x = DPI_Cont, 
        y = Total_Ki67_Pdgfrb)) +
geom_smooth(
  method = "lm", 
  se     = TRUE,
  color  = "black") +
geom_smooth(
  method  = "lm", 
  se      = TRUE, 
  formula = y ~ poly(x, 2), 
  color   = "darkred") +
geom_smooth(
  method  = "lm", 
  se      = TRUE, 
  formula = y ~ poly(x, 3), 
  color   = "darkgreen") +
geom_jitter(
  width = 0.5, 
  shape = 1, 
  size  = 1.5, 
  aes(color = Region)) +
  
  scale_y_continuous(name= expression("Total PDGFRβ/Ki67+ cells")) +
  scale_x_continuous(name="DPI",
                     breaks=c(0, 3, 7)) +
  Plot_theme


Ki67_Coloc_Sctr2 <- 
  ggplot(
    data  = Ki67_Data, 
    aes(x = DPI_Cont, 
        y = Vascular_Ki67_Pdgfrb,
        color = Region)) +
geom_smooth(
  method = "lm", 
  se     = TRUE,
  color  = "black") +
geom_smooth(
  method  = "lm", 
  se      = TRUE, 
  formula = y ~ poly(x, 2), 
  color   = "darkred") +
geom_smooth(
  method  = "lm", 
  se      = TRUE, 
  formula = y ~ poly(x, 3), 
  color   = "darkgreen") +
geom_jitter(
  width = 0.5, 
  shape = 1, 
  size  = 1.5, 
  aes(color = Region)) +
  scale_y_continuous(name= expression("Vascular PDGFRβ/Ki67+ cells")) +
  scale_x_continuous(name="DPI",
                     breaks=c(0, 3, 7)) +
  Plot_theme

Ki67_Coloc_Sctr | Ki67_Coloc_Sctr
```
From the graphs above, we can envisage that most Ki67/PDGFRβ are of perivascular nature. As we only have two time points, we will model them as categorical variables. 

Given that our objective is to find out the proportion of perivascular PDGFRβ colocalizing with Ki67, we'll facilitate the modeling by calculating a ratio between `Vascular_Ki67_Pdgfrb` / `Total_Ki67_Pdgfrb`. This will ensure that when the `Total_Ki67_Pdgfrb` is 0, the value for `Vascular_Ki67_Pdgfrb` is also 0. 

```{r}
#| label: Ki67_Ratio
#| include: true
#| warning: false
#| message: false

Ki67_Data$Ratio_Total_Vascular <- Ki67_Data$Vascular_Ki67_Pdgfrb / Ki67_Data$Total_Ki67_Pdgfrb
Ki67_Data$Ratio_Total_Vascular[is.nan(Ki67_Data$Ratio_Total_Vascular)] <- 0
```

# Statistical modeling

## Total Ki67/PDGFRβ 

First, we'll model the total number of Ki67/PDGFRβ colocalized cells as a support to interpret a second model investigating the proportion of perivascular PDGFRβ+ cells colocalized with Ki67. In this case, we use a negative binomial distribution for count data accounting for overdispersion. The model takes the following notation:

$$
Total\_Ki67\_Pdgfrb \sim NegBinomial(\mu, k) \\
\log(\mu) = \beta_0 + \beta_1 \cdot DPI + \beta_2 \cdot Region + \beta_3 \cdot (DPI \times Region) 
$$
Where:¡ $\mu$ is the mean count varying by the interaction between DPI and Region, and $k$ is the dispersion parameter of the negative binomial distribution.

The model uses the default `brms` flat priors.
```{r}
#| label: Ki67_Coloc_Modeling
#| include: true
#| warning: false
#| message: false
#| cache: true

Ki67_Mdl1 <- bf(Total_Ki67_Pdgfrb ~ DPI * Region)

get_prior(Ki67_Mdl1 , Ki67_Data, family = negbinomial())

# Fit model 1
Ki67_Fit1 <- 
  brm(
    data    = Ki67_Data,
    family  = negbinomial(), 
    formula = Ki67_Mdl1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc/Confocal_20x_ROIs_Ki67-Pdgfrb_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Ki67_Fit1 <- 
  add_criterion(Ki67_Fit1, c("loo", "waic", "bayes_R2"))
```

## Perivascular PDGFRβ/Ki67

Now, we model `Vascular_Ki67_Pdgfrb` relative to `Total_Ki67_Pdgfrb` with DPI and Region as predictors for analyzing the proportion of PDGFRβ colocalizing with Ki67. We use a binomial distribution noted as:

$$
Vascular\_Ki67\_Pdgfrb | trials(Total\_Ki67\_Pdgfrb) \sim Binomial(n, p)
$$

where $n$ represents the number of total trials (Ki67/PDGFRβ colocalized cells) and $p$ is the probability of observing perivascular PDGFRβ/Ki67 cells. The model takes the following notation:
$$
\log\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 \cdot DPI + \beta_2 \cdot Region + \beta_3 \cdot (DPI \times Region)
$$
The model uses the default `brms` flat priors.

To facilitate the regression and its interpretation, we exclude the cases wich do not have any PDGFRβ/Ki67 colocalization.

```{r}
#| label: Ki67_Coloc_Excl
#| include: true
#| warning: false
#| message: false
#| cache: true

Ki67_Data_Sub <- Ki67_Data[Ki67_Data$Total_Ki67_Pdgfrb != 0,]
```

New, we fit the model:

```{r}
#| label: Ki67_Coloc_Modeling2
#| include: true
#| warning: false
#| message: false
#| cache: true

Ki67_Mdl2 <- bf(Vascular_Ki67_Pdgfrb | trials(Total_Ki67_Pdgfrb) ~ DPI * Region)

get_prior(Ki67_Mdl2, Ki67_Data_Sub, family = binomial())

# Fit model 1
Ki67_Fit2 <- 
  brm(
    data    = Ki67_Data_Sub,
    family  = binomial(), 
    formula = Ki67_Mdl2,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc/Confocal_20x_ROIs_Ki67-Pdgfrb_Fit2.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Ki67_Fit2 <- 
  add_criterion(Ki67_Fit2, c("loo", "waic", "bayes_R2"))
```

# Model diagnostics

We check the model fitting using posterior predictive checks

```{r}
#| label: fig-Ki67_Pdgfrb_Diagnostics
#| include: true
#| warning: false
#| message: false
#| cache: true
#| fig-cap: Model dianostics using pp_checks
#| fig-height: 5
#| fig-width: 12
 
set.seed(8807)

Ki67_Fit1_pp <- 
  brms::pp_check(Ki67_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Cells ~ DPI * Region") +
  scale_x_continuous(limits = c(0, 10)) +
  Plot_theme  

Ki67_Fit2_pp <- 
  brms::pp_check(Ki67_Fit2, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Vascular| trials(Total) ~ DPI * Region") +
  Plot_theme  
  
Ki67_Fit1_pp | Ki67_Fit2_pp
```
We observe no significant deviation in both cases. We can explore further the model using `shinystan`.

```{r}
#| label: Ki67_Pdgfrb_Shiny
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true

#launch_shinystan(Ki67_Fit1)
#launch_shinystan(Ki67_Fit2)
```

# Model results

We plot the models results using the `conditional_effects` function from `brms`. 

## Visualization of conditional effects

We plot the conditional effects for our first model depicting the total number of Ki67/PDGFRβ 

```{r}
#| label: fig-Ki67_Pdgfrb_CE1
#| include: true
#| warning: false
#| message: false
#| fig-cap: Conditional effects for the number of ki67/PDGFRβ
#| fig-height: 5
#| fig-width: 5

set.seed(8807)

Ki67_Pdgfrb_DPI <- 
  conditional_effects(Ki67_Fit1)

Ki67_Pdgfrb_DPI <- plot(Ki67_Pdgfrb_DPI, 
       plot = FALSE)[[3]]

Ki67_Pdgfrb_pred <- 
  Ki67_Data %>%
  data_grid(DPI,Region, n = 100) %>%
  add_predicted_draws(
    Ki67_Fit1, 
    ndraws = 100)

Ki67_Pdgfrb_fig <- Ki67_Pdgfrb_DPI + 
  scale_y_continuous(name = expression ("PDGFRβ/Ki67+ cells")) +
  scale_x_discrete(name="DPI") +
  
  geom_point(data=Ki67_Data, 
               aes(y = Total_Ki67_Pdgfrb, 
                   x = DPI, colour=Region),
             inherit.aes=FALSE, 
             alpha=0.5,
             size = 1,
             position=position_jitter(h=0, w=0.07)) +
  
  scale_color_manual(
    values = c("#0048BA", "red", "darkgreen"),
    labels = c("Perilesion", "Striatum", "Cortex"),
    name="Region"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red", "darkgreen"),
    labels = c("Perilesion", "Striatum", "Cortex"),
    name="Region"
    ) +
  
  Plot_theme +
  theme(legend.position = c(0.2, 0.8), 
        legend.direction = "vertical")
  
  ggsave(
  plot     = Ki67_Pdgfrb_fig, 
  filename = "Plots/Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc/Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Ki67_Pdgfrb_fig
```

@fig-Ki67_Pdgfrb_CE1 shows an increasing trend in the colocalization of Ki67/PDGFRβ for striatum and cortex. Specifically, the cortex experience an increasing trend in the colocalization at 7 DPI. However, the number of colocalized cells is small compared to the total number of Ki67 or PDGFRβ. Next, we generate the graph for the second model:

```{r}
#| label: fig-Ki67_Pdgfrb_CE2
#| include: true
#| warning: false
#| message: false
#| fig-cap: Conditional effects for the proportion of vascular ki67/PDGFRβ
#| fig-height: 5
#| fig-width: 5

set.seed(8807)

Ki67_Pdgfrb_Vascular <- 
  conditional_effects(Ki67_Fit2)

Ki67_Pdgfrb_Vascular <- plot(Ki67_Pdgfrb_Vascular, 
       plot = FALSE)[[3]]

Ki67_Pdgfrb_Vascular_pred2 <- 
  Ki67_Data %>%
  data_grid(DPI,Region, Total_Ki67_Pdgfrb, n = 100) %>%
  add_predicted_draws(
    Ki67_Fit2, 
    ndraws = 100)

Ki67_Pdgfrb_Vascular_fig <- Ki67_Pdgfrb_Vascular + 
  scale_y_continuous(name = expression ("Vas. PDGFRβ/Ki67+ cells")) +
  scale_x_discrete(name="DPI") +
  
  scale_color_manual(
    values = c("#0048BA", "red", "darkgreen"),
    labels = c("Perilesion", "Striatum", "Cortex"),
    name="Region"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red", "darkgreen"),
    labels = c("Perilesion", "Striatum", "Cortex"),
    name="Region"
    ) +
  
  Plot_theme +
  theme(legend.position = "none")
  
  ggsave(
  plot     = Ki67_Pdgfrb_Vascular_fig, 
  filename = "Plots/Confocal_20x_ROIs_Ki67-Pdgfrb_Coloc/Confocal_20x_ROIs_Ki67-Pdgfrb_Vascular.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Ki67_Pdgfrb_Vascular_fig
```

## Posterior summary

Next, we plot the posterior summary for both models using the `describe_posterior` function:

```{r}
#| label: Ki67_Pdgfrb_Coloc_Posterior
#| include: true
#| warning: false
#| message: false

describe_posterior(
  Ki67_Fit1,
  effects = "all",
  test = c("p_direction", "rope"),
  component = "all",
  centrality = "median")

modelsummary(Ki67_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "Ki67/PDGFRβ+ cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Confocal_20x_ROIs_Ki67-Pdgfrb_Fit1_Table.html",
             )

Ki67_Fit1_Table <- modelsummary(Ki67_Fit1, 
             shape = term ~ model + statistic,
             centrality = "mean", 
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "gt")
gt::gtsave (Ki67_Fit1_Table, 
            filename = "Tables/tex/Confocal_20x_ROIs_Ki67-Pdgfrb_Fit1_Table.tex")


describe_posterior(
  Ki67_Fit2,
  effects = "all",
  test = c("p_direction", "rope"),
  component = "all",
  centrality = "median")

modelsummary(Ki67_Fit2, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "Vascular Ki67/PDGFRβ+ cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Confocal_20x_ROIs_Ki67-Pdgfrb_Fit2_Table.html",
             )

Ki67_Fit2_Table <- modelsummary(Ki67_Fit2, 
             shape = term ~ model + statistic,
             centrality = "mean", 
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "gt")
gt::gtsave (Ki67_Fit2_Table, 
            filename = "Tables/tex/Confocal_20x_ROIs_Ki67-Pdgfrb_Fit2_Table.tex")


```
The tables show the regression estimates. For the first model, estimating the total number of Ki67/PDGFRβ colocalized cells, the output is shown in the log scale. For the second model, estimating the proportion of parenchymal PDGFRβ/Ki67+ cells, the output is sown on the log odds scale.  


# References

::: {#refs}
:::


```{r}
sessionInfo()
```