Geneset enrichment score time trajectory plots
Geneset enrichment gsva score with TSO
Generate Figure 3 D; Figure 4 J,K; Figure 5 E; Figure S3 D
Scores were calculated based on LE genes for each comparison group (see score generation scripts)
### use the score of LE
# input is the files from output folder of sample_gsva_dftoelist.R
PC1_onset_union_gsva_esetlist <- readRDS("output/module_score_gsva/union_PC1_days_since_onset_module_score_gsva_filtered_samples_genes_esetlist.rds")
PC1mid_gsva_esetlist <- readRDS("output/module_score_gsva/PC1High-low_in_mid_module_score_gsva_filtered_samples_genes_esetlist.rds")# plot time since onset geneset enrichment score
plot_onset_module_score_line <- function(input_gsva_esetlist,
LE_score_used = "GSVA PC1 LE score",
fig_out_pdf = "modulescore_time_allcelltype.pdf",
celltype = celltypes,
width = 8.5, height = 6.5) {
p <- list()
pdf(paste(FIG_OUT_PATH,fig_out_pdf,sep = ""), width = width, height = height)
for(cell in celltype){
if(grepl("dblt", cell, ignore.case = T) | cell == "gated_out" | cell == "Unknown"){
next()
}
# print(cell)
eset <- input_gsva_esetlist[[cell]]
dat <- t(exprs(eset))
dat <- as.data.frame(dat) %>% select(intersect(combined_genesets, colnames(dat)))
if(nrow(dat) < 7){
next()
}
dat$days_since_onset <- eset$days_since_onset
dat$pc1_group <- eset$PC1_cat
dat$Subject <- eset$Subject
dat$severity_outcome <- eset$severity.outcome2
dat <- dat %>%
gather(key = module, value = score,
-c(days_since_onset, pc1_group, severity_outcome, Subject)) %>%
filter(!is.na(pc1_group))
tmp.HC <- filter(dat, pc1_group == "HC" & Subject != "CHI014")
means <- tmp.HC %>% dplyr::group_by(module) %>%
dplyr::summarise(median = median(score),
quantile.25 = quantile(score, .25, na.rm = TRUE),
quantile.75 = quantile(score, .75, na.rm = TRUE))
tmp.covid <- filter(dat, pc1_group %in% c("PC1_low", "PC1_high"))
severity.color <- c("Critical-Alive"="#374E55FF","Critical-Deceased"="#DF8F44FF","Moderate-Alive"="#00A1D5FF" ,"Severe-Alive"="#B24745FF","HC" = "#79AF97FF")
p1 <- ggplot(tmp.covid, aes(x = days_since_onset, y = score)) +
# geom_rect(aes(xmin = 17, xmax = 23, ymin = -Inf, ymax = Inf), fill = "#EADCFA", alpha = 0.1)+
geom_point(alpha=.4,shape=21,aes(fill=pc1_group),size=2) +
scale_fill_manual(name="Severity",values = c("deepskyblue1","red")) +
geom_line(aes(group = Subject), alpha = 0.1)+
ggtitle(paste(cell, LE_score_used, sep = " ")) +
scale_color_manual(name="PC1 Class",values=c("deepskyblue1","red")) +
geom_smooth(se = T,aes(color=pc1_group),alpha=0.1) +
geom_hline(aes(yintercept = median), data = means, alpha = 0.7, linetype = "dashed", color = "#00BA38")+
facet_wrap(~module, scales = "free", ncol = 4) +
geom_vline(xintercept = 20, color = "grey60", linetype='dashed') +
theme_bw()
print(p1)
p[[cell]] <- p1
}
dev.off()
return(p)
}
# functions for highlight the TSO d17-23 period
plot_onset_module_score_highlight <- function(input_gsva_esetlist, LE_score_used = "GSVA PC1 LE score",
fig_out_pdf = "modulescore_time_allcelltype.pdf",
celltype = celltypes,
width = 8.5, height = 6.5) {
p <- list()
pdf(paste(FIG_OUT_PATH,fig_out_pdf,sep = ""), width = width, height = height)
for(cell in celltype){
if(grepl("dblt", cell, ignore.case = T) | cell == "gated_out" | cell == "Unknown"){
next()
}
# print(cell)
eset <- input_gsva_esetlist[[cell]]
dat <- t(exprs(eset))
dat <- as.data.frame(dat) %>% select(intersect(combined_genesets, colnames(dat)))
if(nrow(dat) < 7){
next()
}
dat$days_since_onset <- eset$days_since_onset
dat$pc1_group <- eset$PC1_cat
dat$Subject <- eset$Subject
dat$severity_outcome <- eset$severity.outcome2
dat <- dat %>%
gather(key = module, value = score,
-c(days_since_onset, pc1_group, severity_outcome, Subject)) %>%
filter(!is.na(pc1_group))
tmp.HC <- filter(dat, pc1_group == "HC" & Subject != "CHI014")
means <- tmp.HC %>% dplyr::group_by(module) %>%
dplyr::summarise(median = median(score),
quantile.25 = quantile(score, .25, na.rm = TRUE),
quantile.75 = quantile(score, .75, na.rm = TRUE))
tmp.covid <- filter(dat, pc1_group %in% c("PC1_low", "PC1_high"))
severity.color <- c("Critical-Alive"="#374E55FF","Critical-Deceased"="#DF8F44FF","Moderate-Alive"="#00A1D5FF" ,"Severe-Alive"="#B24745FF","HC" = "#79AF97FF")
p1 <- ggplot(tmp.covid, aes(x = days_since_onset, y = score)) +
geom_rect(aes(xmin = 17, xmax = 23, ymin = -Inf, ymax = Inf), fill = "#EADCFA", alpha = 0.1)+
geom_point(alpha=.4,shape=21,aes(fill=pc1_group),size=2) +
scale_fill_manual(name="Severity",values = c("deepskyblue1","red")) +
geom_line(aes(group = Subject), alpha = 0.1)+
ggtitle(paste(cell, LE_score_used, sep = " ")) +
scale_color_manual(name="PC1 Class",values=c("deepskyblue1","red")) +
geom_smooth(se = T,aes(color=pc1_group),alpha=0.1) +
geom_hline(aes(yintercept = median), data = means, alpha = 0.7, linetype = "dashed", color = "#00BA38")+
facet_wrap(~module, scales = "free", ncol = 4) +
# geom_vline(xintercept = 20, color = "grey60", linetype='dashed') +
theme_bw()
print(p1)
p[[cell]] <- p1
}
dev.off()
return(p)
}set interested pathways to plot
using gsva score from union LE of PC1 and TSO to show the effect of severity and time
celltypes <- names(PC1_onset_union_gsva_esetlist)
p <- plot_onset_module_score_line(PC1_onset_union_gsva_esetlist,
LE_score_used = "GSVA PC1 days_onset union LE score",
fig_out_pdf = "modulescore_time_allcelltype_PC1_onsetunionLE.pdf",
width = 10.5, height = 5)
p$Mono_Classical
p$NK_CD16hi
using gsva score from LE of mid(d17-23 juncture) PC1high vs low comparison
highlight the difference at d17-23 juncture
combined_genesets <- c("HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"HALLMARK_INFLAMMATORY_RESPONSE")
celltypes <- names(PC1mid_gsva_esetlist)
p <- plot_onset_module_score_highlight(PC1mid_gsva_esetlist,
LE_score_used = "GSVA PC1_at_mid LE score",
celltype = celltypes,
fig_out_pdf = "modulescore_time_allcelltype_PC1midLE.pdf",
width = 5, height = 2.5)
p$Mono_Classical
p$NK_CD16hi
sI <- sessionInfo()
utils:::print.sessionInfo(sI[-c(10,11)])## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
##
## 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] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] reshape2_1.4.3 Biobase_2.46.0 BiocGenerics_0.32.0
## [4] GSVA_1.34.0 edgeR_3.28.1 limma_3.42.2
## [7] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.4
## [10] purrr_0.3.3 readr_1.3.1 tidyr_1.0.2
## [13] tibble_3.0.3 ggplot2_3.3.2 tidyverse_1.3.0
## [16] plyr_1.8.5 matrixStats_0.55.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-144 bitops_1.0-6 fs_1.3.1
## [4] lubridate_1.7.4 bit64_0.9-7 RColorBrewer_1.1-2
## [7] httr_1.4.1 tools_3.6.2 backports_1.1.5
## [10] R6_2.4.1 DBI_1.1.0 mgcv_1.8-31
## [13] colorspace_1.4-1 withr_2.1.2 tidyselect_1.0.0
## [16] bit_1.1-15.2 compiler_3.6.2 graph_1.64.0
## [19] cli_2.0.2 rvest_0.3.5 xml2_1.2.2
## [22] labeling_0.3 scales_1.1.0 digest_0.6.25
## [25] rmarkdown_2.1 pkgconfig_2.0.3 htmltools_0.4.0
## [28] dbplyr_1.4.2 fastmap_1.0.1 rlang_0.4.7
## [31] readxl_1.3.1 rstudioapi_0.11 RSQLite_2.2.0
## [34] shiny_1.4.0 farver_2.0.3 generics_0.0.2
## [37] jsonlite_1.6.1 RCurl_1.98-1.1 magrittr_1.5
## [40] Matrix_1.2-18 Rcpp_1.0.3 munsell_0.5.0
## [43] S4Vectors_0.24.3 fansi_0.4.1 lifecycle_0.2.0
## [46] stringi_1.4.6 yaml_2.2.1 grid_3.6.2
## [49] blob_1.2.1 promises_1.1.0 crayon_1.3.4
## [52] lattice_0.20-40 haven_2.2.0 splines_3.6.2
## [55] annotate_1.64.0 hms_0.5.3 locfit_1.5-9.1
## [58] knitr_1.28 pillar_1.4.3 geneplotter_1.64.0
## [61] stats4_3.6.2 reprex_0.3.0 XML_3.99-0.3
## [64] glue_1.3.1 evaluate_0.14 modelr_0.1.6
## [67] vctrs_0.3.4 httpuv_1.5.2 cellranger_1.1.0
## [70] gtable_0.3.0 assertthat_0.2.1 xfun_0.12
## [73] mime_0.9 xtable_1.8-4 broom_0.7.0
## [76] later_1.0.0 shinythemes_1.1.2 AnnotationDbi_1.48.0
## [79] memoise_1.1.0 IRanges_2.20.2 ellipsis_0.3.0
## [82] GSEABase_1.48.0