Integrative analyses identify susceptibility genes underlying COVID-19 hospitalization

TWAS identifies genes for COVID-19 related hospitalization

To identify genes underlying COVID-19 related hospitalization, we tested the predicted expression of 22,207 genes across 49 tissues for association with COVID-19 related hospitalization (see Methods). We identified 123 associations representing 21 genes across 45 tissues at 8 independent genomic regions (p-value < 2.3E-6; Fig. 1, 2; Tables S1-S3). Next, to improve statistical power, we tested for association between predicted gene expression levels from multiple tissues simultaneously with COVID-19 related hospitalization GWAS. Of the 22,207 tested genes, we identified 14 genes across 10 genomic regions, which consisted of 2 additional genes – XCR1 and DNAH3 (p-value < 1.4E-06; Fig. S1.; Table S4.). Overall, we found 23 TWAS-based gene associations across 10 genomic regions.

• Download figure
• Open in new tab

Figure 1: Significant genes.

The integrative analyses identified 27 genes (labeled) across 13 regions color coded) shown in the ideogram.

• Download figure
• Open in new tab

Figure 2: TWAS:

The top panel is a Manhattan plot of genes associated via multiple-tissue WAS. Each data point represents a gene grouped by chromosome (x-axis) and lowest p-value (y-axis) of the gene across significant tissues. The significant genes are shown as pink triangles, wherein triangles facing up and down represent positive and negative z-scores, respectively. The bottom panel show distribution of z-scores across significant gene-tissue airs. The genes are grouped based on chromosome (y-axis) and respective tissues (x-axis).

To find additional support for genetic regulation of identified susceptibility genes by SNPs at risk regions, we tested the 123 gene/tissue pairs identified in single-tissue scans for allelic imbalance within 1Mb GWAS regions (see Methods). We identified 9 genes (ABO, CCR2, CXCR6, FYCO1, IFNAR2, IL10RB, LZTFL1, OAS1, OAS3) with evidence of allelic imbalance at COVID-19 GWAS risk variants, with 3 genes (ABO, OAS3, IL10RB; see Fig. S2., Table S5.) when restricted to leading GWAS index variants (p-value < 0.05 / 21). Together, these results further support a model where risk is conferred through genetic dysregulation at susceptibility genes.

Next, we focused on the impact of genetic regulation of alternative splicing for COVID-19 severity and performed a multi-tissue splicing transcriptome-wide association study (spTWAS; see Methods). Overall, we tested 131,376 splice sites of predicted alternative-splicing expression across 49 tissues for association with COVID19 related hospitalization and identified 420 associations representing 43 splice variants for 11 genes across 49 tissues and 5 genomic regions (see Fig. 1, 3.; Tables S6.-S8.). Next, we performed a multi-tissue analysis (see Methods) and identified 34 splice variants for 12 genes (2 genes – IFIT3 and GNL3 not identified in single-tissue scans) across 40 tissues (p-value < 3.7E-07; see Fig. S3.; Table S9.).

• Download figure
• Open in new tab

Figure 3: Splicing TWAS.

The top panel is a Manhattan plot of genes associated via multiple tissue spTWAS. Each data point represents splice site grouped by chromosome (x-axis) and lowest p-value (y-axis) of the splice site across significant tissues. The annotated genes to splice site are labeled. The significant splice sites are shown as pink diamonds. The bottom panel show distribution of splice sites across significant site-tissue pairs. The genes annotated to splice sites are grouped based on chromosome (y-axis) and respective tissues (x-axis).

Comparing genes identified from TWAS (23 genes) and spTWAS (13 genes), 9 genes were implicated by both approaches - LZTFL1, DPP9, IL10RB, IFNAR2, OAS3, FYCO1, ABO, OAS1 and XCR1. Alternative splicing had stronger overall association signals at the 9 genes in common (p-value = 2.17E-09), with 5/9 genes showing greater signals on average (p-value < 0.05 / 9; see Table S10.).

Next, we interrogated the role of genetic regulation of protein abundances and performed a proteome-wide association study (PWAS) using 1,031 predictive models of plasma proteins fitted from population data in the INTERVAL study (N=3,301; see Methods) (10). Of the 1,031 tests performed, 2 genes (ABO, and OAS1) were significantly associated with COVID-19 related hospitalization (p-value < 4.85E-5; Fig. 1, Table S11.).

Phenome- and laboratory-wide association scans highlight functional role for the 27 genes

We investigated the potential functional role of the 27 TWAS genes using data of 1,404 clinical phenotypes for N=70,439 individuals of European ancestry using the Vanderbilt Biobank, BioVU (see Methods; Fig. 4., Table S12.). Overall, 40 clinical phenotypes were significantly associated with genetically-predicted ABO, IFNAR2, and CCR1 expression levels; ABO accounted for the majority (30 out of 40) of the associations. Across the 17 phenotype categories, we found circulatory system-related phenotypes were enriched for association (see Methods; 7.23-fold enrichment, p-value = 8.62E-22, Table S13). We repeated an enrichment analysis using all association data, regardless of statistical significance, and observed circulatory- and infectious disease-related phenotypes strongly enriched for association signal on average (p-value < 3.9E-44; Table S14). Top associations with genetically-derived ABO gene expression were driven by circulatory system phenotypes, including acute pulmonary heart disease, deep vein thrombosis, other venous embolism and thrombosis, pulmonary heart disease, and acute pulmonary thrombosis and infarction (OR=1.47; p-value = 3.97E-11). IFNAR2 was associated with migraine (OR = 1.35; p-value = 4.10E-06) and with throat pain (OR = 2.05; p-value = 2.62E-05; Table S12).

• Download figure
• Open in new tab

Figure 4: PheWAS Manhattan Plot.

Each data point represents phenotypic associations with genetically-regulated expression of gene-tissue pairs. The data points are grouped and color-coded by phenotype groups x-axis) and -log10(p-value) (y-axis). The dashed line represents Bonferroni threshold, and most significant gene-phenotype associations across all significant tissues are text-labeled.

Next, we focused on laboratory results for N=70,337 individuals of European ancestry using the Vanderbilt Biobank, BioVU (see Methods). For the 323 laboratory traits tested, we found 32 labs significantly associated with four genes (ABO, IFNAR2, KEAP1, and SLC6A20; p-value < 1.55E-04; see Fig. 5., Table S15.). Of these, ABO captured 27/32 significant associations (mean OR = 1.33; p-value = 8.02E-14 < 1.40E-04). Across the 12 broad lab definitions in our data, our significant LabWAS findings were enriched for blood-related lab measurements (see Methods; 4.21-fold enrichment, p-value = 1.23E-11; Table S16.). When extending our enrichment analyses to all associations, we found blood-(see Methods; p-value = 9.23E-22; Table S17.) and immune-related labs (p-value = 2.81E-14) displayed the greatest enrichment, with toxicology-, urinary-, and cancer-related labs exhibiting significant depletion of signal (Table S17.). Genetically-predicted ABO expression was associated with various measures of blood and platelet count, coagulation factors, and ferritin in blood, as well as labs measuring immune and metabolic function (Table S15.). Genetically-predicted IFNAR2 expression was negatively associated with creatine kinase (OR=0.89; p=5.85E-05; Table S15.). Genetically-predicted KEAP1 expression was positively associated with total cholesterol, and non-high-density lipoprotein levels (beta= 0.50, p=3.72E-05). SLC6A20 genetically-predicted expression was negatively associated with basophil volume in blood (beta = -0.19, p= 4.04e-5) and magnesium volume in serum/plasma (beta= -0.26, p= 1.04e-4).

• Download figure
• Open in new tab

Figure 5: LabWAS Manhattan Plot.

Each data point represents laboratory-trait associations with genetically-regulated expression of gene-tissue pairs. The data points are grouped and color-coded by clinical laboratory-test groups (x-axis) and -log10(p-value) (y-axis). The dashed line represents Bonferroni threshold, and most significant gene-laboratory trait associations across all significant tissues are text-labeled.

Cross-ancestry phenotypic comparisons

To perform cross-ancestry validation of the clinical and laboratory phenotypes implicated in the European-based results, we performed PheWAS and LabWAS in the N=15,123 individuals of African ancestry in the BioVU records. Of the 32 identified laboratory measures, we found 22 ABO-associated labs replicated at nominal levels (p-value < 0.05) with 5 replicating after adjusting for the number of tests performed (p-value < 0.05 / 32). We attribute lack of statistical power to the 17 phenotypes that did not replicate after multiple testing correction. Effect sizes across ancestries were highly concordant (slope=0.87, 95CI [0.73, 1.01], p-value = 1.66-13; see Methods; Table S15), with no individual gene/lab pair demonstrating evidence of significant heterogeneity (p-value < 0.05 / 32). Of the 40 clinical phenotypes identified in participants of European ancestry, we found none that replicated at Bonferroni adjusted levels. This was largely due to reduced statistical power from smaller sample sizes, as estimated effect-sizes were similar across ancestries (slope=0.33, 95CI[0.25, 0.41], p-value = 6.72E-10; see Methods; Table S12.). These cross-ancestry analyses suggest similar effects of predicted expression on relevant clinical and laboratory phenotypes.

We next sought to test cross-ancestry replication of identified phenotypes in the Pan-UK Biobank and Biobank Japan. Briefly, we performed PheWAS using LD-independent eQTL/sQTL SNPs (see methods) of TWAS-identified genes. We identified 233 FDR significant SNP-phenotype results dominated by the associations between ABO (214/233, 91.8%) SNPs and blood differential tests such as basophils and monocytes. A subset of 80/233 (34.3%) FDR significant associations also were nominally significant in at least one population of non-European ancestry (p-value < 0.05; Tables S18.-S19.). There were six instances of FDR significant effect estimate heterogeneity across ancestries, all of which involved ABO SNPs and the biomarker alkaline phosphatase or erythrocyte properties. We next tested how significant EUR effect estimates reflect SNP effects across ancestries. We found that EUR SNP effects significantly predicted SNP effects in six global ancestry groups (maximum prediction in AMR; Table S20.).

COVID19-HGI genome-wide association summary statistics

We downloaded GWAS summary statistics for severe COVID-19 outcomes meta-analyzed across 21 studies (hospitalized N=7,885; population N=961,804). A detailed description of the contributing studies, meta analyses, and primary GWAS results for several COVID-19-related phenotypes are presented at https://www.covid19hg.org/results/, specifically the Freeze 4-October 2020 results File: COVID19_HGI_B2_ALL_leave_23andme_20201020.txt.gz. These summary statistics do not include 23&Me cohort results, and their sample size was removed from final sample reported. Genome-wide association statistics consisted of inverse-variance meta-analyzed log-odds ratios and their standard errors to compute a final Wald statistic and p-value. Most of the individuals of the contributing studies to meta-analysis genomic study were of European descent (93%). We performed strict quality control on GWAS data, by filtering statistics at palindromic variants, harmonizing variants with GTEx v8 European-panel genotypes (29). Our quality control procedure resulted in a final count of 10,340,768 autosomal genetic variants with summary statistics.

TWAS and spTWAS using models of predicted gene and alternative splicing expression

To perform TWAS and spTWAS we leveraged pre-trained prediction models fitted in GTEx v8 data (49 tissues, N=838) using the fine-mapping software DAP-G with a biologically-informed prior, Multivariate Adaptive Shrinkage in R (MASHR). For detailed information regarding molecular, genetic, and phenotypic data in GTEx v8, please see (29). Prediction models for each tissue were integrated with COVID19-HGI GWAS data using the software S-PrediXcan (30). In total, we tested 655,563 and 1,728,429 models of total expression and alternative splicing, respectively for association with severe COVID19 outcomes, however, owing to the significant amount of correlation across tissues, we used a per-tissue Bonferroni correction threshold in our multi-tissue analyses (see Table S1). To combine association statistics across all tissues while adjusting for tissue–tissue correlation, we used S-MultiXcan (31). For 22,206 genes, we also performed a joint multi-tissue approach using S-MultiXcan which accounts for tissue correlation and boosts statistical power. Here, we applied a single Bonferroni correction of 0.05 / 22,206.

PWAS using models of predicted expression

To perform a protein-wide association study (PWAS) using predicted protein expression, we fitted predictive models using genetic and plasma proteins from European-ancestry individuals in the INTERVAL study (N=3,301) (10). We performed quality control on genotype data and kept only bi-allelic SNPs with MAF ≧0.01, HWE p > 5e-5, imputation quality INFO > 0.6 and were annotated in HapMap3. Plasma proteins had undergone strict quality control and adjustment in the original study (10). We fit predictive 3,222 predictive models for 3,170 proteins using genotypes within 1Mb flanking the gene body (i.e., ±500kb gene start and stop). For measured proteins consisting of multiple monomers (i.e., dimer, trimer, etc.), we fit multiple predictors, each restricted contributing gene’s region. We included the following covariates into all downstream models of protein abundance: age, sex, duration of blood processing, the top 3 genotyping PCs, contributing cohort, and top 4 protein PCs). To reduce the number of tests and increase statistical power, we restricted to genes whose protein levels exhibited evidence of genetic control by testing for non-zero cis-heritability (p-value < 0.05) using GCTA. Our final set of local/cis-based predictors resulted in 1,031 models of protein with significant cis-SNP heritability (p-value < 0.05). We fit penalized linear models using SuSiE (32) and performed downstream PWAS using the tool FUSION (33) with the same quality controlled genome-wide association statistics as above.

Allelic Specific Expression

To determine the allelic effect of GWAS SNPs (p-value < 5e-5) on identified susceptibility genes, we used haplotype-level ASE data with WASP filtering (34) from the GTEx v8, containing 15,253 samples spanning over 49 human tissues and 838 individuals (29, 35). We used haplotype-aggregated allelic expression generated by phASER (36). To assess cis-acting regulatory effect of expression imbalance between the alleles in heterozygous individuals, we compared allelic imbalance between the individuals homozygous and heterozygous for each SNP. All individuals with minimum coverage of 8 reads (with one pseudocount added) were included. Allelic imbalance was quantified as the log ratio between the two allelic counts or log allelic Fold Change (log aFC) (37) and to ensure robustness to rare variant effects and phasing errors the absolute value of log aFCs are compared, using a one-sided ranksum test. We used a gene-level Bonferonni correction for the total number of genes tested (p-value < 0.05 / 21).

EHR-based PheWAS

To better understand the phenotypic consequences of dysregulated mRNA expression across our genes of interest, we performed a Phenome-Wide Association Study (PheWAS) (38) including patients in the Vanderbilt EHR and linked biobank, BioVU. Phenotypes in BioVU are represented as phecodes, which are assigned as a dichotomous trait and are a hierarchical clustering of the International Classification of Diseases (ICD9/ICD10) codes. For each phenotype we required a minimum number of 100 cases for inclusion in our PheWAS analyses, which resulted in testing 1,404 phecodes in 70,439 individuals of European ancestry and 740 phecodes in 15,174 individuals of African ancestry. We used the PheWAS package in R to perform logistic regressions to identify the phecodes that are significantly associated with imputed gene expression after adjusting for sex, age, and the top ten principal components from genetic data to control for population stratification (Denny et al. 2010, 2013). We corrected for the number of tests (i.e. 0.05/1404 = 3.56e-05) to determine statistical significance.

EHR Biomarker LabWAS

The Lab-Wide Association Scan (LabWAS) (39) allows us to screen clinical lab tests from the Vanderbilt University Medical Center EHR. For each gene identified in the TWAS analyses, we tested the association between its predicted gene expression and all clinical labs. We applied the QualityLab cleaning pipeline (39) with settings to yield median age-adjusted (residual taken after regressing the cubic splines of age with 4 knots) inverse normal quantile transformed lab values (to control for skewness and non-normality). We screened across all labs with measurements for at least 100 individuals, which resulted in testing 323 labs in 70,337 individuals of European ancestry and 241 labs in 15,123 individuals of African ancestry. The lab tests are divided into 12 sub-categories; blood, metabolic, endocrine, kidney, immune, liver, urinary, OB/gyn, toxicology, cardiovascular and cancer. Our analyses included the covariates age, sex, and top ten principal components from genetic data to adjust for genetic ancestry. We used a Bonferroni-corrected threshold accounting for the number of labs present in the associations tested (i.e. 0.05/323 = 1.55e-04).

PheWAS and LabWAS category enrichment analyses

We tested for enrichment of association signals across the clinical categories of phenotypes (laboratories) in PheWAS (LabWAS) in two ways. First, we performed a hypergeometric test using phenotypes/labs that were labelled as significant/not-significant per category. Second, we performed a relaxed test that considers for the average magnitude of association signal in a given category. Here, we computed the mean χ² association statistic per phenotype (laboratory) category and bootstrapped its standard error using 2000 bootstraps. Our enrichment (depletion) statistics for a phenotype (laboratory) category were the difference between its mean χ² from 1 (the expected χ² under the null) divided by the bootstrapped standard error. We used a Bonferroni adjusted p-value < 0.05 / k to determine enrichment or depletion for either aproach, where k=17 for PheWAS category enrichment tests and k=12 for LabWAS category enrichment tests.

Cross-ancestry gene effects

We compared effect-sizes of predicted expression on phenotypes and laboratories estimated in European-ancestry patients from the BioVU with estimates obtained from of African-ancestry records. To do so, we performed a weighted linear regression, , where ., t are the effect-sizes estimated for the ith gene/phenotype or gene/laboratory pair in African-(A) or European-ancestry (E) patients, α is the cross-ancestry relationship, and is the gene/phenotype-specific (gene/laboratory) noise parameterized by the overall variance σ² and squared standard-error around the African-ancestry-based estimate. We report estimates of and its 95 confidence intervals assuming normality.

Cross-ancestry SNP effects

We investigated BioVU phenotypes and laboratory measures significantly associated with TWAS-associated loci for heterogeneous effects in the trans-ethnic Pan-UK Biobank (Pan-UKB). The Pan-UKB represents a multi-ancestry analysis of 7,221 phenotypes in six continental ancestry groups: African (AFR N=6,636), Admixed American (AMR N=980), Centra/South Asian (CSA N=8,876), East Asian (EAS N=2,709), European (EUR N=420,531), and Middle Eastern (MID N=1,599). The Pan-UKB consists of 16,119 genome-wide association studies of biological assays, health status, behavioral information, and lifestyle factors. We clumped SNPs in PLINK using eQTL/sQTL p-values (q-value ≦ 0.05) reported in GTEx v8 and limited pair-wise SNP correlations to r²=0.1 over 250kb windows. Per-tissue clumping resulted in 27 SNPs (ABO, IFNAR2, CCR1, and SLC6A20) that were tested with respect to 1,571 Pan-UKB phenotypes from the same trait domains detected by PheWAS in BioVU (e.g., circulatory system, digestive disorder, neurological). We used linear models to test for consistency of EUR effect estimates at FDR significant and high confidence SNPs with those estimated in AFR, AMR, CSA, EAS (Pan-UKB and Japan Biobank), and MID populations. Additionally, we used Biobank Japan to verify associations between EUR and EAS-based SNP effects in the Pan-UKB due to increased sample size in the latter. Biobank Japan consists of genetic data for over 200,000 participants and ∼120 disease states and quantitative measures (cell type percentage, body mass index, etc.).