Large-scale Multi-omic Analysis of COVID-19 Severity

Sample cohort and experimental design

From 6 April 2020 through 1 May 2020, we collected blood samples from 128 adult patients admitted to Albany Medical Center in Albany, NY for moderate to severe respiratory issues presumably related to infection with SARS-CoV-2. In addition to collection of various clinical data (Table 1), a blood sample was obtained at the time of enrollment. Patients who tested positive (n = 102) and negative (n = 26) for the virus were assigned to the COVID-19 and non-COVID-19 groups, respectively (see Methods for enrollment details). Females comprised 37.3% and 50.0% of the COVID-19 and non-COVID-19 groups, respectively. The average age of patients was similar – 61.3 and 63.1 years in the COVID-19 group (females and males, respectively; p-value = 0.56) – with some statistically insignificant differences between females and males in the non-COVID-19 group (59.5 and 67.0 years, respectively; p-value = 0.25) (Figure 1a). The average number of days hospitalized before study enrollment was 3 and 1 for the COVID-19 and non-COVID-19 groups, respectively (Table 1). The COVID-19 group was more racially diverse than the non-COVID-19 group, with white individuals comprising only 46% of the total (vs. 80% of the non-COVID-19 control group); these demographics are consistent with the reported racial and ethnic health disparities of COVID-19 (Webb Hooper et al., 2020).

• Download figure
• Open in new tab

Figure 1. Overview of sample cohort and experimental design.

a Age and sex distributions of COVID-19 (n = 102) and non-COVID-19 (n = 26) groups. b Distributions of hospital-free days over a continuous 45-day period aggregated with survival (HFD-45, see outcomes selection in the methods section) among COVID-19 and non-COVID-19 groups. c The proportion (%) of female and male patients who were admitted to the intensive care unit (ICU) and required support of a mechanical ventilator. d Overview of the study design, experimental approaches, and primary outcomes. Notice that the leukocytes were separated by filtering (see methods for details).

Severe COVID-19 infections require hospitalization and can lead to death due to respiratory deterioration. Thus, patients with the most severe cases requiring ventilatory support in the ICU entail longer admissions, while those with the most extreme cases die during hospitalization. To integrate length of hospitalization stay with mortality in one single outcome measure reflecting disease burden, we scored hospital-free days at day 45 (HFD-45). This composite outcome variable assigns a zero value (0-free days) to patients who remain admitted longer than 45 days or die while admitted and higher values of HFD-45 to patients with shorter hospitalizations and milder disease severity. Similar composite outcome variables are frequently used in the ICU setting (Jaitovich et al., 2019; Young et al., 2015). Other clinically relevant information was obtained, including: Acute Physiologic Assessment and Chronic Health Evaluation (APACHE II) score, Sequential Organ Failure Assessment (SOFA) score, Charlson Comorbidity Index score, the number of days spent on mechanical ventilators, need for admission into the intensive care unit (ICU), and laboratory measurements of C-reactive protein (CRP), D-dimer, ferritin, lactate, procalcitonin, fibrinogen, and others (Table 1). Members of the non-COVID-19 group had on average considerably more HFD-45 than those of the COVID-19 group (32.3 vs. 22.0, respectively, p-value = 0.004; Figure 1b). APACHE II and SOFA severity scores, used to stratify severity of critical illness, were assigned to ICU-requiring patients, and both metrics exhibited similar distributions between the groups (Table 1). Consistent with previous reports, we observed male sex predominance in the ICU-requiring COVID-19 group (53.1 vs. 44.7%) and the group requiring mechanical ventilatory support (46.9 vs. 34.2%) (Figure 1c).

This comprehensive clinical characterization of the cohort provided substantial data to be integrated with multi-omic measurements. Following collection, blood samples were centrifuged to separate their components for various omic analyses (Figure 1d). We performed four individual mass spectrometry (MS)-based assays on the plasma components using high-resolution, high mass accuracy MS coupled with either gas or liquid chromatography: shotgun proteomics, discovery lipidomics, discovery metabolomics, and targeted metabolomics. RNA-Seq was used to characterize the transcriptomes of leukocytes extracted from the patient samples.

Collectively, these data enabled a comprehensive systems analysis of COVID-19 blood samples – containing abundance measurements of over 17,000 diverse biomolecules (517 proteins, 13,263 transcripts, 646 lipids, 110 metabolites, and 2,786 unidentified small molecule features). We compiled all biomolecule abundance measurements and de-identified patient information with clinical data into a highly curated relational (SQLite) database, which is publicly available for exploration and further analysis (see Data Availability, Supplementary Figure 1). The database includes additional metadata for each molecular measurement, such as gene ontology (GO) terms, alternative identifiers, and analytical metrics used to filter data, among others (see Methods). To facilitate exploration of this dataset by the broader scientific community, we also created a companion web tool (https://covid-omics.app) that accesses data and various analyses and performs on-demand data visualization (Figure 1d). Detailed description of this tool is presented in a later section.

Systems analysis reveals strong biomolecule associations with COVID-19 status and severity

We hypothesized that measurement of a large number of biomolecules across different molecular classes (i.e., nucleic acids, proteins, lipids, and metabolites) would provide insight into the COVID-19 molecular landscape and facilitate understanding of factors associated with higher severity. To test this hypothesis in an unsupervised manner, we initially performed a principal component analysis (PCA) (Figure 2a). Here we note a strong grouping of patient samples based on severity (HFD-45 values) with samples of high and low severity diverging along the diagonal axis of principal components 1 and 2. This, together with a subtler grouping based on status (COVID-19 vs. non-COVID-19), merited further supervised exploration.

• Download figure
• Open in new tab

Figure 2. Multi-omics analysis reveals strong molecular signatures associated with COVID-19 status and severity.

a Principal component analysis using quantitative values from all omics data (leukocyte transcripts, and plasma proteins, lipids, and small molecules) shows principal components 1 and 2 capture 16% and 10% of the variance between patient samples. Plotting samples by these two components show a linear tread with hospital free days at 45 days (HFD-45). b Associations of biomolecules with COVID-19 status was determined using differential expression analysis (EBseq) for transcripts, and linear regression log-likelihood tests for plasma biomolecules, the adjusted p-values (1 - posterior probability or Benjamini Hochberg-adjusted p _values, respectively) are plotted relative to the log2 fold change of mean values between COVID and non-COVID samples. In total, 2,537 leukocyte transcripts, 146 plasma proteins, 168 plasma lipids, and 13 plasma metabolites had adjusted p-values < 0.05. c Associations between biomolecules and HFD-45 was estimated using a univariate linear regression (HFD-45 ∼ biomolecule abundance + age + sex) resulting in 7,408 biomolecules. A multivariate linear regression with elastic net penalty was applied to each omics dataset separately to further refine features of interest and resulted in 946 features. In total 219 features were determined as most important for distinguishing COVID status and severity. d The 219 features abundances were visualized via a heat map and clustered with hierarchical clustering. Features that were elevated (e) or reduced (f) with COVID status and severity were used for GO-term and molecular class enrichment analysis.

To gain biological insight into the host’s response to SARS-CoV-2 and pathways influencing its severity, we integrated our biomolecule measurements with clinical outcome variables. For this supervised analysis, we used univariate and multivariate regression to identify features which associate with (1) COVID-19 status and (2) HFD-45. Significant changes in plasma proteins, lipids, and metabolites associated with COVID-19 were determined by ANOVA and log-likelihood ratio tests, which incorporated potentially confounding variables, such as sex, age, and ICU status. Significant changes in leukocyte transcripts associated with COVID-19 diagnosis were analyzed separately by EBSeq. Figure 2b presents the statistical significance relative to the log₂ of the fold changes of means for COVID-19 vs. non-COVID-19 groups for each biomolecule. In total, 2,537 leukocyte transcripts, 146 plasma proteins, 168 plasma lipids, and 13 plasma metabolites were significantly associated with COVID-19 status. Next, we applied GO and molecular class enrichment analysis to these features (Supplementary Table 1). Among the biological processes that were enriched in biomolecules whose abundance was altered with COVID-19 (combined proteins and transcripts), many support recent findings (Messner et al., 2020) (Shen et al., 2020). These GO terms include mitotic cell cycle, phagocytosis recognition, positive recognition of B cell activation, complement activation (classical pathway), and innate immune response. Some of these enriched GO terms were driven by transcripts, for example mitotic cell cycle (38 transcripts vs. 0 proteins up with COVID-19), while other GO terms, like phagocytosis recognition or positive regulation of B cell activation, were driven by both transcripts and proteins (11 transcripts/18 proteins and 10 transcripts/18 proteins, respectively, Supplementary Table 1). Triacylglycerides (TGs) were also enriched and more abundant in COVID-19 samples, while plasmenyl–phosphatidylcholines (PCs) were less abundant in COVID-19 samples compared to non-COVID-19 samples. These results are in line with prior reports of significant changes in lipid composition in COVID-19 sera (Shen et al., 2020).

In addition to the annotated features, 511 unidentified metabolites and lipids were also significantly associated with COVID-19 diagnosis. We performed additional spectral database searching and manually confirmed the identity of several of these features (Supplementary Table 2). Initially, we explored the features with the largest fold change between COVID-19 and non-COVID-19 (Figure 2b). The tandem mass spectrum from the unidentified feature with a mass-to-charge ratio (m/z) of 350.1625, and the largest fold-change in COVID-19 vs. non-COVID-19, contained a fragment ion diagnostic of chloroquine/hydroxychloroquine (m/z 247.10, C₁₄H₁₆ClN₂) suggesting this top feature is a metabolite of the experimental therapeutic hydroxychloroquine (Supplementary Figure 2). The feature with the second largest increase in abundance in COVID-19-positive plasma (Figure 2b) having m/z value of 749.5155 was matched to the antibiotic Azithromycin (Supplementary Figure 2). By referencing patient medical records (Table 1), we determined that 85% of COVID-19 patients received hydroxychloroquine and 80% of COVID-19 patients received azithromycin treatment prior to sample collection (vs. 0% and 8% non-COVID-19 patients, respectively), thus explaining the strong association with COVID-19 status. This ability to detect drug metabolites highlights both the power of this discovery approach and the complexity of these samples, as these biomolecular data provide an archival-quality detailed snapshot of both the patients’ physiological response and exposures to therapies.

One ongoing challenge with COVID-19 disease is the broad and unpredictable range of patient outcomes, i.e., disease severity. To find biomolecules associated with severity, we regressed HFD-45 against abundance of each biomolecule. Using this univariate approach, we found 6,202 transcripts, 189 plasma proteins, 218 plasma lipids, and 35 plasma small molecules which were associated with disease severity. To further refine features of interest, for each molecule-type we performed multivariate linear regression on HFD-45 using the elastic net penalty (Zou and Hastie, 2005) (Methods). The elastic net simultaneously performs regression and feature selection; the elastic net penalty induces sparsity in the fit coefficients, which leads to selection of key features that best predict the response variable. In contrast to other approaches that also incorporate sparsity, such as LASSO, elastic net is known to have a “grouping effect” where it selects groups of correlated, predictive features enabling more informative results in downstream enrichment analyses. With the elastic net approach, 497 transcripts, 382 proteins, 140 lipids, and 60 metabolites were retained as predictive features for the outcome HFD-45. In total, we used the following criteria: (1) significance with COVID-19 status, (2) significance with HFD-45, and (3) elastic net feature selection, and generated a list of 219 features most significantly associated with COVID-19 status and severity (Figure 2c and Supplementary Table 3). Levels of these biomolecules in each sample are visualized in a heatmap in Figure 2d. Hierarchical clustering in the heatmap shows clear grouping by COVID-19 status and severity and reveals clusters of molecules (across omes) with similar trends across patient samples.

Applying GO term and molecular class enrichment analysis to these 219 most differentially abundant features (Figure 2e-f), we found that biomolecules in GO terms for complement activation (classical pathway), antimicrobial humoral response, and neutrophil degranulation were enriched with COVID-19 severity (Supplementary Table 1). The “complement activation” term encompassed many variable segments of immunoglobulins encoded by genes that undergo recombination to allow expression of light or heavy chains; twelve of them (8 IGHVs, 2 IGHVs, and 2 IGLVs) correlated with and predicted lower HFD-45 (Supplementary Table 3). In general, the complement system is involved in immune response to virus infection (Stoermer and Morrison, 2011) and has been postulated as a target for therapeutics (Risitano et al., 2020). Our observations about immunoglobulins offer a glimpse into the early humoral response to SARS-CoV-2 and the relationship between preferential variable segment use and individual response trajectories. To better understand production of protective or possibly injurious antibodies, serial studies and characterization of the properties of individual antibodies is required (Ju et al., 2020; Seydoux et al., 2020). Another top feature was pulmonary surfactant-associated protein B (SFTPB), notable because circulating SFTPB correlates with decreased lung function in smokers (Leung et al., 2015), and may be a surrogate marker of lung deterioration in COVID-19 individuals.

Features that were reduced in COVID-19 samples and whose abundance was lowest in the most severe cases were enriched in the following categories: response to folic acid, PCs, plasmenyl-PCs, and high density lipoprotein particle remodeling. These features suggest a significant change in plasma lipid regulation. Similar to previous reports, we observed lower levels of APOA1, APOA2, and APOM in COVID-19 patients relative to non-COVID-19 patients (Messner et al., 2020); (Shen et al., 2020).

Cross-ome correlation analysis

A major strength of multi-omic approaches is the ability to uncover functional connections between different classes of biomolecules. This effort is particularly important for small molecules, especially lipids, whose functions are often less well annotated. To reveal these connections, we integrated all omic measurements acquired on plasma samples by performing cross-ome correlation analysis. We calculated Kendall Tau coefficients in the pairwise fashion for all features. The heat map in Figure 3a reveals significant hierarchical clustering of Kendall Tau coefficients between proteins (rows) and small molecules (lipids and metabolites; columns) with statistical significance for their correlation with each other. Significance for COVID-19 status and HFD-45 is also colored as a binary metric where q-value <0.05. The approach identified many dense clusters with strong associations across biomolecule classes.

• Download figure
• Open in new tab

Figure 3: Leveraging the value of multi-omic data through cross-ome correlation analysis.

a Hierarchical clustering of Kendall Tau coefficients calculated for correlations between abundances of proteins (rows) and small molecules (lipids and metabolites; columns) in the pairwise fashion. Significance of their association with HFD-45 and COVID-19 status is indicated above the biomolecule clusters. b Re-clustering of biomolecules found in the clusters highlighted in panel a with molecule annotations. c Enrichment analyses of protein GO terms (purple) and small molecule classes (green) present in the cluster in panel b. d A schematic of a high-density lipoprotein (HDL) particle containing APOA1 and APOA2 proteins surrounded by various lipids, specifically plasmalogens. SAA2, also detected in the cluster in panel b, can replaced APOA1 within the particle. e Relative abundance measurements of plasma gelsolin (pGSN), cellular gelsolin (cGSN), and total gelsolin obtained using parallel reaction monitoring (PRM) on representative peptide sequences. * and ** indicate p-values < 0.05 and 0.001, respectively. f Regression analysis of plasma gelsolin levels and SOFA scores (R² = 0.267, p = 4.53 × 10-5).

One cluster with a diverse mix of biomolecules and multiple associations between each other and with COVID-19 status and severity is highlighted in Figure 3b. This cluster captured the lipids PC and plasmenyl-PCs and high density lipoproteins (Figure 3c). These categories were found to be significantly associated with COVID-19 status and severity (vide supra) and importantly this cross-ome analysis links these molecules as possible molecular interactors. We recently reported that plasma plasmenyl-PCs (containing 16:0) had strong genetic associations to the APOA2 locus in mice (Linke et al., in. press, http://LipidGenie.com), and now, with this cross-ome correlation analysis, we provide further evidence that these specific plasmalogen species are associated with APOA2 and APOA1 in humans. APOA1 and APOA2 are among the most abundant proteins in HDL particles (Figure 3d) (Davidson et al., 2009). During proinflammatory response SAA2 can displace APOA1 within HDL particles and change HDL particles functions, and lead to increased clearance of HDL (Feingold and Grunfeld, 2019; Khovidhunkit et al., 2004). Reduced levels of HDL are associated in increased risk of hospitalization for infections (Trinder et al., 2020) and poorer outcomes during infections (Feingold and Grunfeld, 2019). The simultaneous reduction in apolipoproteins and plasmenyl PCs suggest that loss of these plasmenyl species also leads to poorer outcomes for COVID-19 patients. We suspect this may be due to the reduced capacity to handle oxidative damage, as plasmalogens are potent antioxidants (Messias et al., 2018). These data highlight the potential of mitigating COVID-19 severity by using therapeutic strategies aimed at restoring HDL particles, potentially through supplementing plasmalogen lipids or balancing HDL and LDLs through statins (Feingold and Grunfeld, 2019).

The notable cluster in Figure 3c also includes the protein gelsolin (GSN) and the metabolite citrate. Citrate is a calcium chelator that is used to prevent coagulation; intriguingly, GSN is a Ca²⁺-activated, actin-severing protein. Our cross-ome analysis offers the possibility that GSN and citrate are co-regulated in COVID-19 patients. Both features were among the 219 most significant features relating to COVID-19 infection and severity, and their reduced abundances were associated with worse outcomes. A possible confounding variable is that some patients received citrate as part of renal replacement therapy prior to sample collection; however, we find this confounder has no significant effect on the overall observation that citrate was lower in COVID-19 vs. non-COVID-19 (Supplementary Figure 3). The decrease in GSN in sera was recently reported by Messner et al., and confirming the reduced GSN levels in plasma samples (vs. sera) helps mitigate the confounding issue of variable loss of GSN into the fibrin clot (Piktel et al., 2018).

Two proteoforms of GSN arise from the use of alternative transcriptional start sites: the cytoplasmic form, involved in remodeling of intracellular actin (cGSN) (Feldt et al., 2019), and the plasma form, secreted by multiple cell types, including skeletal muscle (pGSN). Pools of circulating GSN, the aggregation of the two proteoforms, scavenge circulating filamentous actin and possess anti-inflammatory activities (Piktel et al., 2018). As cGSN may enter plasma upon cellular injury, we used parallel reaction monitoring to independently quantify levels of pGSN and cGSN via targeting five distinguishing peptides (Supplementary Table 4) (Peterson et al., 2012).

Abundance of pGSN correlated strongly with COVID-19 and ICU status, while cGSN levels were not reduced and did not correlate with COVID-19 severity (Figure 3e). These results indicate that cGSN from damaged cells is a relatively small contributor to the pool of circulating GSN and suggest that reduction of pGSN in COVID-19 may result from decreased biosynthesis and/or increased clearance from plasma.

Further, in the sub-group of ICU patients, decreased pGSN significantly correlated with SOFA scores (R² = 0.27, p-value = 4.5 ⨯10⁻⁵, Figure 3f), demonstrating its potential utility as a prognostic biomarker of COVID-19 severity. These data also indicate that pGSN supplementation has excellent potential as a COVID-19 therapeutic. Indeed, pGSN is known to decrease in acute lung injury, and in animal models of lung injury repletion of pGSN has favorable effects(Piktel et al., 2018). Recombinant pGSN has been given safely in phase I/II trials to patients with community-acquired pneumonia (NCT03466073, ClinicalTrials.gov). As of this writing, a trial of the recombinant pGSN in COVID-19 is being planned (NCT04358406).

Biological processes dysregulated in COVID-19 patients

Based on recent literature reports on COVID-19 pathophysiology and our systems analyses, we next focused on specific biological processes evidently dysregulated in COVID-19 patients: neutrophil degranulation, vessel damage, platelet activation and degranulation, blood coagulation, and acute phase response.

Neutrophil degranulation

Our list of 219 features strongly correlated to COVID-19 status and severity contained multiple transcripts and proteins involved in neutrophil degranulation (GO term 0043312; Supplementary Table 1). The volcano plot in Figure 4a illustrates the mean fold-change in abundance of GO-associated proteins (pink) and transcripts (purple), plotted against adjusted p-values of significance with COVID-19. We noticed an increased expression of multiple genes related to neutrophil function, including PRTN3, LCN2, CD24, BPI, CTSG, DEFA1, DEFA4, MMP8, and MPO (Supplementary Table 3). The latter encodes neutrophil myeloperoxidase, a protein instrumental in complexing extracellular DNA for the development of neutrophil extracellular traps (NETs) (Jiao et al., 2020). NETs are extracellular aggregates of DNA, histones, toxic proteins, and oxidative enzymes released by neutrophils to control infections, and their overdrive can amplify tissue injury by inflammation and thrombosis (Twaddell et al., 2019). NETs have been associated with the pathogenesis of ARDS (Lefrançais et al., 2018; Middleton et al., 2020; Narasaraju et al., 2011) and thrombosis (Ali et al., 2019), both phenotypes observed in severe COVID-19 patients, and elevated neutrophil counts predict worse outcomes in COVID-19 patients (Wang et al., 2020a). Interestingly, several NETs proteins that are also linked with thrombosis (Healy et al., 2006; Wang et al., 2017), including calprotectin (S100A8/9), ferritin, CRP, and histone H3, were also increased in our cohort (Table 1 and Supplementary Table 3). Thus our data suggest that combating dysregulated NETs may present an avenue towards mitigation of disease severity in COVID-19. For example, dipyridamole, an FDA-approved drug that can inhibit NET formation by activation of adenosine A_2A receptors (Ali et al., 2019), has recently been shown to improve outcomes of COVID-19 patients with respiratory failure (Liu et al., 2020). Similarly, the platelet P2Y12 receptor antagonist ticagrelor has been proposed to attenuate NET formation in COVID-19 (Omarjee et al., 2020) and may simultaneously inhibit platelet activation, a process also involved in COVID-19 pathophysiology (Manne et al., 2020a). Drugs that antagonize the NET– IL1β loop, such as anakinra, are currently being tested (NCT04324021), and colchicine, a potent neutrophil recruiter and IL1β-secretion inhibitor, is also under evaluation (NCT04326790).

• Download figure
• Open in new tab

Figure 4. Biological processes dysregulated in COVID-19.

a Volcano plots highlighting proteins (pink) and transcripts (purple) assigned with the GO term 0043312 “Neutrophil Degranulation.” Increased point size signified the inclusion of the biomolecule in the list of 219 features most significantly associated with COVID-19 status and severity (Figure 2e). b Linear regressions of protein abundance vs. HFD-45 for the indicated proteins as measured in COVID-19 (left) and non-COVID-19 patients (right). Resulting R2 values and their associated +/- slope indicate the goodness of fit and change in abundance of a given protein with severity (HFD-45). Proteins that are more decreased in severe cases appear blue, while proteins that are increased in severe cases appear red. Significance of the protein vs. HFD-45 correlation is denoted by a dot (p-value < 0.01). c Relative abundance measurements of peptides attributed to plasma fibronectin (pFN) and cellular fibronectin (cFN). d Relative abundance measurements of VWF multimer and VWF Antigen-2 (VWF Ag2), as estimated based on relative abundances of its unique peptides. Peptide- and protein-level data are log₂-transformed and grouped into four categories, according to patient status: COVID-19 ICU (red), COVID-19 non-ICU (orange), non-COVID-19 ICU (blue), and non-COVID-19 non-ICU (green). * and ** indicate p-values < 0.05 and 0.001, respectively.

Vessel damage

Numerous proteins involved in the body’s response to blood vessel damage, including AGT, FBLN5, NID1, and SERPINB1, increased in the COVID-19 samples relative to non-COVID-19 group and were particularly higher in abundance in more severe patients (Figure 4b and Supplementary Figure 4, Category #1). Neuropilin-1 (NRP1) is a regulator of the vascular endothelial growth factor (VEGF)-induced angiogenesis, and increase in its abundance with COVID-19 severity is of particular relevance given the recent report describing excessive pulmonary intussusceptive angiogenesis at autopsy in COVID-19 (Ackermann et al., 2020). In contrast, other vessel damage-associated proteins, such as SERPINA4, significantly decreased in abundance in COVID-19 patients (p-value = 2.715 × 10⁻⁴), highlighting the regulatory intricacies of this biological process.

Platelet degranulation

In accord with earlier reports, we observed severe dysregulation of proteins associated with platelet activation and degranulation (Figure 4b, Category #2) (Manne et al., 2020b). For example, abundance of histidine-rich glycoprotein (HRG) - which binds heme, heparin, thrombospondin, and plasminogen (Priebatsch et al., 2017) - was significantly reduced in patients with COVID-19, with the mean level 1.8-fold lower than in non-COVID-19 patients (Supplementary Figure 4, p-value = 8.138 ⨯ 10⁻¹²). Similarly, abundances of GPLD1 and CLEC3B significantly decreased in COVID-19 and correlated with disease severity (Figure 4b Category #2, Supplementary Figure 4). In contrast, the abundances of other platelet-associated proteins were increased in COVID-19 samples vs. non-COVID-19 samples (Supplementary Figure 4), for example, serglycin (SRGN) and the von Willebrand Factor (VWF); VWF has been recently implicated in COVID-19-associated endotheliopathy by the use of an FDA-approved antibody test (Goshua et al., 2020).

VWF is synthesized by endothelial cells as a 2,813 amino acid-long protein that is processed intracellularly to form VWF antigen 2 (VWFAg2) and VWF multimer (Haberichter, 2015). VWFAg2 is released constitutively into circulation, whereas VWF multimer is stored for later release upon stimulation of endothelial cells (Nightingale and Cutler, 2013). To distinguish between these two products, we quantified 19 peptides unique to VWFAg2 and 107 peptides unique to VWF (Figure 4c). Abundances of both VWFAg2 and VWF multimer were increased and correlated with COVID-19 severity (p-values = 4.860 × 10⁻² and 1.715 × 10⁻², respectively). We propose that endothelial cells increase synthesis, package, and release both VWFAg2 and VWF multimers, particularly in severe COVID-19 cases.

Fibronectin (FN1) works alongside VWF and fibrinogen to mediate the interaction of platelets with endothelial surfaces of injured vessels (Wang et al., 2014). Previous studies have reported increased FN1 levels in COVID-19 patients, relative to non-COVID-19 patients (Stukalov et al., 2020), but did not distinguish cellular (cFN) and plasma (pFN) proteoforms. By monitoring unique peptides (IAWESPQGQVSR for cFN and ESVPISDTIIPAVPPPTDLR for pFN), we specifically determined the abundance of each and revealed that cFN was significantly increased (p-value = 2.8 ⨯ 10⁻⁶) in COVID-19 patients, relative to non-COVID-19 patients, whereas pFN was not (Figure 4d). Such a finding is of considerable interest given that forced over-expression of cFN in plasma of mice is associated with increased thrombo-inflammation (Dhanesha et al., 2019).

Coagulation cascade

Consistent with the widely reported hypercoagulative phenotype of COVID-19 patients (Becker, 2020; Goshua et al., 2020; Lax et al., 2020), all our patients exhibited evidence of excessive clotting in vivo as demonstrated by the increased levels of circulating fibrin D-dimer (Table 1). To further examine the coagulative dysregulation in COVID-19 patients at molecular level, we deliberately conducted our proteomic analyses using plasma, which retains clotting factors otherwise depleted in serum. These plasma proteins, working alongside proteins of the vessel wall and platelets, undergo a cascade of regulated proteolytic reactions to generate thrombin, which converts fibrinogen to fibrin (Kattula et al., 2017). We detected significant increases in abundance of fibrinogen alpha (FGA) and beta (FGB) in COVID-19 vs. non-COVID-19 patient plasma, which was consistent with the clinically measured Fibrinogen (Supplementary Figure 4, Table 1), but to our surprise gamma (FGG) chains showed less significant trends in COVID-19 patients. We also observed decreases in abundance of prothrombin (F2) and thrombin-activation factor XIII (F13A1, F13B) in COVID-19 samples compared to non-COVID, and these proteins were further decreased in the most severe patients (Figure 4b and Supplementary Figure 4, Category #3). Significant reductions and correlations with HFD-45 were found for plasminogen (PLG), the zymogen precursor of plasmin (Urano et al., 2018); kallikrein (KLKB1) and kininogen (KNG1), which function in the intrinsic coagulation cascade leading to generation of thrombin (Schmaier, 2016); heparin cofactor 2 (SERPIND1), which forms an inhibitory complex with thrombin (Huntington, 2013); protein C (PROC), which is activated by thrombin to activated protein C (APC) (Griffin et al., 2018); and protein C inhibitor (SERPINA5) (Figure 4b, Supplementary Figure 4). The decreases of abundance of PROC and SERPINA5 with severity assume greater significance when considered alongside the finding that soluble plasma thrombomodulin (THBD) rivals elevated VWF as a predictor of mortality in COVID-19 (Goshua et al., 2020). THBD is normally tethered to the surface of endothelial cells, where it forms a complex with thrombin that efficiently activates PROC to APC (Ikezoe, 2015). APC also serves to modulate inflammatory response, and administration of recombinant APC mutant has been proposed as a therapy to mitigate thromboinflammation which occurs with ARDS and in severe COVID-19 patients (Griffin and Lyden, 2020).

Our data are pertinent to current discourse on prevention and treatment of the thrombotic diathesis in COVID-19. Several serpin-protease interactions, modified by heparin and heparin-like compounds, are being investigated as anti-coagulant therapeutics for COVID-19 patients (Cattaneo et al., 2020). Our data suggest such compounds and alternative antithrombotic therapies (e.g., anti-platelets) may also be useful in the subset of most severely affected COVID-19 patients.

Acute phase response

Next, we evaluated acute phase proteins – i.e., proteins whose plasma concentration changes by at least 25% during inflammatory disorders as a result of reprogramming by cytokines and are often secreted from the liver (Gabay and Kushner, 1999) – including S100A8, S100A9, SAA1 and SAA2. Our findings both support previous reports of significant increase in abundance of these four proteins in COVID-19 patient samples (p-values = 1.673 ⨯ 10⁻³, 1.966 ⨯ 10⁻⁴, 2.830 ⨯ 10⁻², 9.71 ⨯ 10⁻³, respectively) and reveal their correlation with severity as estimated by HFD-45 (p-values = 1.408 ⨯ 10⁻⁸, 1.996 ⨯ 10⁻⁷, 1.808 ⨯ 10⁻⁴, 1.420 ⨯ 10⁻⁵, respectively) (Supplementary Figure 4, Figure 4b Category #4). We also found that attractin (ATRN), a protein involved in chemotactic activity (Duke-Cohan et al., 1998), decreased in abundance COVID-19 samples vs. non-COVID-19 samples and was significantly associated with severity (p-value = 1.124 ⨯ 10⁻³). Four of these proteins were also found in the cross-ome correlation analysis (Figure 3b), where S100A8, S100A9, and SAA2 were negatively correlated with plasmalogens, while ATRN was positively correlated with these lipids.

Our data also revealed the increase in abundance of a liver-derived protein transketolase (TKT). One of the top 219 significant features (Figure 2), TKT was significantly elevated in COVID-19 patient plasma (p-value = 5.592 × 10⁻⁴, Supplementary Table 3). Other acute-phase proteins detected include GSN (vide supra) and SERPINA3, both of which decreased in abundance in COVID-19 patients and correlated with HFD-45 (p-value = 7.823 ⨯ 10⁻⁴ for SERPINA3). In general, abundances of most identified plasma proteins associated with the acute phase response were significantly altered with COVID-19 status and indicative of severity.

An interactive web tool to facilitate COVID-19 multi-omic data analysis

A primary goal of this study was to create a multi-omic compendium of biomolecules and their abundances in COVID-19 that could aid in the elucidation of disease pathophysiology and therapeutic development. To this end we have developed a web-based tool with interactive visualizations that allows for simple and quick navigation of our resource and performs on-demand PCA, linear regression, differential expression, and hierarchical clustering analysis via a clustered heatmap (https://covid-omics.app, Figure 5). The webtool enables supervised and unsupervised, univariate and multivariate analysis of the omic datasets in the context of highly curated clinical data. While the web resource will enable global access without need for code development, the source code is openly available on GitHub (https://github.com/ijmiller2/COVID-19_Multi-Omics) and the underlying relational database and unprocessed data are freely available (Methods).

• Download figure
• Open in new tab

Figure 5. Overview of the COVID-19 Multi-omics Web Tool.

a The home page provides principal component analysis (PCA) scores and loadings plots. Selected biomolecules are presented in a barplot and a boxplot. Each page provides buttons to navigate to the other web tools. b The differential expression page displays a multi-omics volcano plot with the y-axis representing -log10(p-values) where the p-values derive from the analysis in Figure 2 c The linear regression page allows users to select any combination of biomolecule and clinical measurement to analyze via univariate linear regression. R2 and p-values for the F-statistic are displayed on the plot. d The Clustergrammer page offers an interactive clustered heatmap.

To illustrate how this resource might be leveraged, we focus on a standout in the proteomics dataset, Cartilage acidic protein 1 (CRTAC1), which appears on the diagonal axis of the PCA loadings plot separated in both principal components 1 and 2 (Figure 5a). The biomolecule bar plots and boxplots, found just below the PCA plot on the resource’s homepage, show the per-sample abundance and distribution across sample groups of this protein, respectively (Figure 5a). Collectively, these data visualizations pinpoint CRTAC1 as an important feature in the innate structure of the proteomics data with a distinct abundance distribution for each patient subgroup (based on COVID-19 and ICU status), with the COVID-19 ICU group exhibiting the lowest average levels. CRTAC1 is a secreted matrix-associated protein involved in chondrocyte development and, interestingly, also belongs to the GO process “Olfactory Bulb Development” (GO:0021772), which is notable given previous reports of anosmia in COVID-19 patients (Marinosci et al., 2020).

To learn how CRTAC1 relates to our clinical measurements, such as HFD-45, the user can redirect to the linear regression page - a flexible tool for generation of scatter plots and calculations of R² values and (coefficient) p-values between any biomolecule and clinical covariate pair (Figure 5b). For the levels of CRTAC1 and HFD-45, the linear regression tool reports the R² value as 0.14 and a p-value of 1.704e-05 (n=125). Checkboxes on the control panel of the page allow users to select various patient subgroups and independently calculate the metrics for them. For example, for the COVID-19 subset, the R² is 0.079 and p-value is 4.490e-03 (n=100) (Supplementary Figure 5a), whereas for the non-COVID-19 subset of patient samples, the R² is 0.099 and p-value is 1.254e-01 (n=25) (Supplementary Figure 5b).

The third analysis enabled by the web tool is differential expression. The page depicts a volcano plot, where effect size represents the log₂ of the fold change of biomolecule abundance in COVID-19 vs. non-COVID-19 patient samples, and significance is -log₁₀(p-value). Here the p-value is calculated using a likelihood ratio that compares linear models with and without possible confounders, including age, sex, and ICU status. For any given biomolecule, an accompanying table (not shown) lists the effect size, p-value, q-value (adjusted by false discovery rate [FDR]), type of statistical test, and associated formula used to calculate the p-values. This information provides users with a robust means of examining the implications of potential confounding variables. According to the differential expression analysis, abundance of CRTAC1 is the second most reduced in the comparison of COVID-19 vs. non-COVID-19 patients (Figure 5c). Even after accounting for age, sex, and ICU status, the adjusted p-value of the reduction in abundance of CRTAC1 in COVID-19 vs. non-COVID-19 samples is 2.236e-08.

The final data exploration tool is an interactive heatmap, powered by Clustergrammer (Fernandez et al., 2017) and provided to aid in investigation of how a biomolecule of interest may fit into larger patterns across samples (Figure 5d). Interestingly, CRTAC1 groups with a series of apolipoproteins (vide supra). Thus, as exemplified by our examination of CRTAC1, covid-omics.app permits facile and rapid exploration of any one of the >17,000 quantified biomolecules, providing a potential departure point towards novel biological insights.

Two resource applications: comparative transcriptomics analysis and machine learning

To demonstrate how our resource could be leveraged by the scientific community, we present two complementary analyses: 1) we compare our transcriptomic data to the existing public transcriptomic datasets from patients with classical ARDS to elucidate gene expression signatures unique to COVID-19-induced ARDS, and 2) we utilize our multi-omic data to train machine learning models that accurately predict COVID-19 severity.

Severe COVID-19 is associated with development of ARDS, and although the evidence exists that COVID-19-induced ARDS is distinct from the classical ARDS (Gattinoni et al., 2020b), more data are needed to clarify the differences between the two conditions. We performed differential expression analysis on transcriptomic data obtained on samples of COVID-19 patients fulfilling diagnostic criteria of ARDS in this study, i.e., patients in the ICU (Ferguson et al., 2012), and an existing transcriptomic dataset on samples from patients with non-COVID-19 ARDS (Englert et al., 2019a). After removing differentially expressed genes that were likely affected by batch effects (Methods), we found 436 and 470 genes more highly expressed in the COVID-19-ARDS and the non-COVID-ARDS cohorts, respectively (Supplementary Figure 6). We then examined these gene sets separately via the gene set enrichment analysis. Using an FDR threshold of 0.05, we found 30 gene sets to be enriched in the genes more highly expressed in the COVID-19-ARDS patients (Figure 6a and Supplementary Table 5). Among them was TNF-alpha signaling, involved in the cytokine storm of severe COVID-19 patients (Catanzaro et al., 2020). TNF-alpha inhibitors have been proposed for treatment of COVID-19 (Soy et al., 2020), and a TNF-alpha clinical trial has been registered (ChiCTR2000030089). We also identified 79 gene sets enriched in genes more highly expressed in non-COVID-19-ARDS patients (Figure 6b and Supplementary Table 5). Several of them were related to viral response, e.g., decreased interferon signaling, suggesting that, in agreement with a recent report (Blanco-Melo et al., 2020), COVID-19 may be characterized by a dysregulated leukocyte viral response. These results indicate that leukocytes from ARDS patients with and without COVID-19 possess overlapping, yet distinct transcriptomic profiles.

• Download figure
• Open in new tab

Figure 6. Results from analyses demonstrating use-cases of this multi-omic resource.

a The top-ten enriched gene sets ranked by their adjusted p-value. For the gene set “TFNA signalling via NFKB”, we show a heatmap (right) of the z-score normalized expression data (in units of log transcripts per million) partitioned by whether the data came from the COVID-19-ARDS patients (right) or the non-COVID-19-ARDS (i.e. sepsis ARDS) patients from Englert et al. (left). The first row of each heatmap depicts the hospital-free days of each patient. We note that hospital-free days are not available for the Englert et al. dataset. The gene names labelling each row are colored according to whether the gene was deemed by EBSeq to be more highly expressed in COVID-19 ARDS (orange) or non-COVID-19 ARDS (blue). b Similar to (a); however, we instead analyze DE genes that are more lowly expressed in COVID-19 ICU patients. c Data splitting scheme for training and test sets from the 100 COVID patients with all four omic datasets. A random 20% was held out to be used for model evaluation, and the remaining 80% was used to determine the best hyperparameters with 5-fold cross validation. d Extra trees classifier performance metrics on the test set after hyperparameter optimization using each of the four omic datasets separately for training or all omic data combined. e Macro-averaged receiver-operator characteristic curves for the models trained with multi-omic data, Charlson score, or both multi-omic data and charlson score. f Test set predictions of the extra trees model trained on the combined multi-omic dataset showing correct predictions as a function of the disease severity defined by hospital free days. G Top 5 most important predictive features for each of the models trained on the four omic subsets. Feature importance for each set was normalized to the most important feature.

The ability to predict COVID-19 severity early in admission would help clinicians prioritize patients for treatments. Thus, we tested whether a machine learning model trained with the multi-omic data in our resource could accurately predict COVID-19 severity. The 100 COVID-19 patients with complete multi-omic profiles were assigned to “severe” and “less severe” groups based on HFD-45 using the median day 26 as the threshold. They were then randomly split into two groups: 80 were used for training and the remaining 20 were held out as a test set (Figure 6c). We optimized hyperparameters of ExtraTrees classification models using both the combined data and each of the four omic datasets separately via 5-fold cross-validation on the training set. This comparison shows that each omic dataset possessed different predictive utility (combined model AUROC = 0.93, average precision = 0.96; Figure 6d). Examination of the individual predictions from the combined model show that 2/4 incorrect predictions were for the patients with HFD-45 values near the median (Figure 6e), all of whom were also misclassified by the omic subset models (Supplementary Figure 7). Although preexisting comorbidities are associated with poor COVID-19 outcomes (Guan et al., 2020), the model based on the combined multi-omic data was a substantially better predictor of disease severity than a model based on Charlson score, the clinical estimate of 10-year survival (Richardson et al., 2020b). Adding Charlson score to the multi-omic data did not improve model performance (Figure 6f).

The use of the tree-based model also enabled determination of the most predictive features (Figure 6g, Supplementary Table 6). Some of them overlap with other analyses presented above, such as plasmalogen lipids and proteins S100A8 and S100A9. However, this model-based prioritization also highlighted potentially disease-relevant features not yet discussed. For example, for the metabolomic model, the most important features were tryptophan metabolites kynurenine and quinolinic acid, which are known mediators of immune function and inflammation in various contexts (Heyes et al., 1992; Sofia et al., 2018) and have been recently implicated in COVID-19 severity (Migaud et al., 2020). Comparison of the top 25 predictive features from our proteomic and metabolomic models to those from another COVID-19 severity prediction model (Shen et al., 2020) revealed several overlapping features, including metabolite L-Kynurenine and proteins SAA2, SERPINA3, and LCP1. All in all, these results indicate that our multi-omic data can be utilized to accurately predict COVID-19 severity, and the feature importance from these models might further aid in understanding disease pathophysiology.