Simple Computational Models for Predicting Pertussis Boost Response

Setup: Load libraries and helper functions

suppressPackageStartupMessages({
  library(verification)
  library(pROC)
  library(discretefit)
  library(dplyr)
  library(ggplot2)
  library(tidyr)
  library(tidyverse)
  library(corrplot)
  library(data.table)
  library(readr)
  library(kableExtra)
  library(glmnet)
  library(here)
  library(kableExtra)
  library(formattable)
  library(colorRamp2)
  library(ComplexHeatmap)
  library(GetoptLong)
  library(ellipse)
})

Setup: Variables used through out the script

Tasks:

1. Monocytes at day 1

2. CCL3 at day 3

3. IgG-PT at day 14

DATASET <- c("2020_dataset", "2021_dataset")
TIMEPOINTS <- c(0, 1, 3, 14)
dataFile <- "combined_dataset2020_2021.csv"

META_COLS <- c("specimen_id", "subject_id", "timepoint", "dataset", 
               "biological_sex", "infancy_vac", "age_at_boost")
ABTITER_COLS <- c("IgG_PT")
RNA_COLS <- c("CCL3")
CELL_COLS <- c("Monocytes")

DEMOGRAPHY_COLS <- c("age_at_boost", "biological_sex", "infancy_vac")

TASK_COLS <- c("Monocytes_D1", "CCL3_D3", "IgG_PT_D14")
TASKS_BASELINES <- c("Monocytes_D0", "CCL3_D0", "IgG_PT_D0")
BASE_COLS <- c("Monocytes_D0", "CD4Tcells_D0", "IgG_PT_D0", 'IgG_FHA_D0', 
               'IgG_PRN_D0','IgG1_PT_D0', 'IgG1_FHA_D0', 'IgG4_PT_D0', 
               'IgG4_FHA_D0', "CCL3_D0", "IL6_D0", "NFKBIA_D0")

LOG_TRANS_COLS <- c("CCL3_D0", "IgG_PT_D0", "CCL3_D3",  "IgG_PT_D14")

Setup: Load dataset, combine different experiments

knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
df_source <- readRDS(here("./data/master_normalized_data_challenge2_train.RDS"))

metaDf <- data.frame(df_source[["subject_specimen"]])
metaDf["age_at_boost"] <- as.numeric(round(difftime(metaDf$date_of_boost, metaDf$year_of_birth,units="weeks")/52, 2))
metaDf <- metaDf[, META_COLS]

abtiterDf <- data.frame(df_source[["plasma_antibody_levels_wide"]])
abtiterDf["specimen_id"] <- as.numeric(row.names(abtiterDf))
abtiterDf <- data.frame(abtiterDf[, c("specimen_id", ABTITER_COLS)])

rnaDf <- data.frame(df_source[["pbmc_gene_expression_wide"]])
rnaDf["specimen_id"] <- as.numeric(row.names(rnaDf))
tasks_seq <- c('ENSG00000277632')
for (i in 1:length(tasks_seq)){
  rnaDf <- data.frame(rnaDf %>% rename_at(vars(starts_with(tasks_seq[i])), ~RNA_COLS[i]))
}
rnaDf <- data.frame(rnaDf[, c("specimen_id", RNA_COLS)])

cellDf <- data.frame(df_source[["pbmc_cell_frequency_wide"]])
cellDf["specimen_id"] <- as.numeric(row.names(cellDf))
cellDf <- data.frame(cellDf[, c("specimen_id", CELL_COLS)])

list_df <- list(metaDf, cellDf, abtiterDf, rnaDf)
df_merge <- list_df %>% reduce(full_join, by="specimen_id")
df_merge <- df_merge[df_merge$timepoint %in% TIMEPOINTS, ]

df_pivot <- df_merge[, names(df_merge)!="specimen_id"] %>% 
  pivot_wider(id_cols=c("subject_id", "dataset", "biological_sex", 
                        "infancy_vac", "age_at_boost"), 
              names_from = timepoint, 
              values_from = all_of(c(CELL_COLS, RNA_COLS, ABTITER_COLS)), 
              names_sep = "_D")

df_pivot <- df_pivot[df_pivot$dataset %in% DATASET, ]

Setup: Dataset log transform , Ranking and fold change

Fold change for actual values been generated by dividing the actual value of each task to it’s baseline
Fold change for normalized values been generated by dividing the log transformed value of each task to it’s log transformed baseline

knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
targetX <- c("Monocytes_D0", "CCL3_D0", "IgG_PT_D0")
targetY <- c("Monocytes_D1","CCL3_D3", "IgG_PT_D14")

fc_cols <- paste(targetY, "FC", sep="_")

ranked_cols <- paste(c("age_at_boost", targetX, targetY, fc_cols), "Rank", sep="_")

df <- df_pivot[, c("subject_id", DEMOGRAPHY_COLS, targetX, targetY)]

rankingFunction <- function(x) {
  as.numeric(rank(-x, ties.method = "average", na.last = "keep"))
}

targetX <- c("Monocytes_D0", "CCL3_D0", "IgG_PT_D0")
targetY <- c("Monocytes_D1","CCL3_D3", "IgG_PT_D14")

df[,"Monocytes_D1_FC"] <- df[, "Monocytes_D1"] / df[, "Monocytes_D0"]
df[,"CCL3_D3_FC"] <- df[, "CCL3_D3"] / df[, "CCL3_D0"]
df[,"IgG_PT_D14_FC"] <- df[, "IgG_PT_D14"] / df[, "IgG_PT_D0"]

df[, ranked_cols] <- apply(df[, c("age_at_boost", targetX, targetY, fc_cols)],
                           2, rankingFunction)
df <- data.frame(df)

Model construction

No actual model has been developed. We use baseline values of tasks provided and demographic features such as age of participants as model. Two types of models have been developed:

A. Age of the participants

B. Baseline values of task variables

A. Modeling Approach: Age of the participants

Modeling steps:

1.  Extract age of the participants
2.  Convert actual values into ranks
3.  Prepare submission file

The submission file looks as follows:

Possible rank of predicted values for each task when model is trained based on age at boost.

expCols <- c("1.1) IgG-PT-D14-Rank", "1.2) IgG-PT-D14-FC-Rank", 
             "2.1) Monocytes-D1-Rank", "2.2) Monocytes-D1-FC-Rank", 
             "3.1) CCL3-D3-Rank", "3.2) CCL3-D3-FC-Rank")
ageModel_df <- df[, c("subject_id", DEMOGRAPHY_COLS, 
                      "IgG_PT_D14_Rank", "IgG_PT_D14_FC_Rank", 
                      "Monocytes_D1_Rank", "Monocytes_D1_FC_Rank", 
                      "CCL3_D3_Rank", "CCL3_D3_FC_Rank")]

ageModel_df$Monocytes_D1_Rank <- df$age_at_boost_Rank
ageModel_df$CCL3_D3_Rank <- df$age_at_boost_Rank
ageModel_df$IgG_PT_D14_Rank <- df$age_at_boost_Rank

ageModel_df$Monocytes_D1_FC_Rank <- df$age_at_boost_Rank
ageModel_df$CCL3_D3_FC_Rank <- df$age_at_boost_Rank
ageModel_df$IgG_PT_D14_FC_Rank <- df$age_at_boost_Rank

colnames(ageModel_df)[colnames(ageModel_df)=="subject_id"] <- "Subject ID"
colnames(ageModel_df)[colnames(ageModel_df)=="age_at_boost"] <- "Age" 
colnames(ageModel_df)[colnames(ageModel_df)=="infancy_vac"] <- "Vaccine Priming Status" 
colnames(ageModel_df)[colnames(ageModel_df)=="biological_sex"] <- "Biological Sex at Birth" 
colnames(ageModel_df)[colnames(ageModel_df)=="Monocytes_D1_Rank"] <- "2.1) Monocytes_D1_Rank"
colnames(ageModel_df)[colnames(ageModel_df)=="Monocytes_D1_FC_Rank"] <- "2.2) Monocytes_D1_FC_Rank" 
colnames(ageModel_df)[colnames(ageModel_df)=="CCL3_D3_Rank"] <- "3.1) CCL3_D3_Rank" 
colnames(ageModel_df)[colnames(ageModel_df)=="CCL3_D3_FC_Rank"] <- "3.2) CCL3_D3_FC_Rank" 
colnames(ageModel_df)[colnames(ageModel_df)=="IgG_PT_D14_Rank"] <- "1.1) IgG_PT_D14_titer_Rank" 
colnames(ageModel_df)[colnames(ageModel_df)=="IgG_PT_D14_FC_Rank"] <- "1.2) IgG_PT_D14_FC_Rank"  

knitr::kable(ageModel_df, "html", align = "lccrr", booktabs=TRUE, border_left = T, 
             border_right = T, caption = "Age Based Models")  %>% 
  kable_styling("striped", full_width = T) %>% 
  scroll_box(width = "100%", height = "400px")

Age Based Models
Subject ID	Age	Biological Sex at Birth	Vaccine Priming Status	1.1) IgG_PT_D14_titer_Rank	1.2) IgG_PT_D14_FC_Rank	2.1) Monocytes_D1_Rank	2.2) Monocytes_D1_FC_Rank	3.1) CCL3_D3_Rank	3.2) CCL3_D3_FC_Rank
1	30.80	Female	wP	22.0	22.0	22.0	22.0	22.0	22.0
2	51.25	Female	wP	1.0	1.0	1.0	1.0	1.0	1.0
3	33.89	Female	wP	13.0	13.0	13.0	13.0	13.0	13.0
4	28.76	Male	wP	30.0	30.0	30.0	30.0	30.0	30.0
5	25.75	Male	wP	42.0	42.0	42.0	42.0	42.0	42.0
6	28.87	Female	wP	27.0	27.0	27.0	27.0	27.0	27.0
7	35.97	Female	wP	6.0	6.0	6.0	6.0	6.0	6.0
8	34.27	Female	wP	11.0	11.0	11.0	11.0	11.0	11.0
9	20.63	Male	aP	63.0	63.0	63.0	63.0	63.0	63.0
10	34.68	Female	wP	9.5	9.5	9.5	9.5	9.5	9.5
11	30.76	Female	wP	23.0	23.0	23.0	23.0	23.0	23.0
12	34.68	Male	wP	9.5	9.5	9.5	9.5	9.5	9.5
13	19.63	Male	aP	84.0	84.0	84.0	84.0	84.0	84.0
14	23.70	Male	wP	47.0	47.0	47.0	47.0	47.0	47.0
15	27.71	Male	wP	33.0	33.0	33.0	33.0	33.0	33.0
16	29.66	Female	wP	24.0	24.0	24.0	24.0	24.0	24.0
17	36.82	Female	wP	5.0	5.0	5.0	5.0	5.0	5.0
18	19.73	Female	aP	83.0	83.0	83.0	83.0	83.0	83.0
19	22.81	Male	wP	50.0	50.0	50.0	50.0	50.0	50.0
20	35.78	Female	wP	7.0	7.0	7.0	7.0	7.0	7.0
21	33.77	Male	wP	15.0	15.0	15.0	15.0	15.0	15.0
22	31.77	Female	wP	20.0	20.0	20.0	20.0	20.0	20.0
23	25.82	Female	wP	41.0	41.0	41.0	41.0	41.0	41.0
24	24.79	Female	wP	43.0	43.0	43.0	43.0	43.0	43.0
25	28.80	Female	wP	29.0	29.0	29.0	29.0	29.0	29.0
26	33.85	Female	wP	14.0	14.0	14.0	14.0	14.0	14.0
27	19.80	Female	aP	81.5	81.5	81.5	81.5	81.5	81.5
28	34.85	Male	wP	8.0	8.0	8.0	8.0	8.0	8.0
29	19.80	Male	aP	81.5	81.5	81.5	81.5	81.5	81.5
30	28.84	Female	wP	28.0	28.0	28.0	28.0	28.0	28.0
31	27.83	Female	wP	32.0	32.0	32.0	32.0	32.0	32.0
32	19.88	Male	aP	79.0	79.0	79.0	79.0	79.0	79.0
33	26.87	Male	wP	35.0	35.0	35.0	35.0	35.0	35.0
34	33.93	Female	wP	12.0	12.0	12.0	12.0	12.0	12.0
35	25.86	Male	wP	40.0	40.0	40.0	40.0	40.0	40.0
36	19.88	Female	aP	79.0	79.0	79.0	79.0	79.0	79.0
37	18.91	Female	aP	94.5	94.5	94.5	94.5	94.5	94.5
38	19.88	Female	aP	79.0	79.0	79.0	79.0	79.0	79.0
39	31.92	Female	wP	19.0	19.0	19.0	19.0	19.0	19.0
40	22.89	Female	wP	49.0	49.0	49.0	49.0	49.0	49.0
41	31.96	Male	wP	18.0	18.0	18.0	18.0	18.0	18.0
42	19.92	Female	aP	77.0	77.0	77.0	77.0	77.0	77.0
43	18.91	Female	aP	94.5	94.5	94.5	94.5	94.5	94.5
44	18.91	Female	aP	94.5	94.5	94.5	94.5	94.5	94.5
45	19.98	Female	aP	75.0	75.0	75.0	75.0	75.0	75.0
46	18.91	Female	aP	94.5	94.5	94.5	94.5	94.5	94.5
47	20.98	Female	aP	62.0	62.0	62.0	62.0	62.0	62.0
48	19.11	Female	aP	90.0	90.0	90.0	90.0	90.0	90.0
49	20.11	Female	aP	72.0	72.0	72.0	72.0	72.0	72.0
50	19.98	Female	aP	75.0	75.0	75.0	75.0	75.0	75.0
51	19.98	Male	aP	75.0	75.0	75.0	75.0	75.0	75.0
52	19.07	Male	aP	91.5	91.5	91.5	91.5	91.5	91.5
53	19.07	Female	aP	91.5	91.5	91.5	91.5	91.5	91.5
54	20.11	Female	aP	72.0	72.0	72.0	72.0	72.0	72.0
55	20.11	Female	aP	72.0	72.0	72.0	72.0	72.0	72.0
56	20.15	Female	aP	68.5	68.5	68.5	68.5	68.5	68.5
57	21.15	Female	aP	60.0	60.0	60.0	60.0	60.0	60.0
58	20.15	Female	aP	68.5	68.5	68.5	68.5	68.5	68.5
59	20.15	Female	aP	68.5	68.5	68.5	68.5	68.5	68.5
60	20.15	Male	aP	68.5	68.5	68.5	68.5	68.5	68.5
61	32.38	Female	wP	16.0	16.0	16.0	16.0	16.0	16.0
62	25.99	Female	wP	38.5	38.5	38.5	38.5	38.5	38.5
63	23.98	Female	wP	46.0	46.0	46.0	46.0	46.0	46.0
64	25.99	Male	wP	38.5	38.5	38.5	38.5	38.5	38.5
65	29.02	Male	wP	26.0	26.0	26.0	26.0	26.0	26.0
66	43.07	Female	wP	4.0	4.0	4.0	4.0	4.0	4.0
67	47.24	Female	wP	2.5	2.5	2.5	2.5	2.5	2.5
68	47.24	Male	wP	2.5	2.5	2.5	2.5	2.5	2.5
69	29.17	Female	wP	25.0	25.0	25.0	25.0	25.0	25.0
70	21.15	Male	aP	60.0	60.0	60.0	60.0	60.0	60.0
71	21.15	Female	aP	60.0	60.0	60.0	60.0	60.0	60.0
72	28.25	Female	wP	31.0	31.0	31.0	31.0	31.0	31.0
73	24.23	Female	wP	44.5	44.5	44.5	44.5	44.5	44.5
74	24.23	Female	wP	44.5	44.5	44.5	44.5	44.5	44.5
75	21.22	Female	aP	57.5	57.5	57.5	57.5	57.5	57.5
76	21.22	Female	aP	57.5	57.5	57.5	57.5	57.5	57.5
77	31.32	Male	wP	21.0	21.0	21.0	21.0	21.0	21.0
78	26.30	Female	wP	36.5	36.5	36.5	36.5	36.5	36.5
79	32.32	Male	wP	17.0	17.0	17.0	17.0	17.0	17.0
80	27.30	Female	wP	34.0	34.0	34.0	34.0	34.0	34.0
81	26.30	Male	wP	36.5	36.5	36.5	36.5	36.5	36.5
82	21.28	Female	aP	56.0	56.0	56.0	56.0	56.0	56.0
83	20.34	Female	aP	66.0	66.0	66.0	66.0	66.0	66.0
84	22.34	Female	aP	52.0	52.0	52.0	52.0	52.0	52.0
85	19.39	Female	aP	88.0	88.0	88.0	88.0	88.0	88.0
86	21.40	Female	aP	55.0	55.0	55.0	55.0	55.0	55.0
87	19.39	Male	aP	88.0	88.0	88.0	88.0	88.0	88.0
88	19.39	Male	aP	88.0	88.0	88.0	88.0	88.0	88.0
89	22.49	Female	aP	51.0	51.0	51.0	51.0	51.0	51.0
90	20.49	Female	aP	65.0	65.0	65.0	65.0	65.0	65.0
91	21.49	Male	aP	54.0	54.0	54.0	54.0	54.0	54.0
92	19.54	Female	aP	85.5	85.5	85.5	85.5	85.5	85.5
93	23.56	Female	aP	48.0	48.0	48.0	48.0	48.0	48.0
94	20.55	Male	aP	64.0	64.0	64.0	64.0	64.0	64.0
95	21.55	Female	aP	53.0	53.0	53.0	53.0	53.0	53.0
96	19.54	Male	aP	85.5	85.5	85.5	85.5	85.5	85.5

B. Modeling approach: Baseline values of task variables

Modeling steps:

1.  Extract baseline of each task
2.  Convert actual values into ranks
3.  Prepare submission file

The submission file looks as follow:

Possible rank of predicted values for each task when model is trained based on the baseline of that task.

expCols <- c("1.1) IgG-PT-D14-Rank", "1.2) IgG-PT-D14-FC-Rank",  
             "2.1) Monocytes-D1-Rank", "2.2) Monocytes-D1-FC-Rank", 
             "3.1) CCL3-D3-Rank", "3.2) CCL3-D3-FC-Rank")
baselineModel_df <- df[, c("subject_id", DEMOGRAPHY_COLS, 
                           "IgG_PT_D14_Rank", "IgG_PT_D14_FC_Rank", 
                            "Monocytes_D1_Rank", "Monocytes_D1_FC_Rank", 
                            "CCL3_D3_Rank", "CCL3_D3_FC_Rank")]

baselineModel_df$Monocytes_D1_Rank <- df$Monocytes_D1_Rank
baselineModel_df$CCL3_D3_Rank <- df$CCL3_D3_Rank
baselineModel_df$IgG_PT_D14_Rank <- df$IgG_PT_D14_Rank

baselineModel_df$Monocytes_D1_FC_Rank <- df$Monocytes_D1_Rank
baselineModel_df$CCL3_D3_FC_Rank <- df$CCL3_D3_Rank
baselineModel_df$IgG_PT_D14_FC_Rank <- df$IgG_PT_D14_Rank

colnames(baselineModel_df)[colnames(baselineModel_df)=="subject_id"] <- "Subject ID"
colnames(baselineModel_df)[colnames(baselineModel_df)=="age_at_boost"] <- "Age" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="infancy_vac"] <- "Vaccine Priming Status" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="biological_sex"] <- "Biological Sex at Birth" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="Monocytes_D1_Rank"] <- "2.1) Monocytes_D1_Rank"
colnames(baselineModel_df)[colnames(baselineModel_df)=="Monocytes_D1_FC_Rank"] <- "2.2) Monocytes_D1_FC_Rank" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="CCL3_D3_Rank"] <- "3.1) CCL3_D3_Rank" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="CCL3_D3_FC_Rank"] <- "3.2) CCL3_D3_FC_Rank" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="IgG_PT_D14_Rank"] <- "1.1) IgG_PT_D14_titer_Rank" 
colnames(baselineModel_df)[colnames(baselineModel_df)=="IgG_PT_D14_FC_Rank"] <- "1.2) IgG_PT_D14_FC_Rank"    

knitr::kable(baselineModel_df, "html", align = "lccrr", booktabs=TRUE, border_left = T, 
             border_right = T, caption = "Baseline Based Models")  %>% 
  kable_styling("striped", full_width = T) %>% 
  scroll_box(width = "100%", height = "400px")

Baseline Based Models
Subject ID	Age	Biological Sex at Birth	Vaccine Priming Status	1.1) IgG_PT_D14_titer_Rank	1.2) IgG_PT_D14_FC_Rank	2.1) Monocytes_D1_Rank	2.2) Monocytes_D1_FC_Rank	3.1) CCL3_D3_Rank	3.2) CCL3_D3_FC_Rank
1	30.80	Female	wP	20	20	NA	NA	21.0	21.0
2	51.25	Female	wP	NA	NA	NA	NA	NA	NA
3	33.89	Female	wP	40	40	NA	NA	36.0	36.0
4	28.76	Male	wP	35	35	53.0	53.0	64.0	64.0
5	25.75	Male	wP	48	48	NA	NA	42.0	42.0
6	28.87	Female	wP	29	29	1.0	1.0	45.0	45.0
7	35.97	Female	wP	50	50	NA	NA	NA	NA
8	34.27	Female	wP	NA	NA	NA	NA	NA	NA
9	20.63	Male	aP	87	87	NA	NA	25.0	25.0
10	34.68	Female	wP	89	89	NA	NA	48.5	48.5
11	30.76	Female	wP	1	1	52.0	52.0	50.0	50.0
12	34.68	Male	wP	49	49	NA	NA	NA	NA
13	19.63	Male	aP	68	68	NA	NA	14.0	14.0
14	23.70	Male	wP	11	11	NA	NA	NA	NA
15	27.71	Male	wP	46	46	46.0	46.0	30.0	30.0
16	29.66	Female	wP	53	53	NA	NA	NA	NA
17	36.82	Female	wP	28	28	29.0	29.0	26.0	26.0
18	19.73	Female	aP	23	23	NA	NA	67.0	67.0
19	22.81	Male	wP	27	27	NA	NA	64.0	64.0
20	35.78	Female	wP	38	38	10.0	10.0	23.0	23.0
21	33.77	Male	wP	77	77	5.0	5.0	37.5	37.5
22	31.77	Female	wP	52	52	NA	NA	69.0	69.0
23	25.82	Female	wP	39	39	NA	NA	57.5	57.5
24	24.79	Female	wP	82	82	NA	NA	55.0	55.0
25	28.80	Female	wP	83	83	NA	NA	45.0	45.0
26	33.85	Female	wP	71	71	30.0	30.0	64.0	64.0
27	19.80	Female	aP	72	72	NA	NA	45.0	45.0
28	34.85	Male	wP	43	43	NA	NA	NA	NA
29	19.80	Male	aP	60	60	13.0	13.0	33.0	33.0
30	28.84	Female	wP	66	66	NA	NA	NA	NA
31	27.83	Female	wP	78	78	51.0	51.0	19.0	19.0
32	19.88	Male	aP	10	10	NA	NA	18.0	18.0
33	26.87	Male	wP	88	88	25.0	25.0	16.0	16.0
34	33.93	Female	wP	26	26	NA	NA	NA	NA
35	25.86	Male	wP	58	58	NA	NA	37.5	37.5
36	19.88	Female	aP	73	73	26.0	26.0	48.5	48.5
37	18.91	Female	aP	NA	NA	NA	NA	NA	NA
38	19.88	Female	aP	41	41	NA	NA	72.0	72.0
39	31.92	Female	wP	62	62	NA	NA	NA	NA
40	22.89	Female	wP	34	34	NA	NA	NA	NA
41	31.96	Male	wP	21	21	NA	NA	NA	NA
42	19.92	Female	aP	31	31	NA	NA	17.0	17.0
43	18.91	Female	aP	81	81	NA	NA	29.0	29.0
44	18.91	Female	aP	7	7	4.0	4.0	57.5	57.5
45	19.98	Female	aP	75	75	38.0	38.0	NA	NA
46	18.91	Female	aP	19	19	7.0	7.0	NA	NA
47	20.98	Female	aP	9	9	50.0	50.0	3.0	3.0
48	19.11	Female	aP	90	90	22.0	22.0	6.0	6.0
49	20.11	Female	aP	36	36	47.0	47.0	NA	NA
50	19.98	Female	aP	69	69	NA	NA	1.0	1.0
51	19.98	Male	aP	61	61	NA	NA	NA	NA
52	19.07	Male	aP	44	44	15.0	15.0	2.0	2.0
53	19.07	Female	aP	16	16	NA	NA	4.0	4.0
54	20.11	Female	aP	25	25	NA	NA	NA	NA
55	20.11	Female	aP	42	42	33.0	33.0	NA	NA
56	20.15	Female	aP	18	18	NA	NA	NA	NA
57	21.15	Female	aP	55	55	NA	NA	NA	NA
58	20.15	Female	aP	14	14	NA	NA	NA	NA
59	20.15	Female	aP	2	2	NA	NA	NA	NA
60	20.15	Male	aP	86	86	NA	NA	NA	NA
61	32.38	Female	wP	70	70	NA	NA	39.0	39.0
62	25.99	Female	wP	13	13	NA	NA	13.0	13.0
63	23.98	Female	wP	5	5	40.5	40.5	31.0	31.0
64	25.99	Male	wP	51	51	49.0	49.0	22.0	22.0
65	29.02	Male	wP	6	6	36.0	36.0	10.0	10.0
66	43.07	Female	wP	4	4	48.0	48.0	28.0	28.0
67	47.24	Female	wP	3	3	17.0	17.0	40.5	40.5
68	47.24	Male	wP	12	12	44.0	44.0	40.5	40.5
69	29.17	Female	wP	22	22	39.0	39.0	24.0	24.0
70	21.15	Male	aP	79	79	6.0	6.0	52.0	52.0
71	21.15	Female	aP	33	33	42.0	42.0	20.0	20.0
72	28.25	Female	wP	8	8	20.0	20.0	51.0	51.0
73	24.23	Female	wP	15	15	28.0	28.0	43.0	43.0
74	24.23	Female	wP	65	65	3.0	3.0	59.0	59.0
75	21.22	Female	aP	84	84	NA	NA	53.5	53.5
76	21.22	Female	aP	17	17	32.0	32.0	61.0	61.0
77	31.32	Male	wP	56	56	18.0	18.0	11.0	11.0
78	26.30	Female	wP	37	37	11.0	11.0	56.0	56.0
79	32.32	Male	wP	30	30	2.0	2.0	9.0	9.0
80	27.30	Female	wP	74	74	21.0	21.0	8.0	8.0
81	26.30	Male	wP	54	54	12.0	12.0	12.0	12.0
82	21.28	Female	aP	NA	NA	35.0	35.0	5.0	5.0
83	20.34	Female	aP	63	63	19.0	19.0	70.5	70.5
84	22.34	Female	aP	59	59	24.0	24.0	34.0	34.0
85	19.39	Female	aP	80	80	40.5	40.5	35.0	35.0
86	21.40	Female	aP	47	47	23.0	23.0	32.0	32.0
87	19.39	Male	aP	NA	NA	31.0	31.0	15.0	15.0
88	19.39	Male	aP	NA	NA	37.0	37.0	7.0	7.0
89	22.49	Female	aP	85	85	43.0	43.0	27.0	27.0
90	20.49	Female	aP	76	76	34.0	34.0	47.0	47.0
91	21.49	Male	aP	57	57	8.0	8.0	66.0	66.0
92	19.54	Female	aP	67	67	16.0	16.0	61.0	61.0
93	23.56	Female	aP	45	45	45.0	45.0	53.5	53.5
94	20.55	Male	aP	32	32	9.0	9.0	70.5	70.5
95	21.55	Female	aP	24	24	27.0	27.0	61.0	61.0
96	19.54	Male	aP	64	64	14.0	14.0	68.0	68.0

Model Training and Evaluation

Models have been developed and trained using different predictors. Two models have been constructed for each task and the fold change of each task using Spearman Correlation Coefficient.

1. Age of Participants

2. Baseline values of task variables

knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
x <- c("age_at_boost_Rank", ranked_cols[grepl('D0', ranked_cols)])
y <- ranked_cols[!ranked_cols %in% x ]

corDf <- data.frame()
corP_df <- data.frame()

corrTest <- suppressPackageStartupMessages({t(cor(df[,x], df[,y], 
                                                  method="spearman", 
                                                  use= "pairwise.complete.obs"))})
corrPval <- suppressPackageStartupMessages({cor.mtest(as.matrix(df[,ranked_cols]), 
                                                      method="spearman", 
                                                      conf.level = 0.95, 
                                                      na.rm = TRUE)})$p

corDf["Age", "1.1) IgG-PT-D14-titer-Rank"] <- corrTest["IgG_PT_D14_Rank", "age_at_boost_Rank"]
corDf["Age", "1.2) IgG-PT-D14-FC-Rank"] <- corrTest["IgG_PT_D14_FC_Rank", "age_at_boost_Rank"]
corDf["Age", "2.1) Monocytes-D1-Rank"] <- corrTest["Monocytes_D1_Rank", "age_at_boost_Rank"]
corDf["Age", "2.2) Monocytes-D1-FC-Rank"] <- corrTest["Monocytes_D1_FC_Rank", "age_at_boost_Rank"]
corDf["Age", "3.1) CCL3-D3-Rank"] <- corrTest["CCL3_D3_Rank", "age_at_boost_Rank"]
corDf["Age", "3.2) CCL3-D3-FC-Rank"] <- corrTest["CCL3_D3_FC_Rank", "age_at_boost_Rank"]

corDf["Baseline", "1.1) IgG-PT-D14-titer-Rank"] <- corrTest["IgG_PT_D14_Rank", "IgG_PT_D0_Rank"]
corDf["Baseline", "1.2) IgG-PT-D14-FC-Rank"] <- corrTest["IgG_PT_D14_FC_Rank", "IgG_PT_D0_Rank"]
corDf["Baseline", "2.1) Monocytes-D1-Rank"] <- corrTest["Monocytes_D1_Rank", "Monocytes_D0_Rank"]
corDf["Baseline",  "2.2) Monocytes-D1-FC-Rank"] <- corrTest["Monocytes_D1_FC_Rank", "Monocytes_D0_Rank"]
corDf["Baseline", "3.1) CCL3-D3-Rank"] <- corrTest["CCL3_D3_Rank", "CCL3_D0_Rank"]
corDf["Baseline", "3.2) CCL3-D3-FC-Rank"] <- corrTest["CCL3_D3_FC_Rank", "CCL3_D0_Rank" ]

corP_df["Age", "1.1) IgG-PT-D14-titer-Rank"] <- corrPval["IgG_PT_D14_Rank", "age_at_boost_Rank"]
corP_df["Age", "1.2) IgG-PT-D14-FC-Rank"] <- corrPval["IgG_PT_D14_FC_Rank", "age_at_boost_Rank"]
corP_df["Age", "2.1) Monocytes-D1-Rank"] <- corrPval["Monocytes_D1_Rank", "age_at_boost_Rank"]
corP_df["Age", "2.2) Monocytes-D1-FC-Rank"] <- corrPval["Monocytes_D1_Rank", "age_at_boost_Rank"]
corP_df["Age", "3.1) CCL3-D3-Rank"] <- corrPval["CCL3_D3_Rank", "age_at_boost_Rank"]
corP_df["Age", "3.2) CCL3-D3-FC-Rank"] <- corrPval["CCL3_D3_Rank", "age_at_boost_Rank"]

corP_df["Baseline", "1.1) IgG-PT-D14-titer-Rank"] <- corrPval["IgG_PT_D14_Rank", "IgG_PT_D0_Rank"]
corP_df["Baseline", "1.2) IgG-PT-D14-FC-Rank"] <- corrPval["IgG_PT_D14_FC_Rank", "IgG_PT_D0_Rank" ]
corP_df["Baseline", "2.1) Monocytes-D1-Rank"] <- corrPval["Monocytes_D1_Rank", "Monocytes_D0_Rank"]
corP_df["Baseline", "2.2) Monocytes-D1-FC-Rank"] <- corrPval["Monocytes_D1_FC_Rank", "Monocytes_D0_Rank"]
corP_df["Baseline", "3.1) CCL3-D3-Rank"] <- corrPval["CCL3_D3_Rank", "CCL3_D0_Rank"]
corP_df["Baseline", "3.2) CCL3-D3-FC-Rank"] <- corrPval["CCL3_D3_FC_Rank", "CCL3_D0_Rank"]

Model Evaluation using Spearman Rank Correlation Coefficient

Heatmap of Spearman correlation coefficient for each task against age at boost or the baseline of the task. Columns of Heatmap indicate the task and the task fold change, and rows indicate the model type; age base or baseline base models. Crosses indicate insignificant correlations (p-values > 0.05).

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
maxCorr <- max(abs(as.matrix(corDf)), na.rm = T)
corrplot(as.matrix(corDf), tl.col = 'black', 
         addCoef.col = 1,
         mar=c(0,0,0,2),
         number.cex = 1.2, tl.srt = 45,
         p.mat = as.matrix(corP_df), sig.level = 0.05,
         col.lim = c(-maxCorr, maxCorr),
         cl.ratio = 0.2, col = colorRampPalette(c("blue", "white","red"))(100))

Session Info

sessionInfo()

## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      utils     grDevices stats     graphics  methods   base     
## 
## other attached packages:
##  [1] ellipse_0.4.3         GetoptLong_1.0.5      ComplexHeatmap_2.14.0
##  [4] colorRamp2_0.1.0      formattable_0.2.1     here_1.0.1           
##  [7] glmnet_4.1-6          Matrix_1.5-1          kableExtra_1.3.4     
## [10] data.table_1.14.8     corrplot_0.92         lubridate_1.9.2      
## [13] forcats_1.0.0         stringr_1.5.0         purrr_1.0.1          
## [16] readr_2.1.4           tibble_3.2.0          tidyverse_2.0.0      
## [19] tidyr_1.3.0           ggplot2_3.4.1         dplyr_1.1.0          
## [22] discretefit_0.1.2     pROC_1.18.0           verification_1.42    
## [25] dtw_1.23-1            proxy_0.4-27          CircStats_0.2-6      
## [28] MASS_7.3-58.1         boot_1.3-28           fields_14.1          
## [31] viridis_0.6.2         viridisLite_0.4.1     spam_2.9-1           
## 
## loaded via a namespace (and not attached):
##  [1] matrixStats_0.63.0  RColorBrewer_1.1-3  doParallel_1.0.17  
##  [4] webshot_0.5.4       httr_1.4.5          rprojroot_2.0.3    
##  [7] tools_4.2.2         bslib_0.4.2         utf8_1.2.3         
## [10] R6_2.5.1            BiocGenerics_0.44.0 colorspace_2.1-0   
## [13] withr_2.5.0         tidyselect_1.2.0    gridExtra_2.3      
## [16] compiler_4.2.2      cli_3.6.0           rvest_1.0.3        
## [19] xml2_1.3.3          sass_0.4.5          scales_1.2.1       
## [22] systemfonts_1.0.4   digest_0.6.31       rmarkdown_2.20     
## [25] svglite_2.1.1       pkgconfig_2.0.3     htmltools_0.5.4    
## [28] highr_0.10          fastmap_1.1.1       maps_3.4.1         
## [31] GlobalOptions_0.1.2 htmlwidgets_1.6.1   rlang_1.0.6        
## [34] rstudioapi_0.14     shape_1.4.6         jquerylib_0.1.4    
## [37] generics_0.1.3      jsonlite_1.8.4      magrittr_2.0.3     
## [40] dotCall64_1.0-2     S4Vectors_0.36.2    Rcpp_1.0.10        
## [43] munsell_0.5.0       fansi_1.0.4         lifecycle_1.0.3    
## [46] stringi_1.7.12      yaml_2.3.7          plyr_1.8.8         
## [49] parallel_4.2.2      crayon_1.5.2        lattice_0.20-45    
## [52] splines_4.2.2       circlize_0.4.15     hms_1.1.2          
## [55] knitr_1.42          pillar_1.8.1        rjson_0.2.21       
## [58] stats4_4.2.2        codetools_0.2-18    glue_1.6.2         
## [61] evaluate_0.20       vctrs_0.5.2         png_0.1-8          
## [64] tzdb_0.3.0          foreach_1.5.2       gtable_0.3.1       
## [67] clue_0.3-64         cachem_1.0.7        xfun_0.37          
## [70] survival_3.4-0      iterators_1.0.14    IRanges_2.32.0     
## [73] cluster_2.1.4       timechange_0.2.0    ellipsis_0.3.2

Simple Computational Models for Predicting Pertussis Boost Response

2023-08-21

Setup: Load libraries and helper functions

Setup: Variables used through out the script

Tasks:

Setup: Load dataset, combine different experiments

Setup: Dataset log transform , Ranking and fold change

Model construction

A. Modeling Approach: Age of the participants

Modeling steps:

B. Modeling approach: Baseline values of task variables

Modeling steps:

Model Training and Evaluation

Model Evaluation using Spearman Rank Correlation Coefficient

Session Info