Development of a metabolome-based respiratory infection prognostic during COVID-19 arrival

Between 1 March 2020 and 30 September 2020, 163 patients presenting to the Barnes-Jewish Hospital emergency department with symptoms compatible with COVID-19 or at high risk of developing COVID-19 were enrolled (Fig. 1). Among the 150 SARS-CoV-2-positive patients in our cohorts, 45% (68) reached the pre-defined, severe COVID-19 composite endpoint of intubation or death by 90 days from admission, with all but 2 within 30 days. Of these, 82% (56) was intubated, with 50% (11) mortality post-intubation. Of the 40 COVID-19 patients who died during the study period, death was attributed to COVID-19 in 88%. Cause of death was not available for two COVID-19 patients. Patients with severe COVID-19 were significantly older (median 68.4 vs. 61.2 years, P < 0.001) with non-significantly increased diagnoses of diabetes or chronic kidney disease (36.8% vs. 25.6% and 7.4% vs. 2.4%, respectively; Table 1).

Discovery cohort enrollment ended after 4 months (1 March 2020 to 30 June 2020) with 96 patients. Subsequently, 67 patients were recruited (1 July 2020 to 30 September 2020) to the validation cohort. No significant differences in baseline clinical characteristics were observed between cohorts. When 90-day endpoints were known, we divided each cohort into severe and non-severe patients (39 and 57 for discovery, 30 and 37 for validation). Among severe COVID-19 patients, 57% (39/68) was intubated before urine collection while 50% of their urine specimens was collected between 2.9 days before and 3.5 days after reaching the severity-defining endpoint. Restricting to patients with positive SARS-CoV-2 PCR tests later reduced the discovery cohort to 87 patients (38 severe, 49 non-severe) and the validation cohort to 63 patients (30 severe, 33 non-severe).

To identify metabolomic variation associated with progression to severe COVID-19, we obtained full-scan LC-MS profiles of the 96 discovery cohort samples (including PCR-negative samples) in positive- and negative-ion modes. Both unsupervised hierarchical clustering and principal component analysis (PCA) of low-sparsity metabolite logratios could partially discriminate severe and non-severe COVID-19 patient groups (Fig. S17 and S18). To better identify correlates of progression to severe COVID-19, we next applied multiple feature filtering and bounded feature selection techniques (see Materials and Methods). Models combining both logratios and clinical characteristics generally outperformed models trained with either data type alone (Fig. S15). Models including positive ion and/or negative ion-derived logratios performed comparably. Using positive-ion mode logratios and clinical characteristics with a maximum of 70 metabolites in our bounded feature selection yielded models with a mean area under the curve (AUC) of 0.801 from, on average, 20.6 metabolites. The Hamming distance between optimized feature subsets defined a single most representative set of 24 metabolites and 4 clinical variables (age, presence of stroke/dementia/seizure disorder, confusion, and fatigue) for further biomarker development. When sample collection for the discovery cohort began, patients with both confirmed and suspected SARS-CoV-2 infection were enrolled due to limitations in SARS-CoV-2 diagnostic performance and availability. Early steps of feature selection were performed on this combined discovery cohort, which contained nine PCR-negative patients (with one severe outcome). We retained these PCR-negative patients to prevent outcome information from biasing our already-performed bounded feature selection. Subsequent LC-MS/MS detection and classifier training (below) on the discovery cohort utilized only the 87 PCR-positive discovery cohort patients.

We next sought to identify a simplified metabolomic signature of progression to severe COVID-19. To quantify the 24 candidate metabolites with greater sensitivity, specificity, and precision than the full-scan LC-MS mode used for metabolomic discovery, we re-measured each using LC-MS/MS detection (Table S2). Of the 24 metabolites, 5 (21%) were detected in all specimens and 16 (67%) were detected in more than half. As detailed above, we used the DiCoVarML R package (12, 13) to train a model in terms of logratios between five metabolites, yielding the highest average training set AUCs among models trained with four to eight metabolites (Fig. S16).

Because ensemble models offered only modestly higher training AUCs, we chose to use more-interpretable ridge-regularized logistic regression (Fig. S16), yielding a model defined by 6 logratios between five metabolites (labeled by m/z and retention time as [email protected], [email protected], [email protected], [email protected], and [email protected]). Algebraic simplification from 6 logratios to the five underlying (log-transformed) metabolites yields regression coefficients for our final scoring function (Table 2). Model predictions are shown in Fig. 2A. The final model yielded an AUC of 89.4% (95% CI: 82.3–95.1) on the discovery cohort. The logratio derivation ensures the signature is inherently independent of absolute quantities and remains useful if sample dilution and/or instrument responses differ from those of the present study, provided the relative quantities of these five metabolites are accurately obtained.

To evaluate its prospective external validity, we applied our severe COVID-19 signature model to the validation cohort. Instrumental analyses and data interpretation were conducted by operators blinded to the outcomes associated with each specimen. Upon unblinding, the validation cohort AUC for the severe COVID-19 signature (81.2%; 95% CI: 68.4–91.1) was significantly better than the null hypothesis (50%) (Fig. 2B) and was not significantly lower than in the discovery cohort (P > 0.05, bootstrap test; Fig. 3) (14). We observed similar distributions (P > 0.05, Kolmogorov-Smirnov test) of outcome classification probabilities in discovery and validation cohorts (Fig. 2A). Univariate analysis of individual logratios also revealed similar trends between discovery and validation cohorts (Fig. 4B; Fig. S19). Using a standard 50% probability cutoff to assign regression classifications yielded similar accuracy on the discovery and validation cohorts (81.6% vs. 74.6%; Table S3). Together, these results indicate minimal model overfitting and strong beyond-sample generalizability of the severe COVID-19 metabolomic signature.

A network representation (Fig. 4A) of the logratios between the five validated features distinguishes their roles as numerators (3 features) and denominators (2 features) in the model. The two denominator features were present in all specimens and are negatively associated with severe COVID-19 (Wilcoxon test, P < 0.01; Fig. 4B ), although the relative magnitudes of these associations are small. The denominators thus effectively serve as “housekeeping” features, providing ubiquitous physiological references for metabolite normalization, analogous to creatinine. In contrast, each of the three numerator features exhibits bimodal distributions, consistent with three distinctive severity-associated metabotypes (Fig. 4B) (15). Metabotype 1 is defined by the presence of [email protected], metabotype 2 by [email protected], and metabotype 3 by [email protected]. Consistent with the manner of their selection, these metabotypes were not completely disjointed—subjects may simultaneously test positive for one or more metabotypes (Fig. 4C and D). Interestingly, metabotype 2 was associated with more advanced age (Wilcoxon test, P < 0.01; Fig. S20). No metabotype was associated with known clinical risk factors of diabetes or kidney disease (P > 0.05). These results are consistent with three independent biomarkers of COVID-19 severity with distinctive biochemical origins.

The potential for a metabotype to guide early antiviral interventions and patient disposition is influenced by its temporal relationship to severe disease, but this has not been typically considered in COVID-19 biomarker studies. Among the three metabotypes, the metabotype 3 signal was greater in specimens collected after intubation (Wilcoxon test, P < 0.01; Fig. 5). Metabotype 3 was neither uniformly absent in pre-intubation subjects nor uniformly present in post-intubation subjects, consistent with a physiological process that is closely related to acute respiratory failure. The presence of metabotypes 1 and 2, in contrast, was independent of intubation, consistent with baseline correlates of COVID-19 severity.

COVID-19-associated metabotypes present early in infection, before the full inflammatory phase of illness, could potentially affect decisions about patient disposition and antiviral treatment. To evaluate their potential in this role, we fit predictor models including metabotypes 1 and 2 with or without baseline clinical predictor variables (diabetes, kidney disease, and age). The model using only metabotypes 1 and 2 was not significantly inferior to the model using all three metabotypes (Fig. 3, validation AUC 74.2% vs 81.2%; P > 0.05, bootstrap test with Bonferroni correction), consistent with significant predictive power in the absence of metabotype 3 (14). Adding clinical predictors non-significantly lowered the validation AUC versus metabotypes 1 and 2 alone. Notably, the validation cohort performance of clinical predictors alone was not significantly different from chance. These results suggest that metabotypes 1 and 2 are more consistently associated with severe COVID-19 than baseline clinical predictors.

Additional information about COVID-19 progression may be obtained by chemically identifying the metabolites in the severe COVID-19 signature. Multistage fragmentation (MSⁿ) coupled with high-resolution accurate mass analysis of the five metabolites comprising our signature revealed they are oxyhydrocarbons, three of which possess a single nitrogen and no other heteroatoms. The two denominator metabolites in our model were identified with high confidence as an O-acylcarnitine ([email protected]) and δ-valerolactam ([email protected]) (Table 3). The fragmentation pattern of feature [email protected] contains multiple diagnostic fragments observed in other acylcarnitines (Fig. S2 to S5), with differences attributable to the presence of an unsaturated nine-carbon fatty acyl group. Feature [email protected] was definitively confirmed by comparison to authentic δ-valerolactam, an established urinary metabolite generated from cadaverine and subject to cytochrome P450-dependent metabolism (Fig. S6) (16).

The metabolite-defining metabotype 1, [email protected], was identified with high confidence as an electrospray ionization source-decay fragment of a previously-observed deaminated metabolite of the anticonvulsant drug levetiracetam (Keppra) (Fig. S7 to S10) (17). Its MS2 spectrum exactly matched a reference spectrum of levetiracetam in the mzCloud spectral library, while its molecular weight matches the deaminated metabolite. Concomitant detection of unmetabolized levetiracetam in specimens with this feature further increased confidence in this assignment (Fig. S10). This result is consistent with a positive association between levetiracetam use and progression to severe COVID-19.

The metabotype 2 feature, [email protected], has predicted molecular formula C₁₂H₁₆O (Table 3) and shares multiple MS2 fragments with the fecal metabolite 1-phenyl-2-hexanone (Fig. S11A), including a tropylium ion characteristic of carbon-substituted benzene rings. Moreover, the 131.0850 m/z MS2 fragment ion has a neutral loss of C₆H₆ (consistent with benzene) in its MS3 spectrum (Fig. S11B). An H₂O neutral loss in the MS2 spectrum also suggests the presence of an alcohol. Together, these data are consistent with [email protected] as a monounsaturated aliphatic alcohol-substituted benzene. The specific biochemical origin of this metabolite is unclear, though the presence of a benzene ring may implicate phenylalanine or tyrosine metabolism.

The metabotype 3 feature, [email protected], with predicted molecular formula C₁₈H₂₂O₅, yielded no direct match in online spectral databases. However, it shares distinctive fragments and neutral losses with enterodiol (Fig. S12 to S14), the eponymous member of an established group of human urinary metabolites and a product of lignan metabolism by the intestinal microbiota. The predicted formula (Table 3) matches hydroxyenterodiol, which may derive from P450-mediated hydroxylation of enterodiol in mammalian hosts (18, 19). These data are consistent with [email protected] as a member of the enterodiol family of bacterial-human co-metabolites derived from common dietary compounds. Elevation of this metabolite after respiratory failure in a subset of patients with severe COVID-19 may reflect increased systemic absorption associated with increased gut barrier permeability and/or increased human P450 metabolism of xenometabolites.

Patients presenting to the Barnes-Jewish Hospital emergency department with symptoms compatible with COVID-19 or at high risk of developing COVID-19 were enrolled upon submission of a nasopharyngeal swab for SARS-CoV-2 RNA testing using a Food and Drug Administration-approved clinical PCR test (Fig. S1). The primary endpoint for progression of COVID-19 to severe disease was respiratory failure requiring mechanical ventilation and/or death within 90 days of study enrollment. Subjects enrolled from 1 March 2020 to 30 June 2020 comprised the initial discovery data set (n = 96, SARS-CoV-2 positive n = 87). Additional subjects enrolled between 1 July 2020 and 30 September 2020 comprised the validation data set (n = 67, SARS-CoV-2 positive n = 63) (Table 1). The workflow for discovery and validation cohorts is outlined in Fig. S1. Importantly, no information about the validation cohort was used in identifying and formulating our candidate metabolomic signature. Clinical data were extracted from electronic health records by the Washington University Institute of Informatics and provided in a de-identified manner.

Patient urine samples were collected in a sterile cup by the hospital staff as convenient and transported twice a day to the laboratory. Samples were stored at 4°C until ready for processing on the same day. All samples were processed and stored in less than 24 hours. Urine was filtered through a 40-µm mesh followed by 2,000 × g centrifugation for 12 minutes. Supernatant was immediately aliquoted into cryovials and stored at −80°C until ready for experimental analysis.

Urinary metabolites were measured using an AB Sciex 4000 QTrap triple-quadrupole mass spectrometer interfaced with a Shimadzu HPLC as previously described (36, 37), (see supplemental material for details). Reversed-phase chromatography was conducted with an Ascentis-Express fused core phenyl-hexyl column (Millipore Sigma; 100 mm × 2 mm × 2.7 µm). Metabolite profiling was conducted in full-scan mode across a mass/charge range of 50–1,200 m/z. Multiplexed profiling of selected metabolites was conducted in MS/MS mode. Structural characterization of validated metabolites was performed using an Orbitrap ID-X high-resolution accurate mass spectrometer coupled to a Vanquish ultra-high pressure liquid chromatograph (UHPLC; Thermo Scientific), as described in the supplemental materials (Fig. S2 to S14).

To focus our search on broadly applicable features, we restricted analyses to the 103 negative-ion and 163 positive-ion features present in >50% of samples in at least one clinical outcome group. To control for physiologic differences in urinary dilution, we analyzed urinary metabolomic profiles as log-transformed metabolite peak area ratios (hereafter designated “logratios”) and selected from all possible ratios in the data during subsequent feature selection steps (11, 29). Using metabolite ratios allowed us to identify urinary reference metabolites without making a priori assumptions, for example, we did not assume that any single metabolite, such as creatinine, would provide an optimal physiologic reference.

To select a reduced set of LC-MS metabolite features for targeted LC-MS/MS analysis, we used a feature selection approach based in genetic algorithms (for additional details, see supplemental materials). To estimate the ability of metabolomic and/or clinical data to correctly classify severe outcomes and to minimize overfitting, we performed this process with fivefold cross-validation train/test splits of the discovery cohort, repeated 20 times with different random seeds, for each combination of metabolomic data and clinical metadata. After a subset of metabolites was identified using the genetic algorithm approach, we used differential compositional variation scores to select subsets of logratios for model training through multistage feature selection using network-based, recursive elimination combined with penalized regression techniques (28). This broad approach was used to improve model performance and improve interpretability by eliminating uninformative logratios. Clinical metadata features were optionally selected using the Boruta algorithm on individual training partitions (38). Features selected in 8 of 10 repeats of the Boruta algorithm were retained. The resulting clinical metadata subset was then appended to the selected metabolite subset. Finally, the XGBoost algorithm was used to train tree booster models using the selected features (39).

Across all test splits, average model AUC was superior for models combining metabolite and clinical data. Models using positive-mode data, negative-mode data, or both modes yielded similar average model AUCs, with larger models generally exhibiting higher AUCs (Fig. S15). Because positive mode ionization is more common and sensitive, we proceeded with a model combining positive mode features and clinical metadata. To select a smaller set of representative metabolites and clinical attributes, we repeated the feature selection process 50 times with LC-MS features and clinical attributes, with a maximum model size of 70 variables. We aggregated the importance scores for each selected variable across all models and, representing each subset as a binary vector, selected the subset with minimum total pairwise Hamming distances as the most representative (Fig. S22). The 24 metabolites in this most representative subset were selected for targeted metabolomic analysis.

To simplify the 24-metabolite signature while maintaining classification performance, we used the targeted differential compositional variation machine learning framework implemented in the DiCoVarML R package (13). DiCoVarML is a multilevel nested framework that respects the variable dilution of urinary specimens by leveraging network methods that identify and robustly estimate classification performance of minimal logratio signatures (13). Pairwise logratios between the 24 metabolites were calculated with zeroes imputed using multiplicative replacement (28, 40). We selected ridge regression over ensemble models for easier interpretability and because we did not see substantial decrease in accuracy (Fig. S16). As in the previous stages of feature selection and model building, only discovery cohort data were used as input. Out-of-sample classification performance was estimated using five repeats of twofold cross-validation for each target subset size (four to eight metabolites). Multiplicative zero imputation was performed separately for each cross-validation step to avoid leaking information between folds. The DiCoVarML nested classification performance estimate was then used to select the final subset of five metabolites (Fig. S16) and train the final ridge regression model with the corresponding pairwise logratios. For ease of interpretation, the model in terms of 6 logratios was algebraically simplified to 5 log-scaled metabolite signal values to specify our final prognostic function (coefficients in Table 2). This final prognostic function was applied to LC-MS/MS data to predict outcome severity in the validation cohort.

All statistical analysis and machine learning model building and assessment were done in R version 4.1.0. The R mixOmics package was used for PCA (41). We used the R glmnet package to predict outcome probabilities and labels using the final model trained on the discovery data (42). Discovery and validation set AUCs and 95% confidence intervals were computed using the R pROC package (14). ROC AUCs were compared using bootstrapping (R pROC package, roc.test function). Comparison P values were adjusted for multiple comparisons by the Bonferroni method (R stats package, p.adjust function). The prognostic function scores, using the discovery-trained model, were computed for both validation and discovery data using the R predict function with response output. The Kolmogorov-Smirnov test (R stats package, function ks.test) was used to compare distributions of predicted probabilities of severe outcome between discovery and validation cohorts, stratified by outcome. The two-sample Wilcoxon signed-rank test (R stats package, function wilcox.test) was used to describe differences in means of logratio features between outcome groups.