Generate Figure S5D
# input is the pseudobulk eset_list object
eset_list <- readRDS("output/dge_lists/pbulk_eset_list_normalized_WCTcourse_metafiltered.rds")
combined_genes <- c("TSC22D3")
gene_avg_df <- getgsvascore_list_df(eset_list, 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"
gene_avg_df_T0 <- gene_avg_df %>%
rownames_to_column("sample") %>%
filter(Timepoint %in% c("T0","HC"), pc1_group %in% c("HC", "PC1_low", "PC1_high")) %>%
filter(!str_detect(sample, "^CHI014")) %>%
filter(celltype %in% c("Mono_Classical", "NK_CD16hi"))
gene_avg_df_T0$severity_outcome = as.character(gene_avg_df_T0$severity_outcome)
gene_avg_df_T0$steroid.use = factor(gene_avg_df_T0$steroid.use, levels = c("HC", "FALSE", "TRUE"))
jitter <- position_jitter(width = 0.2, height = 0.1)
# filter out low cell number samples not included in the test
gene_avg_df_T0 <- filter(gene_avg_df_T0, n_barcodes>7)
p <- ggplot(gene_avg_df_T0, aes(x = steroid.use, y = TSC22D3))+
geom_boxplot(outlier.shape=NA, aes(color = steroid.use))+
geom_point(aes(shape = severity_outcome, group = 1, color = steroid.use), position = jitter)+
scale_color_manual(name="steroid.use",values=c("#a2de96", "#3ca59d", "#e79c2a"))+
scale_shape_manual(values = c(15:16,3,17:18))+
facet_wrap(~celltype, scales = "free_y")+
theme(axis.text.x = element_text(angle = 90))+
theme_bw()
p
# ggsave("output/TSC22D3.steroiduse.pdf", plot = p, device = "pdf", width = 7, height = 2.5)
Anova test of the effect of steroid use in COVID-19 patients accounting for severity(PC1/DSM), TSO, Age and experimental batch
gene_avg_df_T0_mono <- gene_avg_df_T0 %>% filter(celltype == "Mono_Classical")
gene_avg_df_T0_NK <- gene_avg_df_T0 %>% filter(celltype == "NK_CD16hi")
gene_avg_df_T0_mono_covid <- gene_avg_df_T0_mono %>% filter(steroid.use != "HC")
gene_avg_df_T0_mono_stat <- car::Anova(aov(TSC22D3 ~ steroid.use+PC1+days_since_onset+Age+Batch, gene_avg_df_T0_mono_covid))
gene_avg_df_T0_mono_stat
## Anova Table (Type II tests)
##
## Response: TSC22D3
## Sum Sq Df F value Pr(>F)
## steroid.use 0.2699 1 1.1513 0.3002
## PC1 0.2824 1 1.2045 0.2897
## days_since_onset 0.0950 1 0.4054 0.5339
## Age 0.1032 1 0.4402 0.5171
## Batch 0.0797 2 0.1699 0.8453
## Residuals 3.5165 15
gene_avg_df_T0_NK_covid <- gene_avg_df_T0_NK %>% filter(steroid.use != "HC")
gene_avg_df_T0_NK_stat <- car::Anova(aov(TSC22D3 ~ steroid.use+PC1+days_since_onset+Age+Batch, gene_avg_df_T0_NK_covid))
gene_avg_df_T0_NK_stat
## Anova Table (Type II tests)
##
## Response: TSC22D3
## Sum Sq Df F value Pr(>F)
## steroid.use 0.0185 1 0.0377 0.8484
## PC1 0.0005 1 0.0010 0.9747
## days_since_onset 0.4168 1 0.8510 0.3692
## Age 0.0355 1 0.0724 0.7911
## Batch 0.6152 2 0.6280 0.5456
## Residuals 8.3275 17
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] car_3.0-6 carData_3.0-3 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] Rcpp_1.0.3 lubridate_1.7.4 assertthat_0.2.1 digest_0.6.25
## [5] R6_2.4.1 cellranger_1.1.0 backports_1.1.5 reprex_0.3.0
## [9] evaluate_0.14 httr_1.4.1 pillar_1.4.3 rlang_0.4.7
## [13] curl_4.3 readxl_1.3.1 rstudioapi_0.11 data.table_1.12.8
## [17] rmarkdown_2.1 labeling_0.3 foreign_0.8-75 munsell_0.5.0
## [21] broom_0.7.0 compiler_3.6.2 modelr_0.1.6 xfun_0.12
## [25] pkgconfig_2.0.3 htmltools_0.4.0 tidyselect_1.0.0 rio_0.5.16
## [29] fansi_0.4.1 crayon_1.3.4 dbplyr_1.4.2 withr_2.1.2
## [33] grid_3.6.2 jsonlite_1.6.1 gtable_0.3.0 lifecycle_0.2.0
## [37] DBI_1.1.0 magrittr_1.5 scales_1.1.0 zip_2.0.4
## [41] cli_2.0.2 stringi_1.4.6 farver_2.0.3 fs_1.3.1
## [45] xml2_1.2.2 ellipsis_0.3.0 generics_0.0.2 vctrs_0.3.4
## [49] openxlsx_4.1.4 tools_3.6.2 glue_1.3.1 hms_0.5.3
## [53] abind_1.4-5 yaml_2.2.1 colorspace_1.4-1 rvest_0.3.5
## [57] knitr_1.28 haven_2.2.0