Exploratory plots

PCAs

255 samples | Different brain regions from the same donor | Multi disease cohort.

expression_dir = "~/ad-omics_hydra/microglia_omics/expression_tables/added_pilot_314s/expr_4brain_regions/"
work_plots = "~/pd-omics/katia/scripts/GitHub_scripts/glia_omics/3rd_pass_mic_255s/"
metadata_path = "~/pd-omics/katia/Microglia/mic_255s/metadata_files/"
plots4paper = "/Users/katia/OneDrive/Documentos/MountSinai/Projects/Microglia/Figures4paper/"

load(paste0(expression_dir, "Expression_filt_255s.Rdata")) 
# dim(genes_counts_exp_3rd) # 19376   255

load(paste0(metadata_path, "metadata_255filt_eng_29jun2020.Rdata"))
# str(metadata3rd_pass)
metadata <- metadata3rd_pass
# To use **cause_of_death** in the model, I got the NAs and change it for "Other" category. 
#To use C1-C4 I calculated the median for the missing samples. 

metadata$cause_of_death_categories[metadata$cause_of_death_categories %in% NA] <- "Other"
#table(metadata$cause_of_death_categories)

metadata$C1 = metadata$C1 %>% replace_na(median(metadata$C1, na.rm = T))
metadata$C2 = metadata$C2 %>% replace_na(median(metadata$C2, na.rm = T))
metadata$C3 = metadata$C3 %>% replace_na(median(metadata$C3, na.rm = T))
metadata$C4 = metadata$C4 %>% replace_na(median(metadata$C4, na.rm = T))

PCAs

PCs 1 and 2. Voom without correction.

res.pca = prcomp(t(genes_counts_voom_3rd)) 

g1 <- autoplot(res.pca, data = metadata, colour = 'tissue') +
  scale_color_futurama() +
  easy_add_legend_title("Region") +
theme_classic()

g2 <- autoplot(res.pca, data = metadata, colour = 'cause_of_death_categories') +
  scale_colour_igv() +
  easy_add_legend_title("Cause of death") +
  theme_classic() 

g3 <- autoplot(res.pca, data = metadata, colour = 'main_diagnosis') +
  scale_colour_manual(values = c("darkviolet", "orange", "green", "darkblue", "lightblue", "red", "pink", "yellow", "#5C88DAFF")) +
  easy_add_legend_title("Diagnosis") +
  theme_classic()

g4 <- autoplot(res.pca, data = metadata, colour = 'pmd_minutes') +
  scale_color_viridis_c() +
  easy_add_legend_title("PMD minutes") +
  theme_classic()

g5 <- autoplot(res.pca, data = metadata, colour = 'batch') +
theme_classic() 

# pdf(paste0(plots4paper, "SUPP_PCAs_voom_255s.pdf"), width=12, height=8)
#grid.arrange(g1,g2,g3,g4, ncol = 2)
ggarrange(g1,g2,g3,g4, ncol=2, nrow = 2, labels = c("a", "b", "c", "d"))
# dev.off()
# plot(g2)

Remove BatchEffects

Region (n = 255)

a = Voom without correction. b = After removeBatchEffect without regressing out tissue.

allResiduals <- removeBatchEffect(x = genes_counts_voom_3rd, 
                                  batch = metadata$sex, 
                                  batch2 = metadata$cause_of_death_categories,
                                  design = model.matrix(~ tissue, data = metadata), #force to not regress tissue. 
                                  covariates = as.matrix(metadata[, c("picard_pct_mrna_bases", "picard_summed_median", "picard_pct_ribosomal_bases", "C1","C2","C3","C4" )]))

res.pca = prcomp(t(allResiduals)) 

# pdf(paste0(plots4paper, "PCA_255s_residuals.pdf"), width=4, height=3)
g6 <- autoplot(res.pca, data = metadata, colour = 'tissue') +
  scale_colour_futurama() +
  easy_add_legend_title("Region") +
 # scale_colour_viridis_c() +
  theme_classic()
# dev.off()

# pdf(paste0(plots4paper, "PCAs_255s_voomxresiduals.pdf"), width=10, height=4)
ggarrange(g1,g6, labels = c("a", "b"))
# dev.off()

Remove BatchEffects

Age (n = 255)

a = Voom without correction. b = After removeBatchEffect without regressing out age -> colored by age.

# MFG = #FF6F00FF
# STG = #C71000FF
# SVZ = #008EA0FF
# THA = #8A4198FF

res.pca = prcomp(t(genes_counts_voom_3rd)) 

k1 <- autoplot(res.pca, data = metadata, colour = 'age') +
 # scale_color_gradient2(low = "#008EA0FF", midpoint = 60, mid = "#FF6F00FF", high = "#8A4198FF", na.value = NA ) +
  scale_color_viridis_c(option = "magma") + 
  easy_add_legend_title("Age") +
  theme_classic()

allResiduals <- removeBatchEffect(x = genes_counts_voom_3rd, 
                                  batch = metadata$sex, 
                                  batch2 = metadata$cause_of_death_categories,
                                  design = model.matrix(~ tissue, data = metadata), #force to not regress tissue. 
                                  covariates = as.matrix(metadata[, c("picard_pct_mrna_bases", "picard_summed_median", "picard_pct_ribosomal_bases", "C1","C2","C3","C4" )]))

res.pca = prcomp(t(allResiduals)) 

k2 <- autoplot(res.pca, data = metadata, colour = 'age') +
  scale_color_viridis_c(option = "magma") + 
  easy_add_legend_title("Age") +
  theme_classic()

# pdf(paste0(plots4paper, "PCAs_255s_voomxresiduals_Age.pdf"), width=10, height=4)
ggarrange(k1,k2, labels = c("a", "b"))
# dev.off()

(n = 112)

Only the donors who share the four main brain regions.
a = Voom without correction. b = After removeBatchEffect without regressing out tissue.

Covariates removed: sex, cause_of_death_categories, picard_pct_mrna_bases, picard_summed_median, picard_pct_ribosomal_bases and C1-C4. Without regressing out tissue.

metadata_numb <- metadata
donors_4tissues <- metadata_numb %>% group_by(donor_id) %>% summarise(n=n()) %>% filter(n==4)
metadata_numb <- metadata[metadata$donor_id %in% donors_4tissues$donor_id ,]
metadata_numb$cause_of_death_categories <- as.factor(as.character(metadata_numb$cause_of_death_categories))

genes_voom_numb <- genes_counts_voom_3rd[, colnames(genes_counts_voom_3rd) %in% metadata_numb$donor_tissue]
# all(colnames(genes_voom_numb) == metadata_numb$donor_tissue) # Check the order 
# dim(genes_voom_numb)

res.pca = prcomp(t(genes_voom_numb)) 

g7 <- autoplot(res.pca, data = metadata_numb, colour = 'tissue') +
  scale_colour_futurama() +
  easy_add_legend_title("Region") +
  theme_classic()

allResiduals <- removeBatchEffect(x = genes_voom_numb, 
                                  batch = metadata_numb$sex, 
                                  batch2 = metadata_numb$cause_of_death_categories,
                                 # design = model.matrix(~ tissue, data = metadata_numb), #force to not regress tissue. 
                                  covariates = as.matrix(metadata_numb[, c("picard_pct_mrna_bases", "picard_summed_median", "picard_pct_ribosomal_bases", "C1","C2","C3","C4" )]))

res.pca = prcomp(t(allResiduals)) 

g8 <- autoplot(res.pca, data = metadata_numb, colour = 'tissue') +
  scale_colour_futurama() +
  easy_add_legend_title("Region") +
  theme_classic()

# pdf(paste0(plots4paper, "PCA_112s_voomxresiduals.pdf"), width=10, height=4)
ggarrange(g7,g8, labels = c("a", "b"))
# dev.off()

Donors are numbers

PCA with donors who share the 4 main brain tissues. 112 samples from 28 donors.

a = Voom without correction. b = After removeBatchEffect without regressing out tissue.

metadata_numb$donor_numb = as.numeric(as.factor(metadata_numb$donor_id)) # Gera números por donor! 

metadata_numb$tissue_sm = as.character(metadata_numb$tissue)
metadata_numb$tissue_sm[metadata_numb$tissue_sm == "MFG"] <- "m"
metadata_numb$tissue_sm[metadata_numb$tissue_sm == "STG"] <- "s"
metadata_numb$tissue_sm[metadata_numb$tissue_sm == "THA"] <- "t"
metadata_numb$tissue_sm[metadata_numb$tissue_sm == "SVZ"] <- "z"
metadata_numb$tissue_sm <- as.factor(metadata_numb$tissue_sm)

metadata_numb$donor_numb2 <- paste0(metadata_numb$donor_numb, metadata_numb$tissue_sm)
rownames(metadata_numb) <- metadata_numb$donor_numb2

genes_voom_numb <- genes_counts_voom_3rd[, colnames(genes_counts_voom_3rd) %in% metadata_numb$donor_tissue]
# all(colnames(genes_voom_numb) == metadata_numb$donor_tissue) # Check the order 
colnames(genes_voom_numb) <- rownames(metadata_numb)

res.pca = prcomp(t(genes_voom_numb)) 

g9 <- autoplot(res.pca, data = metadata_numb, colour = 'tissue', label=T, shape = F) +
  scale_colour_futurama() +
  easy_add_legend_title("Region") +
 # scale_colour_viridis_c() +
  theme_classic()
 
############# Using residuals 
residuals_numb <- allResiduals[, colnames(allResiduals) %in% metadata_numb$donor_tissue ] #Calculated in the later chunck for 112 samples
# all(colnames(residuals_numb) == metadata_numb$donor_tissue)
colnames(residuals_numb) <- rownames(metadata_numb)

res.pca = prcomp(t(allResiduals)) 

g10 <- autoplot(res.pca, data = metadata_numb, colour = 'tissue', label=T, shape = F) +
  scale_colour_futurama() +
  easy_add_legend_title("Region") +
 # scale_colour_viridis_c() +
  theme_classic()

#pdf(paste0(plots4paper, "PCA_112s_voomxresiduals_numbers.pdf"), width=10, height=4)
ggarrange(g9,g10, labels = c("a", "b"))
#dev.off()

sessionInfo()

R version 3.6.2 (2019-12-12) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Catalina 10.15.5

Matrix products: default BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages: [1] grid stats graphics grDevices utils datasets methods
[8] base

other attached packages: [1] edgeR_3.28.0 limma_3.42.0 ggpubr_0.2.4 magrittr_1.5
[5] ggeasy_0.1.0 tidyr_1.1.0 factoextra_1.0.6 ggsci_2.9
[9] ggfortify_0.4.10 gridExtra_2.3 RColorBrewer_1.1-2 ggplot2_3.3.2
[13] kableExtra_1.1.0 dplyr_1.0.0 knitr_1.26

loaded via a namespace (and not attached): [1] Rcpp_1.0.4.6 pillar_1.4.4 compiler_3.6.2 tools_3.6.2
[5] digest_0.6.25 lattice_0.20-38 evaluate_0.14 lifecycle_0.2.0
[9] tibble_3.0.1 gtable_0.3.0 viridisLite_0.3.0 pkgconfig_2.0.3
[13] rlang_0.4.6 rstudioapi_0.11 ggrepel_0.8.2 yaml_2.2.0
[17] xfun_0.11 withr_2.2.0 stringr_1.4.0 httr_1.4.1
[21] xml2_1.2.2 generics_0.0.2 vctrs_0.3.1 hms_0.5.3
[25] cowplot_1.0.0 locfit_1.5-9.1 webshot_0.5.2 tidyselect_1.1.0 [29] glue_1.4.1 R6_2.4.1 rmarkdown_2.0 farver_2.0.3
[33] readr_1.3.1 purrr_0.3.4 scales_1.1.1 ellipsis_0.3.1
[37] htmltools_0.4.0 rvest_0.3.5 colorspace_1.4-1 ggsignif_0.6.0
[41] labeling_0.3 stringi_1.4.6 munsell_0.5.0 crayon_1.3.4