A diagnostic host response biosignature for COVID-19 from RNA profiling of nasal swabs and blood

Abstract

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes coronavirus disease-19 (COVID-19), has emerged as the cause of a global pandemic. We used RNA sequencing to analyze 286 nasopharyngeal (NP) swab and 53 whole-blood (WB) samples from 333 patients with COVID-19 and controls. Overall, a muted immune response was observed in COVID-19 relative to other infections (influenza, other seasonal coronaviruses, and bacterial sepsis), with paradoxical down-regulation of several key differentially expressed genes. Hospitalized patients and outpatients exhibited up-regulation of interferon-associated pathways, although heightened and more robust inflammatory responses were observed in hospitalized patients with more clinically severe illness. Two-layer machine learning-based host classifiers consisting of complete (>1000 genes), medium (<100), and small (<20) gene biomarker panels identified COVID-19 disease with 85.1-86.5% accuracy when benchmarked using an independent test set. SARS-CoV-2 infection has a distinct biosignature that differs between NP swabs and WB and can be leveraged for COVID-19 diagnosis.

Fig. 1. Overview of sample collection and metatranscriptomic analysis.

(A) Flowchart of NP swab and WB sample collection. (B) Box and whisker plots of RT-PCR cycle threshold (Ct) values of SARS-CoV-2–positive individuals who are outpatients (n = 55) were compared with those who are hospitalized, non-ICU (n = 17), or in the ICU (n = 7). There was no difference in viral load, inversely related to the Ct value, regardless of disease severity [P = 0.89 by analysis of variance (ANOVA)]. (C to E) Box and whisker plots of abundance (C), Chao richness (D), and Shannon diversity (E) of the viral metatranscriptome in patients with SARS-CoV-2 (COVID-19) (n = 137), respiratory viruses (“Other virus”) (n = 41), and without respiratory viruses (“No virus”) (n = 108), stratified by the inclusion (“Including respiratory viral reads”) or exclusion (“Exclusion respiratory viral reads”) of respiratory viral reads. (F to H) Box and whisker plots of abundance (F), Chao richness (G), and Shannon diversity (H) of the bacterial metatranscriptome. (I) Distribution of viral families in each group, expressed as log₁₀-normalized RPM. (J) Distribution of the top 10 bacterial families in each group. For box and whisker plots, the median is represented by a dotted line, boxes represent the first to third quartiles, whiskers represent the minimum and maximum values, and jitters represent the distribution of the population. For (C) to (H), statistical analysis was conducted by Kruskal-Wallis test, followed by the Nemenyi test for post hoc analysis.

Fig. 2. Comparison of canonical pathways predicted to be involved in COVID-19.

(A) Comparison of NP swab pathways in COVID-19, influenza, other respiratory viruses (human metapneumovirus, human rhinovirus, and human parainfluenza 2), and seasonal coronaviruses. (B) Comparison of NP swab pathways in COVID-19 outpatients and hospitalized patients (non-ICU and ICU). (C) Comparison of NP swab and WB pathways in COVID-19–positive patients. (D) Comparison of WB pathways in COVID-19, influenza, and bacterial sepsis. Pathway prediction is determined by z-score; a positive value denotes up-regulation, and a negative value denotes down-regulation. All pathways and z-scores were calculated relative to NP swab and WB donor controls.

Fig. 3. Heatmap of pathways predicted to be involved in COVID-19 and other respiratory viruses from NP swabs.

Comparison of pathways for (A) all patients with COVID-19 compared with donor controls, (B) outpatients compared with donor controls, and (C) hospitalized (non-ICU and ICU) patients compared with donor controls. (D) Influenza compared with donor samples. (E) Seasonal coronaviruses (HCoV HKU-1, HCoV NL-63, HCoV OC43, and HCoV 229E) compared with donor samples. (F) Other respiratory viruses (human metapneumovirus and human rhinovirus) compared with donor samples. Dev, development; Fxn, function; GI, gastrointestinal; morph, morphology.

Fig. 4. Heatmaps of NP swab and WB pathways.

(A) Hospitalized patients compared with outpatients with COVID-19 using NP swabs. Comparison of WB pathways in hospitalized patients with COVID-19 compared with (B) donor controls, (C) bacterial sepsis, and (D) influenza.

Fig. 5. Hierarchically clustered heatmaps of DEGs.

(A) NP swab DEGs from patients with infection by SARS-CoV-2 compared with influenza, seasonal coronaviruses, and other respiratory viruses, or virus-negative patients. Samples are stratified by diagnosis (“Infection”) and highest level of care (“Level of care” including outpatient, hospitalized, non-ICU, and ICU). (B) WB DEGs from patients with infection by SARS-CoV-2 compared with influenza, bacterial, and bacterial sepsis. Samples are stratified by diagnosis (“Infection”). (C) Comparison of DEGs among outpatients and hospitalized patients admitted (“ICU”) or not admitted (“non-ICU”) to the ICU. For (A), DEGs with Bonferroni-corrected P value <0.001 were included. For (B) and (C), DEGs with Bonferroni-corrected P value <0.01 were included. All noncoding genes were removed. For (A), DEGs are calculated relative to individuals with nonviral ARIs, and the top 100 DEGs are included; for (B), DEGs are calculated relative to donor controls, and the top 100 DEGs are included; for (C), DEGs are calculated on the basis of pairwise comparisons between outpatients and hospitalized patients, and all DEGs are included. Up-regulated genes are colored in red, and down-regulated genes are colored in blue. The distribution of predicted pathways (“Host response pathways”) is shown by pie graphs above each group in the heatmap; cellular functions are indicated in shades of green, immune functions are indicated in shades of red, and other pathways are indicated in black. DEG groups (A to L) were identified on the basis of hierarchical clustering with complete linkage, after exclusion of noncoding genes. The lettering is presented in alphabetical order, starting with A in the top left corner. Dotted outlines highlight key DEG groups, which are distinctly up-regulated in SARS-CoV-2 infection compared with other disease states.

Fig. 6. Venn diagrams of DEGs.

(A) Comparison of NP swab DEGs in patients with SARS-CoV-2 and influenza. (B) Comparison of NP swab and WB DEGs in patients with SARS-CoV-2. (C) Comparison of NP swab DEGs in COVID-19 outpatients and hospitalized patients. (D) Comparison of NP swab and WB DEGs in patients with influenza. The DEGs are calculated relative to donor controls.

Fig. 7. Diagnostic classifier for COVID-19.

(A) Overview of classifier design and distribution of samples for training and test sets in layers 1 and 2. (B) Pie graphs showing the distribution of predicted pathways (“Host response pathways”) represented by the genes within the full gene panel (left), medium gene panel (middle), or the small gene panel (right).

Fig. 8. Performance characteristics of the two-layered (combined layers 1 and layer 2) classifier for full, medium, and small gene panels.

(A to C) Training set ROC curve (left) and test set violin plot (middle) and confusion matrix (right) for the two-layer classifier, using either the full gene panel (A), medium gene panel (B), or the small gene panel (C).