Common and rare variant analyses combined with single-cell multiomics reveal cell-type-specific molecular mechanisms of COVID-19 severity

RefMap analysis of common variants uncovers cell-type-specific genetic basis of COVID-19 severity

We used the RefMap machine learning model (Methods) to identify the genomic regions and genes associated with severe COVID-19. Briefly, RefMap is a Bayesian network that combines genetic signals (e.g., allele Z-scores) with functional genomic profiling (e.g., ATAC-seq and ChIP-seq) to fine-map risk regions for complex diseases. With RefMap, we can scan the genome for functional regions in which disease-associated genetic variation is significantly shifted from the null distribution. The power of the RefMap model for gene discovery and recovery of missing heritability has been demonstrated in our recent work³². Here, to achieve cell-type-specific resolution within multicellular tissue, we modified RefMap to integrate single-cell multiomic profiling of human lungs with COVID-19 GWAS data (Fig. 1a). In particular, we obtained summary statistics (COVID-19 Host Genetics Initiative, Release 5, phenotype definition A2; 5,101 cases versus 1,383,241 population controls) from the largest GWAS study of COVID-19⁷, where age, sex, and 20 first principal components were included in the analysis as covariates. Severe COVID-19 was defined by the requirement for respiratory support or death attributed to COVID-19. Human lung single-cell multiomic profiling, including snRNA-seq and snATAC-seq, was retrieved from a recent study of healthy individuals³³. There are 19 cell types identified in both snATAC-seq and snRNA-seq profiles, including epithelial (alveolar type 1 (AT1), alveolar type 2 (AT2), club, ciliated, basal, and pulmonary neuroendocrine (PNEC)), mesenchymal (myofibroblast, pericyte, matrix fibroblast 1 (matrix fib. 1), and matrix fibroblast 2 (matrix fib. 2)), endothelial (arterial, lymphatic, capillary 1 (cap1), and capillary 2 (cap2)), and hematopoietic (macrophage, B-cell, T-cell, NK cell, and enucleated erythrocyte) cell types. We adopted these 19 cell types as the reference set within lung tissue throughout our study. Based on snATAC-seq peaks called in one or more of the 19 cell types to annotate functional regions, we used RefMap to identify disease-associated genomic regions from the COVID-19 GWAS data, which resulted in 6,662 1kb regions passing the 5% significance threshold (Q^+/--score>0.95, Methods; referred to as RefMap COVID-19 regions). These identified regions were further intersected with open chromatin in individual cell types based on corresponding snATAC-seq peaks, resulting in cell-type-specific RefMap regions (mean per cell type =1732.8, standard deviation (SD)=623.5; Fig. 1b, Supplementary Table 1). After removing RefMap regions present in more than one cell type (mean per cell type =121.2, SD=142.7), we observed only a weak correlation between the number of unique RefMap regions and the number of snATAC-seq peaks detected per cell type (Spearman ρ=0.40, P>0.05; Fig. 1b), indicating enrichment of genetic signals within certain cell types.

• Download figure
• Open in new tab

Figure 1. Common variant analysis of COVID-19 severity integrated with lung single-cell multiomics.

a, Schematic of the study design for fine-mapping cell-type-specific genes from COVID-19 GWAS (Panel 1). The diagram of the RefMap model is shown in Panel 2, where grey nodes represent observations, green nodes are local hidden variables, and pink nodes indicate global hidden variables (Methods). Cell-type-specific RefMap genes are mapped using single-cell multiomic profiling (Panel 3). Heritability (Panel 4), Mendelian randomization (Panel 5), and transcriptome analysis (Panel 6) validate the functional importance of RefMap genes, particularly for NK cells, in severe COVID-19. b, Total number and number of unique genomic regions containing genetic variation associated with severe COVID-19 for different cell types. c, Total number and number of unique genes implicated by genetic variation associated with severe COVID-19 for different cell types. d, Fraction of unique genomic regions and genes associated with severe COVID-19 for major cell types. e, Similarity between different cell types quantified by the overlap of RefMap genes. Gene set overlapping was calculated by the Jaccard index. f, RefMap regions overlap significantly with COVID-19-associated genetic variation in an independent COVID-19 GWAS study. cCRE: candidate cis-regulatory element. *: P<0.05.

Next, we sought to map the target genes of RefMap COVID-19 regions in a cell-type-specific manner. In particular, we identified the closest genes that are expressed in the corresponding cell type for individual RefMap regions (Methods). In total, we discovered 1,370 genes (referred to RefMap COVID-19 genes; mean per cell type =279.9 and SD=80.3; Fig. 1c, Supplementary Table 1) associated with the severe disease. Interestingly, hematopoietic cells have the largest number of unique RefMap regions and genes among all major cell types (Fig. 1d); for example, there is a significant enrichment of unique RefMap regions observed for hematopoietic cells versus epithelial cells (P=5.2e-03, odds ratio (OR)=1.15, Fisher’s exact test; Fig. 1d). This indicates a critical role of immune cells, which are primarily hematopoietic, in the development of severe COVID-19^16–19. To profile the cell-cell interactions underlying severe COVID-19 from a genetic perspective, we constructed a cell correlation matrix based on the overlap of RefMap genes between cell types (Fig. 1e). We discovered that the correlation is strongest between functionally related cells, demonstrating that the RefMap signal is consistent with known biology³³.

To replicate our findings, we obtained SNPs associated with severe COVID-19 from a GWAS for an entirely independent sample set (the 23andMe cohort, 15,434 COVID-19-positive cases and 1,035,598 population controls)⁵. The total union of RefMap regions is significantly enriched with SNPs associated with multiple COVID-19 phenotypes defined in this new sample set (mean P<5e-03, Fisher’s exact test; Fig. 1f, Supplementary Table 2; Methods). Specifically, the most significant enrichment is with SNPs associated with COVID-19 requiring respiratory support (mean P=4.68e-04, Fisher’s exact test; Fig. 1f). We further performed the enrichment analysis per cell type; only RefMap regions associated specifically with T cells and NK cells are significantly enriched with disease-associated SNPs across all measured COVID-19 phenotypes (mean P<0.05, Fisher’s exact test; Supplementary Table 2).

Heritability analysis and Mendelian randomization link NK cell function to COVID-19 severity

The LD score regression (LDSC)³⁴ has been used to measure the total SNP-based heritability (h²) from the GWAS study of severe COVID-19 (COVID-19 Host Genetics Initiative, Release 5, phenotype definition A2)⁷. Here, we examined the partitioning of SNP-based heritability for severe COVID-19 within RefMap genes (Methods). We discovered that the heritability of severe COVID-19 is significantly enriched for RefMap genes (OR=4.6, standard error (SE)=0.78, P=1.55e-07; Fig. 2a, Supplementary Table 3). We compared the proportion of SNP-based heritability captured by RefMap to other methods (Methods), including naïve GWAS⁷ and MAGMA³⁵. The proportion of heritability within naïve GWAS genes is 0.15 compared to 0.37 within MAGMA genes, but 0.77 within RefMap genes (Fig. 2b, Supplementary Table 3), representing a five-fold improvement in the recovered heritability based on RefMap over traditional methods. The proportion of SNP-based heritability for hospitalized COVID-19 (COVID-19 Host Genetics Initiative, Release 5, phenotype definition B2) within RefMap genes is 0.62 and within COVID-19 independent of severity (COVID-19 Host Genetics Initiative, Release 5, phenotype definition C2) it is 0.52 (Fig. 2b, Supplementary Table 3). In both cases the improvement in captured heritability based on RefMap compared to traditional methods is three-fold. Consistent with the design of our model, the recovered heritability is highest in severe COVID-19.

• Download figure
• Open in new tab

Figure 2. Severe-COVID-19-associated common variants are linked to NK cell function.

a, Heritability enrichment estimated by LDSC for different cell types. Enrichment was calculated as the proportion of total SNP-based heritability adjusted for SNP number. b, Proportion of SNP-based heritability associated with risk genes identified using RefMap or conventional methodology. c, d, e, Significant Mendelian randomization results for three exposures linked to severe COVID-19, including blood counts of (c) CD335+ CD314-, (d) CCR7-CD314-, and (e) CD314+ NK cells. f, Sensitivity analyses and robust tests for MR analyses (Methods). g, Comparative gene expression analysis of NK-cell RefMap genes in NKG2D+ and NKG2D-NK cells. Fold change was calculated as the ratio of gene expression levels in NKG2D+ NK cells to NKG2D-NK cells. The transcriptome was defined by all the expressed genes (with at least one UMI (unique molecular identifier)) in NK cells. Violin plots show the distributions of fold change values within each group, and boxplots indicate the median, interquartile range (IQR), Q1−1.5×IQR, and Q3+1.5×IQR. The red dashed line denotes the median value of fold change distribution for the transcriptome.

Next, we used cell-type-specific RefMap genes to determine which cell types are involved in the development of severe COVID-19. Specifically, we calculated the partitioned heritability per cell type within the severe COVID-19 GWAS (A2) and also within GWAS for hospitalized versus non-hospitalized COVID-19 (B2) and COVID-19 versus population (C2) (Methods). For severe COVID-19, of all 19 cell types tested, NK cells are the most enriched with SNP-based heritability (OR=8.87, SE=3.68, P=0.016; Fig. 2a, Supplementary Table 3). The same is also true for hospitalized COVID-19 (OR=10.57, SE=4.95, P=0.039), but not for COVID-19 irrespective of severity (OR=5.74, SE=3.09, P=0.077; Supplementary Table 3). Thus, we conclude that NK cell function is enriched with severe disease-associated genetic variation.

Two-sample Mendelian randomization (MR) facilitates identification of a causal relationship between an exposure and an outcome³⁶. We examined whether NK cell populations measured in the blood are causally related to severe COVID-19. In total, 46 GWAS measures of NK cell subtypes were identified³⁷ (Methods). After harmonizing exposure and outcome genetic instruments, we excluded tests with less than five SNPs (Methods). With MR, three exposures were shown to be causally related to severe COVID-19 after correcting for multiple testing (P<1e-03, multiplicative random effects (MRE), inverse-variance weighted (IVW)). All three exposures relate to NKG2D/CD314 expression on the cell surface, where a higher number of NKG2D/CD314-cells was linked to severe COVID-19 (A2) (Figs. 2c-d) and a higher number of NKG2D/CD314+ cells is protective (Fig. 2e). Evidence of genetic pleiotropy (MR PRESSO intercept not significantly different from zero, P>0.05; Fig. 2f) or instrument heterogeneity (P>0.05, Cochran’s Q test, and I²_GX >0.95; Fig. 2f) are not evident. Moreover, robust measures are significant for all three exposures (Fig. 2f). We also tested the identical phenotypes with alternative COVID-19 phenotype GWAS; we discovered that CD335+ CD314-cell counts are also causally associated with hospitalized COVID-19 (B2), and with COVID-19 independent of severity (C2) (Supplementary Fig. 1), but in each case the effect size is reduced compared to severe COVID-19 (A2). NKG2D/CD314 is a primary receptor responsible for NK cell activation³⁸ and in light of this, we conclude that severe COVID-19 is associated with a loss of NK cell cytotoxicity rather than a gain of function linked to NKG2D/CD314-cells.

Inspired by our MR analysis, we further tested if the expression of RefMap genes reflects a functional difference between NKG2D/CD314+ and NKG2D/CD314-cells. We examined the expression levels of NK-cell RefMap genes based on scRNA-seq data from healthy lungs³⁹, and discovered that RefMap gene expression is higher within NKG2D/CD314+ cells than NKG2D/CD314-cells (P=0.036, one-tailed Wilcoxon rank-sum test; Fig. 2g). This further supports the functional significance of RefMap genes in COVID-19 severity and associates the genetic risk of severe COVID-19 directly with NK cell activity.

Transcriptome analysis supports the functional significance of RefMap genes in health and severe COVID-19

To link RefMap COVID-19 genes to the underlying biology, we first performed functional enrichment analyses based on gene ontology (GO) and KEGG pathways⁴⁰ (Figs. 3a and 3b, Supplementary Tables 4 and 5). We observed that RefMap NK-cell genes are enriched with pathways and ontologies related to intra- and intercellular signalling important for NK cell activation, including “Phospholipase D signalling pathway”⁴¹, “Antigen processing and presentation”, “regulation of small GTPase mediated signal transduction (GO:0051056)”⁴², and “regulation of intracellular signal transduction (GO:1902531)” (adjusted P<0.1; Figs. 3a and 3b, Supplementary Tables 4 and 5). This is consistent with the hypothesis that COVID-19 severity is determined by failed activation of NK cells. Furthermore, the pathway with the highest enrichment is “human immunodeficiency virus (HIV) 1 infection” (adjusted P=3e-04; Fig. 3b). Since HIV-1 works to suppress NK cell activation⁴³, and NK cell function has been associated with an effective immune response to HIV⁴⁴, this result is also consistent with a role of NK cells in severe COVID-19. Other cell-type-specific RefMap gene lists are also enriched with relevant biological pathways. For example, AT2-cell genes are linked to pathways associated with viral infection such as ‘human papillomavirus infection’ and ‘viral carcinogenesis’ (adjusted P<0.1; Supplementary Table 5), which is consistent with the established role of AT2 cells as the initial site of SARS-CoV-2 entry into host cells⁴⁵. T-cell genes are enriched with ‘IL17-signalling pathway’ (adjusted P=0.021; Supplementary Table 5) which is interesting in light of previous literature highlighting the production of IL-17 by T cells from COVID-19 patients as a potential therapeutic target⁴⁶.

• Download figure
• Open in new tab

Figure 3. Functional enrichment and transcriptome analyses of RefMap COVID-19 genes.

a, Gene Ontology (GO) terms that are significantly enriched in cell-type-specific RefMap gene lists corresponding to hematopoietic cell types; only terms with adjusted P<0.05, OR>3, and character number<60 are visualized. b, KEGG Pathways that are significantly enriched in cell-type-specific RefMap gene lists corresponding to hematopoietic cell types; only terms with adjusted P<0.05, OR>5, and character number<50 are visualized. c, Gene expression analysis of RefMap genes across different cell types in healthy lungs. The transcriptome was defined as the total set of expressed genes for each cell type (Methods). Violin plots show the distributions of log expression levels within each group, and point plots indicate the median and IQR. d, Overlap between cell-type-specific RefMap regions and H3K27ac and H3K4me3 ChIP-seq peaks from ENCODE lung and immune cell samples. Z-scores calculated by regionR⁴⁹ (1,000 permutations) were normalized into the 0-1 range for visualization. e, f, Comparative gene expression analysis of cell-type-specific RefMap genes in severe COVID-19 patients versus moderately affected patients based on scRNA-seq datasets from (e) Ren et al. and (f) Liao et al., respectively. The Z-score of Wilcoxon rank-sum test was used to indicate the gene expression change between severe and moderate patient groups, where a positive value means higher gene expression in severe patients. The Benjamini-Hochberg (BH) procedure was used to calculate FDRs throughout the study. Violin plots show the distribution of gene expression changes within each group, and boxplots indicate the median, IQR, Q1−1.5×IQR, and Q3+1.5×IQR. *: FDR<0.1. +: FDR<0.01.

Next, we investigated the baseline expression pattern of RefMap genes in healthy lungs. In particular, we calculated mean expression levels of genes in different cell types based on the lung snRNA-seq data from Wang et al.³³, and then compared the expression of RefMap genes with the total set of expressed genes in each cell type. Interestingly, although the gene expression level was not an input to the RefMap model, RefMap genes are expressed at a higher level compared to the background transcriptome in all 19 cell types, including immune and epithelial cells (false discovery rate (FDR)<0.1, one-tailed Wilcoxon rank-sum test; Fig. 3c) with the exception of pericytes (FDR=0.11, Z-score=1.25); it is interesting to note that pericytes may be downstream in the pathogenesis of COVID-19 because they are protected by an endothelial barrier⁴⁷. This supports the functional significance of RefMap genes across multiple cell types in healthy human lungs. As a negative control, we performed a similar expression comparison between non-developmental genes and all expressed genes in lungs, which yielded no significant difference (Supplementary Fig. 2; Methods). In summary, our transcriptome analyses indicate that RefMap genes are expressed above background in relevant cell types, supporting their important role in lung function.

To obtain further insights into the function of RefMap COVID-19 regions, we tested whether RefMap regions are enriched with cell-type-specific candidate cis-regulatory elements (cCREs, or enhancers and promoters) defined by H3K27ac and H3K4me3. We obtained cCREs for lung tissues and primary cells from the ENCODE project⁴⁸ (Supplementary Table 6), and examined the overlap between those cCREs and RefMap regions by permutation test⁴⁹. We discovered that RefMap regions specific to immune cells (e.g., T cells, B cells, and NK cells) are significantly enriched with cCREs in corresponding cell types (FDR<0.1; Fig. 3d). For other cell types, RefMap regions are generally enriched with cCREs from lung tissue (FDR<0.1; Fig. 3d). These observations are consistent with an important role of RefMap regions in the regulation of gene expression. Moreover, the enrichment with cCREs across a variety of individuals and datasets supports the generalizability of the genetic architecture defined by RefMap. A similar association between genome-wide snATAC-seq peaks and cCREs was also observed (Supplementary Fig. 3).

We have shown that RefMap COVID-19 regions are enriched within promoters and enhancers responsible for regulating gene expression in healthy lung tissue. On this basis, we hypothesized that genetic variation within RefMap regions would alter the expression of corresponding target genes in the context of severe disease. Specifically, we proposed that RefMap genes would be expressed at a lower level in lung tissue from severe COVID-19 patients than moderately affected patients. To validate this hypothesis, we obtained scRNA-seq data from the respiratory system for a large COVID-19 cohort²³, including 12 bronchoalveolar lavage fluid (BALF) samples, 22 sputum samples, and 1 sample of pleural fluid mononuclear cells (PFMCs) from 27 severely and 8 mildly affected patients. Severity was classified based on the World Health Organization (WHO) guidelines (https://www.who.int/publications/i/item/WHO-2019-nCoV-clinical-2021-1). For individual cell types, we compared the expression level of RefMap genes in severe patients versus moderately affected patients (Methods). Compared to the background transcriptome, we observed that RefMap genes are relatively lower expressed in severe patients in corresponding cell types than in moderate patients (FDR<0.01, one-tailed Wilcoxon rank-sum test; Fig. 3e), supporting the functional significance of RepMap genes in severe COVID-19. As a replication experiment, we carried out a similar analysis based on an independent COVID-19 scRNA-seq dataset²², including 9 BALF samples from 6 severe patients and 3 moderate patients (Methods). The lower expression of RefMap genes in severe patients is consistent across multiple cell types (FDR<0.01, one-tailed Wilcoxon rank-sum test; Fig. 3f). Altogether, these transcriptome-based orthogonal analyses are consistent with the hypothesis that identified cell-type-specific RefMap genes are functionally linked to COVID-19 severity.

PULSE analysis of rare variants enables population-independent prediction of COVID-19 severity

RefMap utilizes common genetic variation profiled in GWAS. Biological dysfunction can also be determined by rare variants and therefore we assessed whether there is a significant burden of severe-COVID-19-associated rare variants within RefMap genes. First, we employed a standard methodology using SKAT⁵⁰ rare-variant burden testing applied to whole-exome sequencing (WES) data from the GEN-COVID cohort⁵¹, including non-elderly patients who suffered severe COVID-19 requiring respiratory support, and individuals who suffered non-severe COVID-19 not requiring hospitalization (Methods). No individual gene is enriched with significant genetic burden after adjusting for multiple testing (Supplementary Fig. 4). This is true whether we tested genome-wide or only for RefMap COVID-19 genes. However, for one subset of cell-type-specific RefMap genes, the median P-value was lower than expected: NK cells (P<0.05, permutation test; Fig. 4a; Methods). This result from the analysis of rare genetic variation in an independent cohort is convergent with our common variant analysis, highlighting NK cell biology as a critical determinant of COVID-19 severity.

• Download figure
• Open in new tab

Figure 4. Rare variant analysis informs individual risk of critical illness of COVID-19.

a, Enrichment analysis of cell-type-specific RefMap COVID-19 genes with rare variants using SKAT burden testing. The red dashed line indicates P=0.05. b, Schematic of the study design for our rare variant analysis based on PULSE. We examine two independent cohorts in which rare variants were profiled by different technologies: whole-exome sequencing (WES) and whole-genome sequencing (WGS) (Panel 1). Variants are annotated using ANNOVAR (Panel 2) and encoded as input for the PULSE model (Panel 3, Methods), where grey nodes are observations and pink nodes represent hidden variables. PULSE is trained to differentiate cases and controls (Panel 4), where the gene weights are useful for gene discovery (Panel 5). Functional characterization of risk genes is performed based on scRNA-seq and PPIs (Panel 6). c, Receiver operating characteristic (ROC) curves of different models, including PULSE, age+sex, and integrative models, in the 5-fold cross-validation. Solid lines represent the mean values, and the grey area indicates the standard errors. d, Summary statistics of the VA COVID-19 cohort. e, AUROC (area under the receiver operating characteristics) scores of predictions in multiple test datasets. Prediction performance is shown for PULSE, age+sex, and integrative models. f, Comparison of prediction sensitivity between PULSE and age+sex models. EHRs: electronic health records.

Traditional rare-variant burden testing failed to identify any enrichment of COVID-19-associated variants within a single gene although there is significant enrichment in the group of 236 NK-cell RefMap genes. This suggests that traditional burden testing is underpowered when applied to the GEN-COVID dataset. We decided to develop a new method with increased sensitivity. Here we present PULSE (probabilistic burden analysis based on functional estimation), a discriminative Baysian network that integrates functional annotations of rare variants to model the relationship between genotype and phenotype (Fig. 4b; Methods, Supplementary Notes). In particular, PULSE combines multiple predictions of functional effects for different types of variants, including missense, nonsense, splicing-site, and small insertion-deletion (indel) mutations (Supplementary Table 7). After aggregating those functional scores for individual genes, PULSE learns the importance of different annotations and genes from the training data and maps the phenotype from the genotype in a bilinear form (Fig. 4b; Methods). With PULSE as a discovery-by-prediction strategy, we are able to (i) predict individual phenotypes from personal genotypes and (ii) discover phenotype-associated genes by model interpretation.

We applied PULSE to study rare genetic variants associated with COVID-19 severity based on the GEN-COVID cohort of whole-exome sequencing from 1,339 COVID-19 patients with 5 severity gradings⁵² (Fig. 4b, Supplementary Table 8; Methods). After quality controls (QCs), we constructed a discovery cohort (training dataset) of non-elderly European (EUR) adults (age >30 and <60 years) who were critically ill (cases, n=109) or not hospitalized (controls, n=269) (Methods). There is no significant age difference between cases and controls after filtration (P=0.29, Wilcoxon rank-sum test; Supplementary Fig. 5). We then performed genome annotation⁵³ and feature engineering (Methods), where only rare variants (i.e., absent from the EUR cohort within the 1000 Genomes Project Phase 3⁵⁴) in autosomes were utilized for downstream analysis. To test the prediction performance of PULSE, we first performed 5-fold cross-validation (CV) based on the GEN-COVID discovery cohort, where a mean AUROC (area under the receiver operating characteristics) of 0.631 was achieved (SE=0.062; Fig. 4c). This demonstrates the predictability of COVID-19 severity from personal genomes. The AUROC scores (0.629±0.073) of a logistic regression model built from patient age and sex information are comparable to the PULSE genetic model (Fig. 4c). However, combining scores (by averaging) of PULSE and age+sex produced a further improvement in prediction performance (AUROC=0.653±0.072; Fig. 4c), demonstrating that host genetics is relatively independent of the effect of age and sex on disease severity. We note that since we removed the age bias in the discovery cohort for genetic concentration, the largest contribution in the age+sex model came from the sex information (model coefficients: 1.33±0.052 for sex versus 0.001±0.004 for age; Supplementary Fig. 6). We also note that in our PULSE model only autosomal variants were considered to remove the effect of sex in genetic modelling.

To further validate the prediction power of PULSE, we analyzed whole-genome sequencing (WGS) data of an independent cohort from the Veterans Health Administration (VA), consisting of 590 COVID-19 patients with variable disease severity (Fig. 4d, Supplementary Table 9; Methods). Extensive QCs (without filtering based on ancestry) resulted in 571 genomes (Methods). Genome annotation and feature engineering were conducted as for the GEN-COVID cohort. Similarly, to remove the effect of age, we focused on non-elderly adults (age >30 and <65 years) who were critically ill or not hospitalized, yielding 243 samples (24 cases and 219 controls). In this analysis, we relaxed the upper threshold of age from 60 to 65 years to include more samples in testing. The PULSE model trained on the whole GEN-COVID cohort was applied to predict severity within the VA EUR samples (14 cases versus 125 controls). We found that PULSE succeeded in predicting severe disease solely from personal genomes for this independent cohort with an AUROC of 0.675 (Fig. 4e).

Next, we asked if the prediction accuracy is generalizable across different populations. We constructed a test set of non-EUR non-elderly adults (age >30 and <60 years) that passed all other QC criteria within the GEN-COVID cohort, resulting in 12 cases (critically ill) and 6 controls (not hospitalized). The PULSE model trained on EUR samples was then applied to this non-EUR dataset, yielding AUROC of 0.667, which is comparable to the prediction solely based on age and sex (AUROC=0.722; Fig. 4e). Combining two scores further increased the prediction performance (AUROC=0.799; Fig. 4e). Furthermore, we applied the same trained PULSE model to predict severe disease for African (AFR) individuals (10 cases versus 92 controls) within the VA cohort. Similarly, we discovered that PULSE succeeded in the cross-population prediction with an AUROC of 0.784 (Fig. 4e). A similar result was observed for the whole VA dataset with mixed populations (AUROC=0.716; Fig. 4e). These results demonstrate the prediction power of PULSE and suggest that the rare-variant genetic architecture of COVID-19 severity is conserved across multiple populations. Importantly, we observed that the prediction in the VA cohort based on just age and sex information trained on the GEN-COVID cohort is inferior to PULSE (AUROC=0.655, 0.474, and 0.577 for EUR, AFR, and all samples, respectively; Fig. 4e). This may be linked to the different sex distribution with fewer females in the VA cohort (Fig. 4d, Supplementary Figs. 7 and 8), but is further evidence of the robustness of host genetic signals in determining COVID-19 severity and demonstrates that the PULSE prediction is independent of age and sex.

We investigated additional performance measures including sensitivity and specificity based on different cutoffs. Importantly, although the specificity scores are comparable between PULSE and age+sex models cross cutoffs (Supplementary Fig. 9), we discovered that PULSE yielded a significantly higher sensitivity (median values: 0.857 versus 0.688, 0.900 versus 0.525, and 0.875 versus 0.560 for EUR, AFR, and all VA samples, respectively; Fig. 4f). High sensitivity is important for the clinical application of severity prediction to guide the identification of at-risk individuals. Predictions for the GEN-COVID non-EUR samples yielded similar sensitivity and specificity (Supplementary Fig. 10).

Common and rare variant analyses of severe COVID-19 converge on NK cell function

The trained PULSE model assigns a weighting to individual genes as a measure of association between gene function and severe COVID-19 (Supplementary Fig. 11; Methods), where a larger weight indicates a higher gene mutation burden in severe patients. To test for the convergence between our common and rare variant analyses, we compared the absolute values of model weights for RefMap genes per cell type with all genes considered by the PULSE model. After correcting for multiple testing, we concluded that for all cell types, RefMap genes tend to have weights with larger absolute values, indicating an association with severe COVID-19 (FDR<0.01, one-tailed Wilcoxon rank-sum test; Fig. 5a). Common and rare genetic variations are largely independent^55–57, and therefore, this convergence of common and rare variant signals indicates shared biology underlying severe disease. Among all cell types, club-cell RefMap genes are the most enriched with PULSE genes (FDR<1e-05, Z-score=5.33; Fig. 5a). Interestingly, of hematopoietic cells, NK cell genes are the most enriched with PULSE genes (FDR=1.1e-03, Z-score=3.16; Fig. 5a), consistent with our previous conclusion that NK cells are an essential component of the immune response against SARS-CoV-2.

• Download figure
• Open in new tab

Figure 5. Transcriptome analysis of PULSE COVID-19 genes.

a, Analysis of convergence between PULSE and RefMap COVID-19 genes. The Z-scores were calculated per cell-type by Wilcoxon rank-sum test of the difference in PULSE weights between RefMap genes and the background transcriptome. Non-zero Z-scores indicate biological overlap between common and rare variant architectures detected by RefMap and PULSE, respectively. b, Gene expression analysis of PULSE genes across different cell types in healthy lungs. The transcriptome was defined as the total set of expressed genes for each cell type (Methods). Violin plots show the distributions of log expression levels within each group, and point plots indicate the median and IQR. c, d, Comparative gene expression analysis of cell-type-specific PULSE genes in severe COVID-19 patients versus moderate patients based on scRNA-seq datasets from (c) Ren et al. and (d) Liao et al., respectively. The Z-score of Wilcoxon rank-sum test was used to indicate the gene expression change between severe and moderate patient groups. Violin plots show the distribution of gene expression changes within each group, and boxplots indicate the median, IQR, Q1−1.5×IQR, and Q3+1.5×IQR. *: FDR<0.1. +: FDR<0.01.

To further validate the importance of genes captured by PULSE for severe COVID-19, we identified 657 genes (referred to as PULSE COVID-19 genes) based on model weights in the top 5% from all genes (Supplementary Table 10). We re-examined our SKAT burden analysis results for the GEN-COVID cohort and observed that PULSE genes are significantly enriched with severe-disease-associated rare variants (median P<1e-5, permutations test). Similar enrichment was confirmed based on an independent rare-variant burden analysis from Regeneron¹², where the PULSE genes are significantly enriched with Regeneron genes implicated in the analysis of severe COVID-19 versus non-hospitalized COVID-19 (n=68 genes; P=0.02, OR=2.9, Fisher’s exact test; Methods).

To gain further insights into the cell-type-specificity of PULSE genes, we investigated their expression levels in healthy lungs per cell type. We confirmed the function of PULSE genes across different cell types by observing their non-random overlapping with lung snATAC-seq peaks³³ (FDR<0.1, permutation test). Furthermore, the expression levels of PULSE genes measured by scRNA-seq were examined, where we found that PULSE genes are higher expressed in B cells, club cells, lymphatics, matrix fibroblast 1, and NK cells (FDR<0.1, one-tailed Wilcoxon rank-sum test; Fig. 5b). This supports the functional importance of PULSE genes in lung function.

PULSE genes carry a higher mutation burden in severe COVID-19 and therefore we hypothesized that loss of function of PULSE genes leads to severe symptoms. To validate this, we analyzed the expression levels of PULSE genes based on the scRNA-seq data of COVID-19 patients²³. Consistent with our hypothesis, we observed a down-regulation of PULSE genes in severe disease compared to moderate disease across B cells, ciliated cells, macrophages, NK cells, and T cells (FDR<0.1, one-tailed Wilcoxon rank-sum test; Fig. 5c). A similar analysis in another cohort²² led to the same conclusion for macrophages, NK cells, and T cells (FDR<0.1, one-tailed Wilcoxon rank-sum test; Fig. 5d). Our transcriptome study demonstrates the functional role of PULSE genes in severe disease across multiple cell types. Notably, among all the cell types we investigated, only NK cells are consistently associated with severe COVID-19 across all observations. This supports the conclusion of our common variant analysis, suggesting that NK cells are vital determinants of COVID-19 severity.

Systems analysis implicates association of NK cell activation with COVID-19 severity

All of our analyses have suggested that NK cell dysfunction is a determinant of COVID-19 severity. To obtain a comprehensive landscape of NK cell biology underlying severe COVID-19, we examined the function of NK-cell genes identified by either RefMap or PULSE (377 genes; Supplementary Table 11). Indeed, genes do not function in isolation^58,59 and therefore, rather than examining individual genes, we mapped NK-cell genes to the global protein-protein interaction (PPI) network and then inspected functional enrichment of COVID-19-associated network modules.

In particular, we extracted high-confidence (combined score >700) PPIs from STRING v11.0⁶⁰, which include 17,161 proteins and 839,522 protein interactions. To eliminate the bias of hub genes⁶¹, we performed the random walk with restart algorithm over the raw PPI network to construct a smoothed network based on edges with weights in the top 5% (Supplementary Table 12; Methods). Next, this smoothed PPI network was decomposed into non-overlapping subnetworks using the Leiden algorithm⁶². This process yielded 1,681 different modules (Supplementary Table 13), in which genes within modules are densely connected but sparsely connected with genes in other modules.

NK-cell COVID-19 genes were mapped to individual modules, and four modules were found to be significantly enriched with NK-cell genes: M237 (n=471 genes; FDR<0.1, hypergeometric test; Fig. 6a), M1164 (n=396 genes; FDR<0.1, hypergeometric test; Fig. 6b), M1311 (n=14 genes; FDR<0.1, hypergeometric test), and M1540 (n=226 genes; FDR<0.1, hypergeometric test; Fig. 6c) (Supplementary Table 13). We excluded M1311 from our downstream analysis due to its limited size and lack of functional enrichment.

• Download figure
• Open in new tab

Figure 6. Network analysis of NK-cell genes identified in common and rare variant analyses.

a, b, c, Three PPI network modules, including (a) M237, (b) M1164, and (c) M1540, are significantly enriched with NK-cell genes identified in either common or rare variant analysis. Blue nodes represent NK-cell genes and yellow nodes indicate other genes within each module. Edge thickness is proportional to STRING confidence score (>700). d, e, f, Gene Ontology (GO) terms that are significantly enriched in modules (d) M237, (e) M1164, and (f) M1540. Selected terms are shown for visualization and the complete lists can be found in Supplementary Tables 14 and 15. g, Gene expression analysis of module genes in NK cells. The transcriptome was defined as the total set of expressed genes in NK cells (Methods). Violin plots show the distributions of log expression levels within each group, and boxplots indicate the median, IQR, Q1−1.5×IQR, and Q3+1.5×IQR. The red dashed line indicates the median expression level of the transcriptome. h, i, Comparative gene expression analysis of module genes in severe COVID-19 patients versus moderate patients based on scRNA-seq datasets from (h) Ren et al. and (i) Liao et al., respectively. The Z-score of Wilcoxon rank-sum test was used to indicate the gene expression change between severe and moderate patient groups. Violin plots show the distribution of gene expression changes within each group, and boxplots indicate the median, IQR, Q1−1.5×IQR, and Q3+1.5×IQR. The red dashed line indicates the median expression change of the transcriptome. +: FDR<0.01. GOBP: gene ontology biological process.

Functionally, M237, M1164, and M1540 are all enriched with gene expression linked to NK cells (P<0.05, Human Gene Atlas), demonstrating their specificity in the NK cell function. Moreover, these three modules relate to different stages of NK cell activation. M237 is enriched with GO/KEGG terms including ‘mRNA processing (GO:0006397)’ and ‘Spliceosome’, which are important for the transcriptional response involved in NK cell activation (adjusted P<0.1; Fig. 6d, Supplementary Tables 14 and 15). M1164 is enriched with GO/KEGG terms linked to intracellular signalling (e.g., ‘regulation of small GTPase mediated signal transduction (GO:0051056)’), including pathways (e.g., ‘Rap1 signalling pathway’) key for NK cell activation (adjusted P<0.1; Fig. 6e, Supplementary Tables 14 and 15). M1540 is highly enriched with GO/KEGG terms linked to type I interferon signalling (e.g., ‘type I interferon signalling pathway (GO:0060337)’ and ‘Antigen processing and presentation’) (adjusted P<0.1; Fig. 6f, Supplementary Tables 14 and 15). In summary, the functional enrichment of M237, M1164, and M1540 genes includes extracellular, cytoplasmic, and nuclear processes necessary for NK cell activation; thus the genetic architecture we have discovered places NK cell activation upstream in determining severe COVID-19.

To further characterize the function of the identified NK-cell modules, we investigated the expression of module genes based on scRNA-seq data from healthy and diseased lung tissues. Genes in all three modules are relatively over-expressed in NK cells of healthy lungs³³ than the background transcriptome (FDR<0.01, one-tailed Wilcoxon rank-sum test; Fig. 6g). In contrast, in lung tissues infected with SARS-CoV-2, we observed a down-regulation of M237 and M1540 genes in NK cells of severe disease²³ (FDR<0.01, one-tailed Wilcoxon rank-sum test; Fig. 6h). M1164 genes are also down-regulated in NK cells from severe COVID-19 patients in another cohort²² (FDR<0.01, one-tailed Wilcoxon rank-sum test; Fig. 6i) along with M237 and M1540 genes. These results are consistent with our previous findings and functionally link the modules we have detected to NK cell biology in the context of severe COVID-19.

The RefMap model

Allele Z-scores were calculated as Z=b/se, where b and se are effect size and standard error, respectively, as reported by the COVID-19 GWAS⁷ (COVID-19 Host Genetics Initiative, Release 5, phenotype definition A2, EUR only) where the sample age, sex, and ancestry information were included as covariates. Given Z-scores and lung snATAC-seq peaks, we aim to identify functional genomic regions in which the Z-score distribution is significantly shifted from the null distribution. Suppose we have K 1Mb linkage disequilibrium (LD) blocks, where each LD block contains J_k(k=1, …, K) 1kb regions and each region harbors I_j,k(j=1, …, J_k, I_j,k>0) SNPs, the Z-scores follow a multivariate normal distribution, i.e., in which the Z-score of the i-th SNP in the j-th region of the k-th block is denoted as z_i,j,k (i=1, …, I_j,k) and u_k are the effect sizes that can be expressed as In addition, in Eq. (1) represents the in-sample LD matrix comprising of the pairwise Pearson correlation coefficients between SNPs within the k-th block, where I_k is the total number of SNPs calculated by . Here, since we have no access to the individual-level data, we used EUR samples from the 1000 Genomes Project (Phase 3) to estimate Σ_k, yielding the out-sample LD matrix. A modified Cholesky algorithm⁷¹ was used to get a symmetric positive definite (SPD) approximation of the LD matrix.

Further, we assume u_i,j,k (i=1, …, I_j,k) are independent and identically distributed (i.i.d.), following a normal distribution given by where the precision λ_j,k follows a Gamma distribution, i.e., Moreover, to characterize the shift of the expectation in Eq. (3) from the null distribution, we model m_j,k by a three-component Gaussian mixture model given by where the precisions follow and v_-1 and v₊₁ are non-negative variables measuring the absolute values of effect size shifts for the negative and positive components, respectively.

To impose non-negativity over v_-1 and v₊₁, we adopt the rectification nonlinearity technique proposed previously⁷². In particular, we assume v_-1 and v₊₁ follow in which the rectified Gaussian distribution is defined via a dumb variable. In particular, we first define v_-1 and v₊₁ by which guarantees that v_-1 and v₊₁ are non-negative. The dump variable r_-1 and r₊₁ follow the Gaussian distributions given by where m_± and λ_± follow the Gaussian-Gamma distributions, i.e., The indicator variables in Eq. (5) denote whether that region is disease-associated or not. Indeed, we define the region to be disease-associated if or , and to be non-associated otherwise. To simplify the analysis, we put a symmetry over and , and define the distribution by Furthermore, the probability parameter π_j,k in Eq. (15) is given by Where σ(.) is the sigmoid function, S_j,k is the vector of epigenetic features for the j-th region in the k-th LD block, and the weight vector w follows a multivariate normal distribution, i.e., and ⋀ follows In this study, the epigenetic feature s_j,k was calculated as the overlapping ratios of that region with the snATAC-seq peaks detected in any of the cell types in healthy human lungs.

Based on the model defined in Eqs. (1) to (18), we are interested in calculating the posterior probability p(T | Z, S), where the mean-field variational inference (MFVI)⁷³ was adopted to solve the intractability. More technical details, including a coordinate ascent-based inference algorithm, can be found in our previous work³².

In this study, we ran the MFVI algorithm per chromosome to accelerate the computation. The Q⁺- and Q^--scores were defined as q(t^{(+ 1)} = 1) and q(t^{(− 1)} = 1), respectively, and we also defined the Q-score as Q=Q⁺+Q^-. RefMap regions were identified by Q⁺- or Q^--score >0.95.

Mapping cell-type-specific genes from RefMap regions

For each cell type within lung tissue, we defined cell-type-specific RefMap regions as the overlap between RefMap regions and the total set of snATAC-seq peaks detected in that cell type (Supplementary Table 1). Cell-type-specific RefMap genes were then identified if the extended gene body (i.e., the region up to 10kb either side of the annotated gene body) overlapped with any of the cell-type-specific regions. To get the final gene lists, RefMap genes were further filtered based on their expression levels. In particular, with the lung snRNA-seq data³³, we defined expressed genes in each cell type as those with Seurat⁷⁴ log-normalized value>0.6931. In addition, we note that there are non-adult samples (∼30 weeks gestation and ∼3 years) sequenced in the single cell profiling data³³. To remove the bias towards lung development, we first calculated the fold change of gene expression levels between the adult sample (∼30 years) and non-adult ones, and defined non-developmental genes (nDG) as those with FC>1.5. Only RefMap genes that were identified as expressed and non-developmental in each cell type were kept for downstream analysis (Supplementary Table 1).

Validation of RefMap COVID-19 regions in the 23andMe dataset

We calculated the overlap of total RefMap regions and of cell-type-specific RefMap regions with genomic regions shown to contain COVID-19-associated SNPs (P<1e-04) based on the GWAS of an independent cohort recruited by 23andMe⁵. To determine whether the observed overlap is statistically significant, we examined the average overlap with ten sets of control regions of equivalent length to RefMap regions. Control regions were +/-1Mb-5Mb distant from the RefMap regions⁷⁵.

Heritability analysis

We used LD score regression (LDSC)³⁴ to calculate overall heritability for severe COVID-19 (A1), hospitalized COVID-19 (B2), and COVID-19 overall (C2), respectively. Heritability partitioning within genes identified by traditional methods and within cell-type-specific RefMap genes was performed as previously described⁷⁶. Briefly, for all gene lists, we examined the proportion of total SNP-based heritability carried by SNPs +/-100kb from the transcription start site (TSS) of each gene in the list. Enrichment was calculated by comparing the ratio of partitioned heritability to the quantity of genetic materials.

Mendelian randomization

In total, 46 GWAS measures of NK cell subtypes were identified from the IEU Open GWAS Project, including “prot-a-180”, “met-b-124”, “met-b-245”, “met-b-242”, “prot-c-5244_12_3”, “met-b-237”, “met-b-258”, “prot-a-1669”, “prot-c-2917_3_2”, “met-b-246”, “prot-a-1671”, “met-b-249”, “met-b-140”, “met-b-240”, “prot-a-3159”, “prot-c-5104_57_3”, “prot-c-3056_11_1”, “prot-a-13”, “prot-a-3160”, “met-b-123”, “met-b-250”, “met-b-239”, “met-b-120”, “met-b-154”, “prot-a-3162”, “met-b-247”, “met-b-251”, “met-b-238”, “met-b-243”, “prot-a-2487”, “met-b-244”, “prot-c-2734_49_4”, “met-b-153”, “prot-a-3161”, “prot-c-3003_29_2”, “met-b-248”, “prot-a-1674”, “prot-a-1675”, “met-b-152”, “met-b-122”, “met-b-121”, “prot-a-1670”, “prot-c-5424_55_3”, “met-b-252”, “prot-a-3233” and “met-b-241”^37,77,78. Exposure SNPs or instrumental variables (IVs) are chosen based on an arbitrary P-value cutoff^79,80. A cutoff that is too low will lose informative instruments, but a cutoff that is too high could introduce non-informative instruments. We chose to set the cutoff at 5e-06 in line with our previous work⁸¹. We employed a series of sensitivity analyses to ensure that our analysis was not confounded by invalid IVs. Identified SNPs were clumped for independence using PLINK clumping in the TwoSampleMR tool⁸². A stringent cutoff of R²≤0.001 and a window of 10,000kb were used for clumping within a European reference panel. Where SNPs were in LD, those with the lowest P-value were retained. SNPs that were not present in the reference panel were excluded. Where an exposure SNP was unavailable in the outcome dataset, a proxy with a high degree of LD (R²≥0.9) was identified in LDlink within a European reference population⁸³. Where a proxy was identified to be present in both datasets, the target SNP was replaced with the proxy in both exposure and outcome datasets in order to avoid phasing issues⁸⁴. Where a SNP was not present in both datasets and no SNP was available in sufficient LD, the SNP was excluded from the analysis. The effects of SNPs on outcomes and exposures were harmonized in order to ensure that the beta values were signed with respect to the same alleles. For palindromic alleles, those with minor allele frequency (MAF) > 0.42 were omitted from the analysis in order to reduce the risk of errors due to strand issues⁸⁴.

The MR measure with the greatest power is the inverse-variance weighted (IVW) method, but this is contingent upon the exposure IV assumptions being satisfied ⁸⁵. With the inclusion of a large number of SNPs within the exposure IV, it is possible that not all variants included are valid instruments and therefore, in the event of a significant result, it is necessary to include a range of robust methods which provide valid results under various violations of MR principles at the expense of power⁸⁶. Robust methods applied in this study include MR-Egger, MR-PRESSO, weighted median, weighted mode, and MR-Lasso.

With respect to the IVW analysis, a fixed-effects (FE) model is indicated in the case of homogeneous data, whilst a multiplicative random effects (MRE) model is more suitable for heterogeneous data. Burgess et al. recommended that an MRE model be implemented when using GWAS summary data to account for heterogeneity in variant-specific causal estimates⁸⁶. In the interest of transparency, we calculated both results but present the MRE in the text.

MR analyses should include evaluation of exposure IV strength. In order to achieve this, we provided the F-statistic, MR-Egger intercept, MR-PRESSO global test, Cochran’s Q test, and I² for our data. The F-statistic is a measure of instrument strength with >10 indicating a sufficiently strong instrument⁸⁷. We provided F-statistics for individual exposure SNPs and the instrument as a whole. Cochran’s Q test is an indicator of heterogeneity in the exposure dataset and serves as a useful indicator that horizontal pleiotropy is present as well as directing decisions to implement FE or MRE IVW approaches⁸⁸. The MR-Egger intercept test determines whether there is directional horizontal pleiotropy. The MR-PRESSO global test determines if there are statistically significant outliers within the exposure-outcome analysis⁸⁹. I² was calculated as a measure of heterogeneity between variant specific causal estimates, with a low I² indicating that Egger is more likely to be biased towards the null⁹⁰. Finally, we performed a leave-one-out analysis using the method of best fit for each exposure SNP within the IV in order to determine if any single variants were exerting a disproportionate effect upon the results of our analysis⁸⁶.

MAGMA analysis of COVID-19 GWAS data

MAGMA (v1.08)⁹¹ was applied using default settings. Input consisted of summary statistics for all SNPs genome-wide as measured in the COVID-19 GWAS⁷. We estimated LD structure using EUR samples from the 1000 Genomes Project (Phase 3). The top 50 MAGMA genes³⁵ were used for downstream analysis.

DNA sequencing in rare variant analysis

Data quality control

Ancestry analysis

We performed population admixture analysis using ADMIXURE⁹⁸ (v1.3.0) referencing five super populations, including AFR (African), AMR (Ad Mixed American), EAS (East Asian), EUR (European), and SAS (South Asian), in the 1000 Genomes Project (Phase 3) and inferred the ancestry for each genome. For the GEN-COVID cohort, samples with >90% EUR ancestry fraction were kept in the discovery cohort. For the VA COVID-19 cohort, we relaxed the ancestry fraction cutoff to 70% for including more samples in testing. Inferred sample ancestry can be found in Supplementary Tables 8 and 9.

Variant- and gene-level annotations

Genome annotation was performed by Annovar⁵³ integrating multiple databases. Variant frequency was estimated using the 1000 Genomes Project (Phase 3). Nonsynonymous (missense and nonsense) variants were annotated using dbNSFP⁹⁹ (v3.5). The mutation effect of splicing-site variants was predicted by dbscSNV¹⁰⁰ (v1.1) and regSNP-intron¹⁰¹.

Rare-variant burden testing

Rare-variant burden testing was performed to determine whether any genes were differentially enriched with rare variants between severe COVID-19 patients and non-severe COVID-19-positive controls. We utilized whole-exome sequencing data from the GEN-COVID cohort⁵¹, including 122 individuals aged ≤60 years who suffered severe COVID-19 requiring respiratory support, and 465 individuals aged ≥20 years who suffered non-severe COVID-19 not requiring hospitalisation. Variants were included if they altered an amino acid, were rare (MAF<1%) and absent from the EUR cohort of the 1000 Genomes Project (Phase 3). Burden was calculated using SKAT⁵⁰ adjusted for sample imbalance using a saddlepoint approximation¹⁰². Sex and the first ten principal components were included as covariates. Genetic burden was compared with the complete set of coding genes; genes caring <10 variants were removed because of insufficient data. After filtering a total set of 4,280 genes were tested for severe-COVID-19-associated rare genetic variation of which 625 were also RefMap COVID-19 genes. A QQ-plot confirmed that there was no significant genomic inflation (λGC=1.1; Supplementary Fig. 4).

Regeneron’s burden testing results were obtained from the Regeneron results browser (https://rgc-covid19.regeneron.com/results), where only semi-significant genes (P<1e-03, REGENIE¹⁰³) were available. Data consists of exome-wide association studies of various COVID-19 outcomes across 662,403 individuals (11,356 with COVID-19) aggregated from four studies: UK Biobank (UKB; n=455,838), AncestryDNA COVID-19 Research Study (n=83,930), Geisinger Health System (GHS; n=113,731), and Penn Medicine BioBank (PMBB; n=8,904). For the Regeneron study of severe COVID-19 versus non-hospitalized COVID-19, we obtained a list of 68 genes harboring disease-associated missense mutations at a significance cutoff of P<1e-03 in EUR samples. Overlap between PULSE genes and Regeneron gene lists was tested by Fisher’s exact test, assuming a background of 19,396 coding genes in the genome which is the total number profiled by Regeneron¹².

The PULSE model

Network analysis

We first downloaded the human PPIs from STRING v11, including 19,567 proteins and 11,759,455 protein interactions. To eliminate the bias caused by hub proteins, we first carried out the random walk with restart algorithm¹⁰⁸ over the PPI network, wherein the restart probability was set to 0.5, resulting in a smoothed network after retaining the top 5% predicted edges. To decompose the network into different subnetworks/modules, we performed the Leiden algorithm⁶², a community detection algorithm that searches for densely connected modules by optimizing the modularity. After the algorithm converged, we obtained 1,681 modules with an average size of 9.98 nodes (SD=53.35; Supplementary Table 13).

Transcriptome analysis

Four single-cell RNA-seq datasets were used in the transcriptome analyses, including human healthy lungs³³,39 and COVID-19 patients^22,23. Data after QC was acquired for each study. Only samples from the respiratory system were considered in the analyses. For the healthy lung data, a cutoff of 0.6931 was used to define expressed genes in the transcriptome throughout the study if not specified. For the disease samples, we removed the overlap of severe patients between the two cohorts^22,23. In the comparative expression analysis of severe versus moderate patients, to stabilize the analysis we estimated the change of gene expression levels using the Z-score estimated from Wilcoxon rank-sum test, wherein a positive Z-score indicates a higher expression level in severe patients and a negative value suggests the lower expression. The Benjamini-Hochberg (BH) procedure was used for multiple testing correction throughout the study.