Effect size (NES) of DSM_high and DSM_low with TSO association respectively

using type-I IFN response as an example

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")

plot IFN NES scatter

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