High-throughput mass spectrometry maps the sepsis plasma proteome and differences in response

Abbreviation table

The UK Genomic Advances in Sepsis (GAinS) Study

Other cohorts

Vasopressin vs Norepinephrine as Initial Therapy in Septic Shock (VANISH; REC reference 12/SC/0014) was a randomized clinical trial enrolling adult patients who had septic shock requiring vasopressors despite fluid resuscitation. Details of the study design and patient recruitment have been previously described (Gordon et al., 2016). The study reported no significant difference in outcome between vasopressin and norepinephrine groups, or between hydrocortisone and placebo groups. Up to four timepoints were included for each of the 45 VANISH patients in this study. The baseline timepoint (TP0) was mostly taken within 6 hours after the onset of shock and before the study drugs were given. Time from diagnosis of shock to study drug1 administration was within 6hr for 43 of the patients, and had a median of 4hr. Blood samples were collected into 10mL EDTA vacutainers on ice. The tube was inverted 5-8 times and then centrifuged at 1000RCF for 10min before the plasma layer was aliquoted.

The Oxford BioBank (OBB; REC reference 18/SC/0588) is a random, population-based biobank of healthy participants initially recruited between the ages of 30 and 50 years from the Oxfordshire general population (Karpe et al., 2018). Individuals with: previous diagnosis of myocardial infarction or heart failure currently on treatment; untreated malignancies; or other systemic ongoing disease, and pregnant women were excluded. One sample was included for each of the 152 OBB participants in this study. To match the demographics of other cohorts as closely as possible, older OBB participants were selected. Plasma samples were collected with EDTA as the anti-coagulant.

X-MINS (Acquired loss of cardiac vagal activity is associated with myocardial injury in patients undergoing non-cardiac surgery; REC reference 16/LO/0635) was an observational mechanistic cohort study at the University College London Hospital, recruiting adult patients undergoing major elective noncardiac surgeries including orthopedic, upper GI and colorectal surgeries. The X-MINS cohort has been utilized to examine whether serial measures of cardiac vagal dysfunction were associated with perioperative myocardial injury and noncardiac morbidity (May et al., 2019), and had measured microRNA concentrations in extracellular vesicles (May et al., 2020). Two samples for each of the 106 X-MINS patients were included in this study, one taken pre-operation (PreOp (disc)) and one taken within 24hr post-operation (PostOp (disc)). Blood samples were collected into EDTA vacutainers and centrifuged at 3500RCF for 10min before the plasma layer was aliquoted sterile filtered.

BIONIC (Biomarker based Identification Of Nosocomial Infective Complications; REC reference 14/EM/1223) was a prospective observational study at the Royal London Hospital that recruited patients undergoing elective major abdominal surgeries encompassing colorectal, upper GI (gastrointestinal) and HPB (Hepato-Pancreato-Biliary) surgeries. For each of the 43 patients included in this study, three samples were included, taken either immediately before induction of anesthesia (PreOp (vali)), or 24hr post-operation (PostOp (vali)), or 48hr post-operation. Ten patients also had a 2-6hr post-operation sample included. Blood samples were collected into 4.5mL citrate vacutainers and centrifuged at 1860RCF for 10min before the plasma layer was aliquoted.

The MOTION (Methylnaltrexone for the Treatment of Opioid Induced Constipation and Gastrointestinal Stasis in Intensive Care Patients; REC reference 14/LO/2004) trial recruited adult ICU patients who were mechanically ventilated, receiving opioids and were constipated (Patel et al., 2020). For the 50 patients included in this study, reasons for ICU admission included non-operative medical causes (n=34) and emergency(n=13) or elective (n=3) operative procedures. Patients admitted to ICU for infection or respiratory causes were not included. One sample from each patient, taken before the study drug was given, was included for this study. Blood samples were collected into citrate vacutainers on ice and centrifuged at 1500RCF for 15min before the plasma layer was aliquoted.

The MONOGRAM (REC reference 15/LO/0933) study was set up for diagnosing and monitoring infection in critically ill patients using metabolic and immunological signatures. Baseline samples taken within 48hr of ICU admission for 14 mechanically ventilated patients with systemic inflammatory response syndrome (SIRS) but without identified infection were included in this study. Blood samples were collected into 10mL EDTA vacutainers on ice. The tube was inverted 5-8 times and then centrifuged at 1000RCF for 10min before the plasma layer was aliquoted.

TACE (REC reference 10/H0709/77) is a general abbreviation for observational studies on the mechanisms of monocyte priming and tolerance in vitro and in vivo to ascertain TNF-α converting enzyme (TACE) activity and metabolic signatures of patients with direct and indirect acute lung injury (ALI). Adult patients with or at risk of ALI, admitted to ICU within Imperial College Healthcare

NHS Trust, were recruited within 48hr of onset of ALI or intubation. More details were described in associated publications (Antcliffe et al., 2017, 2018; O’Callaghan et al., 2015). The baseline samples on study entry for twelve patients with noninfectious SIRS (ten with brain injury, one with cardiac arrest, one with motor neurone disease) were included in this study. Plasma samples were collected with EDTA or heparin as the anti-coagulant.

Sample preparation of non-depleted plasma for high-throughput LC-MS/MS platform

Plasma samples were thawed at 4◦C. 50μL aliquots of each of the 2622 clinical plasma samples from the 8 cohorts were pipetted into 28 96-well plates in randomized positions. Three microliters of plasma were then diluted in 50mM ammonium bicarbonate at 1:30 ratio (sample: buffer) in 96-well plates. Plates were centrifuged at 2000g for five minutes before transferring thirty microliters of diluted plasma into clean 96-well plates for subsequent in-solution digest using the BravoAssaymap liquid handler robot (Agilent). In brief, samples were reduced with DTT (15mM final concentration) for 30 minutes at room temperature, alkylated with iodoacetamide (70mM final concentration) for 30 min at room temperature followed by the addition of 100 μl of 50mM ammonium bicarbonate containing 10mM DTT to quench the reaction. Samples were incubated with 1 μg of Trypsin (Worthington TPCK) at 37C overnight. Trypsin digestion was stopped by adding trifluoroacetic acid at 1% final concentration.

Pools

A pool was built using thirty microliters of the diluted plasma for each patient. This pool was included in each digestion plate as a “QC Plate pool” to monitor trypsin digestion. An additional pool was created post digestion (labelled as QC system) and run every twenty-four samples to monitor the performance of the LC system across the whole dataset.

Generation of a plasma library combining top 64 depletion with High pH fractionation

Plasma has a high dynamic range at protein level. To increase the detection range of proteins in plasma and to generate a protein library we created a super-depleted plasma by combining two-consecutive antibody-based depletions using first the SepproIgY14 column (SEP030 column Sigma-Aldrich) to deplete for the top 14 most abundant plasma proteins followed by a Seppro Supermix column (SPE050 column, Sigma-Aldrich). 200μL of the pooled plasma was diluted to 1ml with ‘Dilution buffer’ containing 0.1%N-Octyl-B-D-glucopyranoside then applied to a prewashed Costar 45uM spin filter and cleared by centrifugation 14,000 x g for 10mins at room temperature. 1ml of the filtrate was injected into a BioRad NGC medium pressure system following the method described (Keshishian et al., 2017), with some additional changes to accommodate the use of the two columns. In brief, the two depletion columns were attached in tandem, first passing the Seppro IgY14 column then the Supermix column resulting in the unbound depleted plasma directly flowing into the second column. The columns were pre-equilibrated for 3.33mins at 1.5ml/min in dilution buffer followed by sample application at 1ml/min for 1min then 0.5ml/min for a further 3mins. Flow rate was then increased to 1.5ml/min and dilution buffer isocratically for 15mins. The column was then stripped of the proteins which had been depleted, using stripping buffer at a flow rate 2mls/min for 6mins. Depleted proteins were collected in 0.5ml fractions between min 12.33 and 22.33 and in 1ml fractions until min 28.33. The column was then neutralized and re-equilibrated into dilution buffer for 9.33mins at 1.5ml/min. 50μL of each fraction was evaluated using Ag stain and the depleted fractions pooled together according to the Abs280nm. The plasma elution profile showed that the depleted peak appeared between min 3.5 and 13.5 and accounted for approximately 40% of the total protein present in the sample. The fractions were pooled into 15 pools containing approximately 3% total proteins present and labelled ‘A’ to ‘O’. Each fraction ‘A’ to ‘O’ was chloroform/methanol precipitated and the pellets resuspended into 50μL of 6M Urea in 0.1M Tris pH 7.8. To ensure protein pellets were solubilized, samples were mixed at room temperature for 20min and sonicated for 10min in a sonicating water bath. Urea concentration was diluted down to 1M by adding 250μL H₂O. Samples were incubated with 1μg of smart trypsin (ThermoFisher) at 70°C for an 1h shaking at 300rpm. 150μL of the sample was then HpH fractionated into 12 fractions on RP-W cartridges using the Bravo Assaymap liquid handler (Agilent technologies). Each individual fraction was then loaded onto C18 EvoTips for subsequent LC-MS/MS analysis.

LC_MS/MS using the high-throughput Evosep One - Bruker timsTOF Pro platform

Samples were analyzed using an Evosep One LC system connected to the timsTOF Pro mass spectrometer (Bruker Daltonics) (COvid-19 Multi-omics Blood ATlas (COMBAT) Consortium, 2022). Briefly, Peptides were analyzed using the pre-built 100 samples / per day method (EvosepOne) with an 11.5 min gradient (total cycle time of 14.4min) at a 1.2 µl/min flow rat (Bache et al.). In brief, tryptic peptides were transferred from the pre-loaded C18 Evotips with a pre-built gradient to a sample loop and separated on an 8cm C18 analytical column (Evosep Pepsep, 3um beads, 100 um ID) with an overall gradient from 3 to 40% acetonitrile.

Mass spectrometry data were acquired in PASEF mode (oTOF control v6.0.0.12). The ion mobility window was set to 1/k0 start = 0.85 Vs/cm2 to 1/k0 end = 1.3 Vs/cm2, ramp time 100 ms with locked duty cycle, mass range 100 - 1700 m/z. MS/MS were acquired in 4 PASEF frames (3 cycles overlap). Target intensity was set to 6000 and threshold intensity 200. The data acquisition method has been deposited with the raw data in the PRIDE data repository (deposition in progress, accession ID TBC). To assess the technical reproducibility, QC system pools were run every twenty-four samples. QC plate pools were run to monitor the performance of the trypsin digestion in each plate.

QE-HF mass spectrometry

Luminex assay

Sixty-five inflammatory mediators were measured in 204 samples from 146 patients from the GAinS study using the ProcartaPlex™ Luminex platform (ThermoFisher Scientific, Waltham, MA, USA), as has been described previously (Antcliffe et al., 2022). Briefly, plasma samples were measured across three 96-well plates with data acquired on Luminex 100 system. For each analyte, background fluorescence was corrected for and the lower/upper limits of quantification (LOQs) were determined. Absolute concentrations were derived from concentration-response standard curves fitted with a five-parameter logistic model for each analyte on each plate. Sample concentrations beyond the LOQs were censored at the LOQ values. Outlying samples of poor quality were identified and removed.

Leukocyte whole-genome expression profiling

For the GAinS patient-timepoints included in this study, 1354 (1010 patients) have paired leukocyte gene expression data available through either microarray profiling (518 patients, 642 samples), RNA sequencing (649 patients, 837 samples), or both platforms.

The microarray datasets were acquired using Illumina Human-HT-12 v4 Expression BeadChips, full details of which have been described previously (Davenport et al. 2016; Burnham et al. 2017). Four batches of microarray data have been cleaned-up, combined and the batches normalized. Gene expression was represented at gene-level for 22130 genes. The RNAseq dataset was generated on NovaSeq with 100bp paired end reads across 12 lanes. After quality control, gene expression was represented as TMM (trimmed mean of M-value)-normalized log-transformed counts per million for 20416 features.

The SRS transcriptomic endotypes have been assigned based on these datasets as well as for 7 samples with qPCR data using a 7-gene set random forest model previously described (Cano-Gamez et al., 2022). For the small number (11) of border-line samples that had different assignments from the platforms, assignments were taken in the order of priority as microarray before RNAseq before qPCR. Since there was only a small number (48) of patient-timepoints assigned to SRS3 (transcriptionally closer to healthy volunteers) and that SRS3 patients are more similar to SRS2, SRS3 was combined with SRS2 as a ”non-SRS1” category in the analysis of this study.

Protein and sample filtering

In order to retain proteins that reach a detectable level in only certain biological conditions, protein detection was counted within each of the biological groups. Among the cohorts, the GAinS study is the majority group constituting >70% samples. Filtering based on the heterogeneous syndrome of sepsis as one group could potentially miss out the proteins with distinct patterns between unknown biological subgroups. Therefore, the GAinS samples were further separated by unsupervised consensus clustering on a matrix of binary states detected (1) or missing (0), which identified four GAinS protein detection subgroups, the smallest of which comprised 292 samples. Protein detection was then counted within nine groups including: the healthy volunteers of OBB; the pre-operation samples of X-MINS and BIONIC; the post-operation samples of X-MINS and BIONIC; the no-infection ICU samples of MOTION, MONOGRAM and TACE; the septic shock samples of VANISH; and the four GAinS detection groups (Figure S1B). We kept 291 proteins that were detected in at least 50% sample injections in at least one of the nine groups.

Sample quality represented by sample-wise protein detection numbers showed variation corresponding to the changes of the chromatographic column. We excluded 32 sample injections for having outlying low detections (<125 proteins), the threshold of which was determined by plotting the detection numbers (Figure S1C). The raw intensity data matrix post-filtering contained quantifications for 291 proteins in 2615 sample injections.

Correction for cell residue contaminations

After filtering, the data was normalized, missing values imputed and batches corrected to derive an “interim processed data” version, as shown in the flowchart (Figure S1A). Although all plasma samples were extracted by the standard approach of centrifuging whole blood taken with anticoagulants, samples from different studies recruited in a variety of clinical settings could potentially have differences in pre-analytical handling. Even within a single study following the same protocol, in the interim processed data the GAinS samples, recruited over 14 years in 34 ICUs, showed a small subset of ∼12% samples that were significantly higher in proteins including actins and actin-binding proteins that point to higher platelet residues in plasma. Therefore, it is necessary to assess and account for the impact of the technical effects.

We calculated three contamination indices (Figure S1D) to assess plasma sample contamination by erythrocyte lysis, platelet contamination, or partial coagulation, based on the protein marker panels proposed by Geyer et al. (Geyer et al., 2019). We calculated and considered the contamination indices using both the log-transformed raw intensities and the interim processed data as both have uncorrected biases. X-MINS samples were collected using the strictest protocol (higher speed in centrifugation and filtered through 0.8μm membrane filters) and were among the samples with the lowest level and least variation in the cell residues (erythrocytes index and platelets index) (Figure S1D). As the coagulation index was less variable across the cohorts, and as sepsis pathophysiology involves dysregulated thrombosis, the difference in coagulation index was considered to be biological signal and not corrected for.

We corrected for cell residue contamination by identifying and removing the most affected proteins species from the dataset. In addition to proteins on the marker panels, we also removed proteins that correlated with either the platelet or erythrocyte index in both the raw and interim processed data (Pearson’s r > 0.20). In total 22 proteins were removed, leaving 269 to be further processed (Figure S1E). The potentially contaminated samples did not need to be excluded, since these were not outlying in protein detection numbers or on PCA after excluding the contamination proteins, and the removed proteins did not compromise the detection of other proteins.

Normalization and imputation

We used Variance Stabilizing Normalization (VSN) (Huber et al., 2002) to normalize the raw intensities and account for systematic bias. This method was selected based on running a smaller sepsis pilot dataset through Normalyzer (Chawade et al., 2014) which compared different normalization algorithms, as well as based on others’ work showing the advantage of VSN in addressing the variance heterogeneity in mass spectrometry-based proteomics (Karp et al., 2010; Välikangas et al., 2016).

Assuming proteins are missing not completely at random, we applied a hybrid imputation approach. For 170 proteins detected in ≥60% of the samples, we applied the K-Nearest Neighbors (KNN) algorithm which is more tolerant on the assumption of the causes of missingness. For each protein with missing values to be imputed, KNN finds the k (k=6) most similar proteins (neighbors) based on the non-missing values, and fills in the missing value by the mean observed for the neighbor proteins in this sample (Hastie et al., 2021). When there is a higher proportion of missingness in a protein, KNN is not suitable since neighbor proteins cannot be accurately located. For the 99 proteins detected in <60% of the samples, missing values were imputed by random draw from a protein-specific normal distribution with a standard deviation (s.d.) of 30% s.d. of the measured values and a down-shift of the mean by 1.8 s.d.

Batch correction

The 28 randomized sample plates were profiled in a single mass spectrometry experiment but the long acquisition duration introduced technical variation. We performed batch correction across the plates using the ComBat function from the sva package (Leek et al., 2021). Seventeen batches were determined based on the pattern of PC1 scores of the plates in uncorrected data and the spectral yields and detection numbers in spiked-in pool injections repeated along the plates. Seven plates that were the most consistent between each other were used as the reference batch where measurements were not altered.

Discovery and validation cohorts

For downstream analysis, we defined the discovery and validation cohort sets from the sepsis and control cohorts (Figure 1A). The GAinS and OBB patients were each separated into a discovery and a validation cohort with a split of 2:1 by random draws after pre-processing. Serial samples from the same patient were kept together.

Exploratory analysis

We performed principal component analysis on all samples in the pre-processed data using the prcomp R function without scaling. We explored the relationship between the first four PCs or individual proteins with sepsis severity by testing Spearman correlation with the clinical variables, restricted to the first available samples of GAinS patients. Clinical variables detected in less than 30% of patients were excluded. Variables only available in certain groups of the patients (e.g. estimated days from CAP/FP onset) were retained. All statistical analysis were performed in R (R Core Team, 2015).

Statistical tests

Differential abundance analysis for proteins and differential expression analysis for genes were performed by fitting the intensities in linear models using the limma R package (Ritchie et al., 2015), using only the first available sample of each patient and including age and sex as covariates. The Benjamini-Hochberg procedure was used to adjust for multiple testing. Significance for downstream analysis was defined as FDR < 0.05 and fold change (FC) > 1.5 unless otherwise specified. Comparisons between the pre- and post- operative samples are paired and with no additional covariates. All tests are two-sided.

Welch’s t-test, also known as unequal variances t-test, is more conservative when the group size and variance differs by a large amount between the two groups. The differential abundance analysis between GAinS and MONOGRAM+TACE was also performed using Welch’s t-test and highly consistent results were obtained as when using limma (Pearson’s r for log₂FC > 0.99).

Cytokine concentrations measured by Luminex assay were compared by Wilcoxon rank-sum tests (i.e. Mann-Whitney tests) using the first available sample of each patient.

Categorical clinical variables were compared using Chi-squared tests without Yate’s correction. Numerical clinical variables were compared between two groups using Wilcoxon rank-sum tests. For comparing numerical clinical variables between three groups, Kruskal–Wallis test (i.e. one-away analysis of variance on ranks) was used to determine whether there was any significant difference (FDR<0.05) between the groups. For variables where the null hypothesis in Kruskal–Wallis tests were rejected, Dunn’s post-hoc tests were performed to compare between each pair of the three groups, the significance levels of which were labelled on boxplots after adjusting for multiple testing within the variable. SOFA scores and ARDS levels were considered as numerical variables. Only the first available samples following ICU admission of each patients were included in comparing the total of 66 clinical variables.

Survival differences were assessed and Kaplan-Meier curves plotted using the R packages survival and survminer (Kassambara et al., 2021; Therneau and Grambsch, 2000). The input data is a data frame specifying the time to event from the sampling day, the event (death or end of 28-day or 6-month observation) and the patient subgroups. For patients with multiple timepoints sampled, cluster assignment from the last available sample within five days of ICU admission was used. The p values were given by log-rank tests. Hazard ratios and the confidence intervals were calculated using univariate (SPC) or multivariate (SPC and SRS) Cox proportional hazard models.

Data visualization

On all boxplots, hinges show the first and third quartiles on both sides of the median, and the whiskers extend to the highest or lowest values within 1.5*IQR of the hinges. Data points are plotted on horizontally jittered positions on the boxplots. For sepsis patients with serial sample points, only the first available samples following ICU admission were plotted.

On all volcano plots, proteins more abundant in Group A of the contrast “Group A vs B” are plotted on the right-hand side of the plot. Differentially abundant proteins (FDR<0.05 and |FC|>1.5) are plotted in red points.

Trajectories of patients’ movements between the proteomic clusters were visualized using the ggalluvial package (Brunson, 2020), restricted to patients with a baseline (Day1) sample available. Widths of the flows on the y axis are scaled to the number of patients in each movement type. Figure 5F and a subset of the icons in Figure 1A “Downstream analysis” panel were generated with BioRender.com.

Unsupervised clustering

We applied consensus clustering to generate unsupervised clusters in 1236 discovery cohort GAinS samples, using the ConsensusClusterPlus R package (Wilkerson and Hayes, 2010). Compared with a single-run hierarchical clustering performed on all available samples and features, the consensus clustering approach has the advantage of cluster results being robust to small variations in sample or feature composition. For each cluster number k tested from two to ten, 500 iterations were performed, randomly sampling 80% of samples and 90% of proteins in each iteration. In each iteration, unsupervised hierarchical clustering was performed based on Euclidean distance for dissimilarity between samples and Ward’s method as linkage for cluster agglomeration.

The cumulative distribution function (CDF) curves (Figure 4C) were inspected to assess cluster stability along an increasing number of clusters (k). The increases in area under the CDF were the largest when k increased from 2 to 3, and reached an elbow point at k=4 (Figure 4D). We also assessed how model accuracy improved as cluster number increased by calculating the total within-cluster variation in 10-fold cross-validation with k-means clustering. The fit of test points to the nearest centers improved gradually as the cluster numbers increased from 1 to 10, and the step widths were the largest from k=1 to k=3. Therefore, k=3 was chosen as the optimal number of partitioning for this dataset, balancing the need for explaining most of the variability and for the simplicity of describing the major patient subgroups.

Cluster prediction models

To validate the three patient clusters (SPC) in independent samples, cluster prediction models were built based on SPC assignments in GAinS discovery cohort samples, and the model with the highest accuracy was applied to the GAinS validation cohort samples. In evaluating the prediction model performances, accuracy is the sum of the diagonal of the contingency table divided by the total number. A random 3-cluster classification keeping the same proportions in each cluster as in the discovery cohort would have an expected accuracy of (21%)²+(24%)²+(55%)² = 40%. Sensitivity for a particular event is the proportion of events correctly predicted, out of the number of true events. The best-performing model was selected based on the test set accuracy.

We tested eight models in total for SPC prediction using three statistical learning methods: partial least squares discriminative analysis (PLS-DA), generalized linear models (GLM), and class prediction by random forest through the R packages ropls, glmnet, caret and randomForest (Friedman et al., 2010; Kuhn, 2021; Liaw and Wiener, 2002; Thévenot et al., 2015). Input candidates are either the total 269 proteins in pre-processed data, or restricted to the 138 proteins that were differentially abundant between any pair of the 3 discovery cohort clusters at a less stringent threshold (FDR<0.05 and |FC|>1.2), or the 41 proteins at a larger fold change (FDR<0.05 and |FC|>1.5).

The numbers of predictive components to include in PLS-DA models were determined by calculating the proportion of variation explained in 7-fold cross-validation and evaluating model significance with label permutations GLM models were tested with lasso, ridge, and elastic net regressions to include the overfitting penalty and enable variable selection. Symmetric multinomial models were fit so that each of the 3 SPC classes were represented by a GLM. Two parameters were optimized: the tuning parameter λ which controls the overall strength of the penalty by defining the amount of coefficient shrinkage; and the elastic net mixing parameter α which controls whether the penalty behaves more towards a lasso (α=1) or ridge (α=0) regression for correlated predictors.

For lasso and ridge regressions, the value of λ was determined by minimizing the prediction error rate in 10-fold cross validation (CV) within the training set. When there was not a great increase in CV misclassification rate between the λ that gives the smallest error, and the largest λ such that the error is within the smallest error plus one standard error, the latter value of λ was chosen so that the model was less heavily based on the training set. For elastic net regression (0<α<1), the best values for λ and α were selected by testing a grid of 20*20 combinations, and selecting the combination that gave the highest prediction accuracy in cross validation.

Since the elastic net model with 269 protein candidates had the best performance in GLM models, 269 proteins were used as candidate inputs to train the random forest model. The 3-cluster accuracy from 10-fold cross validation was used as the optimizing metric for tuning the parameters. In the optimized random forest model, 35 variables were randomly sampled for each iteration; the maximum number of terminal nodes of each tree was 42; 350 trees were trained; the other parameters were kept as default values in the package.

The 8 protein predictors in the minimal elastic net model were selected by restricting to 3 proteins with the strongest signals in the differential abundance analyses between each cluster pair, whilst also ensuring at least one higher- and one lower- abundance protein were retained for each pairwise contrast.

Matrix decomposition

We performed matrix decomposition on the 837 samples with both RNAseq (for 20,416 genes) and plasma proteome data (for 269 proteins), using the SDA method (Hore et al., 2016), with the same number of iterations (10), clustering and cutoffs to end up with 284 latent components. Within the components, we also used the same cutoff for the posterior inclusion probability (0.5) as (Hore et al., 2016) for the inclusion of the gene or protein loading scores. This resulted in 76 components with protein loading scores. Similarly to our work in (COvid-19 Multi-omics Blood ATlas (COMBAT) Consortium, 2022) we matched these to clinical features as well as the SRS or SPC using either Spearman’s correlation coefficient or Mann-Whitney U rank test, adjusting for multiple testing using the Benjamini/Hochberg (non-negative) method. Correlation of the latent components with SRS were tested with the quantitative SRS scores (Cano-Gamez et al., 2022). Z-scores of representative proteins shown in Figure 6 were calculated within the first available samples used in matrix decomposition analysis.