Multiomic analysis reveals cell-type-specific molecular determinants of COVID-19 severity

Abstract

The determinants of severe COVID-19 in healthy adults are poorly understood, which limits the opportunity for early intervention. We present a multiomic analysis using machine learning to characterize the genomic basis of COVID-19 severity. We use single-cell multiome profiling of human lungs to link genetic signals to cell-type-specific functions. We discover >1,000 risk genes across 19 cell types, which account for 77% of the SNP-based heritability for severe disease. Genetic risk is particularly focused within natural killer (NK) cells and T cells, placing the dysfunction of these cells upstream of severe disease. Mendelian randomization and single-cell profiling of human NK cells support the role of NK cells and further localize genetic risk to CD56^bright NK cells, which are key cytokine producers during the innate immune response. Rare variant analysis confirms the enrichment of severe-disease-associated genetic variation within NK-cell risk genes. Our study provides insights into the pathogenesis of severe COVID-19 with potential therapeutic targets.

Keywords: COVID-19; GWAS; Mendelian randomization; NK cell; gene discovery; genome-wide association study; machine learning; network analysis; rare variant analysis; single-cell multiomic profiling.

Figure 1

Schematic of the study design (A-H) The COVID-19 GWAS and human lung single-cell multiome (A) are integrated by the RefMap model shown in (B), where gray nodes represent observations, green nodes are local hidden variables, and pink nodes indicate global hidden variables (STAR Methods). Cell-type-specific risk genes are mapped using single-cell multiome profiling (C). Heritability analysis (D), Mendelian randomization (E), transcriptome analysis (F), and network analysis (G) together characterize the functional importance of RefMap genes, particularly for NK cells, in severe COVID-19. Rare variant analysis (H) orthogonally supports the role of NK cells in severe disease. cCRE, candidate cis-regulatory element. See also Table S1.

Figure 2

RefMap identifies cell-type-specific risk genes associated with severe COVID-19 (A) Total number and number of unique genomic regions containing genetic variation associated with severe COVID-19 for 19 different cell types. (B) Total number and number of unique genes implicated by genetic variation associated with severe COVID-19 for 19 different cell types. (C) Fraction of unique genomic regions and genes associated with severe COVID-19 for major cell types; statistical comparison of enrichment was determined by Fisher’s exact test. (D) Similarity between different cell types quantified by the overlap of RefMap genes. Gene set overlap was calculated by the Jaccard similarity index. See also Table S1.

Figure 3

Severe-COVID-19-associated common variants are linked to multiple cell types (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. The Benjamini-Hochberg (BH) procedure was used to calculate FDRs throughout the study. ^∗: FDR < 0.1. (B) Heritability enrichment for different subsets of T cells and NK cells. (C–E) Significant Mendelian randomization results for three exposures linked to severe COVID-19, including blood counts of CD335+ CD314− (C), CCR7− CD314− (D), and CD314+ (E), NK cells; points indicate effect size (β) and standard errors for each SNP-outcome relationship. (F) Sensitivity analyses and robust tests for MR analyses (STAR Methods). (G) Comparative gene expression analysis of RefMap NK-cell 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 unique molecular identifier [UMI]) 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. Distributions were compared by one-tailed Wilcoxon rank-sum test. See also Figures S1 and S2 and Table S3.

Figure 4

Transcriptomic signature of RefMap COVID-19 genes in different cell types (A) UMAP of iterative latent semantic indexing (LSI) for combined gene expression and open chromatin over the gene body and promoter. Cells are colored by relative expression of NKG2D, CD56, CD16, and NKG2A, respectively. Expression is quantified as log₂(normalized gene counts + 1); yellow/orange cells have relatively high expression of each marker. (B) Enriched TF motifs of RefMap COVID-19 regions across 19 cell types. Relative enrichment per cell type is indicated by circle size and significant enrichment (HOMER, Q value < 0.1) is annotated with black circle; TF expression is indicated by color according to log(normalized gene counts + 1). Only highly expressed (expression level in the top 95% in corresponding cell types) TFs were considered. (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. Violin plots show the distributions of log expression levels within each group, and point plots indicate the median and IQR. (D and E) 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 Ren et al. (D) and Liao et al. (E), respectively. The Z score of Wilcoxon rank-sum test was used to indicate the change of gene expression between severe and moderate patient groups, where a positive value means higher gene expression in severe patients. 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. ^∗: 0.01 ≤ FDR < 0.1. +: FDR < 0.01. Distributions were compared by one-tailed Wilcoxon rank-sum test. See also Figure S3.

Figure 5

PPI modules enriched with COVID-19 genes and their functional characterization (A–C) Three PPI network modules, including M546 (A), M1164 (B), and M1540 (C), are significantly enriched with ciliated-cell gene, epithelial-cell genes, and immune-cell genes, respectively. Blue nodes represent RefMap COVID-19 genes and yellow nodes indicate other genes within each module. Edge thickness is proportional to STRING confidence score (>400). (D–F) Gene functions that are significantly enriched (Fisher’s exact test, adjusted p < 0.1) in modules M546 (D), M1164 (E), and M1540 (F). GOBP, GO biological process. (G) Gene expression analysis of module genes in NK and T cells. The transcriptome was defined as the total set of expressed genes in NK and T cells, respectively. (H and I) Comparative gene expression analysis of module genes in severe COVID-19 patients versus moderate patients based on scRNA-seq datasets from Ren et al. (H) and Liao et al. (I), respectively. The Z score of Wilcoxon rank-sum test was used to indicate the change of gene expression 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. Distributions were compared by one-tailed Wilcoxon rank-sum test. See also Table S6.

Figure 6

Rare variant analysis supports the association of NK cells with severe COVID-19 (A) Enrichment analysis of cell-type-specific RefMap COVID-19 genes based on rare variant associations. ^∗: Q value < 0.1, permutation test. (B) Overlapping genes between the RefMap NK-cell gene set and the rare-variant-associated gene set (p < 0.05, REGENIE). MAC, minor allele count. (C) Q-Q plot of rare variant association test for RefMap NK-cell genes and all RefMap genes, including expected p values based on the null hypothesis and observed p values by REGENIE. See also Figure S4 and Table S7.