Geneset enrichment gsva scores and cytokine correlations

Generate Figure 4 B,C,G,I,L

NK_CD16hi fatty acid metabolism correlations with inflammatory pathways and circulating IL15 level

PC1t0score_list_df_NK <- filter(PC1t0score_list_df, celltype == "NK_CD16hi") %>%
  melt(id = colnames(.)[-c(2,3,5,6,7,19)])

p <- ggplot(PC1t0score_list_df_NK, aes(x = value, y = `reactome_Fatty acid metabolism`)) +
  geom_point(shape=21,aes(fill=severity_outcome),size=3, color="white")+
  scale_fill_manual(name="Severity_outcome",values = c("#374E55FF","#DF8F44FF","#00A1D5FF","#B24745FF")) +
  geom_smooth(se = F,method = "lm") + 
  stat_cor(method = "pearson") +
  facet_wrap(~variable, scales = "free") +
  theme_bw()

p

# ggsave(filename = "FAvs.score.IL15.NK.pdf", plot = p, device = "pdf", width = 9, height = 5)

NK_CD16hi fatty acid metabolism correlations with IFG mRNA level

input is the pseudobulk eset_list object from pseudobulk.normalize.esetlist.R script

# input is the pseudobulk eset_list object from pseudobulk.normalize.esetlist.R script
eset_list <- readRDS("output/dge_lists/pbulk_eset_list_normalized_WCTcourse_metafiltered.rds")

# take T0 only and remove technical CHI014 control
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]
})

combined_genes <- c("IFNG")
exprs_list_df <- getgsvascore_list_df(eset_list_T0, combined_genesets = combined_genes)
## [1] "B_Mem"
## [1] "B_Naive"
## [1] "CD4_Mem"
## [1] "CD4_Naive"
## [1] "CD8_Mem"
## [1] "CD8_Naive"
## [1] "cDC"
## [1] "DNT"
## [1] "DPT"
## [1] "gammadeltaT"
## [1] "Granulocytes"
## [1] "MAIT"
## [1] "Mono_Classical"
## [1] "Mono_Intermediate"
## [1] "Mono_NonClassical"
## [1] "NK_CD16hi"
## [1] "NK_CD56hiCD16lo"
## [1] "NK_CD56loCD16lo"
## [1] "PB_Plasmablasts"
## [1] "pDC"
## [1] "Platelets"
## [1] "RBC"
## [1] "Tcell"
## [1] "TissueResMemT"
## [1] "Treg"
exprs_list_df <- exprs_list_df %>% filter(Timepoint == "T0",pc1_group != "HC") %>%
  mutate(sample_id = paste(Subject, Timepoint, sep = "_"))

exprs_list_df_NK <- exprs_list_df %>%
  filter(celltype == "NK_CD16hi") %>%
  left_join(PC1t0score_list_df, by = c("sample_id","celltype"))

p <- ggplot(exprs_list_df_NK, aes(x = IFNG, y = `reactome_Fatty acid metabolism`)) +
  geom_point(shape=21,aes(fill=severity_outcome.x),size=3, color="white")+
  scale_fill_manual(name="Severity_outcome",values = c("#374E55FF","#DF8F44FF","#00A1D5FF","#B24745FF")) +
  geom_smooth(se = F,method = "lm") + 
  stat_cor(method = "pearson") +
  theme_bw()

p

# ggsave(filename = "FAvs.score.IL15.NK.pdf", plot = p, device = "pdf", width = 5, height = 3.3)

circulating IL15 correlation with severity at T0

#extract only CITEseq T0 samples
IL15_T0 <- filter(IL15_T0, !is.na(subject_id)) %>%
  filter(sample_id %in% unique(PC1t0score_list_df$sample_id))
IL15_T0$severity.outcome <- paste(IL15_T0$severity, IL15_T0$outcome, sep="-")
p <- ggplot(IL15_T0, aes(x = test_value_log10, y = PC1)) +
  geom_point(shape=21,aes(fill=severity.outcome),size=3, color="white")+
  scale_fill_manual(name="Severity_outcome",values = c("#374E55FF","#DF8F44FF","#00A1D5FF","#B24745FF"))+
  geom_smooth(se = F,method = "lm") + 
  stat_cor(method = "pearson") +
  theme_bw()

p

# ggsave(filename = "output/IL15vsPC1.pdf", plot = p, device = "pdf", width = 5, height = 3.3)

circulating IL15 change overtime – figure 4L

covid19.samples.pbmc <- covid19.samples[covid19.samples$material_type == "PBMC",] %>%
  mutate(sample_id = paste(subject_id, visit, sep = "_"))

IL15_alltime <- IL15 %>% 
  filter(subject_id %in% unique(PC1$subject_id)) %>%
  group_by(sample_id) %>% 
  dplyr::summarise_at(vars(test_value_log10), mean) %>%
  left_join(covid19.samples.pbmc, by = "sample_id") %>%
  filter(!is.na(days_from_symptom_onset_to_sample_drawn)) %>%
  left_join(PC1[,which(colnames(PC1) %in% c("subject_id", "PC1", "PC1class"))], by = "subject_id") %>%
  mutate(PC1class = factor(PC1class, levels = c("PC1_low","PC1_high")))

p <- ggplot(IL15_alltime)+
  geom_point(alpha=.4, shape=21,
             aes(x = days_from_symptom_onset_to_sample_drawn, y = test_value_log10, fill=PC1class), size = 3.5)+
  scale_fill_manual(name="Severity", values = c("deepskyblue1","red")) +
  geom_line(aes(x = days_from_symptom_onset_to_sample_drawn, y = test_value_log10, group = subject_id), alpha = 0.2)+
  stat_summary(aes(x = days_from_symptom_onset_to_sample_drawn, y = test_value_log10, group = PC1class), fun = median, alpha=0)+
  stat_smooth(aes(x = days_from_symptom_onset_to_sample_drawn, y = test_value_log10, group = PC1class, color=PC1class), se = TRUE, alpha=0.1, size = 1.5)+# group = 1 overwrite group id
  scale_color_manual(name="PC1 Class",values=c("deepskyblue1","red"))+
  geom_vline(xintercept = 20, color = "grey60", linetype='dashed')+
  theme(axis.text.x = element_text(angle = 90))+
  facet_wrap(~condition)+
  theme_bw()
p

# ggsave(filename = "output/IL15.TSO.pdf", plot = p, device = "pdf", width = 5, height = 3.3)
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] gridExtra_2.3       ggpubr_0.4.0.999    reshape2_1.4.3     
##  [4] Biobase_2.46.0      BiocGenerics_0.32.0 forcats_0.4.0      
##  [7] stringr_1.4.0       dplyr_0.8.4         purrr_0.3.3        
## [10] readr_1.3.1         tidyr_1.0.2         tibble_3.0.3       
## [13] ggplot2_3.3.2       tidyverse_1.3.0     plyr_1.8.5         
## [16] matrixStats_0.55.0 
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.1        jsonlite_1.6.1    splines_3.6.2     carData_3.0-3    
##  [5] modelr_0.1.6      assertthat_0.2.1  cellranger_1.1.0  yaml_2.2.1       
##  [9] pillar_1.4.3      backports_1.1.5   lattice_0.20-40   glue_1.3.1       
## [13] digest_0.6.25     ggsignif_0.6.0    rvest_0.3.5       colorspace_1.4-1 
## [17] htmltools_0.4.0   Matrix_1.2-18     pkgconfig_2.0.3   broom_0.7.0      
## [21] haven_2.2.0       scales_1.1.0      openxlsx_4.1.4    rio_0.5.16       
## [25] mgcv_1.8-31       generics_0.0.2    farver_2.0.3      car_3.0-6        
## [29] ellipsis_0.3.0    withr_2.1.2       cli_2.0.2         magrittr_1.5     
## [33] crayon_1.3.4      readxl_1.3.1      evaluate_0.14     fs_1.3.1         
## [37] fansi_0.4.1       nlme_3.1-144      rstatix_0.6.0     xml2_1.2.2       
## [41] foreign_0.8-75    tools_3.6.2       data.table_1.12.8 hms_0.5.3        
## [45] lifecycle_0.2.0   munsell_0.5.0     reprex_0.3.0      zip_2.0.4        
## [49] compiler_3.6.2    rlang_0.4.7       grid_3.6.2        rstudioapi_0.11  
## [53] labeling_0.3      rmarkdown_2.1     gtable_0.3.0      abind_1.4-5      
## [57] DBI_1.1.0         curl_4.3          R6_2.4.1          lubridate_1.7.4  
## [61] knitr_1.28        stringi_1.4.6     Rcpp_1.0.3        vctrs_0.3.4      
## [65] dbplyr_1.4.2      tidyselect_1.0.0  xfun_0.12