Fig2.regional_astrocyte_heterogeneity.Rmd

---
title: "Regional astrocyte heterogeneity"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Import packages and set directories
```{r}
library ("data.table")
library("dplyr")
library("textshape")
require(ComplexHeatmap)
require(Seurat)
require(data.table)
require(tidyverse)
require(MAST)
library (ggVennDiagram)
library (ggpubr)
```

Differential analysis was run GPU/cluster on the rPCA integrated astrocytes. Input data will be available upon request.
```{r}
load("../Data/ast_ctrl_rPCA")
# @ ast_ctrl_rPCA - Suerat object which has the astrocyte filtered data
region_de = FindAllMarkers(ast_ctrl_rPCA,
                    min.pct = 0, logfc.threshold = 0, 
                    min.cells.feature = 0, return.thresh = 1, 
                    test.use = "LR", latent.vars = "Donor.ID")
region_de = region_de %>% rownames_to_column("gene_name") %>% setDT()
fwrite(region_de, file.path("../Data/ast_ctrl_de.csv"))
```

```{r}
#load the file - Differential expression results
path1_de = fread(file.path("Example_Data/ast_ctrl_de.csv"))
```

Figure 2a
Venn diagrams of the differentially expressed genes (DEGs)
```{r}
path1_de[, region := factor(region, levels = c("EC", "ITG", "PFC", "V2", "V1"))]
# up
# FC > 1.2
up_gs = path1_de[p_val_adj < 0.05 & avg_log2FC > 0.26 & pct.1 > 0.5, ]
up_gs = split(up_gs, up_gs$region) %>% 
  map(.f = ~(.x$gene))
fig2a_up = ggVennDiagram(
  up_gs, label_alpha=0,
  label = "none",  #set this to counts to view the count, removed for aesthetics
  color = "black", 
  lwd = 0.8, 
  lty = 1,
  set_size = 6
) +
  scale_fill_gradient(
    high = "#FF1700", 
    low = "white",
    name="# of DEGs",
    guide = guide_colorbar(
      frame.colour = "black", 
      ticks.colour = "black",
      frame.linewidth = 1,
      font.size = 5
  )) +
  scale_color_manual(values = rep("black", 5)) +
  theme(legend.text=element_text(size=12),
        legend.title=element_text(size=14))
#ggsave(fig2a_up, file = file.path(p_dir, "fig2c_ast_ctrl_venn_up.pdf"), height = 3.5, width = 4)

# down
# -0.16: log2(1/1.12)
down_gs = path1_de[p_val_adj < 0.05 & avg_log2FC < -0.16 & pct.1 > 0.5, ]
down_gs = split(down_gs, down_gs$region) %>% 
  map(.f = ~(.x$gene))

fig2a_down = ggVennDiagram(
  down_gs, 
  label_alpha=0,
  label = "none", #set this to counts to view the count, removed for aesthetics
  color = "black", 
  lwd = 0.8, lty = 1,
  set_size = 6
) +
  scale_fill_gradient(
    high = "#4B8D4F", 
    low = "white",
    name="# of DEGs",
    guide = guide_colorbar(
      frame.colour = "black", 
      ticks.colour = "black",
      frame.linewidth = 1,
      font_size = 5
  )) +
  scale_color_manual(values = rep("black", 5))  +
  theme(legend.text=element_text(size=12),
        legend.title=element_text(size=14))

ggarrange (fig2a_up, fig2a_down)
```

Figure 2b
Heatmap of DEGs in each brain region (upregulated in red and downregulated in green) ranked by average log fold-change. 
```{r}
require(openxlsx)
path1_region_up_gs = read.xlsx(
  "Example_Data/ast_ctrl_region_de_summary.xlsx",
  sheet = "up genes") %>%
  setDT() %>%
  .[, .(cluster, gene)]
path1_region_down_gs = read.xlsx(
  "Example_Data/ast_ctrl_region_de_summary.xlsx",
  sheet = "down genes") %>%
  setDT() %>%
  .[, .(cluster, gene)]

path1_region_up_gs[, cluster := gsub("SVC", "V2", cluster)]
path1_region_up_gs[, cluster := gsub("PVC", "V1", cluster)]
path1_region_up_gs[, cluster := factor(cluster, levels = c("EC", "ITG", "PFC", "V2", "V1"))]
path1_region_down_gs[, cluster := gsub("PVC", "V1", cluster)]
path1_region_down_gs[, cluster := gsub("SVC", "V2", cluster)]
path1_region_down_gs[, cluster := factor(cluster, levels = c("EC", "ITG", "PFC", "V2", "V1"))]
```

Up-regulated genes heatmap
```{r}
plot_fig2b = function(de, gs, annot_de, cluster, side = "left", 
                      alpha = 0.5, col_fun = NULL){
  # define color
  # require(circlize)
  # col_fun = colorRamp2(c(-1, 0, 1), c("purple", "black", "yellow"))
  
  # prepare data
  mtx = de[gene %in% gs$gene, ] %>%
    dcast(., gene ~ region, value.var = "avg_log2FC") %>%
    as.data.frame() %>%
    column_to_rownames("gene") %>%
    as.matrix()
  mtx = mtx[gs$gene, ]
  
  require(ComplexHeatmap)

  
  ha = rowAnnotation(foo = anno_mark(at = c(30,21,90,14,86,16,28,8,18,65,54,13,157,153,186,195,215,213,217,221,210,211,216,233,260,312,276,256,231,240,282,301,292,300,230,245,302,306,261,193,79,51,102,80,50,47,35,38), labels = up_gs$gene[c(30,21,90,14,86,16,28,8,18,65,54,13,157,153,186,195,215,213,217,221,210,211,216,233,260,312,276,256,231,240,282,301,292,300,230,245,302,306,261,193,79,51,102,80,50,47,35,38)], which = "row", side = "left", labels_gp = grid::gpar(fontsize = 12)))
  
  library(circlize)
  col_fun = colorRamp2(c(-1, 0, 1), c("#5CC93A", "#F5ECEC", "#9F1F1D"))
  col_fun(seq(-3, 3))

  #F5ECEC 
  #require(ComplexHeatmap)
  require(ComplexHeatmap)
  
  p = Heatmap(
    mtx,
    col = col_fun, 
    #row_split = gs$region,
    cluster_columns = F,
    left_annotation = ha,
    cluster_rows = F,
    cluster_row_slices = F,
    show_row_names = F,
    column_names_side = "top",
    column_names_rot = 0,
    column_names_centered = T,
    row_split = cluster,
    row_gap = unit(3, "mm"), border = FALSE,
    name = "average log2 FC",
    heatmap_legend_param = list(
     title_position = "leftcenter-rot",
      legend_height = unit(4, "cm")
    ),
    #legend = F,
    column_names_gp = grid::gpar(fontsize = 12),
    row_names_gp = grid::gpar(fontsize = 12),
    row_title_side = "left"
  )
  
  return(p)
}

up_gs = path1_de[p_val_adj < 0.05 & avg_log2FC > 0.26 & pct.1 > 0.5, ] %>%
  .[, n_sig_region := .N, by = .(gene)] %>%
  .[n_sig_region == 1, ]
up_gs = up_gs[order(region, -avg_log2FC), ]


path1_de[up_gs, on = .(gene), marker_region := i.region]
path1_de[, marker_region := factor(marker_region, levels = c("EC", "ITG", "PFC", "V2", "V1"))]

fig2b_up = plot_fig2b(
  de = path1_de,
  gs = up_gs,
  cluster = up_gs$marker_region,
  annot_de = path1_region_up_gs
)
```

Down-regulated genes heatmap 
```{r}
plot_fig2b = function(de, gs, annot_de, cluster, side = "left", 
                      alpha = 0.5, col_fun = NULL){

  # prepare data
  mtx = de[gene %in% gs$gene, ] %>%
    dcast(., gene ~ region, value.var = "avg_log2FC") %>%
    as.data.frame() %>%
    column_to_rownames("gene") %>%
    as.matrix()
  mtx = mtx[gs$gene, ]
  
  
ha = rowAnnotation(foo = anno_mark(at = c(  152,10,74,51,82,43,106,146,52,75,13,6,25,140,29,65,31,155,177,193,181,203,220,175,179,180,221,187,252,253,230,228,249,237,244,282,311,338,301,325,348,270,344,284,261,296,305,376,281,280,310,75,60,11,55), labels = down_gs$gene[c(  152,10,74,51,82,43,106,146,52,75,13,6,25,140,29,65,31,155,177,193,181,203,220,175,179,180,221,187,252,253,230,228,249,237,244,282,311,338,301,325,348,270,344,284,261,296,305,376,281,280,310,75,60,11,55)],labels_gp = grid::gpar(fontsize = 12))) 
  

  library(circlize)
  col_fun = colorRamp2(c(-1, 0, 1), c("#5CC93A", "#F5ECEC", "#9F1F1D"))
  col_fun(seq(-3, 3))

   
  #require(ComplexHeatmap)
  require(ComplexHeatmap)
  
  p = Heatmap(
    mtx,
    col = col_fun, 
    #row_split = gs$region,
    cluster_columns = F,
    right_annotation = ha,
    cluster_rows = F,
    cluster_row_slices = F,
    show_row_names = F,
    column_names_side = "top",
    column_names_rot = 0,
    column_names_centered = T,
    row_split = cluster,
    row_gap = unit(3, "mm"), border = FALSE,
    #column_split = cluster,
    name = "average log2 FC",
    heatmap_legend_param = list(
      title_position = "leftcenter-rot",
      legend_height = unit(4, "cm")
    ),
    column_names_gp = grid::gpar(fontsize = 12),
    row_names_gp = grid::gpar(fontsize = 12),
    row_title_side = "left"
  )
  
  return(p)
}

down_gs = path1_de[p_val_adj < 0.05 & avg_log2FC < -0.16 & pct.1 > 0.5, ] %>%
  .[, n_sig_region := .N, by = .(gene)] %>%
  .[n_sig_region == 1, ]
down_gs = down_gs[order(region, avg_log2FC), ]

path1_de[down_gs, on = .(gene), marker_region := i.region]

fig2b_down = plot_fig2b(
  de = path1_de,
  gs = down_gs,
  cluster = down_gs$region,
  annot_de = path1_region_down_gs
)
```

The next 9 chunks were on GPU cluster. Input data will be available upon request.

Figure 2d:
Correlation analysis between pTau/Tau and gene expression levels in pathology stage 1 astrocytes. 
```{r setup, include=FALSE}
knitr::opts_knit$set(
  root.dir = "/autofs/space/mindds_001/projects/AbbvieSnRNASeq/scripts/all_brain_regions/MAST"
)
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE, 
  message = FALSE
) 

# directories
data_dir = file.path("..", "..", "..", "data")
rdata_dir = file.path(data_dir, "Rdata")
rdir = file.path("..", "..", "..", "results")
p_dir = file.path(rdir, "plots", "all_brain_regions")

systime = function(){
  paste0("[", Sys.time(),"]")
}

rdata_f = file.path(rdata_dir, "MAST", "seuratobj_ast_sampled_allpath_ptau_04222022.Rdata")
sca_f = file.path(rdata_dir, "MAST", "scaobj_allregion_ast_sampled_path1_ptau_06022022.Rdata")
geneinfo_f = file.path(rdata_dir, "MAST", "geneinfo_allregion_ast_sampled_path1_ptau_06022022.csv")
de_f_r = file.path(rdata_dir, "MAST", "mast_allregion_ast_sampled_path1_ptau_06022022_modelde_")
dt_f_r = file.path(rdir, "MAST", "mast_allregion_ast_sampled_path1_ptau_06022022_modeldt_")
summary_f_r = file.path(rdata_dir, "MAST", "mast_allregion_ast_sampled_path1_ptau_06042022_modelsummary_")

pct_cut = 0.2
var_cut = 0.2

source("cell_level_model_MAST.R")
source("sample_validation.R")
```

- use sampled astrocytes from file: `r rdata_f`
- sca file: `r sca_f`
- use cells from brain region: BA20 & BA46 and path group 1,2,3
- variables:
  - cngeneson: the cellular detection rate of a cell is a number between 0 and 1 representing the fraction of the total number of features measured that have a count > 0 in that cell
  - Braak_stage: 
    - B1: 0, I, II
    - B2: III (no IV from this data)
    - B3: V, VI
  - percent.mito: percentage of mito genes
```{r}
load(rdata_f)
ast_sampled$sample_id = paste(ast_sampled$Donor.ID, ast_sampled$Region, sep = "_")
allregion_path1_sampled = ast_sampled[, ast_sampled$Path..Group. == "1"]
allregion_path1_sampled = prepare_variable(allregion_path1_sampled)
```

Filtering of genes
- filtering genes based on following cut off:
  - gene has count > 0 in at least `r pct_cut * 100`% of the cells
  - gene with variance.standardized > `r var_cut`
Percentage of expression
- expressed genes were defined as if gene has count > 0 in at least `r pct_cut * 100`% of the cells
```{r}
p_dist_allregion_path1 = p_distribution_of_gene(
  allregion_path1_sampled,
  pct_cut = pct_cut
)
p_dist_allregion_path1$p
```

Variance of expression
- cut off for the variant genes:
  - variance.standardized > `r var_cut`
```{r}
allregion_path1_sampled = FindVariableFeatures(allregion_path1_sampled, nfeatures = 10000)
allregion_hvf = HVFInfo(allregion_path1_sampled) %>%
  rownames_to_column("gene")
setDT(allregion_hvf)
allregion_hvf[p_dist_allregion_path1$pct_exp_dt, on = .(gene), pct_exp := i.pct_exp]
```

Distribution of the standardized variance
```{r}
ggplot(
    allregion_hvf,
    aes(x = variance.standardized)
  ) +
    geom_density() +
    geom_vline(aes(xintercept=var_cut),
            color="red", linetype="dashed", size=1) +
    labs(
      x = "variance.standardized",
      title = "dirstribution of variance.standardized of genes across all genes"
    ) +
    ggthemes::theme_few() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )
```

Visualization of selected genes
```{r}
allregion_hvf[, is_model_genes := pct_exp > pct_cut & variance.standardized > var_cut]
n_genes = sum(allregion_hvf$is_model_genes)

ggplot(
  allregion_hvf,
  aes(x = pct_exp, y = variance.standardized,
      color = is_model_genes)
) +
  geom_point(size = 0.4) +
  geom_vline(aes(xintercept=pct_cut),
            color="red", linetype="dashed", size=0.5) +
  geom_hline(aes(yintercept=var_cut),
            color="red", linetype="dashed", size=0.5) +
  ggthemes::theme_few() +
  scale_color_manual(values = c("gray", "red")) +
  labs(title = paste("selected", n_genes, "genes"))
```

Summary
```{r}
allregion_path1_sca = seurat_to_sca(
  allregion_path1_sampled,
  use_genes = allregion_hvf[is_model_genes == T, gene]
)
save(allregion_path1_sca, file = sca_f)
fwrite(allregion_hvf, file = geneinfo_f)
```

MAST analysis
```{r}
test_genes = c("APOE", "HSPB1", "HSP90AA1", "APP")
p_cut = 0.05
logFC_cut = 0
```

f12: ~ (1|Donor.ID/sample_id) + ptau + cngeneson + ratio.mito
```{r}
fit_f = paste0(de_f_r, "ctrl_f12.Rdata")
dt_f = paste0(dt_f_r, "ctrl_f12.csv")
summary_f = paste0(summary_f_r, "ctrl_f12.Rdata")
f12 = "~ (1|Donor.ID/Region) + ptau + cngeneson + ratio.mito"

ptau_de_f12_path1 = mast_zlm(
  sca = allregion_path1_sca,
  f = f12
)
save(ptau_de_f12_path1, file = fit_f)
load(fit_f)

ptau_sum_f12_path1 = sum_zlm(
  zlm_de = ptau_de_f12_path1,
  cont = "ptau"
)
fwrite(ptau_sum_f12_path1$dt, file = dt_f)
save(ptau_sum_f12_path1, file = summary_f)

```

Figure 2d: 
Volcano plot
```{r}
p_volcano = function(
  dt, choose_level, level_col = "levels",
  abslogFC_cut = 0, p_cut = 0.05, pcol = "p"
){
  dt2 <- dt[get(level_col) == choose_level,]
  top10_up = dt2[direction == "up", ] %>%
    .[order(coef, decreasing = T), primerid] %>%
    head(10)
  top10_down = dt2[direction == "down", ] %>%
    .[order(coef, decreasing = T), primerid] %>%
    tail(10)
  dt2[primerid %in% c(top10_down, top10_up), label_genes := primerid]
  dt2[, y := get(pcol)]
  ggplot(
    dt2,
    aes(x = coef, y = -log10(y), color = direction)
  ) +
    geom_point() +
    geom_text_repel(aes(label = label_genes),
                    colour = "black") +
    geom_vline(aes(xintercept = abslogFC_cut), linetype="dashed") +
    geom_vline(aes(xintercept = -abslogFC_cut), linetype="dashed") +
    geom_hline(aes(yintercept = -log10(p_cut)), linetype="dashed") +
    scale_color_manual(
      values = c(
        "up" = "#FF1818",
        "down" = "#08872e",
        "ns" = "gray"
      )
    ) +
    ggthemes::theme_few()
}

dt = fread(file.path("Example_Data/mast_results.csv"))

library(ggrepel)
setnames(dt, "Pr(>Chisq)", "p")
dt[adj.pvalue < 0.05 & coef > 0, direction := "up"]
dt[adj.pvalue < 0.05 & coef < 0, direction := "down"]

fig2d_P1 = p_volcano(
  dt, choose_level = "ptau", 
  level_col = "contrast", pcol = "adj.pvalue"
) +
  labs(x = expression(paste("Coefficient (", ~beta,")")), y = expression("-log"[10]*"FDR")) +
  labs(title = "pTau") +
  theme(plot.title = element_text(hjust = 0.5))



ggsave(fig2d_P1,filename = "../Results/Fig2/fig2d.png", height = 4)
```

Figure 2e: 
Venn diagrams show the number of genes upregulated in EC (EC high) and/or positively correlated with the pTau/Tau ratio (top) and the number of genes downregulated in EC (EC low) and/or genes negatively correlated with pTau/Tau ratio (bottom). 
```{r}
path1_de = fread("Example_Data/ast_ctrl_de.csv")
# load pTau genes
# venn diagram
#' @param a [vector] a vector of category 1
#' @param b [vector] a vector of category 2
#' @param a_name [string] label of category 1
#' @param b_name [string] label of category 2
#' @param colors [vector] a vector of colors for cate 1 & 2
do_fig2e = function(a, b, a_name, b_name, colors){
  require(ggvenn)
  p = ggvenn(
    list(a, b) %>% set_names(a_name, b_name),
    fill_color = colors,
    show_percentage = FALSE,
    text_size = 6,
    columns = c(a_name, b_name)
  ) %>% suppressWarnings()
  return(p)
}

# EC high genes
path1_high_gs = split(path1_de, path1_de$region) %>%
  lapply(., function(de){
    de[p_val_adj < 0.05 & avg_log2FC > 0.26 & pct.1 > 0.5, gene]
  })
ec_high_gs = setdiff(
  path1_high_gs$EC,
  Reduce(union, path1_high_gs[c("ITG", "PFC", "V2", "V1")])
)
# EC low genes
path1_low_gs = split(path1_de, path1_de$region) %>%
  lapply(., function(de){
    de[p_val_adj < 0.05 & avg_log2FC < -0.16 & pct.1 > 0.5, gene]
    # de[p_val_adj < 0.05 & avg_log2FC < 0 & pct.1 > 0.5, gene]
  })
ec_low_gs = setdiff(
  path1_low_gs$EC,
  Reduce(union, path1_low_gs[c("ITG", "PFC", "V2", "V1")])
)


# load 
dt_f = file.path("Example_Data/mast_results.csv")
path1_ptau_dt = fread(dt_f)
path1_ptau_dt[coef > 0 & adj.pvalue < 0.05, direction := "pos"]
path1_ptau_dt[coef < 0 & adj.pvalue < 0.05, direction := "neg"]
path1_ptau_dt[is.na(direction), direction := "ns"]


# venn of EC high & pTau/Tau pos corr genes

figure_2ea<- do_fig2e(
  a = path1_ptau_dt[direction == "pos", primerid],
  b = ec_high_gs,
  a_name = "pTau/Tau pos corr",
  b_name = "EC high",
  colors = c("#EE9595", "#FFCDA3")
)

figure_2eb<- do_fig2e(
  a = path1_ptau_dt[direction == "neg", primerid],
  b = ec_low_gs,
  a_name = "pTau/Tau neg corr",
  b_name = "EC low",
  colors = c("#90E0EF", "#C1F4C5")
)

```