Compare cytokine profiles between deceased patients and those with good outcome. Specific focus on 17-23 days since symptom onset.

Create sub-cohort

Identify all deceased patients and match them with critical but recovered patients by age and gender (excluding those staying hospitalized for more than 50 days). Patients in the citiseq cohort can optionally be excluded.

load(file.path(data.folder,"covid19.metadata.paper1.RData"),verbose=F)
covid19.lab.results$severity.outcome <- paste0(covid19.lab.results$severity,"-",covid19.lab.results$outcome)

severity.color <- c("Critical-Alive"="#374E55FF","Critical-Deceased"="#DF8F44FF","Moderate-Alive"="#00A1D5FF" ,"Severe-Alive"="#B24745FF","HC" = "#79AF97FF")

deceased.patients.cytokines <- subset(covid19.lab.results,outcome == "Deceased" & test_type == "ELISA" & !is.na(days_from_symptom_onset_to_test))
deceased.patients <- unique(deceased.patients.cytokines$subject_id)

cat("Map of sampling TSO time points for deceased patients:\n")
## Map of sampling TSO time points for deceased patients:
deceased.time.points <- as.data.frame.matrix(table(deceased.patients.cytokines[,c("subject_id","days_from_symptom_onset_to_test")])) 
colSums(deceased.time.points > 0)
##  3  8 10 11 12 13 14 16 17 18 19 20 21 22 23 24 25 26 28 29 30 32 35 36 39 42 
##  1  1  1  3  1  2  1  2  1  1  2  3  1  4  1  3  3  2  1  2  1  1  2  2  2  1 
## 45 47 53 
##  1  2  1
cat("Age and sex distribution:\n")
## Age and sex distribution:
table(subset(covid19.patient.demo,subject_id %in% deceased.patients)$sex)
## 
##  F  M 
##  2 15
summary(subset(covid19.patient.demo,subject_id %in% deceased.patients)$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   46.00   61.00   66.00   66.24   73.00   78.00
# recovered patients with similar age and sex
age.tol <- 3
recovered.patients <- c()
for (i in deceased.patients) {
  target.age <- covid19.patient.demo[i,"age"]
  target.sex <- covid19.patient.demo[i,"sex"]
  cat("Target age:",target.age,"sex:",target.sex,"\n")
  tmp <- subset(covid19.lab.results,outcome != "Deceased" & test_type == "ELISA" & !is.na(days_from_symptom_onset_to_test) & 
           age >= target.age - age.tol & age <= target.age + age.tol & severity == "Critical" & sex == target.sex & 
             days_from_symptom_onset_to_test > 15 & days_from_symptom_onset_to_test <= 25 &
           number_of_days_hospitalized < 50)

  matched.subject <- unique(tmp$subject_id)
  cat("Matched subject:",matched.subject,"\n")
  recovered.patients <- c(recovered.patients,matched.subject)
}
## Target age: 73 sex: M 
## Matched subject: HGR0000015 HGR0000088 HGR0000093 HGR0000101 
## Target age: 65 sex: M 
## Matched subject: HGR0000039 HGR0000077 
## Target age: 66 sex: M 
## Matched subject: HGR0000063 HGR0000085 HGR0000149 
## Target age: 59 sex: M 
## Matched subject: HGR0000039 HGR0000049 HGR0000077 HGR0000102 HGR0000108 HGR0000135 
## Target age: 46 sex: M 
## Matched subject: HGR0000032 HGR0000092 HGR0000142 HGR0000145 
## Target age: 66 sex: M 
## Matched subject: HGR0000063 HGR0000085 HGR0000149 
## Target age: 78 sex: M 
## Matched subject: HGR0000088 HGR0000089 
## Target age: 61 sex: M 
## Matched subject: HGR0000039 HGR0000077 HGR0000102 HGR0000135 
## Target age: 68 sex: M 
## Matched subject: HGR0000063 HGR0000085 HGR0000093 HGR0000149 
## Target age: 71 sex: F 
## Matched subject:  
## Target age: 66 sex: M 
## Matched subject: HGR0000063 HGR0000085 HGR0000149 
## Target age: 49 sex: F 
## Matched subject: HGR0000044 HGR0000075 HGR0000139 
## Target age: 73 sex: M 
## Matched subject: HGR0000015 HGR0000088 HGR0000093 HGR0000101 
## Target age: 76 sex: M 
## Matched subject: HGR0000015 HGR0000088 HGR0000089 HGR0000101 
## Target age: 60 sex: M 
## Matched subject: HGR0000039 HGR0000049 HGR0000077 HGR0000102 HGR0000135 
## Target age: 72 sex: M 
## Matched subject: HGR0000015 HGR0000063 HGR0000085 HGR0000093 HGR0000101 HGR0000149 
## Target age: 77 sex: M 
## Matched subject: HGR0000015 HGR0000088 HGR0000089
recovered.patients <- unique(recovered.patients[!is.na(recovered.patients)])
recovered.patients.cytokines <- subset(covid19.lab.results,test_type == "ELISA" & subject_id %in% recovered.patients)
cat("Map of samping TSO time points for recovered patients:\n")
## Map of samping TSO time points for recovered patients:
recovered.time.points <- as.data.frame.matrix(table(recovered.patients.cytokines[,c("subject_id","days_from_symptom_onset_to_test")])) 
colSums(recovered.time.points > 0)
##  3  8  9 10 11 12 13 14 15 16 17 18 19 20 22 23 24 25 26 27 29 30 31 32 33 34 
##  1  1  1  3  2  2  3  3  4  4  2  5  5  6  1  5  2  3  2  3  2  3  2  1  2  1 
## 35 37 38 39 40 44 45 48 51 56 
##  2  1  1  1  1  2  1  1  1  1
cat("Age and sex distribution:\n")
## Age and sex distribution:
table(subset(covid19.patient.demo,subject_id %in% recovered.patients)$sex)
## 
##  F  M 
##  3 18
summary(subset(covid19.patient.demo,subject_id %in% recovered.patients)$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   45.00   49.00   59.00   60.33   69.00   78.00
selected.cytokine.profiles <- rbind(deceased.patients.cytokines,recovered.patients.cytokines)
# remove patients in citeseq
#citeseq.patients <- subset(covid19.samples,!is.na(batch))$subject_id
#intersect(citeseq.patients,c(recovered.patients,deceased.patients))
#selected.cytokine.profiles <- subset(selected.cytokine.profiles,!(subject_id %in% citeseq.patients)) 
parameter.labels <- unique(selected.cytokine.profiles[,c("test_id","test_name")])
rownames(parameter.labels) <- parameter.labels$test_id
parameter.labels$test_name <- gsub("/ ","/",parameter.labels$test_name,fixed = T)

Grouping of samples

Separate samples into interval groups relative to juncture: early (7-16 days), mid (juncture/17-23 days), and late (24-30 days)

selected.cytokine.profiles$unique.id <- make.unique(apply(selected.cytokine.profiles[,c("subject_id","days_from_admission_to_test","test_id")],1,
                                                          paste0,collapse="."),sep = "|")
selected.cytokine.profiles$unique.id <- sapply(selected.cytokine.profiles$unique.id,function(x){unlist(strsplit(x,"|",fixed=T))[2]})
selected.cytokine.profiles.matrix <- reshape2::dcast(selected.cytokine.profiles,
                                                     subject_id + age + sex + outcome + severity + days_from_symptom_onset_to_test + unique.id ~ test_id,
                                                     value.var = "test_value")
selected.cytokine.profiles.matrix[,colnames(selected.cytokine.profiles.matrix) %in% selected.cytokine.profiles$test_id] <-
  apply(selected.cytokine.profiles.matrix[,colnames(selected.cytokine.profiles.matrix) %in% selected.cytokine.profiles$test_id],2,function(x){scale(log10(x+1))})
selected.cytokine.profiles.matrix$stage <- NA
tmp <- which(selected.cytokine.profiles.matrix$days_from_symptom_onset_to_test >= 7 & selected.cytokine.profiles.matrix$days_from_symptom_onset_to_test <= 16)
selected.cytokine.profiles.matrix[tmp,"stage"] <- "early"
tmp <- which(selected.cytokine.profiles.matrix$days_from_symptom_onset_to_test >= 17 & selected.cytokine.profiles.matrix$days_from_symptom_onset_to_test <= 23)
selected.cytokine.profiles.matrix[tmp,"stage"] <- "mid"
tmp <- which(selected.cytokine.profiles.matrix$days_from_symptom_onset_to_test >= 24 & selected.cytokine.profiles.matrix$days_from_symptom_onset_to_test <= 30)
selected.cytokine.profiles.matrix[tmp,"stage"] <- "late"

cat("Sample distribution:\n")
## Sample distribution:
table(selected.cytokine.profiles.matrix[,c("outcome","stage")])
##           stage
## outcome    early late mid
##   Alive       24   16  25
##   Deceased    11   12  13


Time course of selected cytokines

display.cytokine.groups <- unique(unlist(bin1.cytokines))

display.cytokine.profiles <- subset(selected.cytokine.profiles,test_id %in% display.cytokine.groups)
#display.cytokine.profiles$test_name <- factor(display.cytokine.profiles$test_name,levels=unlist(display.cytokine.groups)) 

display.cytokine.profiles.avg <- aggregate(test_value ~ subject_id + days_from_symptom_onset_to_test + test_name + test_id + outcome,
                                            display.cytokine.profiles,mean)
subjects_with_early_tp <- unique(subset(display.cytokine.profiles.avg,days_from_symptom_onset_to_test <= 23)$subject_id)
ggplot(display.cytokine.profiles.avg,aes(days_from_symptom_onset_to_test,log10(test_value+1),group=outcome)) + 
  #geom_vline(xintercept = 20,color="grey60",linetype="dashed") +
  geom_rect(aes(xmin=17,xmax=23,ymax=Inf,ymin=-Inf),alpha=0.05,fill = "#EADCFA") +
  geom_point(alpha=.25,pch=21,aes(fill=outcome),size=.5) + 
  geom_smooth(se = T,aes(color=outcome,fill=outcome),alpha=0.15,size=.8) + #,data = subset(display.cytokine.profiles.avg,subject_id %in% subjects_with_early_tp)) + 
  geom_line(aes(group=subject_id),alpha=0.1) +
  scale_color_manual(name="Outcome",values=c("#374E55FF","#DF8F44FF")) + scale_fill_manual(name="Outcome",values=c("#374E55FF","#DF8F44FF")) +
  facet_wrap(~test_name,scales = "free",ncol=5,dir = "v") + theme_bw() + xlim(0,60) + 
  xlab("Days from Symptom Onset") + ylab("Cytokine Level (log10 pg/mL)") + 
  theme(strip.background = element_rect(fill="white"),strip.text = element_text(size=7))
<b>Fig</b> Time course of selected circulating protein levels of critical ill patients with recovery or deceased outcomes. Longitudinal samples from the same individual are connected by gray lines. Trajectories were fitted to the two age- and gender-matched recovered versus deceased patient groups separately and are shown with the shaded areas representing standard error.

Fig Time course of selected circulating protein levels of critical ill patients with recovery or deceased outcomes. Longitudinal samples from the same individual are connected by gray lines. Trajectories were fitted to the two age- and gender-matched recovered versus deceased patient groups separately and are shown with the shaded areas representing standard error.


Serology time course differences

SARS-Covid2 Spike and nucleocapsid antibodies

selected.sero.profiles <- subset(covid19.lab.results,test_type == "Serology" & subject_id %in% selected.cytokine.profiles$subject_id)
selected.sero.profiles$test_value <- log10(selected.sero.profiles$test_value+1)

ggplot(selected.sero.profiles,aes(days_from_symptom_onset_to_test,test_value)) +
  geom_vline(xintercept = 20,color="grey60",linetype="dashed") +
  geom_point(alpha=.25,pch=21,aes(fill=outcome),size=1) + 
  facet_wrap(~test_name,scales = "free",ncol = 2) + theme_bw() + 
  geom_smooth(se = T,aes(color=outcome,fill=outcome),alpha=0.15,size=1.25) + 
  geom_line(aes(group=subject_id),alpha=0.1) +
  scale_fill_manual(name="Outcome",values=c("#374E55FF","#DF8F44FF")) + 
  scale_color_manual(name="Outcome",values=c("#374E55FF","#DF8F44FF")) + 
  ylab("Antibody Level (log10 LU)") + xlab("Days from Symptom Onset") 
<b>Fig</b> Antibody time course for deceased and recovered patients.

Fig Antibody time course for deceased and recovered patients.

CBC cell count time course

selected.cbc.profiles <- subset(covid19.lab.results,test_type == "Clinical" & subject_id %in% c(recovered.patients,deceased.patients))
selected.cbc.profiles <- selected.cbc.profiles[grep("((T|C|s|\\.)_x10.[36].uL$)|[sCT]_pct$",selected.cbc.profiles$test_id),]
selected.cbc.profiles <- subset(selected.cbc.profiles,test_id %in% names(which(table(selected.cbc.profiles$test_id) > 100)))

# reference range
cell.counts.range <- unique(selected.cbc.profiles[,c("test_name","test_id","normal_range")])
cell.counts.range[,c("lower.bound","upper.bound")] <- t(sapply(cell.counts.range$normal_range,function(x){as.numeric(unlist(strsplit(x,"-")))}))
cell.counts.range <- cbind(aggregate(lower.bound ~ test_id + test_name,cell.counts.range,min),
                           upper.bound=aggregate(upper.bound ~ test_id + test_name,cell.counts.range,max)[,"upper.bound"])

# remove possible erroneous values
tolerance <- 10
to.remove.rows <- c()
for (i in unique(selected.cbc.profiles$test_id)) {
  cat(i,"\n")
  test.values <- subset(selected.cbc.profiles,test_id == i)$test_value
  median.value <- median(test.values,na.rm = T)
  range <- quantile(test.values,.75,na.rm = T) - quantile(test.values,.25,na.rm = T)
  to.remove.rows <- c(to.remove.rows,rownames(subset(selected.cbc.profiles,test_id == i & 
                                                       (test_value > (median.value + tolerance*range) | 
                                                       test_value < (median.value - tolerance*range)))))
}
## White.Blood.Cells.WBC_x10.3.uL 
## Hematocrit.HCT_pct 
## Basophils_pct 
## Monocytes_x10.3.uL 
## Eosinophils_x10.3.uL 
## Red.Blood.Cells.RBC_x10.6.uL 
## Neutrophils_pct 
## Eosinophils_pct 
## Monocytes_pct 
## Basophils_x10.3.uL 
## Lymphocytes_x10.3.uL 
## Lymphocytes_pct 
## Platelets.PLT_x10.3.uL 
## Neutrophils_x10.3.uL
#View(selected.cbc.profiles[to.remove.rows,])
selected.cbc.profiles <- subset(selected.cbc.profiles,!(rownames(selected.cbc.profiles) %in% to.remove.rows))

selected.cbc.profiles.avg <- aggregate(test_value ~ subject_id + days_from_symptom_onset_to_test + test_name + test_id + severity.outcome,
                                            selected.cbc.profiles,mean)
to.display <- selected.cbc.profiles.avg[grep("_pct$",selected.cbc.profiles.avg$test_id),]
ggplot(to.display) + 
  geom_hline(data=subset(cell.counts.range,test_id %in% to.display$test_id),aes(yintercept=lower.bound),alpha=1,color="darkseagreen3",linetype="dashed") +   
  geom_hline(data=subset(cell.counts.range,test_id %in% to.display$test_id),aes(yintercept=upper.bound),alpha=1,color="darkseagreen3",linetype="dashed") +
  geom_rect(aes(xmin=17,xmax=23,ymax=Inf,ymin=-Inf),alpha=0.05,fill = "#EADCFA") +
  #geom_vline(xintercept = 20,color="grey60",linetype="dashed") +
  geom_line(aes(x=days_from_symptom_onset_to_test,y=test_value,group=subject_id),alpha=0.1) + 
  geom_point(alpha=.1,pch=21,aes(fill=severity.outcome,x=days_from_symptom_onset_to_test,y=test_value)) +
  geom_smooth(se = T,aes(x=days_from_symptom_onset_to_test,y=test_value,group=severity.outcome,color=severity.outcome),alpha=0.3) + xlim(0,40) +
  scale_color_manual(name="Severity-Outcome",values=severity.color) + scale_fill_manual(name="Severity-Outcome",values=severity.color) +
    facet_wrap(~test_name,scales = "free",ncol=5,dir = "v") + theme_bw() + ylab("Cell Counts (% of Parent Pop.)") + xlab("Days from Symptom Onset") +
    theme(strip.background = element_rect(fill="white"),strip.text = element_text(size=7))
<b>Fig</b> Cell populations as % parent for deceased and recovered patients.

Fig Cell populations as % parent for deceased and recovered patients.

to.display <- selected.cbc.profiles.avg[grep("_pct$",selected.cbc.profiles.avg$test_id,invert = T),]
ggplot(to.display) + 
 geom_hline(data=subset(cell.counts.range,test_id %in% to.display$test_id),aes(yintercept=lower.bound),alpha=1,color="darkseagreen3",linetype="dashed") +   
  geom_hline(data=subset(cell.counts.range,test_id %in% to.display$test_id),aes(yintercept=upper.bound),alpha=1,color="darkseagreen3",linetype="dashed") +
  geom_rect(aes(xmin=17,xmax=23,ymax=Inf,ymin=-Inf),alpha=0.05,fill = "#EADCFA") +
  #geom_vline(xintercept = 20,color="grey60",linetype="dashed") +
  geom_line(aes(x=days_from_symptom_onset_to_test,y=test_value,group=subject_id),alpha=0.1) + 
  geom_point(alpha=.1,pch=21,aes(fill=severity.outcome,x=days_from_symptom_onset_to_test,y=test_value)) +
  geom_smooth(se = T,aes(x=days_from_symptom_onset_to_test,y=test_value,group=severity.outcome,color=severity.outcome),alpha=0.3) + xlim(0,40) +
  scale_color_manual(name="Severity-Outcome",values=severity.color) + scale_fill_manual(name="Severity-Outcome",values=severity.color) +
    facet_wrap(~test_id,scales = "free",ncol=5,dir = "v") + theme_bw() + ylab("Cell Counts") + xlab("Days from Symptom Onset") + 
    theme(strip.background = element_rect(fill="white"),strip.text = element_text(size=7))
<b>Fig</b> Absolute cell counts for deceased and recovered patients.

Fig Absolute cell counts for deceased and recovered patients.

Session Info

# save all results
save(merged.de.df,mid.did.de.df,late.did.de.df,file = file.path(output.folder,"deceased.v.recovered.cytokine.test.results.RData"))

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] precrec_0.11.2   plsRglm_1.2.5    lmerTest_3.1-2   lme4_1.1-23     
##  [5] Matrix_1.2-18    tidymv_3.0.0     mgcv_1.8-31      nlme_3.1-144    
##  [9] ggfortify_0.4.10 ggpubr_0.4.0     ggsci_2.9        ggplot2_3.3.2   
## [13] knitr_1.30      
## 
## loaded via a namespace (and not attached):
##  [1] numDeriv_2016.8-1.1  tools_3.6.3          backports_1.1.10    
##  [4] R6_2.4.1             vegan_2.5-6          colorspace_1.4-1    
##  [7] permute_0.9-5        withr_2.3.0          tidyselect_1.1.0    
## [10] gridExtra_2.3        curl_4.3             compiler_3.6.3      
## [13] network_1.16.0       labeling_0.3         scales_1.1.1        
## [16] mvtnorm_1.1-1        stringr_1.4.0        digest_0.6.25       
## [19] foreign_0.8-75       minqa_1.2.4          rmarkdown_2.4       
## [22] rio_0.5.16           pkgconfig_2.0.3      htmltools_0.5.0     
## [25] highr_0.8            maps_3.3.0           rlang_0.4.7         
## [28] readxl_1.3.1         farver_2.0.3         generics_0.0.2      
## [31] statnet.common_4.3.0 dplyr_1.0.2          zip_2.1.1           
## [34] car_3.0-10           magrittr_1.5         dotCall64_1.0-0     
## [37] bipartite_2.15       Rcpp_1.0.5           munsell_0.5.0       
## [40] abind_1.4-5          lifecycle_0.2.0      stringi_1.4.6       
## [43] yaml_2.2.1           carData_3.0-4        MASS_7.3-51.5       
## [46] plyr_1.8.6           grid_3.6.3           parallel_3.6.3      
## [49] ggrepel_0.8.2        forcats_0.5.0        crayon_1.3.4        
## [52] lattice_0.20-38      haven_2.3.1          splines_3.6.3       
## [55] hms_0.5.3            sna_2.5              pillar_1.4.6        
## [58] igraph_1.2.5         boot_1.3-24          ggsignif_0.6.0      
## [61] reshape2_1.4.4       glue_1.4.2           evaluate_0.14       
## [64] data.table_1.13.0    vctrs_0.3.4          spam_2.5-1          
## [67] nloptr_1.2.2.2       cellranger_1.1.0     gtable_0.3.0        
## [70] purrr_0.3.4          tidyr_1.1.2          xfun_0.17           
## [73] openxlsx_4.2.2       broom_0.7.0          coda_0.19-3         
## [76] rstatix_0.6.0        tibble_3.0.3         fields_11.4         
## [79] cluster_2.1.0        statmod_1.4.34       ellipsis_0.3.1