Upper airway gene expression differentiates COVID-19 from other acute respiratory illnesses and reveals suppression of innate immune responses by SARS-CoV-2
============================================================================================================================================================

* Eran Mick
* Jack Kamm
* Angela Oliveira Pisco
* Kalani Ratnasiri
* Jennifer M. Babik
* Carolyn S. Calfee
* Gloria Castañeda
* Joseph L. DeRisi
* Angela Detweiler
* Samantha Hao
* Kirsten N. Kangelaris
* G. Renuka Kumar
* Lucy M. Li
* Sabrina A. Mann
* Norma Neff
* Priya A. Prasad
* Paula Hayakawa Serpa
* Sachin J. Shah
* Natasha Spottiswoode
* Michelle Tan
* Stephanie A. Christenson
* Amy Kistler
* Charles Langelier

## Abstract

We studied the host transcriptional response to SARS-CoV-2 by performing metagenomic sequencing of upper airway samples in 238 patients with COVID-19, other viral or non-viral acute respiratory illnesses (ARIs). Compared to other viral ARIs, COVID-19 was characterized by a diminished innate immune response, with reduced expression of genes involved in toll-like receptor and interleukin signaling, chemokine binding, neutrophil degranulation and interactions with lymphoid cells. Patients with COVID-19 also exhibited significantly reduced proportions of neutrophils and macrophages, and increased proportions of goblet, dendritic and B-cells, compared to other viral ARIs. Using machine learning, we built 26-, 10- and 3-gene classifiers that differentiated COVID-19 from other acute respiratory illnesses with AUCs of 0.980, 0.950 and 0.871, respectively. Classifier performance was stable at low viral loads, suggesting utility in settings where direct detection of viral nucleic acid may be unsuccessful. Taken together, our results illuminate unique aspects of the host transcriptional response to SARS-CoV-2 in comparison to other respiratory viruses and demonstrate the feasibility of COVID-19 diagnostics based on patient gene expression.

## Introduction

The emergence of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in December 2019 has precipitated a global pandemic with over 4.5 million cases and 300,000 deaths1. Coronavirus disease 2019 (COVID-19), the clinical syndrome caused by SARS-CoV-2, varies from asymptomatic infection to critical illness, with dysregulated inflammatory response to infection a hallmark of severe cases2. Defining the host response to SARS-CoV-2, as compared to other respiratory viruses, is fundamental to identifying mechanisms of pathogenicity and potential therapeutic targets.

Metagenomic next generation RNA sequencing (mNGS) is a powerful tool for assessing host/pathogen dynamics3,4 and a promising modality for developing novel respiratory diagnostics that integrate host transcriptional signatures of infection3,5. While proven for diagnosis of other acute respiratory infections3,5, transcriptional profiling has not yet been evaluated as a diagnostic tool for COVID-19, despite its potential to mitigate the risk of false negatives associated with standard naso/oropharyngeal (NP/OP) swab-based PCR testing6–8.

### Results and Discussion

To interrogate the molecular pathogenesis of SARS-CoV-2 and evaluate the feasibility of a COVID-19 diagnostic based on host gene expression, we conducted a multicenter observational study of 238 patients with acute respiratory illnesses (ARIs) who were tested for SARS-CoV-2 by NP/OP swab PCR, and performed host/viral mNGS on the same specimens. The cohort (**Table S1**) included 94 patients who tested positive for SARS-CoV-2 by PCR, 41 who tested negative but had other pathogenic respiratory viruses detected by mNGS (**Methods, Figure S1A**), and 103 with no virus detected (non-viral ARIs).

We began by performing pairwise differential expression analyses between the three patient groups (**Methods, Supp. File 1**). Hierarchical clustering of the union of the 50 most significant genes in each of the comparisons revealed that the transcriptional response to SARS-CoV-2 was distinct from the response to other viruses (**Figure 1A**). We detected gene clusters that were up- (cluster I) or down-regulated (cluster II) by other viruses as compared to non-viral ARIs, but relatively unaffected by SARS-CoV-2. Importantly, we also identified a small number of genes that were upregulated by SARS-CoV-2 more than by other viruses (cluster III). And many genes upregulated in all viral ARIs (cluster IV) appeared to respond to SARS-CoV-2 proportionally to viral load, as measured by the relative abundance of sequencing reads mapped to the virus (**Methods, Figure S1B**).

![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/19/2020.05.18.20105171/F1.medium.gif)

[Figure 1.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/F1)

Figure 1. Host Transcriptional Signatures of SARS-CoV-2 Infection as Compared to Other Respiratory Viruses
**A**. Hierarchical clustering of 120 genes comprising the union of the top 50 DE genes by significance in each of the pairwise comparisons between patients with COVID-19, other viral ARIs and non-viral ARIs. Patient group labels and viral abundance of SARS-CoV-2 are shown in the annotation bars. rpM, reads-per-million. **B**. Normalized enrichment scores of selected REACTOME pathways that achieved statistical significance (P < 0.05) in at least one of the gene set enrichment analyses, comparing either COVID-19 vs. non-viral ARIs or other viral ARIs vs. non-viral ARIs. If a pathway could not be tested in one of the comparisons since it had less than 10 members in the input gene set, the enrichment score was set to 0. **C**. *In silico* estimation of cell type proportions. All pairwise comparisons were performed with a two-sided Mann-Whitney-Wilcoxon test followed by Bonferroni’s correction. **D**. Scatter plots of normalized gene counts (log2 scale) as a function of SARS-CoV-2 abundance. Robust regression was performed on COVID-19 patients with log10(rpM) > 0 to highlight the relationship to viral load. Shown are inflammasome related genes selected from among the genes most depressed in expression in COVID-19 compared to other viral ARIs. Statistical results for each gene (from top to bottom) refer to the regression analysis, the DE analysis between COVID-19 and non-viral ARIs, and the DE analysis between COVID-19 and other viral ARIs.

To investigate the pathways driving these distinctions, we performed gene set enrichment analyses9 (GSEA) on the genes differentially expressed (DE) between SARS-CoV-2 and non-viral ARIs, and separately, those DE between other viral ARIs and non-viral ARIs (**Methods, Supp. File 2**). We found that both SARS-CoV-2 and other viruses elicited an interferon response in the upper airway (**Figure 1B**). The most significant genes upregulated by SARS-CoV-2 were interferon inducible, including *IFI6, IFI44L, IFI27* and *OAS2* (**Figure S2A**), in agreement with previous reports10,11. *IFI27* was induced by SARS-CoV-2 significantly more than by other viruses, even at low viral load. Most other top DE genes, however, did not distinguish COVID-19 from other viral ARIs. *ACE2*, which encodes the cellular receptor for SARS-CoV-2, was also non-specifically induced, consistent with its recent identification as a general interferon stimulated gene12.

Notably, GSEA of DE genes in the direct comparison of SARS-CoV-2 and other viruses suggested elements of the interferon response to SARS-CoV-2 were attenuated (**Figure S2B, Supp. File 2**). Indeed, numerous interferon response genes, such as *IRF7* and *OASL*, were more potently induced by other viruses, and high SARS-CoV-2 abundance was required to achieve comparable induction (**Figure S2C**). These results may be related to observations of a blunted interferon response in cellular models of SARS-CoV-2 infection13, though the effects in patients appear more nuanced.

A striking contrast between SARS-CoV-2 and other viruses emerged in the activation of additional innate immune signaling pathways (**Figure 1B, S2B**). Other viral ARIs caused significant upregulation of gene expression associated with toll-like receptors, interleukin signaling, chemokine binding, neutrophil degranulation and interactions with lymphoid cells, yet the response of these pathways to SARS-CoV-2 was markedly attenuated (**Figure 1B, S2B**). While other viral ARIs appeared to depress expression of genes involved in cilia functions and antioxidant responses, this was not observed for SARS-CoV-2 (**Figure 1B, S2B**).

*In silico* estimation of cell type proportions revealed significant differences between the groups (**Figure 1C, S3**). Compared to patients with other viral and non-viral ARIs, those infected with SARS-CoV-2 exhibited significantly reduced fractions of monocytes/macrophages and neutrophils, and significantly increased proportions of goblet, dendritic and B-cells. Patients with other viral ARIs exhibited decreased ciliated cell and ionocyte fractions, and increased macrophage, neutrophil and T-cell fractions, compared to those with non-viral ARIs. These results closely aligned with the GSEA findings and suggested that the diminished innate immune responses in COVID-19 patients were coupled to differences in the cellular composition of the airway microenvironment.

The gene that was most decreased in expression in COVID-19 patients compared to those with other viral ARIs was *IL1B*, which encodes a pro-inflammatory cytokine produced by the inflammasome complex, particularly in macrophages14 (**Figure 1D, Supp. File 1**). Among the top 100 differentially decreased genes were those involved in inflammasome activation and activity (*NLRP3, CASP5, IL1A, IL1B, IL18RAP, IL1R2*) and in chemokine signaling for recruiting innate immune cells to the epithelium (*CCL2, CCL3, CCL4*). Given that IL1-β and other pro-inflammatory cytokines are primary targets of monoclonal antibody therapeutics under investigation15, these results raise the question of whether further suppression may be detrimental in the setting of an already suppressed inflammatory response to SARS-CoV-2.

Relatively few genes were specifically upregulated in COVID-19 patients compared to both other viral and non-viral ARIs. These included *TRO*, which encodes a membrane-bound cell adhesion molecule; *WDR74*, which plays a role in rRNA processing and associates with the RNA helicase MTR416; *EIF4A2*, a translation initiation factor that has been shown to interact with other coronaviruses as well as HIV17,18; and *FAM83A*, which is involved in epidermal growth factor receptor (EGFR) signaling19.

We next asked whether host gene expression data could be used to construct a classifier capable of accurately differentiating COVID-19 from other ARIs (viral or non-viral). By employing a combination of lasso regularized regression and random forest (**Methods**), we first identified a 26-gene signature that performed with an area under the receiver operating characteristic curve (AUC) of 0.980 (range of 0.951-1.000), as estimated by 5-fold cross validation (**Figure 2A, Tables S2, S3**). Even though many patients undergoing testing for COVID-19 may not be infected with other respiratory viruses, we recognized the need for classifier specificity in this circumstance and examined how well the classifier performed at distinguishing SARS-CoV-2 from other respiratory viruses. We found that it achieved an AUC of 0.966 (range 0.895-1.000) when tested only on patients with other viral ARIs, indicating robust specificity for SARS-CoV-2 (**Tables S2, S3**). Using a cut-off of 40% predicted out-of-fold probability for COVID-19 to call a patient positive, this translated into a sensitivity of 97% and a specificity of 96% for patients with non-viral ARIs and 83% for patients with other viral ARIs (**Figure 2B**).

![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/19/2020.05.18.20105171/F2.medium.gif)

[Figure 2.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/F2)

Figure 2. Performance of COVID-19 Diagnostic Classifiers Based on Patient Gene Expression
**A**. Receiver operating characteristic (ROC) curve for a 26-gene classifier that differentiates COVID-19 from other acute respiratory illnesses (viral and non-viral). **B**. Accuracy of the 26-gene classifier within each patient group, using a cut-off of 40% out-of-fold predicted probability for COVID-19. **C**. ROC curve for a 10-gene classifier. **D**. ROC curve for a 3-gene classifier. **E**. Out-of-fold predicted probability of COVID-19 derived from the 26-gene classifier plotted as a function of SARS-CoV-2 viral abundance. Dashed lines indicate 40%, our chosen cut-off, and 50%.

Given that a parsimonious gene set could enable practical incorporation into a clinical PCR assay, we implemented a more restrictive regression penalty and identified a 10-gene classifier that could differentiate SARS-CoV-2 from other respiratory illnesses with an AUC of 0.950 (range 0.918-0.974) (**Figure 2C, Tables S2, S3**). Classification performance specifically against other viral ARIs suffered slightly but still achieved an AUC of 0.905 (range 0.842-0.959). Existing SARS-CoV-2 PCR assays typically employ 3 gene targets and thus we tested the potential to further reduce host classifier gene size. We found that a sparse 3-gene (*IL1B, IFI6, IL1R2*) classifier still achieved an AUC of 0.893 (range 0.844-0.933) (**Figure 2D, Tables S2, S3**).

A host-based diagnostic might have particular utility if it could increase the sensitivity of standard NP/OP swab PCR testing, which may return falsely negative in a significant proportion of patients6–8. Presumably, false negatives are in large part due to insufficient viral abundance in the collected specimen. While our cohort did not include PCR-negative samples from patients with clinically confirmed COVID-19, we evaluated whether classifier performance was affected by viral load. The predicted probability of SARS-CoV-2 infection had little apparent relationship to the abundance of SARS-CoV-2, suggesting host gene expression has the potential to provide an orthogonal indication of COVID-19 status even when viral abundance is low (**Figure 2E**).

In summary, we studied 238 patients with acute respiratory illnesses to define the human upper respiratory tract gene expression signature of COVID-19. Our study is limited by sample size, incomplete demographic data and the need for an independent validation cohort. Notwithstanding, our results illuminate unique aspects of the host transcriptional response to SARS-CoV-2 in comparison to other respiratory viruses and provide insight regarding molecular pathogenesis. We also leveraged these data to develop an accurate, clinically practical, COVID-19 diagnostic classifier that may help overcome the limitations of direct detection of viral nucleic acid. Future prospective studies in a larger cohort will be needed to validate these findings, determine the prognostic value of host signatures, and assess the performance of integrated host/viral diagnostic assays.

## Data Availability

Gene counts, sample metadata, IDSeq reports and code for the DE analysis, cell type proportions analysis and COVID-19 gene expression classifiers are available on github.

[https://github.com/czbiohub/covid19-transcriptomics-pathogenesis-diagnostics-results](https://github.com/czbiohub/covid19-transcriptomics-pathogenesis-diagnostics-results)

## Funding

This study was supported by the Chan Zuckerberg Biohub, the Chan Zuckerberg Initiative, and the National Heart, Lung, and Blood Institute (1K23HL138461-01A1).

## Materials and Methods

### Study design, clinical cohort and ethics statement

We conducted an observational cohort study of patients with acute respiratory illnesses suspected to be COVID-19 at the University of California, San Francisco (UCSF) and Zuckerberg San Francisco General Hospital between 03/10/2020 and 04/07/2020. Through UCSF IRB protocol 17-24056, a waiver of consent was granted to evaluate unused clinical specimens in the UCSF Clinical Microbiology Laboratory and assess demographics and basic clinical features from the Epic-based electronic health record.

### SARS-CoV-2 detection by clinical PCR

Testing for COVID-19 was carried out in the UCSF Clinical Microbiology Laboratory using polymerase chain reaction (PCR) of NP swab or pooled NP + OP swab specimens using primers targeting either two regions of the SARS-CoV-2 N gene (n=156, 66%), or the E and RNA-dependent RNA polymerase genes (n=82, 34%), plus human RNAse P as a positive control. In all our analyses, we defined patients with COVID-19 as those with a positive SARS-CoV-2 result by PCR.

### Metagenomic sequencing

To evaluate host gene expression and detect the presence of other viruses, metagenomic next generation sequencing (mNGS) of RNA was performed on the same specimens subjected to SARS-CoV-2 PCR testing. Following DNase treatment, human cytosolic and mitochondrial ribosomal RNA was depleted using FastSelect (Qiagen). To control for background contamination, we included negative controls (water and HeLa cell RNA) as well as positive controls (spike-in RNA standards from the External RNA Controls Consortium (ERCC))1. RNA was then fragmented and subjected to a modified metagenomic spiked sequencing primer enrichment (MSSPE) library preparation method2. Briefly, a 1:1 mixture of the NEBNext Ultra II RNAseq Library Prep (New England Biolabs) random primers and a pool of SARS-CoV-2 primers at 100 μM was used at the first strand synthesis step of the standard RNAseq library preparation protocol to enrich for reads spanning the length of the SARS-CoV-2 genome. RNA-seq libraries underwent 146 nucleotide paired-end Illumina sequencing on an Illumina Novaseq 6000 instrument.

### Quantification of SARS-CoV-2 abundance by mNGS

All samples were processed through a SARS-CoV-2 reference-based assembly pipeline that involved removing non-SARS-CoV-2 reads with Kraken23, and aligning to the SARS-CoV-2 reference genome MN908947.3 using minimap24. We calculated SARS-CoV-2 reads-per-million (rpM) using the number of reads that aligned with mapq >= 20. When calculating log10(rpM), 0.1 was added to the rpM values of all samples to avoid taking the log of 0.

### Detection of other respiratory pathogenic viruses by mNGS

All samples were processed through the IDSeq pipeline5,6, which performs reference based alignment at both the nucleotide and amino acid level against sequences in the National Center for Biotechnology Information (NCBI) nucleotide (NT) and non-redundant (NR) databases, followed by assembly of the reads matching each taxon detected. We further processed the results for viruses with established pathogenicity in the respiratory tract3. We evaluated whether one of these viruses was present in a patient sample if it met the following three initial criteria: i) at least 10 counts mapped to NT sequences, ii) at least 1 count mapped to NR sequences, iii) average assembly nucleotide alignment length of at least 70bp.

Negative control (water and HeLa cell RNA) samples enabled estimation of the number of background reads expected for each virus, which were normalized by input mass as determined by the ratio of sample reads to spike-in positive control ERCC RNA standards7. Viruses were then additionally tested for whether the number of sequencing reads aligned to them in the NT database was significantly greater compared to negative controls. This was done by modeling the number of background reads as a negative binomial distribution, with mean and dispersion fitted on the negative controls. For each batch (sequencing run) and taxon (virus), we estimated the mean parameter of the negative binomial by averaging the read counts across all negative controls after normalizing by ERCCs, slightly regularizing this estimate by including the global average (across all batches) as an additional sample. We estimated a single dispersion parameter across all taxa and batches, using the functions glm.nb() and theta.md() from the R package MASS8. We considered a patient to have a respiratory pathogenic virus detected by mNGS if the virus achieved an adjusted p-value < 0.05 after Holm-Bonferroni correction for all tests performed in the same sample.

### Host differential expression (DE) analysis

Following demultiplexing, sequencing reads were pseudo-aligned with kallisto9 (v. 0.46.1; including bias correction) to an index consisting of all transcripts associated with human protein coding genes (ENSEMBL v. 99), cytosolic and mitochondrial ribosomal RNA sequences, and the sequences of ERCC RNA standards. Samples retained in the dataset had a total of at least 400,000 estimated counts associated with transcripts of protein coding genes, and the average across all samples was 5.79 million. Gene-level counts were generated from the transcript-level abundance estimates using the R package tximport10, with the lengthScaledTPM method.

Genes were retained for differential expression analysis if they had at least 10 counts in at least 20% of samples (n=15,900). The analysis was performed with the R package limma11 using quantile normalization and the design: ~0 + viral status + gender + age + sequencing batch, where viral status was either “SARS-CoV-2”, “other virus” or “no virus”. We note that the gender of patients for whom we lacked this information was inferred based on chromosome Y gene expression, and the age of patients for whom we lacked this information was taken as the mean age of samples with age reported in the respective viral status group.

To generate the gene expression heatmap, hierarchical clustering was performed on the union of the top 50 genes (by p-value) in each of the pairwise comparisons among the three groups (n=120 genes). Gene counts were subjected to the variance stabilizing transformation, as implemented in the R package DESeq212, centered and scaled prior to clustering. For both rows and columns, Euclidean distance was the distance measure and Ward’s criterion (ward.D2) was the agglomeration method.

### Gene set enrichment analysis

Gene set enrichment analyses13 were performed using the fgseaMultilevel function in the R package fgsea14 on REACTOME15 pathways with a minimum size of 10 genes and a maximum size of 1,000. The genes included in each pairwise comparison were those with Benjamini-Hochberg adjusted p-value < 0.1 and |log2(FC)| > log(1.5) in the respective DE analysis, preranked by fold change.

The gene sets shown in Figure 1B were manually selected to reduce redundancy and highlight diverse biological functions from among those with a Benjamini-Hochberg adjusted p-value < 0.05 in at least one of the comparisons i) SARS-CoV-2 vs. no virus, and ii) other virus vs. no virus. And the gene sets shown in Figure S2B were similarly selected from among those with an adjusted p-value < 0.05 in the direct comparison of SARS-CoV-2 vs. other virus. Full results of all analyses are provided as supplementary.

### Regression of gene counts against viral abundance

We performed robust regression of the limma-generated quantile normalized gene counts against log10 (rpM) of SARS-CoV-2 for all genes with a Benjamini-Hochberg adjusted p-value < 0.001 in either the DE analysis of SARS-CoV-2 vs. no virus, or SARS-CoV-2 vs. other virus (n=2,920). The samples included were those in the SARS-CoV-2 patient group with log10 (rpM) > 0. Robust regression was used to better account for outlier data points.

The analysis was performed using the R package robustbase16, which implements MM-type estimators for linear regression17,18, using the KS2014 setting and the model: quantile normalized counts (log2 scale) ~ gender + age + sequencing batch + log10 (rpM). Model predictions for the log10 (rpM) co-variate were used for display in the individual gene plots. Reported p-values for significance of the difference of the regression coefficient from 0 were Benjamini-Hochberg adjusted, and reported R2 values represent the adjusted robust coefficient of determination19.

### *In silico* analysis of cell type proportions

Cell-type proportions were estimated from bulk host transcriptome data using the CIBERSORT X algorithm20. We used the human lung cell atlas dataset21 to derive the single cell signatures. The cell types estimated with this reference cover all expected cell types in the airway samples. The estimated proportions were compared between the three patient groups using a two-sided Mann-Whitney-Wilcoxon test with Bonferroni correction.

### Classifier construction

We built sparse classifiers for COVID-19 status using a combined lasso and random forest approach. For feature selection, we used the logistic lasso (as implemented in the R package glmnet22), and then trained random forests on the selected features (using the R package randomForest23). We used 5-fold cross-validation to evaluate model error. For each train-test split, we used a nested cross-validation within the training set to select the lasso tuning parameter. For the random forest, we used 10,000 trees, and left all tuning parameters at their defaults. For the initial input features (before selection), we used gene counts with a variance-stabilizing transform derived from the training set only (using the R package DESeq212). Classifiers were built using a gold standard of COVID-19 diagnosis based on SARS-CoV-2 PCR positivity.

### Data availability

Gene counts, sample metadata, IDSeq reports and code for the DE analysis, cell type proportions analysis and COVID-19 gene expression classifiers are available at: [https://github.com/czbiohub/covid19-transcriptomics-pathogenesis-diagnostics-results](https://github.com/czbiohub/covid19-transcriptomics-pathogenesis-diagnostics-results)

## Supplemetary Tables

View this table:
[Table S1.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/T1)

Table S1. 
Cohort Clinical and Demographic Characteristics.

View this table:
[Table S2.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/T2)

Table S2. 
A. Performance of classifier models.

View this table:
[Table S3.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/T3)

Table S3. 
Lasso-selected features and coefficients of classifier models.

## Supplementary Files

**Supplementary File 1. Differential expression analyses**.

**Supplementary File 2. Gene set enrichment analyses**.

## Figure Legends

![Supp. Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/19/2020.05.18.20105171/F3.medium.gif)

[Supp. Figure 1.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/F3)

Supp. Figure 1. 
**A**. Breakdown of subjects with other pathogenic respiratory viruses identified by mNGS. Three patients had viral/viral co-infections: SARS-CoV-2/HRV (n=1) and RSV/HRV (n=2). CoV=Coronavirus, HRV=Human Rhinovirus, Flu=Influenza Virus, HMPV=Human Metapneumovirus, RSV=Respiratory Syncytial Virus, PIV=Parainfluenza Virus. **B**. Correlation of SARS-CoV-2 PCR Crossing Threshold (Ct) and mNGS Reads Per Million (rpM) Sequenced. Ct represents an average across the SARS-CoV-2 genomic loci assessed.

![Supp. Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/19/2020.05.18.20105171/F4.medium.gif)

[Supp. Figure 2.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/F4)

Supp. Figure 2. 
**A**. Gene expression scatter plots for the most significant interferon response genes induced by SARS-CoV-2, and the SARS-CoV-2 receptor gene ACE2. **B**. Gene set enrichment analysis for the direct comparison between COVID-19 and other viral ARIs. **C**. Gene expression scatter plots for selected interferon response genes in the leading edge of the interferon signaling gene set, showing lagging expression in SARS-CoV-2 compared to other viral ARIs.

![Supp. Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/19/2020.05.18.20105171/F5.medium.gif)

[Supp. Figure 3.](http://medrxiv.org/content/early/2020/05/19/2020.05.18.20105171/F5)

Supp. Figure 3. 
*In silico* estimation of cell type proportions. All pairwise comparisons were performed with a two-sided Mann-Whitney-Wilcoxon test followed by Bonferroni’s correction.

*   Received May 18, 2020.
*   Revision received May 18, 2020.
*   Accepted May 19, 2020.


*   © 2020, Posted by Cold Spring Harbor Laboratory

This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/)

## References

1.  1.Dong, E., Du, H. & Gardner, L. An interactive web-based dashboard to track COVID-19 in real time. Lancet Infect. Dis. 20, 533–534 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1473-3099(20)30120-1&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

2.  2.Wu, C. et al. Risk Factors Associated With Acute Respiratory Distress Syndrome and Death in Patients With Coronavirus Disease 2019 Pneumonia in Wuhan, China. JAMA Intern. Med. (2020) doi:10.1001/jamainternmed.2020.0994.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/jamainternmed.2020.0994&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32167524&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

3.  3.Langelier, C. et al. Integrating host response and unbiased microbe detection for lower respiratory tract infection diagnosis in critically ill adults. Proc. Natl. Acad. Sci. 115, E12353 LP-E12362 (2018).
    
    

4.  4.Plasschaert, L. W. et al. A single-cell atlas of the airway epithelium reveals the CFTR-rich pulmonary ionocyte. Nature 560, 377–381 (2018).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41586-018-0394-6&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30069046&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

5.  5.Tsalik, E. L. et al. Host gene expression classifiers diagnose acute respiratory illness etiology. Sci. Transl. Med. 8, 322ra11 LP-322ra11 (2016).
    
    

6.  6.Yang, Y. et al. Evaluating the accuracy of different respiratory specimens in the laboratory diagnosis and monitoring the viral shedding of 2019-nCoV infections. medRxiv 2020.02.11.20021493 (2020) doi:10.1101/2020.02.11.20021493.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoibWVkcnhpdiI7czo1OiJyZXNpZCI7czoyMToiMjAyMC4wMi4xMS4yMDAyMTQ5M3YyIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjAvMDUvMTkvMjAyMC4wNS4xOC4yMDEwNTE3MS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 

7.  7.Wölfel, R. et al. Virological assessment of hospitalized patients with COVID-2019. Nature (2020) doi:10.1038/s41586-020-2196-x.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41586-020-2196-x&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32235945&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

8.  8.Wang, W. et al. Detection of SARS-CoV-2 in Different Types of Clinical Specimens. JAMA 323, 1843–1844 (2020).
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

9.  9.Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. 102, 15545 LP–15550 (2005).
    
    

10. 10.Butler, D. J. et al. Shotgun Transcriptome and Isothermal Profiling of SARS-CoV-2 Infection Reveals Unique Host Responses, Viral Diversification, and Drug Interactions. bioRxiv 2020.04.20.048066 (2020) doi:10.1101/2020.04.20.048066.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYmlvcnhpdiI7czo1OiJyZXNpZCI7czoxOToiMjAyMC4wNC4yMC4wNDgwNjZ2NSI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzA1LzE5LzIwMjAuMDUuMTguMjAxMDUxNzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

11. 11.Huang, L. et al. Blood single cell immune profiling reveals the interferon-MAPK pathway mediated adaptive immune response for COVID-19. medRxiv 2020.03.15.20033472 (2020) doi:10.1101/2020.03.15.20033472.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoibWVkcnhpdiI7czo1OiJyZXNpZCI7czoyMToiMjAyMC4wMy4xNS4yMDAzMzQ3MnYxIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjAvMDUvMTkvMjAyMC4wNS4xOC4yMDEwNTE3MS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 

12. 12.Ziegler, C. G. K. et al. SARS-CoV-2 Receptor ACE2 Is an Interferon-Stimulated Gene in Human Airway Epithelial Cells and Is Detected in Specific Cell Subsets across Tissues. Cell (2020) doi:10.1016/j.cell.2020.04.035.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2020.04.035&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32413319&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

13. 13.Blanco-Melo, D. et al. Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19. Cell (2020) doi:10.1016/j.cell.2020.04.026.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2020.04.026&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32416070&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

14. 14.Lopez-Castejon, G. & Brough, D. Understanding the mechanism of IL-1β secretion. Cytokine Growth Factor Rev. 22, 189–195 (2011).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cytogfr.2011.10.001&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22019906&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000298202500002&link_type=ISI) 

15. 15.Cao, X. COVID-19: immunopathology and its implications for therapy. Nat. Rev. Immunol. 20, 269–270 (2020).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41577-020-0308-3&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

16. 16.Hiraishi, N., Ishida, Y.-I., Sudo, H. & Nagahama, M. WDR74 participates in an early cleavage of the pre-rRNA processing pathway in cooperation with the nucleolar AAA-ATPase NVL2. Biochem. Biophys. Res. Commun. 495, 116–123 (2018).
    
    

17. 17.Ndzinu, J. K., Takeuchi, H., Saito, H., Yoshida, T. & Yamaoka, S. eIF4A2 is a host factor required for efficient HIV-1 replication. Microbes Infect. 20, 346–352 (2018).
    
    

18. 18.Song, Z. et al. EIF4A2 interacts with the membrane protein of transmissible gastroenteritis coronavirus and plays a role in virus replication. Res. Vet. Sci. 123, 39–46 (2019).
    
    

19. 19.Lee, S.-Y. et al. FAM83A confers EGFR-TKI resistance in breast cancer cells and in mice. J. Clin. Invest. 122, 3211–3220 (2012).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1172/JCI60498&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22886303&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000308513100021&link_type=ISI) 

## References

1.  1.Pine, P. S. et al. Evaluation of the External RNA Controls Consortium (ERCC) reference material using a modified Latin square design. BMC Biotechnol. 16, 54 (2016).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12896-016-0281-x&link_type=DOI) 

2.  2.Deng, X. et al. Metagenomic sequencing with spiked primer enrichment for viral diagnostics and genomic surveillance. Nat. Microbiol. 5, 443–454 (2020).
    
    

3.  3.Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken 2. Genome Biol. 20, 257 (2019).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s13059-019-1891-0&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31779668&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

4.  4.Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100 (2018).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/bty191&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29750242&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

5.  5.Kalantar, K. L. et al. IDseq - An Open Source Cloud-based Pipeline and Analysis Service for Metagenomic Pathogen Detection and Monitoring. bioRxiv 2020.04.07.030551 (2020) doi:10.1101/2020.04.07.030551.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYmlvcnhpdiI7czo1OiJyZXNpZCI7czoxOToiMjAyMC4wNC4wNy4wMzA1NTF2MyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzA1LzE5LzIwMjAuMDUuMTguMjAxMDUxNzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

6.  6.Ramesh, A. et al. Metagenomic next-generation sequencing of samples from pediatric febrile illness in Tororo, Uganda. PLoS One 14, e0218318 (2019).
    
    

7.  7.Mayday, M. Y., Khan, L. M., Chow, E. D., Zinter, M. S. & DeRisi, J. L. Miniaturization and optimization of 384-well compatible RNA sequencing library preparation. PLoS One 14, e0206194 (2019).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0206194&link_type=DOI) 

8.  8.Venables, W. N., & Ripley, B. D. Modern Applied Statistics with S. (Springer, 2002).
    
    

9.  9.Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34, 525 (2016).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nbt.3519&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27043002&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

10. 10.Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4, 1521 (2015).
    
    

11. 11.Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47-e47 (2015).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gkv007&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25605792&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

12. 12.Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s13059-014-0550-8&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25516281&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

13. 13.Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. 102, 15545 LP–15550 (2005).
    
    

14. 14.Korotkevich, G., Sukhov, V. & Sergushichev, A. Fast gene set enrichment analysis. *bioRxiv* 60012 (2019) doi:10.1101/060012.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYmlvcnhpdiI7czo1OiJyZXNpZCI7czo4OiIwNjAwMTJ2MiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzA1LzE5LzIwMjAuMDUuMTguMjAxMDUxNzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

15. 15.Croft, D. et al. The Reactome pathway knowledgebase. Nucleic Acids Res. 42, D472-D477 (2014).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gkt1102&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24243840&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000331139800070&link_type=ISI) 

16. 16.Maechler, M. et al. robustbase: Basic Robust Statistics. (2020).
    
    

17. 17.Yohai, V. J. High Breakdown-Point and High Efficiency Robust Estimates for Regression. Ann. Stat. 15, 642–656 (1987).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1214/aos/1176350366&link_type=DOI) 

18. 18.Koller, M. & Stahel, W. A. Sharpening Wald-type inference in robust regression for small samples. Comput. Stat. Data Anal. 55, 2504–2515 (2011).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.csda.2011.02.014&link_type=DOI) 

19. 19.Renaud, O. & Victoria-Feser, M.-P. A robust coefficient of determination for regression. J. Stat. Plan. Inference 140, 1852–1862 (2010).
    
    

20. 20.Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457 (2015).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nmeth.3337&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25822800&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom) 

21. 21.Travaglini, K. J. et al. A molecular cell atlas of the human lung from single cell RNA sequencing. bioRxiv 742320 (2020) doi:10.1101/742320.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYmlvcnhpdiI7czo1OiJyZXNpZCI7czo4OiI3NDIzMjB2MiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzA1LzE5LzIwMjAuMDUuMTguMjAxMDUxNzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

22. 22.Friedman, J., Hastie, T. & Tibshirani, R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J. Stat. Software, Artic. 33, 1–22 (2010).
    
    

23. 23.Liaw, A. & Wiener, M. Classification and Regression by randomForest. R News 2, 18–22 (2002).
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1177/154405910408300516&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12973787&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F19%2F2020.05.18.20105171.atom)