Load required libraries
library(omicade4)## Loading required package: ade4library(mogsa)
library(RSpectra)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(glmnet)## Loading required package: Matrix## Loaded glmnet 4.1-7library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:mogsa':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(cowplot)##
## Attaching package: 'cowplot'## The following object is masked from 'package:lubridate':
##
## stamplibrary(ggplot2)Source code from MCIA model
Requires additional MCIA source code from [https://github.com/akonstodata/mcia_mbpca]
source('https://raw.githubusercontent.com/akonstodata/mcia_mbpca/main/R/MCIA_mbpca_extra.R')Load metadata and task matrices
Data source: [https://github.com/joreynajr/cmi-pb-multiomics/tree/main/results/main/cmi_pb_datasets/processed/harmonized] Obtain files below from data source and save to local drive
# data_dir is path to local drive where data is saved
setwd("/Users/rnili/Desktop/repo/gitLab/cmi-pb-multiomics-main/results/main/cmi_pb_datasets/processed/harmonized")
# Read in tasks (2021)
task_mat<-read.table('task_matrix.common_names.mfi_normalised.tsv',sep='\t',header=TRUE,stringsAsFactors=TRUE,row.names=1)
# log normalize certain tasks
# task_mat[,1:7] = log2(task_mat[,1:7]+1)
# task_mat[,11:14] = log2(task_mat[,11:14]+1)
# Read in metadata
meta.2020<-read.table('clinical_metadata.2020.tsv',sep='\t',header=TRUE,stringsAsFactors=TRUE,row.names=1)
meta.2021<-read.table('clinical_metadata.2021.tsv',sep='\t',header=TRUE,stringsAsFactors=TRUE,row.names=1)
# Get age at boost
library(lubridate)
meta.2020$date_of_boost<-parse_date_time(meta.2020$date_of_boost,"ymd")
meta.2020$year_of_birth<-parse_date_time(meta.2020$year_of_birth,"ymd")
meta.2020$age_at_boost<- as.numeric(round(difftime(meta.2020$date_of_boost,
meta.2020$year_of_birth,units="weeks")/52,2))
meta.2021$date_of_boost<-parse_date_time(meta.2021$date_of_boost,"ymd")
meta.2021$year_of_birth<-parse_date_time(meta.2021$year_of_birth,"ymd")
meta.2021$age_at_boost<- as.numeric(round(difftime(meta.2021$date_of_boost,
meta.2021$year_of_birth,units="weeks")/52,2))Load imputed data
Data source: [https://github.com/joreynajr/cmi-pb-multiomics/tree/main/results/main/cmi_pb_datasets/processed/imputed]
Data was imputed using code in: [https://github.com/CMI-PB/2020_2021_data_merge/tree/main/scripts/CMIPB_2020_2021_imputation.Rmd]
Obtain files below from data source and save to local drive
##plot distributions
distPlot <- function(col, df){
p <-
ggplot(df) +
aes_string(col)
if(is.numeric(df[[col]])) {
p <- p + geom_density()
} else {
p <- p + geom_bar()
}
}
# is.vector(names(dataDf))
# length(dataDf)
# names(task_mat)[2:length(task_mat)]
typeof(task_mat)## [1] "list"#
distPlots <- lapply(names(task_mat), distPlot, df=task_mat)## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plot_grid(plotlist = distPlots)## Warning: Removed 37 rows containing non-finite values (`stat_density()`).## Warning: Removed 40 rows containing non-finite values (`stat_density()`).## Warning: Removed 37 rows containing non-finite values (`stat_density()`).## Warning: Removed 21 rows containing non-finite values (`stat_density()`).
## Removed 21 rows containing non-finite values (`stat_density()`).
## Removed 21 rows containing non-finite values (`stat_density()`).
## Removed 21 rows containing non-finite values (`stat_density()`).
# imputed_dir is path to local drive where data is saved
setwd("/Users/rnili/Desktop/repo/gitLab/cmi-pb-multiomics-main/results/main/cmi_pb_datasets/processed/imputed")
# Import imputed datasets
rnaseq_baseline_mat_imputed_20 <- read.csv('rnaseq_baseline_mat_imputed_20_051022.csv',row.names=1)
cytof_baseline_mat_imputed_20 <- read.csv('cytof_baseline_mat_imputed_20_051022.csv',row.names=1)
olink_baseline_mat_imputed_20 <- read.csv('olink_baseline_mat_imputed_20_051022.csv',row.names=1)
abtiters_baseline_mat_imputed_20 <- read.csv('abtiters_baseline_mat_imputed_20_051022.csv',row.names=1)
rnaseq_baseline_mat_imputed_21 <- read.csv('rnaseq_baseline_mat_imputed_21_051022.csv',row.names=1)
cytof_baseline_mat_imputed_21 <- read.csv('cytof_baseline_mat_imputed_21_051022.csv',row.names=1)
olink_baseline_mat_imputed_21 <- read.csv('olink_baseline_mat_imputed_21_051022.csv',row.names=1)
abtiters_baseline_mat_imputed_21 <- read.csv('abtiters_baseline_mat_imputed_21_051022.csv',row.names=1)Check whether features are the same between 2020 and 2021 data
# Check features the same
identical(colnames(rnaseq_baseline_mat_imputed_20),colnames(rnaseq_baseline_mat_imputed_21))## [1] TRUEidentical(colnames(cytof_baseline_mat_imputed_20),colnames(cytof_baseline_mat_imputed_21))## [1] TRUEidentical(colnames(olink_baseline_mat_imputed_20),colnames(olink_baseline_mat_imputed_21))## [1] TRUEidentical(colnames(abtiters_baseline_mat_imputed_20),colnames(abtiters_baseline_mat_imputed_21))## [1] TRUEPrepare data matrix for MCIA model on 2020 data
data_2020<-list(cytof=t(cytof_baseline_mat_imputed_20),
seq=t(rnaseq_baseline_mat_imputed_20),
olink=t(olink_baseline_mat_imputed_20),
abtiters=t(abtiters_baseline_mat_imputed_20))
lapply(data_2020,dim)## $cytof
## [1] 22 36
##
## $seq
## [1] 11589 36
##
## $olink
## [1] 27 36
##
## $abtiters
## [1] 19 36Run MCIA with 10 factors
set.seed(0)
num_comps=10
mcia_out_gs<-mcia_mbpca(data_2020,num_comps=num_comps,preprocess='nsc',block_prep='lambda_all',
deflat_method="globalScore") ## Function 'svd' used, all singular vectors will be computed
## Function 'svd' used, all singular vectors will be computed
## Function 'svd' used, all singular vectors will be computed
## Function 'svd' used, all singular vectors will be computed## calculating component 1 ...
## calculating component 2 ...
## calculating component 3 ...
## calculating component 4 ...
## calculating component 5 ...
## calculating component 6 ...
## calculating component 7 ...
## calculating component 8 ...
## calculating component 9 ...
## calculating component 10 ...Create lists of associated prediction tasks
tasks_ab<-c('IgG.PT','IgG.FHA','IgG.PRN','IgG1.PT','IgG1.FHA','IgG1.PRN','IgG4.PT','IgG4.FHA','IgG4.PRN')
tasks_cytof<-c('Monocytes','ASCs..Plasmablasts.','CD4Tcells')
tasks_seq<-c('ENSG00000277632','ENSG00000136244','ENSG00000100906','ENSG00000229807')Showing required tasks being added from imputed files are already normalized
names(rnaseq_baseline_mat_imputed_20[tasks_seq])## [1] "ENSG00000277632" "ENSG00000136244" "ENSG00000100906" "ENSG00000229807"distPlot <- function(col, df){
p <-
ggplot(df) +
aes_string(col)
if(is.numeric(df[[col]])) {
p <- p + geom_density()
} else {
p <- p + geom_bar()
}
}
#
distPlots <- lapply(tasks_seq, distPlot, df=rnaseq_baseline_mat_imputed_20)
plot_grid(plotlist = distPlots)
Add demographic and baseline features to MCIA factor model to create MCIAplus model
gs_2020_temp = mcia_out_gs$mcia_result$t
gs_2020<-data.frame(gs_2020_temp)
gs_2020$age<-meta.2020[rownames(gs_2020),"age_at_boost"]
gs_2020$infancy_vac<-as.numeric(meta.2020[rownames(gs_2020),'infancy_vac'])
gs_2020$biological_sex<-as.numeric(meta.2020[rownames(gs_2020),'biological_sex'])
n_cols = dim(gs_2020)[2]
gs_2020[,(n_cols+1):(n_cols+length(tasks_ab))]<-abtiters_baseline_mat_imputed_20[rownames(gs_2020),tasks_ab]
colnames(gs_2020)[(n_cols+1):(n_cols+length(tasks_ab))]<-tasks_ab
n_cols = dim(gs_2020)[2]
gs_2020[,(n_cols+1):(n_cols+length(tasks_seq))]<-rnaseq_baseline_mat_imputed_20[rownames(gs_2020),tasks_seq]
colnames(gs_2020)[(n_cols+1):(n_cols+length(tasks_seq))]<-tasks_seq
n_cols = dim(gs_2020)[2]
gs_2020[,(n_cols+1):(n_cols+length(tasks_cytof))]<-cytof_baseline_mat_imputed_20[rownames(gs_2020),tasks_cytof]
colnames(gs_2020)[(n_cols+1):(n_cols+length(tasks_cytof))]<-tasks_cytofShowing everything except demographic data and MCIA results are normalized. Regardless of
names(gs_2020)## [1] "X1" "X2" "X3"
## [4] "X4" "X5" "X6"
## [7] "X7" "X8" "X9"
## [10] "X10" "age" "infancy_vac"
## [13] "biological_sex" "IgG.PT" "IgG.FHA"
## [16] "IgG.PRN" "IgG1.PT" "IgG1.FHA"
## [19] "IgG1.PRN" "IgG4.PT" "IgG4.FHA"
## [22] "IgG4.PRN" "ENSG00000277632" "ENSG00000136244"
## [25] "ENSG00000100906" "ENSG00000229807" "Monocytes"
## [28] "ASCs..Plasmablasts." "CD4Tcells"distPlot <- function(col, df){
p <-
ggplot(df) +
aes_string(col)
if(is.numeric(df[[col]])) {
p <- p + geom_density()
} else {
p <- p + geom_bar()
}
}
# is.vector(names(dataDf))
# length(dataDf)
# names(gs_2020)[2:length(gs_2020)]
typeof(gs_2020)## [1] "list"#
distPlots <- lapply(names(gs_2020), distPlot, df=gs_2020)
plot_grid(plotlist = distPlots)
Test model quality
Analyze correlation between predicted and true values in leave-one-out cross-validation assay
pred_cor<-data.frame(matrix(nrow=ncol(task_mat),ncol=3))
rownames(pred_cor)<-colnames(task_mat)
colnames(pred_cor)<-c('pearson.cor.pred.true', 'spearman.cor.pred.true', 'ranked.spearman.cor.pred.true')
# Loop through all tasks
for (i in 1:ncol(task_mat)){
all_preds<-c()
all_true<-c()
set.seed(1)
x_out <-data.frame(task_mat[,i])
rownames(x_out)<-rownames(task_mat)
x_out$temp<-'temp'
names(x_out)<-c('Y','temp')
x_out_r<-na.omit(x_out)
row_int<-intersect(rownames(x_out_r),rownames(gs_2020))
x_out_r<-x_out_r[row_int,]
gs_2020_filt<-gs_2020[rownames(x_out_r),]
#train = 1:nrow(gs_2020_filt)
# train all leave-one-out models
for (j in 1:nrow(gs_2020_filt)){
train = 1:nrow(gs_2020_filt)
train = train[-c(j)]
a <- colnames(gs_2020_filt)
# create lasso model
suppressWarnings(cvfit_out<-cv.glmnet(x=as.matrix(gs_2020_filt[train,]), x_out_r[train,'Y'], family='gaussian',
alpha=1,nfolds=nrow(gs_2020_filt[train,])))
preds<-predict(cvfit_out,newx=as.matrix(data.frame(gs_2020_filt[-train,])),s='lambda.min')
all_preds<-c(all_preds,preds)
all_true<-c(all_true,x_out_r[-train,'Y'])
}
# pred_cor[i,'cor.pred.true']<-cor(all_preds,all_true) # Can also use spearman correlation
pred_cor[i,'pearson.cor.pred.true'] <- cor(all_preds,all_true)
pred_cor[i,'spearman.cor.pred.true'] <- cor(all_preds,all_true, method="spearman")
pred_cor[i,'ranked.spearman.cor.pred.true'] <- cor(rank(all_preds,na.last="keep",ties.method="min"),rank(all_true,na.last="keep",ties.method="min"), method="spearman")
}
# print(pred_cor)For each model, can assess which features contribute to non-zero coefficients
Consider only choosing models for follow-on analysis that show good correlation scores
all_models_coef<-vector(mode='list',length=14)
all_models_names<-vector(mode='list',length=14)
all_models<-vector(mode='list',length=14)
for (i in 1:14){
set.seed(1)
x_out <-data.frame(task_mat[,i])
rownames(x_out)<-rownames(task_mat)
x_out$temp<-'temp'
names(x_out)<-c('Y','temp')
x_out_r<-na.omit(x_out)
row_int<-intersect(rownames(x_out_r),rownames(gs_2020))
x_out_r<-x_out_r[row_int,]
gs_2020_filt<-gs_2020[rownames(x_out_r),]
# create lasso model
suppressWarnings(cvfit_out<-cv.glmnet(x=as.matrix(gs_2020_filt), y=as.matrix(x_out_r[,'Y']), family='gaussian',
alpha=1,nfolds=nrow(gs_2020_filt-1)))
all_models_coef[i]=list(coef(cvfit_out, s = 'lambda.min')[coef(cvfit_out, s = 'lambda.min')[,1]!= 0])
all_models_names[i]=list(rownames(coef(cvfit_out, s = 'lambda.min'))[coef(cvfit_out, s = 'lambda.min')[,1]!= 0])
}
names(all_models_coef)<-colnames(task_mat)
names(all_models_names)<-colnames(task_mat)
for (i in 1:14){
all_models[[i]] = data.frame(cbind(all_models_names[[i]],all_models_coef[[i]]))
colnames(all_models[[i]])<-c("Variable","Coefficient")
all_models[[i]]$Coefficient<-as.numeric(all_models[[i]]$Coefficient)
all_models[[i]]$Coefficient=round(all_models[[i]]$Coefficient,3)
}
names(all_models)<-colnames(task_mat)
for (i in 1:14){
all_models[[i]]$Variable[all_models[[i]]$Variable==tasks_seq[1]]="CCL3"
all_models[[i]]$Variable[all_models[[i]]$Variable==tasks_seq[2]]="IL6"
all_models[[i]]$Variable[all_models[[i]]$Variable==tasks_seq[3]]="NFKBIA"
all_models[[i]]$Variable[all_models[[i]]$Variable==tasks_seq[4]]="XIST"
all_models[[i]]<-all_models[[i]] %>% arrange(desc(abs(Coefficient)))
}
all_models[1]## $IgG.PT_day14
## Variable Coefficient
## 1 X4 19.118
## 2 (Intercept) 5.326
## 3 IgG1.PT 0.9652021
Generate global scores on 2021 data using 2020 model
#preprocess2021 data
num_comps=10
block_prep_use = 'lambda_all'
data_2021<-list(cytof=t(cytof_baseline_mat_imputed_21),
seq=t(rnaseq_baseline_mat_imputed_21),
olink=t(olink_baseline_mat_imputed_21),
abtiters=t(abtiters_baseline_mat_imputed_21))
table_out<-nsc_prep(data_2021, num_comps)
final_out<-processOpt(table_out,scale=FALSE,center=FALSE,num_comps,option=block_prep_use)## Function 'svd' used, all singular vectors will be computed
## Function 'svd' used, all singular vectors will be computed
## Function 'svd' used, all singular vectors will be computed
## Function 'svd' used, all singular vectors will be computed# generate 2021 global scores using 2020 models
gs_2021_temp=new_gs(final_out,mcia_out_gs$mcia_result)
rownames(gs_2021_temp)<-rownames(meta.2021)
gs_2021<-data.frame(gs_2021_temp)
# Add demographic variables to 2021 MCIA global scores
gs_2021$age<-meta.2021[rownames(gs_2021),"age_at_boost"]
gs_2021$infancy_vac<-as.numeric(meta.2021[rownames(gs_2021),'infancy_vac'])
gs_2021$biological_sex<-as.numeric(meta.2021[rownames(gs_2021),'biological_sex'])
# Add day 0 values of tasks to 2021 MCIA global scores to make 2021 'MCIAplus' model
n_cols = dim(gs_2021)[2]
gs_2021[,(n_cols+1):(n_cols+length(tasks_ab))]<-abtiters_baseline_mat_imputed_21[rownames(gs_2021),tasks_ab]
colnames(gs_2021)[(n_cols+1):(n_cols+length(tasks_ab))]<-tasks_ab
n_cols = dim(gs_2021)[2]
gs_2021[,(n_cols+1):(n_cols+length(tasks_seq))]<-rnaseq_baseline_mat_imputed_21[rownames(gs_2021),tasks_seq]
colnames(gs_2021)[(n_cols+1):(n_cols+length(tasks_seq))]<-tasks_seq
n_cols = dim(gs_2021)[2]
gs_2021[,(n_cols+1):(n_cols+length(tasks_cytof))]<-cytof_baseline_mat_imputed_21[rownames(gs_2021),tasks_cytof]
colnames(gs_2021)[(n_cols+1):(n_cols+length(tasks_cytof))]<-tasks_cytofGenerate predictions for 2021 tasks
# initiate matrices
val_predict<-data.frame(matrix(nrow = nrow(gs_2021), ncol=ncol(task_mat)))
rank_predict<-data.frame(matrix(nrow=nrow(gs_2021),ncol=ncol(task_mat)))
rownames(rank_predict)<-rownames(val_predict)<-rownames(gs_2021)
colnames(rank_predict)<-colnames(val_predict)<-colnames(task_mat)
# If results of 2020 CV testing were poor, make predictions based on day 0 values
val_predict$IgG.FHA_day14<-gs_2021$IgG.FHA
val_predict$IgG.PRN_day14<-gs_2021$IgG.PRN
val_predict$IgG1.FHA_day14<-gs_2021$IgG1.FHA
val_predict$Monocytes_day1<-gs_2021$Monocytes
val_predict$ASCs..Plasmablasts._day7<-gs_2021$ASCs..Plasmablasts.
val_predict$CD4Tcells_day3<-gs_2021$CD4Tcells
# If results of 2020 CV testing were good, use model to make predictions
set.seed(1)
for (i in c(1,4,6,7,11,12,13,14)){ # choose variables here which to predict with model vs. with day 0 values
x_out <-data.frame(task_mat[,i])
rownames(x_out)<-rownames(task_mat)
x_out$temp<-'temp'
names(x_out)<-c('Y','temp')
x_out_r<-na.omit(x_out)
row_int<-intersect(rownames(x_out_r),rownames(gs_2020))
x_out_r<-x_out_r[row_int,]
gs_2020_filt<-gs_2020[rownames(x_out_r),]
# create lasso model using 2020 data
cvfit_out<-cv.glmnet(x=as.matrix(gs_2020_filt), x_out_r$Y, family='gaussian',
alpha=1,nfolds=nrow(gs_2020_filt),type.measure="mse")
# make predictions on 2021 data
preds<-data.frame(predict(cvfit_out,newx=as.matrix(gs_2021),s='lambda.min'))
val_predict[,i]<-preds
}## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldModify predictions to ‘NA’ where original data had ‘NA’
val_predict[c("61","62","75"),c("Monocytes_day1","ASCs..Plasmablasts._day7","CD4Tcells_day3")]<-NaN
val_predict[c("82","87","88"),1:7]<-NaN
for (i in 1:14){
rank_predict[,i]<-rank(-val_predict[,i],na.last="keep",ties.method="min")
}Clean up presentation to fit template
rank_predict_reord<-colnames(rank_predict)
rank_predict_reord[8:10]<-c(rank_predict_reord[9],rank_predict_reord[10],rank_predict_reord[8])
rank_predict_out<-rank_predict[,rank_predict_reord]
rank_predict_out$subjectID = row.names(rank_predict_out)
print(row.names(rank_predict_out))## [1] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70" "71" "72" "73" "74" "75"
## [16] "76" "77" "78" "79" "80" "81" "82" "83" "84" "85" "86" "87" "88" "89" "90"
## [31] "91" "92" "93" "94" "95" "96"new_rank_predict_out <- rank_predict_out %>% select(subjectID, rank_predict_reord)## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(rank_predict_reord)
##
## # Now:
## data %>% select(all_of(rank_predict_reord))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.write.table(new_rank_predict_out,'Anna_cmipb_2020_preds_05132022.csv',sep=',', row.names = FALSE)setwd("/Users/rnili/Desktop/repo/gitLab/cmi-pb-multiomics-main/rasteh/models")
a= append(list(pred_cor),all_models)
names(a)[1] <- 'pred_cor'
sink("allModels_predCor_annas_notNormalised.txt")
print(a)## $pred_cor
## pearson.cor.pred.true spearman.cor.pred.true
## IgG.PT_day14 0.05005005 0.25199485
## IgG.FHA_day14 -0.97222212 -0.99903179
## IgG.PRN_day14 0.00249668 0.03140283
## IgG1.PT_day14 0.58118535 0.54337194
## IgG1.FHA_day14 0.18925577 0.22651223
## IgG4.PT_day14 0.65688921 0.71636669
## IgG4.FHA_day14 0.73200440 0.80746461
## Monocytes_day1 0.23070491 0.23235294
## ASCs..Plasmablasts._day7 -0.43629061 -0.76968380
## CD4Tcells_day3 0.29858463 0.44705882
## CCL3_day3 0.68393955 0.58249678
## IL.6_day3 0.39187182 0.40489060
## NFKBIA_day7 0.37273291 0.44555985
## XIST_day14 0.85969487 0.63552124
## ranked.spearman.cor.pred.true
## IgG.PT_day14 0.25199485
## IgG.FHA_day14 -0.99903179
## IgG.PRN_day14 0.03140283
## IgG1.PT_day14 0.54337194
## IgG1.FHA_day14 0.22651223
## IgG4.PT_day14 0.71636669
## IgG4.FHA_day14 0.80746461
## Monocytes_day1 0.23235294
## ASCs..Plasmablasts._day7 -0.76968380
## CD4Tcells_day3 0.44705882
## CCL3_day3 0.58249678
## IL.6_day3 0.40489060
## NFKBIA_day7 0.44555985
## XIST_day14 0.63552124
##
## $IgG.PT_day14
## Variable Coefficient
## 1 X4 19.118
## 2 (Intercept) 5.326
## 3 IgG1.PT 0.965
##
## $IgG.FHA_day14
## Variable Coefficient
## 1 (Intercept) 6.055
##
## $IgG.PRN_day14
## Variable Coefficient
## 1 X3 61.510
## 2 X5 34.396
## 3 X8 6.916
## 4 (Intercept) 4.895
## 5 IgG1.FHA 4.249
## 6 X4 2.755
## 7 IgG1.PT -1.329
## 8 IgG4.FHA -1.329
## 9 IgG.PT -0.221
## 10 biological_sex 0.089
##
## $IgG1.PT_day14
## Variable Coefficient
## 1 X4 20.327
## 2 (Intercept) 9.396
## 3 IgG1.PT 3.201
## 4 IgG4.PT -0.597
##
## $IgG1.FHA_day14
## Variable Coefficient
## 1 (Intercept) 2.389
## 2 X1 -1.080
## 3 IgG1.FHA 0.278
## 4 IgG.FHA 0.022
##
## $IgG4.PT_day14
## Variable Coefficient
## 1 X5 -187.127
## 2 (Intercept) -132.456
## 3 IgG4.PT 65.695
## 4 NFKBIA 3.739
## 5 Monocytes 3.346
## 6 ASCs..Plasmablasts. 0.047
##
## $IgG4.FHA_day14
## Variable Coefficient
## 1 X9 140.643
## 2 IgG4.FHA 7.474
## 3 IgG1.FHA -4.271
## 4 IL6 1.854
## 5 (Intercept) -1.102
## 6 ASCs..Plasmablasts. 0.657
## 7 infancy_vac -0.534
## 8 IgG1.PT -0.337
## 9 IgG4.PRN 0.312
##
## $Monocytes_day1
## Variable Coefficient
## 1 X7 -128.695
## 2 X10 -105.372
## 3 X2 78.080
## 4 (Intercept) 20.358
## 5 IL6 -1.872
## 6 Monocytes 0.319
## 7 X9 0.149
##
## $ASCs..Plasmablasts._day7
## Variable Coefficient
## 1 (Intercept) 1.642
##
## $CD4Tcells_day3
## Variable Coefficient
## 1 (Intercept) 28.882
## 2 biological_sex -6.052
## 3 CD4Tcells 0.321
##
## $CCL3_day3
## Variable Coefficient
## 1 (Intercept) -2441.556
## 2 X10 611.878
## 3 NFKBIA 267.649
## 4 infancy_vac -57.397
## 5 CCL3 2.808
##
## $IL.6_day3
## Variable Coefficient
## 1 (Intercept) -745.467
## 2 X10 485.221
## 3 NFKBIA 79.557
## 4 infancy_vac -2.439
##
## $NFKBIA_day7
## Variable Coefficient
## 1 IL6 431.859
## 2 infancy_vac -295.524
## 3 CCL3 172.327
## 4 (Intercept) -141.283
##
## $XIST_day14
## Variable Coefficient
## 1 X10 196.115
## 2 (Intercept) 146.176
## 3 X4 138.266
## 4 biological_sex -99.696
## 5 X1 -69.147
## 6 X6 13.201
## 7 NFKBIA 5.279
## 8 IgG.PRN -3.186
## 9 CCL3 1.807
## 10 IgG1.PT 1.446
## 11 ASCs..Plasmablasts. -0.183sink()