Study design
This current study is part of the Predict Sepsis study, which is a prospective cohort study in the ambulance setting, with patient inclusion between 2017 and 2018, in Stockholm County, Sweden (for details see Wallgren et al.18). The study received approval from the Stockholm Regional Ethical Review Board (reference number 2016/2001–31/2 and 2018/2202). Written informed consent was obtained from all participants. The study was registered at ClinicalTrials.gov, identifier: NCT03249597. The outline of the current study is illustrated in Fig. 1.
Study population
The study included a total of 551 adult, non-trauma patients assessed to have a new onset infection according to clinical judgment made by the ambulance personnel. The inclusion and exclusion criteria have been published elsewhere18. A selection of candidate molecular markers reflecting immune responses were performed in a smaller group of consecutively included patients, i.e., the screening cohort (n = 96). The selected candidate molecular markers from the screening cohort were analyzed in the remaining patients (n = 455), i.e., the prediction cohort, and used as predictors in combination with the available clinical variables measured in the ambulance in the final prediction modeling (details of cohorts, see Table 2).
Blood sampling
Blood was drawn in EDTA-tubes in the ambulance, and at arrival to the hospital bound emergency departments (EDs), tubes were centrifuged, aliquoted, and frozen in − 70 °C in biobank. Furthermore, blood was drawn directly into PAXgene tubes (PreAnalytix, GmbH, Hombrechtikon, Switzerland), with immediate stabilization of intracellular RNA, before being frozen in − 70 °C at arrival to the ED.
Quantification of circulating inflammatory mediators
Initially, a total of 71 circulating proteins were analyzed within the screening cohort using U-PLEX Biomarker Group 1 kits (Meso Scale Discovery, Rockville, MD) detected by electrochemiluminescence in Meso QuickPlex SQ 120 (Meso Scale Discovery), according to the manufacturer’s instructions. Three additional proteins, CXCL6, HGF, and TGF-α were measured by Human Magnetic Luminex Assay (R&D systems, Inc. Minneapolis, MN), according to the manufacturer’s instructions. The samples were analyzed on a Luminex®200™ instrument (Invitrogen, Merelbeke, Belgium), and the data were collected using the xPONENT 3.1™ software (Luminex Corporation, Austin, TX). Later in the prediction cohort, nine selected mediators, i.e., CCL24, CX3CL1, CCL27, CCL11, IL-17AF, IL-17A, IL-1Ra, TNF, and CCL19, were analyzed using customized U-PLEX kits (Meso Scale Discovery, Rockville, MD).
All values were expressed as pg/mL deduced from the standard curve, using a 5-parameter logistic algorithm. Values below the detection limit were given half the value of the detection limit. All samples were run in duplicates and a coefficient of variation (CV) below 20% was considered acceptable. In Supplementary Table 1, the average CVs and detection limits of the nine proteins analyzed in the prediction cohort are listed.
RNA extraction and cDNA extraction
All samples were arranged in random order prior to RNA extraction. For the screening cohort, RNA extraction was done using PAXgene Blood RNA Kit (Cat. No. 172021754, Qiagen, GmbH, Hilden, Germany), according to the manufacturer’s instructions. RNA quality and concentration was measured with NanoDrop 2000 (Thermo Fisher Scientific, MA, USA) and 2100 Bionalyzer (Agilent, CA, USA). The A260/A280 ratios were above 1.7 and the RNA integrity number (RIN) values were above 7. cDNA was synthesized using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, CA, USA) in a LifePro Thermal Cycler (Bioer, Hangzhou, P.R. China), using 200 ng RNA per 20 μL reaction. Gene expression was performed in a Quantstudio 7 Flex Real-Time PCR system (Applied Biosystems, CA, USA), using TaqMan Gene Expression Assays and TaqMan Fast Universal PCR Master Mix (Applied Biosystems, CA, USA) in a 20 μL reaction, according to the manufacturer’s instructions.
For the prediction cohort, the samples together with a negative extraction control (consisting of RNase-free water) were extracted using the QIAsymphony extraction robot (Qiagen GmbH, Hilden, Germany). RNA was eluted in a total volume of 80 μL and immediately denatured at 65 °C for 10 min using a thermal cycler (T-100, BioRad, CA, UAS). Sample concentration and purity were determined by spectrophotometry on the Lunatic instrument (Unchained Labs, CA, USA) and RNA integrity was analyzed on capillary gel electrophoresis, Fragment Analyzer (Agilent, CA, USA) using RNA Standard Sensitivity Fragment Analyzer kit (Cat. No. DNF-471, Agilent, CA, USA). The A260/A280 ratios were above 1.5 except for 21 samples. These samples did not turn out as outliers in neither univariate nor multivariate analyses, thus were not excluded. None of the negative control samples (ENTCs) showed cross-contamination. All samples were reversed transcribed into cDNA using the TATAA GrandScript cDNA Synthesis Kit (Cat. No. A103, TATAA Biocenter AB, Gothenburg, Sweden). The reverse transcription was performed using 450 ng RNA per 20 μL reaction.
Real-time PCR
In the screening cohort, a total of 15 genes encoding inflammatory mediators, inflammasome components, and transcription factors (PYCARD, CASP1, NLRP3, IL1B, IL18, TNF, IL6, IL10, IL1RN, HLA-DRA, HIF1A, SPI1, EPAS1, SIRT1, NFKBIA) were analyzed using qPCR. Gene expression was performed in a Quantstudio 7 Flex Real-Time PCR system (Applied Biosystems, CA, USA), using TaqMan Gene Expression Assays and TaqMan Fast Universal PCR Master Mix (Applied Biosystems, CA, USA), according to the manufacturer’s instructions (TaqMan assay IDs are listed in Supplementary Table 2). HPRT1 was used as a reference gene, determined by NormFinder R package (MOMA, Aarhus University Hospital, Denmark) for normalization among a total of three candidate reference genes. All samples of the study were analyzed in duplicates, and the mean quantity values were used in further data analysis. The accepted CV of technical sample replicates was ≤ 15%. Samples with a CV > 15% for each specific assay were re-analyzed. Cycle threshold (CT) cut-off value was set to 35 and all reactions had an efficiency between 90 and 110%. In all cases, gene expression levels were obtained from a six-point serially four-fold diluted calibration curve. The calibration curve was developed from cDNA of PBMCs stimulated by 1 μg/mL LPS.
In the prediction cohort, three genes, EPAS1, HIF1A, and NLRP3 were analyzed using assays designed and validated by TATAA Biocenter AB. qPCR was performed with TATAA SYBR®GrandMaster Mix Low Rox (Cat. No. TA01, TATAA Biocenter AB, Gothenburg, Sweden) in 10 μL reaction volume. Human ValidPrimeTM (Cat. No. A105P10, TATAA Biocenter AB, Gothenburg, Sweden) was used to monitor and correct for contaminating gDNA25. An inter-plate calibrator (Cat. No. IPC250S, TATAA Biocenter, Gothenburg, Sweden) was run on each plate to be able to correct for inter-run differences.
All samples were run in duplicates in 384-well plate format using QuantStudio™ 7 Pro Real-Time PCR system (384-well, ThermoFisher Scientific). The pipetting was performed by a pipetting robot OT-2 (Opentrons, NY, USA). qPCR raw data were pre-processed and analyzed with GenEx software v.7 (MultiD Analyses AB, Gothenburg, Sweden). The limit of quantification of the assays were determined using standard dilution series for which the relative standard deviation of a replicate was < 35%. The accepted standard deviation of technical sample replicates was ≤ 0.5, whereas the accepted standard deviation of the IPC (Inter Plate Calibrator) replicates was ≤ 0.2. Samples with a standard deviation > 0.5 for each specific assay were re-analyzed. Three reference genes, beta-glucuronidase (GUSB), peptidyl-propyl isomerase A, cyclophilinA (PPIA), and ubiquitin C (UBC) were selected from a list of 12 reference gene candidates using the geNorm and NormFinder functions in GenEx software v.7 (MultiD Analyses AB). The relative gene expression was calculated using the delta CT method.
Statistical modeling and data analysis
In the data analysis pipeline, clinical and molecular variables were assessed with regard to their differences between sepsis and non-sepsis cases and regarding their quality as predictors, using univariate and multi-variable models. For all models, prediction performance was measured in a nested cross-validation approach as AUCs for the hold-out testing set. Finally, variable importance measures, as eligible for the different analysis methods, were applied.
This data analysis pipeline was first run on data from the screening cohort, followed by a consensus variable selection for further analysis in the prediction cohort. Then, the pipeline was run again, this time only involving the selected molecular variables, on the data from the prediction cohort, with prediction performances and variables importance reported as before.
Selection of molecular variables in the screening cohort
Based on screening cohort data, molecular variables as synergistic predictors with the clinical parameters were selected through a stepwise process of (i) univariate analyses, (ii) multivariate analyses and (iii) literature review. From this process, a weighted curation was performed for the final selection of molecular variables, which were used for further analysis in the prediction cohort.
The univariate variable selection of the most relevant molecular variables was performed by fitting individual mixed effect models (lmerTest package in R; mixed effects model with the sepsis/non-sepsis as fixed effect and sex as random effect) of the 74 inflammatory mediators as well as the expression levels of 15 genes, to differentiate between non-septic and septic patients, followed by the false discovery rate (FDR) estimation for multiple comparisons. Molecular variables with fold-changes (FC) above the thresholds, set to FC ˃ 1.2 for proteins and FC > 2.0 for mRNAs, and a Benjamini–Hochberg FDR ˂ 0.05, were selected as candidates for analysis in the prediction cohort.
The multivariate variable selection of the most relevant molecular variables was performed by machine learning implemented with a nested cross-validation workflow assessing variables as classifiers of non-septic and septic patients. A set of seven different machine learning algorithms were trained in parallel and tested to evaluate different algorithms with regards to classification performance on 7 different variable sets. The variable sets were; (a) all molecular variables, (b–d) previously reported Predict Sepsis screening tools 1, 2 and 3 (using un-categorized original values of the clinical parameters presented by previous study18, summary of parameters in the screening tools see Table 1), and (e–g) combining Predict Sepsis screening tools 1, 2 and 3 with the molecular variables. In addition, two-way interaction between all variables were created by multiplying variables for evaluation of interaction effects.
Briefly, all variables included had missing values below 20%. Data were partitioned into training (75%) and testing (25%) sets. Each set of data was standardized, and knn-imputation of missing values and class balanced with Synthetic minority over-sampling technique (SMOTE; themis package in R) was performed on the training set. Thereafter, penalized regularized logistic Lasso regressions (Lasso; glmnet package in R), Random forests (ranger package in R), XGBoosted trees (xgboost package in R), Neural network (nnet package in R), Naïve Bayes (klaR package in R) and lightGBM (bonsai package in R), were in parallel trained with hyperparameter tuning. Hyperparameters for Lasso (penalty); Random Forest (number of variables included in each random tree and minimum n for split); XGBoosted trees (number of variables included in each tree, tree depth, loss reduction, learning rate, and minimum n for split), Neural network (number of hidden layers and penalty), Naïve Bayes (smoothness and laplace), and lightGBM (number of variables included in each tree, tree depth, loss reduction and learn rate) were tuned with a Latin Hypercube search approach with internal validation on 25 bootstraps of the training data with classification performance evaluation scored as AUC. Optimal hyperparameters were used in the final models and performance was assessed on the hold-out testing data. To estimate the robustness of predictions with regards to random effects of partitioning, this workflow procedure was repeated iteratively 20 times as randomized nested cross-validation with random partitioning of samples to training and testing set at each iteration. The variability of model performances from the nested cross-validations was estimated by fitting Bayesian models and Markov Chain Monte Carlo26 via the tidyposterior and rstanarm packages in R, with 5000 iterations, four chains and a prior normal distribution for the random (nest) intercepts. Posterior probability distributions of mean AUC and their contrasted differences between all models were evaluated for practical equivalence. Model variable importance was assessed with permutation of variable values (iml package in R) with loss of classification error as metric, and with model specific weights when applicable (vip package in R). Molecular variables with the highest permuted importance and model specific weights were selected as candidates for analysis in the prediction cohort.
Finally, an evaluation of the literature of the molecular variables with the highest permuted importance and model specific weights was performed in PubMed by using the search terms: cytokine/gene name (both official and alternative); “Sepsis Predicting sepsis using a combination of clinical information and molecular immune markers sampled in the ambulance – Scientific Reports”; “Prediction”. The inclusion criteria for the search were: maximum of 10 years old; species “Human”; language: English; original papers. No exclusion criteria were employed. This evaluation aimed to select molecular variables that had been both previously studied in association with sepsis and novel predictive markers of sepsis.
All candidate molecular variables with especially high importance from either univariate, multivariate in combination with literature mining approaches were selected for further analysis in the prediction cohort, as well as variables that were selected by multiple approaches.
The selected molecular variables from the screening cohort were evaluated in the prediction cohort with respect to their univariate and multivariate discrimination of sepsis and non-sepsis patients. Both the univariate and multivariate analyzes were performed according to the workflows described above. A summary of the workflow is shown in Fig. 1.
Ethics approval and consent to participate
The study received approval from the Stockholm Regional Ethical Review Board (reference number 2016/2001–31/2 and 2018/2202). Written informed consent was obtained from all participants. This study was conducted in accordance with the Declaration of Helsinki27.