Plotting the Effect Sizes

work_plots = "/Users/katia/OneDrive/Documentos/MountSinai/Projects/Microglia/Figures4paper/Fig_age/age_effect_size/"
lists2compare = "~/pd-omics/katia/scripts/GitHub_scripts/glia_omics/3rd_pass_mic_255s/age_related/age_by_region/"
inter_terms = "~/pd-omics/katia/scripts/GitHub_scripts/glia_omics/3rd_pass_mic_255s/age_interaction/"

deg_list_1 = read.table(file = paste0(lists2compare, "List_age_limma_MFG.txt"), header = T)
deg_list_1 = deg_list_1[order(deg_list_1$adj.P.Val) ,]
#length(which(deg_list_1$adj.P.Val<0.05)) 
deg_list_1$study = "age_MFG"
deg_list_1$ensembl = rownames(deg_list_1)

deg_list_2 = read.table(file = paste0(lists2compare, "List_age_limma_STG.txt"), header = T)
deg_list_2 = deg_list_2[order(deg_list_2$adj.P.Val) ,]
# length(which(deg_list_2$adj.P.Val<0.05)) 
deg_list_2$study = "age_STG"
deg_list_2$ensembl = rownames(deg_list_2)

deg_list_3 = read.table(file = paste0(lists2compare, "List_age_limma_THA.txt"), header = T)
deg_list_3 = deg_list_3[order(deg_list_3$adj.P.Val) ,]
# length(which(deg_list_3$adj.P.Val<0.05)) 
deg_list_3$study = "age_THA"
deg_list_3$ensembl = rownames(deg_list_3)

deg_list_4 = read.table(file = paste0(lists2compare, "List_age_limma_SVZ.txt"), header = T)
deg_list_4 = deg_list_4[order(deg_list_4$adj.P.Val) ,]
# length(which(deg_list_4$adj.P.Val<0.05)) 
deg_list_4$study = "age_SVZ"
deg_list_4$ensembl = rownames(deg_list_4)

inter_terms_list = read.table(file = paste0(inter_terms, "interaction_list_77g.txt"), header = T)

Age - all genes

Beta values of all genes.

deg_list_all = deg_list_1[,c("ensembl","symbol","logFC","study")] %>% 
  left_join(deg_list_2[,c("ensembl","logFC")], by = "ensembl") %>%
  left_join(deg_list_3[,c("ensembl","logFC")], by = "ensembl") %>%
  left_join(deg_list_4[,c("ensembl","logFC")], by = "ensembl")
colnames(deg_list_all) = c("ensembl","symbol","logFC_MFG","MFG","logFC_STG","logFC_THA","logFC_SVZ")

genes_m = melt(deg_list_all, id.vars = c("ensembl","symbol","logFC_MFG","MFG"))

genes_m$variable = gsub("(.*)_(.*)", "\\2", genes_m$variable)
genes_m$variable = factor(genes_m$variable, levels = c("STG","THA","SVZ"))

#pdf(file = paste0(work_plots, "beta_age_all.pdf"), width = 6, height = 4)
ggplot2::ggplot(genes_m, aes(x = logFC_MFG, y = value, color = variable)) + 
  geom_point(alpha = .5) +
  scale_color_manual(values = c("#C71000FF", "#8A4198FF", "#008EA0FF")) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "black") + 
  stat_smooth(method = "lm", se=F) + # Add Regression Line 
  stat_regline_equation(aes(label = ..adj.rr.label..), show.legend = F)  + # Add R-Square
#  stat_regline_equation(aes(label = ..rr.label..))  +
  easy_labs(x = expression(paste("MFG age-related (", beta,")")), y = expression(paste("age-related (", beta, ")"))) +
  easy_add_legend_title("Region") +
  theme_classic()
#dev.off()

Age - filtered 1

Shows only the genes with FDR < 0.05 in MFG.

X = beta in MFG. Y = beta in the other brain regions (by color).

genes_m_filt = genes_m[genes_m$ensembl %in% deg_list_1$ensembl[deg_list_1$adj.P.Val < 0.05], ]

#pdf(file = paste0(work_plots, "beta_age_005_MFG.pdf"), width = 6, height = 4)
ggplot2::ggplot(genes_m_filt, aes(x = logFC_MFG, y = value, color = variable)) + 
  geom_point() +
  scale_color_manual(values = c("#C71000FF", "#8A4198FF", "#008EA0FF")) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "black") + 
  stat_smooth(method = "lm", se=F) + # Add Regression Line 
  stat_regline_equation(aes(label = ..adj.rr.label..), show.legend = F)  + # Add R-Square
#  stat_regline_equation(aes(label = ..rr.label..))  +
  easy_labs(x = expression(paste("MFG age-related (", beta,")")), y = expression(paste("age-related (", beta, ")"))) +
  easy_add_legend_title("Region") +
  theme_classic()
#dev.off()

Age - filtered 2

Shows only the genes with FDR < 0.05 in at least one region + the interaction terms genes (77 genes from MFGxSVZ).

X = beta in MFG. Y = beta in the other brain regions (by color).

genes_m_filt = genes_m[genes_m$ensembl %in% unique(c(deg_list_1$ensembl[deg_list_1$adj.P.Val < 0.05],
                                                     deg_list_2$ensembl[deg_list_2$adj.P.Val < 0.05],
                                                     deg_list_3$ensembl[deg_list_3$adj.P.Val < 0.05],
                                                     deg_list_4$ensembl[deg_list_4$adj.P.Val < 0.05],
                                                   as.character(inter_terms_list$ensembl))), ]

symbol4label = c("MS4A6A", "MERTK", "ANXA1", "MRC1", "LYVE1") 
genes_m_filt4lab = genes_m_filt
genes_m_filt4lab$overlap = genes_m_filt4lab$symbol
genes_m_filt4lab$overlap[(!genes_m_filt4lab$overlap %in% symbol4label) | genes_m_filt4lab$variable!="SVZ"] = ""

# pdf(file = paste0(work_plots, "beta_age_sign_005_oneregion_label_interac.pdf"), width = 6, height =4)
ggplot2::ggplot(genes_m_filt4lab, aes(x = logFC_MFG, y = value, color = variable)) + 
  geom_point() +
  scale_color_manual(values = c("#C71000FF", "#8A4198FF", "#008EA0FF")) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "lightgrey") +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "lightgrey") + 
  stat_smooth(method = "lm", se=F) + # Add Regression Line +
  stat_poly_eq(formula = y ~ x, aes(label = paste(..adj.rr.label..,..p.value.label..,sep = "*`,`~")), parse=TRUE) +
 # stat_regline_equation(aes(label = paste(..adj.rr.label..)), show.legend = F)  + # Add R-Square
#  stat_regline_equation(aes(label = ..rr.label..))  +
  #geom_text_repel(aes(label = overlap),color="black",force = 4) + # Descoment to show the symbol
  geom_label_repel(aes(label = overlap), size = 3, color="black", box.padding = 0.4, label.size = NA, fill = alpha(c("white"),0.5)) + # Descoment to show the symbol
  easy_labs(x = expression(paste("MFG age-related (", beta,")")), y = expression(paste("age-related (", beta, ")"))) +
  easy_add_legend_title("Region") +
  theme_classic()
#dev.off()

Age - filtered 3

Shows the top 100 DE genes in each region (lists might overlap).

X = beta in MFG. Y = beta in the other brain regions (by color).

genes_m_filt = genes_m[genes_m$ensembl %in% unique(c(deg_list_1$ensembl[1:100],
                                                     deg_list_2$ensembl[1:100],
                                                     deg_list_3$ensembl[1:100],
                                                     deg_list_4$ensembl[1:100])), ]

#pdf(file = paste0(work_plots, "beta_age_sign_top100_oneregion.pdf"), width = 6, height = 4)
ggplot2::ggplot(genes_m_filt, aes(x = logFC_MFG, y = value, color = variable)) + 
  geom_point() +
  scale_color_manual(values = c("#C71000FF", "#8A4198FF", "#008EA0FF")) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "black") + 
  geom_text_repel(data = genes_m_filt[genes_m_filt$logFC_MFG*genes_m_filt$value < 0, ], aes(label = symbol)) + # Descoment to show the symbol
  stat_smooth(method = "lm", se=F) + # Add Regression Line 
  stat_regline_equation(aes(label = ..adj.rr.label..), show.legend = F)  + # Add R-Square
#  stat_regline_equation(aes(label = ..rr.label..))  +
  easy_labs(x = expression(paste("MFG age-related (", beta,")")), y = expression(paste("age-related (", beta, ")"))) +
  easy_add_legend_title("Region") +
  theme_classic()
#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] stats graphics grDevices utils datasets methods base

other attached packages: [1] forcats_0.4.0 stringr_1.4.0 purrr_0.3.4 readr_1.3.1
[5] tidyr_1.1.0 tibble_3.0.1 tidyverse_1.3.0 ggpmisc_0.3.5
[9] dplyr_1.0.0 ggeasy_0.1.0 readxl_1.3.1 ggsci_2.9
[13] venn_1.8 ggrepel_0.8.2 ggpubr_0.2.4 magrittr_1.5
[17] reshape2_1.4.4 ggplot2_3.3.2 gplots_3.0.3 RColorBrewer_1.1-2

loaded via a namespace (and not attached): [1] httr_1.4.1 jsonlite_1.6.1 splines_3.6.2
[4] modelr_0.1.5 gtools_3.8.1 assertthat_0.2.1
[7] BiocManager_1.30.10 cellranger_1.1.0 yaml_2.2.0
[10] pillar_1.4.4 backports_1.1.8 lattice_0.20-38
[13] glue_1.4.1 digest_0.6.25 ggsignif_0.6.0
[16] rvest_0.3.5 colorspace_1.4-1 htmltools_0.4.0
[19] Matrix_1.2-18 plyr_1.8.6 pkgconfig_2.0.3
[22] broom_0.5.6 haven_2.2.0 scales_1.1.1
[25] gdata_2.18.0 mgcv_1.8-31 farver_2.0.3
[28] generics_0.0.2 admisc_0.5 ellipsis_0.3.1
[31] withr_2.2.0 cli_2.0.2 crayon_1.3.4
[34] evaluate_0.14 fs_1.3.1 fansi_0.4.1
[37] nlme_3.1-142 xml2_1.2.2 tools_3.6.2
[40] hms_0.5.3 lifecycle_0.2.0 munsell_0.5.0
[43] reprex_0.3.0 compiler_3.6.2 caTools_1.18.0
[46] rlang_0.4.6 grid_3.6.2 rstudioapi_0.11
[49] bitops_1.0-6 labeling_0.3 rmarkdown_2.0
[52] gtable_0.3.0 DBI_1.1.0 polynom_1.4-0
[55] R6_2.4.1 lubridate_1.7.9 knitr_1.26
[58] KernSmooth_2.23-16 stringi_1.4.6 Rcpp_1.0.4.6
[61] vctrs_0.3.1 dbplyr_1.4.2 tidyselect_1.1.0
[64] xfun_0.11