Generate Figure 5C
read in enrichment results tables and filter selected pathways for each celltype Enrichment results from PC1/DSM and days_since_onset group interaction model reflecting the difference of DSM high and low group in early and mid time periods
FGSEA_IN_PATH_EARLY <- "input/differential_expression/fgsea_tables/PC1High-low_in_early"
FGSEA_IN_PATH_MID <- "input/differential_expression/fgsea_tables/PC1High-low_in_mid"
remove_celltype_early <- paste(FGSEA_IN_PATH_EARLY, c("/Dblt", "/DPT", "/gated_out", "/PB_Plasmablasts", "/TissueResMemT", "/Unknown", "/Granulocytes"), "--model@PC1_onset_group_interaction--coef@PC1High-low_in_early--fgsea.tsv", sep = "")
remove_celltype_mid <- paste(FGSEA_IN_PATH_MID, c("/Dblt", "/DPT", "/gated_out", "/PB_Plasmablasts", "/TissueResMemT", "/Unknown", "/Granulocytes"), "--model@PC1_onset_group_interaction--coef@PC1High-low_in_mid--fgsea.tsv", sep = "")
files_early <- list.files(FGSEA_IN_PATH_EARLY, full.names = TRUE)
files_early <- files_early[-which(files_early %in% remove_celltype_early)]
files_mid <- list.files(FGSEA_IN_PATH_MID, full.names = TRUE)
files_mid <- files_mid[-which(files_mid %in% remove_celltype_mid)]
fgsea_list_early <- lapply(files_early, function(x) read.delim(x, header = TRUE, sep = "\t"))
names(fgsea_list_early) <- gsub("--model@PC1_onset_group_interaction--coef@PC1High-low_in_early--fgsea\\.tsv", "", basename(files_early))
fgsea_list_mid <- lapply(files_mid, function(x) read.delim(x, header = TRUE, sep = "\t"))
names(fgsea_list_mid) <- gsub("--model@PC1_onset_group_interaction--coef@PC1High-low_in_mid--fgsea\\.tsv", "", basename(files_mid))
# read in selected pathways
pathways <- read.xlsx("input/selected_pathway_summary.xlsx", sheet = "reduced", startRow = 1)
fgsea_list_early_selected <- lapply(fgsea_list_early, function(x) {filter(x, pathway %in% pathways$Primary_gene_sets)})
fgsea_list_mid_selected <- lapply(fgsea_list_mid, function(x) {filter(x, pathway %in% pathways$Primary_gene_sets)})
fgsea_list_selected <- list()
for (i in names(fgsea_list_early_selected)) {
fgsea_list_selected[[i]] <- bind_rows(fgsea_list_early_selected[[i]], fgsea_list_mid_selected[[i]], .id = "TSO_group")
}
fgsea_selected_df <- bind_rows(fgsea_list_selected, .id = "celltype")
fgsea_selected_df$TSO_group <- replace(fgsea_selected_df$TSO_group, fgsea_selected_df$TSO_group == "1", "TSO < d17")
fgsea_selected_df$TSO_group <- replace(fgsea_selected_df$TSO_group, fgsea_selected_df$TSO_group == "2", "TSO d17-d23")
fgsea_selected_df$direction <- ifelse(fgsea_selected_df$NES>0, "up", "dn")
fgsea_selected_df_NK_mono <- filter(fgsea_selected_df, celltype %in% c("NK_CD16hi","Mono_Classical"))
# selected genesets bar graph
fgsea_selected_df_NES_inflam <- filter(fgsea_selected_df_NK_mono,
pathway %in% c("HALLMARK_INFLAMMATORY_RESPONSE",
"HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"reactome_Fatty acid metabolism",
"KEGG_OXIDATIVE_PHOSPHORYLATION"))
fgsea_selected_df_NES_inflam$direction = factor(fgsea_selected_df_NES_inflam$direction, levels = c("up","dn"))
fgsea_selected_df_NES_inflam$TSO_group = factor(fgsea_selected_df_NES_inflam$TSO_group, levels = c("TSO d17-d23","TSO < d17"))
p <- ggplot(fgsea_selected_df_NES_inflam, aes(x = celltype, y = NES, fill = TSO_group, group = TSO_group)) +
geom_bar(stat = "identity", position = "dodge", alpha=0.7) +
coord_flip()+
scale_fill_manual(name="TSO_group", values = c("#bcace5","#4bc26c"))+
facet_wrap(~pathway, scales = "free_x") +
geom_vline(xintercept = 0)+
theme_bw()
p
# ggsave(filename = "output/mono_NK_TSOgroup_NES.pdf", plot = p, device = "pdf", height = 4, width = 7)
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.1.1 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