pseudobulk heatmaps

Figure 3 C,F; Figure 4 D,E,F,H; Figure S3 C; Figure S4 G,H
input data from pseudobulk eset object – output from pseudobulk.normalize.esetlist.R script

DGE_IN_PATH <- "input/differential_expression/"
eset_list <- readRDS("output/dge_lists/pbulk_eset_list_normalized_WCTcourse_metafiltered.rds")
eset_list_sorted <- readRDS("output/dge_lists/pbulk_eset_list_normalized_specific_gating_metafiltered.rds")

# get only T0 samples and remove technical control CHI014
eset_list_T0 <- lapply(eset_list, function(dge){
  filter <- colnames(dge)[dge@phenoData@data$Timepoint %in% c("T0", "HC") & 
                            !is.na(dge@phenoData@data$PC1_cat) & 
                            dge@phenoData@data$Donor != "CHI014"]
  dge <- dge[,filter]
})

eset_list_sorted_T0 <- lapply(eset_list_sorted, function(dge){
  filter <- colnames(dge)[dge@phenoData@data$Timepoint %in% c("T0", "HC") & 
                            !is.na(dge@phenoData@data$PC1_cat) & 
                            dge@phenoData@data$Donor != "CHI014"]
  dge <- dge[,filter]
})

# read in selected_pathway_summary.xlsx
pathways <- read.xlsx("input/selected_pathway_summary.xlsx", sheet = "reduced", startRow = 1)
IFNsets <- filter(pathways, Category %in% c("Type I IFN response", 
                                            "Vaccine_infection induced early blood transcriptomic signature"))
Agpresets <- filter(pathways, Primary_gene_sets %in% c("KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION"))
translationsets <- filter(pathways, Category %in% c("Translation/ribosome"))
NFkBsets <- filter(pathways, Primary_gene_sets %in% c("HALLMARK_TNFA_SIGNALING_VIA_NFKB"))
apoptsets <- filter(pathways, Category %in% c("Apoptosis and cell death"))
inflamsets <- filter(pathways, Primary_gene_sets %in% c("HALLMARK_INFLAMMATORY_RESPONSE"))
fattysets <- filter(pathways, Primary_gene_sets %in% c("reactome_Fatty acid metabolism"))
chemoksets <- filter(pathways, Primary_gene_sets %in% c("KEGG_CHEMOKINE_SIGNALING_PATHWAY","GO_CELLULAR_RESPONSE_TO_INTERLEUKIN_1"))

UnSorted populations

Comparison groups: PC1 & COVIDvsHC T0 only ### Mono_Classical #### IFN showing the shared leading edge genes of PC1 & COVIDvsHC

Mono_Classical <- eset_list_T0$Mono_Classical
# read in fgsea results to get LE genes
Mono_Classical_PC1 <- read.delim(paste0(DGE_IN_PATH, "fgsea_tables/PC1/Mono_Classical--model@PC1--coef@PC1--fgsea.tsv"), header = TRUE)
Mono_Classical_covid <- read.delim(paste0(DGE_IN_PATH, "fgsea_tables/COVID-Healthy/Mono_Classical--model@healthy_vs_covid--coef@COVID-Healthy--fgsea.tsv"), header = TRUE)

Mono_Classical_PC1_IFN <- Mono_Classical_PC1 %>% filter(pathway %in% IFNsets$Primary_gene_sets)
Mono_Classical_covid_IFN <- Mono_Classical_covid %>% filter(pathway %in% IFNsets$Primary_gene_sets)
Mono_Classical_PC1_IFN_LE <- unique(unlist(sapply(Mono_Classical_PC1_IFN$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_covid_IFN_LE <- unique(unlist(sapply(Mono_Classical_covid_IFN$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_IFN_LE <- intersect(Mono_Classical_PC1_IFN_LE,Mono_Classical_covid_IFN_LE)

Mono_Classical_mtx_co <- exprs(Mono_Classical)[Mono_Classical_IFN_LE, ]
meta <- pData(Mono_Classical) %>% dplyr::arrange(Class,PC1)
Mono_Classical_mtx_co <- Mono_Classical_mtx_co[,rownames(meta)]

# pdf("output/Subjecthm_Mono_Classical_coLE_IFN.pdf", width = 12, height = 3.6)
subject.hm(exprs_mtx = Mono_Classical_mtx_co, meta = meta, celltype = "Mono_Classical", module = "IFN")

# dev.off()

translation

showing the shared leading edge genes of PC1 & COVIDvsHC

Mono_Classical_PC1_translation <- Mono_Classical_PC1 %>% filter(pathway %in% translationsets$Primary_gene_sets)
Mono_Classical_covid_translation <- Mono_Classical_covid %>% filter(pathway %in% translationsets$Primary_gene_sets)
Mono_Classical_PC1_translation_LE <- unique(unlist(sapply(Mono_Classical_PC1_translation$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_covid_translation_LE <- unique(unlist(sapply(Mono_Classical_covid_translation$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_translation_LE <- intersect(Mono_Classical_PC1_translation_LE,Mono_Classical_covid_translation_LE)

Mono_Classical_mtx_co <- exprs(Mono_Classical)[Mono_Classical_translation_LE, ]
meta <- pData(Mono_Classical) %>% dplyr::arrange(Class,PC1)
Mono_Classical_mtx_co <- Mono_Classical_mtx_co[,rownames(meta)]

# pdf("output/Subjecthm_Mono_Classical_coLE_translation.pdf", width = 12, height = 4.5)
subject.hm(exprs_mtx = Mono_Classical_mtx_co, meta = meta, celltype = "Mono_Classical", module = "translation")

# dev.off()

NFkB

showing the leading edege gene of PC1 association (severity)

Mono_Classical_PC1_NFkB <- Mono_Classical_PC1 %>% filter(pathway %in% NFkBsets$Primary_gene_sets)
Mono_Classical_covid_NFkB <- Mono_Classical_covid %>% filter(pathway %in% NFkBsets$Primary_gene_sets)
Mono_Classical_PC1_NFkB_LE <- unique(unlist(sapply(Mono_Classical_PC1_NFkB$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_covid_NFkB_LE <- unique(unlist(sapply(Mono_Classical_covid_NFkB$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_NFkB_LE <- intersect(Mono_Classical_PC1_NFkB_LE,Mono_Classical_covid_NFkB_LE)

Mono_Classical_mtx_PC1 <- exprs(Mono_Classical)[Mono_Classical_PC1_NFkB_LE, ]
meta <- pData(Mono_Classical) %>% dplyr::arrange(Class,PC1)
Mono_Classical_mtx_PC1 <- Mono_Classical_mtx_PC1[,rownames(meta)]

# pdf("output/Subjecthm_Mono_Classical_PC1LE_NFkB.pdf", width = 12, height = 7)
subject.hm.PC1(exprs_mtx = Mono_Classical_mtx_PC1, meta = meta, 
               celltype = "Mono_Classical", module = "NFkB",
               PC1LEset = Mono_Classical_PC1_NFkB_LE, 
               covidLEset = Mono_Classical_covid_NFkB_LE)

# dev.off()

HALLMARK Inflammation

showing the leading edege gene of PC1 association (severity)

Mono_Classical_PC1_inflam <- Mono_Classical_PC1 %>% filter(pathway %in% inflamsets$Primary_gene_sets)
Mono_Classical_covid_inflam <- Mono_Classical_covid %>% filter(pathway %in% inflamsets$Primary_gene_sets)
Mono_Classical_PC1_inflam_LE <- unique(unlist(sapply(Mono_Classical_PC1_inflam$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_covid_inflam_LE <- unique(unlist(sapply(Mono_Classical_covid_inflam$leadingEdge, function(x) str_split(x, pattern = " "))))
Mono_Classical_inflam_LE <- intersect(Mono_Classical_PC1_inflam_LE,Mono_Classical_covid_inflam_LE)

Mono_Classical_mtx_PC1 <- exprs(Mono_Classical)[Mono_Classical_PC1_inflam_LE, ]
meta <- pData(Mono_Classical) %>% dplyr::arrange(Class,PC1)
Mono_Classical_mtx_PC1 <- Mono_Classical_mtx_PC1[,rownames(meta)]

# pdf("output/Subjecthm_Mono_Classical_PC1LE_inflam.pdf", width = 12, height = 6.5)
subject.hm.PC1(exprs_mtx = Mono_Classical_mtx_PC1, meta = meta, 
               celltype = "Mono_Classical", module = "inflam",
               PC1LEset = Mono_Classical_PC1_inflam_LE, 
               covidLEset = Mono_Classical_covid_inflam_LE)

# dev.off()

NK_CD16hi

NFkB

showing the leading edege gene of PC1 association (severity)

NK_CD16hi <- eset_list_T0$NK_CD16hi
# read in fgsea results to get LE genes
NK_CD16hi_PC1 <- read.delim(paste0(DGE_IN_PATH, "fgsea_tables/PC1/NK_CD16hi--model@PC1--coef@PC1--fgsea.tsv"), header = TRUE)
NK_CD16hi_covid <- read.delim(paste0(DGE_IN_PATH, "fgsea_tables/COVID-Healthy/NK_CD16hi--model@healthy_vs_covid--coef@COVID-Healthy--fgsea.tsv"), header = TRUE)

NK_CD16hi_PC1_NFkB <- NK_CD16hi_PC1 %>% filter(NK_CD16hi_PC1$pathway %in% NFkBsets$Primary_gene_sets)
NK_CD16hi_covid_NFkB <- NK_CD16hi_covid %>% filter(NK_CD16hi_covid$pathway %in% NFkBsets$Primary_gene_sets)
NK_CD16hi_PC1_NFkB_LE <- unique(unlist(sapply(NK_CD16hi_PC1_NFkB$leadingEdge, function(x) str_split(x, pattern = " "))))
NK_CD16hi_covid_NFkB_LE <- unique(unlist(sapply(NK_CD16hi_covid_NFkB$leadingEdge, function(x) str_split(x, pattern = " "))))
NK_CD16hi_NFkB_LE <- intersect(NK_CD16hi_PC1_NFkB_LE,NK_CD16hi_covid_NFkB_LE)

NK_CD16hi_mtx_PC1 <- exprs(NK_CD16hi)[NK_CD16hi_PC1_NFkB_LE, ]
meta <- pData(NK_CD16hi) %>% dplyr::arrange(Class,PC1)  
NK_CD16hi_mtx_PC1 <- NK_CD16hi_mtx_PC1[,rownames(meta)]

# pdf("output/Subjecthm_NK_CD16hi_PC1LE_NFkB.pdf", width = 11.5, height = 5.5)
subject.hm.PC1(exprs_mtx = NK_CD16hi_mtx_PC1, meta = meta, celltype = "NK_CD16hi", module = "NFkB",
               PC1LEset = NK_CD16hi_PC1_NFkB_LE, covidLEset = NK_CD16hi_covid_NFkB_LE)

# dev.off()

chemokine and IL1 response in severity network

showing the leading edege gene of PC1 association (severity)

NK_CD16hi_PC1_chemok <- NK_CD16hi_PC1 %>% filter(NK_CD16hi_PC1$pathway %in% chemoksets$Primary_gene_sets)
NK_CD16hi_covid_chemok <- NK_CD16hi_covid %>% filter(NK_CD16hi_covid$pathway %in% chemoksets$Primary_gene_sets)
NK_CD16hi_PC1_chemok_LE <- unique(unlist(sapply(NK_CD16hi_PC1_chemok$leadingEdge, function(x) str_split(x, pattern = " "))))
NK_CD16hi_covid_chemok_LE <- unique(unlist(sapply(NK_CD16hi_covid_chemok$leadingEdge, function(x) str_split(x, pattern = " "))))
NK_CD16hi_chemok_LE <- intersect(NK_CD16hi_PC1_chemok_LE,NK_CD16hi_covid_chemok_LE)

NK_CD16hi_mtx_PC1 <- exprs(NK_CD16hi)[NK_CD16hi_PC1_chemok_LE, ]
meta <- pData(NK_CD16hi) %>% dplyr::arrange(Class,PC1)  
NK_CD16hi_mtx_PC1 <- NK_CD16hi_mtx_PC1[,rownames(meta)]

# pdf("pbulk/Subjecthm_NK_CD16hi_PC1LE_chemok.pdf", width = 11.5, height = 5)
subject.hm.PC1(exprs_mtx = NK_CD16hi_mtx_PC1, meta = meta, 
               celltype = "NK_CD16hi", module = "chemok",
               PC1LEset = NK_CD16hi_PC1_chemok_LE, covidLEset = NK_CD16hi_covid_chemok_LE)

# dev.off()

Fatty acid metabolism

showing the leading edege gene of PC1 association (severity)

NK_CD16hi_PC1_fatty <- NK_CD16hi_PC1 %>% filter(NK_CD16hi_PC1$pathway %in% fattysets$Primary_gene_sets)
NK_CD16hi_covid_fatty <- NK_CD16hi_covid %>% filter(NK_CD16hi_covid$pathway %in% fattysets$Primary_gene_sets)
# add IFNg manually as contrast
NK_CD16hi_PC1_fatty_LE <- c(unique(unlist(sapply(NK_CD16hi_PC1_fatty$leadingEdge, function(x) str_split(x, pattern = " ")))), "IFNG")
NK_CD16hi_covid_fatty_LE <- c(unique(unlist(sapply(NK_CD16hi_covid_fatty$leadingEdge, function(x) str_split(x, pattern = " ")))), "IFNG")
NK_CD16hi_fatty_LE <- intersect(NK_CD16hi_PC1_fatty_LE,NK_CD16hi_covid_fatty_LE)

NK_CD16hi_mtx_PC1 <- exprs(NK_CD16hi)[NK_CD16hi_PC1_fatty_LE, ]
meta <- pData(NK_CD16hi) %>% dplyr::arrange(Class,PC1)  
NK_CD16hi_mtx_PC1 <- NK_CD16hi_mtx_PC1[,rownames(meta)]

# set cluster genesets
catabolicgene <- read.xlsx("input/fatty acid LE NK16hi gene categories (biosynthetic vs catabolic).xlsx", sheet = "catabolic", rowNames = TRUE, colNames = TRUE)$Gene.Symbol
biosyngene <- read.xlsx("input/fatty acid LE NK16hi gene categories (biosynthetic vs catabolic).xlsx", sheet = "biosynthetic", rowNames = TRUE, colNames = TRUE)$Gene.Symbol
fatty_NK <- vector(length = length(NK_CD16hi_PC1_fatty_LE))
fatty_NK <- replace(fatty_NK, which(rownames(NK_CD16hi_mtx_PC1) %in% biosyngene), "biosynthetic")
fatty_NK <- replace(fatty_NK, which(rownames(NK_CD16hi_mtx_PC1) %in% catabolicgene), "catabolic")
fatty_NK <- replace(fatty_NK, which(rownames(NK_CD16hi_mtx_PC1) == "IFNG"), "IFNG")
fatty_NK <- replace(fatty_NK, which(fatty_NK == "FALSE"), "others")

# pdf("output/Subjecthm_NK_CD16hi_PC1LE_fatty.pdf", width = 11.5, height = 5)
subject.hm.wsplit(exprs_mtx = NK_CD16hi_mtx_PC1, meta = meta, celltype = "NK_CD16hi", module = "fatty", rowsplit = fatty_NK)

# dev.off()

pDC

apoptosis/cell death

pDC <- eset_list_T0$pDC
pDC_PC1 <- read.delim(paste0(DGE_IN_PATH, "fgsea_tables/PC1/pDC--model@PC1--coef@PC1--fgsea.tsv"), header = TRUE)
pDC_covid <- read.delim(paste0(DGE_IN_PATH, "fgsea_tables/COVID-Healthy/pDC--model@healthy_vs_covid--coef@COVID-Healthy--fgsea.tsv"), header = TRUE)

pDC_PC1_apopt <- pDC_PC1 %>% filter(pDC_PC1$pathway %in% apoptsets$Primary_gene_sets)
pDC_covid_apopt <- pDC_covid %>% filter(pDC_covid$pathway %in% apoptsets$Primary_gene_sets)
pDC_PC1_apopt_LE <- unique(unlist(sapply(pDC_PC1_apopt$leadingEdge, function(x) str_split(x, pattern = " "))))
pDC_covid_apopt_LE <- unique(unlist(sapply(pDC_covid_apopt$leadingEdge, function(x) str_split(x, pattern = " "))))
pDC_apopt_LE <- intersect(pDC_PC1_apopt_LE,pDC_covid_apopt_LE)

pDC_mtx_union <- exprs(pDC)[union(pDC_PC1_apopt_LE,pDC_covid_apopt_LE), ]
meta <- pData(pDC) %>% dplyr::arrange(Class,PC1)
pDC_mtx_union <- pDC_mtx_union[,rownames(meta)]

# customize pDC LE heatmap by highlighting the reactome_Oxidative Stress Induced Senescence geneset
# set marked genesets
PC1_COVID <- rownames(pDC_mtx_union) %in% pDC_apopt_LE
oxi_stress_LE <- filter(pDC_PC1, pathway %in% "reactome_Oxidative Stress Induced Senescence")$leadingEdge
oxi_stress <- rownames(pDC_mtx_union) %in% unlist(str_split(oxi_stress_LE, pattern = " "))
names(oxi_stress) <- rownames(pDC_mtx_union)

label <- which(rownames(pDC_mtx_union) %in% union(pDC_PC1_apopt_LE,pDC_covid_apopt_LE)[1:20])
label_genes <- rownames(pDC_mtx_union)[label]

# set colors
col_fun = colorRamp2(c(-1.5, 0, 1.5), c("#FF00FF", "#000000", "#FFFF00"))
severity.color <- c("Critical-Alive"="#374E55FF","Critical-Deceased"="#DF8F44FF","Moderate-Alive"="#00A1D5FF" ,"Severe-Alive"="#B24745FF","HC" = "#79AF97FF")
PC1.color <- colorRamp2(c(-3, 3.6), c("white", "#e97171"))
Class.color <- c("HC" = "#F8766D", "COVID" = "#eebb4d")
PC1class.color <- c("HC"="#00BA38", "PC1_low"="#619CFF", "PC1_high"="#F8766D")

# set heatmap annotation
ha = HeatmapAnnotation(
  PC1class = meta$PC1_cat,
  PC1 = meta$PC1, 
  Severity_outcome = meta$severity.outcome2,
  col = list(PC1 = PC1.color,
             PC1class = PC1class.color,
             Severity_outcome = severity.color)
)

# pdf("pbulk/Subjecthm_pDC_unionLE_apopt2.pdf", width = 10, height = 4)
p <- Heatmap(pheatmap:::scale_rows(pDC_mtx_union), name = paste("pDC", "apopt", sep = "_"), 
             show_column_names = FALSE,
             top_annotation = ha, cluster_columns = FALSE,
             column_split = meta$PC1_cat,
             row_split = PC1_COVID,
             col = col_fun,
             # row_names_gp = gpar(col = c("blue", "black"), fontsize = c(10, 10)),
             column_gap = unit(4, "mm")
)+
  Heatmap(PC1_COVID + 0, name = "coLE", col = c("0" = "white", "1" = "blue"), 
          show_heatmap_legend = FALSE, width = unit(5, "mm")) +
  Heatmap(oxi_stress + 0, name = "OxiSenes", col = c("0" = "white", "1" = "green"), 
          show_row_names = TRUE, row_names_gp = gpar(fontsize = 10),
          show_heatmap_legend = FALSE, width = unit(5, "mm")) +
  rowAnnotation(foo = anno_mark(at = label,
                                labels = label_genes,
                                labels_gp = gpar(fontsize = 10)))
p

# dev.off()
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] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] circlize_0.4.10      pheatmap_1.0.12      ComplexHeatmap_2.5.4
##  [4] openxlsx_4.1.5       forcats_0.5.0        stringr_1.4.0       
##  [7] dplyr_1.0.0          purrr_0.3.4          readr_1.3.1         
## [10] tidyr_0.8.3          tibble_3.0.3         ggplot2_3.1.1       
## [13] tidyverse_1.2.1      Biobase_2.46.0       BiocGenerics_0.32.0 
## [16] edgeR_3.28.1         limma_3.42.2        
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5          locfit_1.5-9.4      lubridate_1.7.9    
##  [4] lattice_0.20-41     png_0.1-7           assertthat_0.2.1   
##  [7] digest_0.6.25       R6_2.4.1            cellranger_1.1.0   
## [10] plyr_1.8.6          backports_1.1.8     stats4_3.6.1       
## [13] evaluate_0.14       httr_1.4.2          pillar_1.4.6       
## [16] GlobalOptions_0.1.2 rlang_0.4.7         lazyeval_0.2.2     
## [19] readxl_1.3.1        rstudioapi_0.11     S4Vectors_0.24.4   
## [22] GetoptLong_1.0.2    rmarkdown_2.3       munsell_0.5.0      
## [25] broom_0.5.2         compiler_3.6.1      modelr_0.1.4       
## [28] xfun_0.15           pkgconfig_2.0.3     shape_1.4.4        
## [31] htmltools_0.5.0     tidyselect_1.1.0    IRanges_2.20.2     
## [34] fansi_0.4.1         crayon_1.3.4        withr_2.2.0        
## [37] nlme_3.1-148        jsonlite_1.7.0      gtable_0.3.0       
## [40] lifecycle_0.2.0     magrittr_1.5        scales_1.1.1       
## [43] zip_2.0.4           cli_2.0.2           stringi_1.4.6      
## [46] xml2_1.3.2          ellipsis_0.3.1      generics_0.0.2     
## [49] vctrs_0.3.2         RColorBrewer_1.1-2  rjson_0.2.20       
## [52] tools_3.6.1         glue_1.4.1          hms_0.5.3          
## [55] yaml_2.2.1          clue_0.3-57         colorspace_1.4-1   
## [58] cluster_2.1.0       rvest_0.3.5         knitr_1.29         
## [61] haven_2.1.0