Generate Figure 3E
read in enrichment results tables and filter selected pathways for each celltype Enrichment results from PC1/DSM group and continuous days_since_onset (TSO) interaction model reflecting the association with TSO in DSM high and DSM low groups respectively
# filter selected pathways for each celltype
FGSEA_IN_PATH_HIGH <- "input/differential_expression/fgsea_tables/days_onset_in_PC1high"
FGSEA_IN_PATH_LOW <- "input/differential_expression/fgsea_tables/days_onset_in_PC1low"
FGSEA_IN_PATH_DIFF <- "input/differential_expression/fgsea_tables/diff_days_onset_PC1high-low"
remove_celltype_high <- paste(FGSEA_IN_PATH_HIGH, c("/Dblt", "/DPT", "/gated_out", "/PB_Plasmablasts", "/TissueResMemT", "/Unknown", "/Granulocytes"), "--model@PC1group_onsetcontinous_interaction--coef@days_onset_in_PC1high--fgsea.tsv", sep = "")
remove_celltype_low <- paste(FGSEA_IN_PATH_LOW, c("/Dblt", "/DPT", "/gated_out", "/PB_Plasmablasts", "/TissueResMemT", "/Unknown", "/Granulocytes"), "--model@PC1group_onsetcontinous_interaction--coef@days_onset_in_PC1low--fgsea.tsv", sep = "")
remove_celltype_diff <- paste(FGSEA_IN_PATH_DIFF, c("/Dblt", "/DPT", "/gated_out", "/PB_Plasmablasts", "/TissueResMemT", "/Unknown", "/Granulocytes"), "--model@PC1group_onsetcontinous_interaction--coef@diff_days_onset_PC1High-low--fgsea.tsv", sep = "")
files_high <- list.files(FGSEA_IN_PATH_HIGH, full.names = TRUE)
files_high <- files_high[-which(files_high %in% remove_celltype_high)]
files_low <- list.files(FGSEA_IN_PATH_LOW, full.names = TRUE)
files_low <- files_low[-which(files_low %in% remove_celltype_low)]
files_diff <- list.files(FGSEA_IN_PATH_DIFF, full.names = TRUE)
files_diff <- files_diff[-which(files_diff %in% remove_celltype_diff)]
fgsea_list_high <- lapply(files_high, function(x) read.delim(x, header = TRUE, sep = "\t"))
names(fgsea_list_high) <- gsub("--model@PC1group_onsetcontinous_interaction--coef@days_onset_in_PC1high--fgsea\\.tsv", "", basename(files_high))
fgsea_list_low <- lapply(files_low, function(x) read.delim(x, header = TRUE, sep = "\t"))
names(fgsea_list_low) <- gsub("--model@PC1group_onsetcontinous_interaction--coef@days_onset_in_PC1low--fgsea\\.tsv", "", basename(files_low))
fgsea_list_diff <- lapply(files_diff, function(x) read.delim(x, header = TRUE, sep = "\t"))
names(fgsea_list_diff) <- gsub("--model@PC1group_onsetcontinous_interaction--coef@diff_days_onset_PC1High-low--fgsea\\.tsv", "", basename(files_diff))
# read in selected pathways
pathways <- read.xlsx("input/selected_pathway_summary.xlsx", sheet = "reduced", startRow = 1)
fgsea_list_high_selected <- lapply(fgsea_list_high, function(x) {filter(x, pathway %in% pathways$Primary_gene_sets)})
fgsea_list_low_selected <- lapply(fgsea_list_low, function(x) {filter(x, pathway %in% pathways$Primary_gene_sets)})
fgsea_list_diff_selected <- lapply(fgsea_list_diff, function(x) {filter(x, pathway %in% pathways$Primary_gene_sets)})
fgsea_list_selected <- list()
for (i in names(fgsea_list_high_selected)) {
fgsea_list_selected[[i]] <- bind_rows(fgsea_list_high_selected[[i]], fgsea_list_low_selected[[i]], .id = "PC1_group")
}
fgsea_selected_df <- bind_rows(fgsea_list_selected, .id = "celltype")
fgsea_diff_selected_df <- bind_rows(fgsea_list_diff_selected, .id = "celltype")
fgsea_selected_df$PC1_group <- replace(fgsea_selected_df$PC1_group, fgsea_selected_df$PC1_group == "1", "PC1_high")
fgsea_selected_df$PC1_group <- replace(fgsea_selected_df$PC1_group, fgsea_selected_df$PC1_group == "2", "PC1_low")
fgsea_selected_df$direction <- ifelse(fgsea_selected_df$NES>0, "up", "dn")
fgsea_selected_df_IFN <- filter(fgsea_selected_df, pathway %in% c("GO_RESPONSE_TO_TYPE_I_INTERFERON"))
fgsea_diff_selected_df_IFN <- filter(fgsea_diff_selected_df, pathway %in% c("GO_RESPONSE_TO_TYPE_I_INTERFERON"))
fgsea_selected_df_IFN_wide <- dcast(fgsea_selected_df_IFN, celltype+pathway ~ PC1_group,value.var = "NES") %>%
left_join(select(fgsea_diff_selected_df_IFN, "celltype", "padj"), by = "celltype") %>%
mutate(diff_padj_sig = padj < 0.05)
p <- ggplot(fgsea_selected_df_IFN_wide, aes(x = PC1_low, y = PC1_high, label = celltype)) +
geom_point(aes(color = diff_padj_sig), size = 3) +
scale_color_manual(values=c("#495464", "#f05454")) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "blue")+
geom_hline(yintercept = 0, linetype = "dashed", color = "#9ba4b4")+
geom_vline(xintercept = 0, linetype = "dashed", color = "#9ba4b4")+
xlim(-3.2,1.5)+
ylim(-3.2,1.5)+
theme_bw() +
geom_text_repel(
data = fgsea_selected_df_IFN_wide,
aes(label = celltype),
color = "#495464",
size = 5,
box.padding = unit(0.35, "lines"),
point.padding = unit(0.35, "lines")
)
p
# ggsave(filename = "output/TSO_PC1group_NES.pdf", plot = p, device = "pdf", width=6, height = 5)
sI <- sessionInfo()
utils:::print.sessionInfo(sI[-c(10,11)])
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server 7.2 (Maipo)
##
## Matrix products: default
##
## locale:
## [1] en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.8.1 reshape2_1.4.4 openxlsx_4.1.5 forcats_0.5.0
## [5] stringr_1.4.0 dplyr_1.0.0 purrr_0.3.4 readr_1.3.1
## [9] tidyr_0.8.3 tibble_3.0.3 ggplot2_3.1.1 tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.15 haven_2.1.0 lattice_0.20-41
## [5] colorspace_1.4-1 vctrs_0.3.2 generics_0.0.2 htmltools_0.5.0
## [9] yaml_2.2.1 rlang_0.4.7 pillar_1.4.6 glue_1.4.1
## [13] withr_2.2.0 modelr_0.1.4 readxl_1.3.1 lifecycle_0.2.0
## [17] plyr_1.8.6 munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0
## [21] rvest_0.3.5 zip_2.0.4 evaluate_0.14 labeling_0.3
## [25] knitr_1.29 fansi_0.4.1 broom_0.5.2 Rcpp_1.0.5
## [29] scales_1.1.1 backports_1.1.8 jsonlite_1.7.0 farver_2.0.3
## [33] hms_0.5.3 digest_0.6.25 stringi_1.4.6 grid_3.6.1
## [37] cli_2.0.2 tools_3.6.1 magrittr_1.5 lazyeval_0.2.2
## [41] crayon_1.3.4 pkgconfig_2.0.3 ellipsis_0.3.1 xml2_1.3.2
## [45] lubridate_1.7.9 assertthat_0.2.1 rmarkdown_2.3 httr_1.4.2
## [49] rstudioapi_0.11 R6_2.4.1 nlme_3.1-148 compiler_3.6.1