An atlas connecting shared genetic architecture of human diseases and molecular phenotypes provides insight into COVID-19 susceptibility

iCPAGdb: An atlas for discovery of cross-phenotype associations

We created iCPAGdb to facilitate exploration of cross-phenotype associations of human phenotypes and discovery of shared genetics connecting traits that were previously not known to be related. iCPAGdb utilizes 85639 SNP-trait associations (p < 5 x 10⁻⁸) across 3793 traits from the NHGRI-EBI GWAS catalog, incorporates additional GWAS datasets (see below and Table 1), and allows for uploading and analysis of user GWAS summary statistics (Fig. 1A). In contrast, the original CPAG (published in 2015) used only 14198 SNP-trait associations for 887 traits from the NHGRI-EBI GWAS catalog.

• Download figure
• Open in new tab

Figure 1. An improved method for finding shared genetic architecture of human traits.

(A) The overall framework of the iCPAGdb pipeline. GWAS summary statistics (from published GWAS datasets or from user-uploaded GWAS) undergo LD clumping to obtain a lead variant for each signal below a specified p-value threshold. These SNPs are queried against an LD proxy database generated from 1000 Genomes African, Asian, or European population to identify cross-phenotype associations through direct overlap or LD proxy at R² > 0.4. Significance of overlap for each trait pair is calculated using Fisher’s exact test. Outputs can be visualized/downloaded from the iCPAGdb web browser.

(B) Comparison of the number of shared SNPs for each NHGRI-EBI GWAS catalog trait pair identified through direct overlap vs. both direct and indirect (LD-proxy) overlap.

(C) iCPAGdb detected more significant cross-phenotypes associations than CPAG1 at FDR < 0.1. Expansion of the NHGRI-EBI GWAS catalog and improvements in capturing by LD proxy in iCPAGdb fueled a large increase in detected cross-phenotype associations across human traits. Comparisons between CPAG1 and iCPAGdb on the same 2013 dataset are in Fig. S3.

(D) Circle plot of cross-phenotype associations detected by iCPAGdb in the NHGRI-EBI GWAS catalog. After excluding compound phenotypes (phenotypes described by NHGRI-EBI GWAS catalog as > 1 comma-separated phenotype in their ontology), a total of 1709 traits involved in a total of 53314 cross-phenotype associations were left. These were categorized into 17 EFO Parental groups. Inner ribbons link phenotypes connected by cross-phenotype associations with the width of ribbon corresponding to the number of cross-phenotype associations. The axis outside the circle represents the cumulative number of associations for each group vs all other groups.

(E) Comparison of genetic correlation from LD score regression (LDSC) and the Chao-Sorensen similarity index implemented in iCPAG demonstrates significant correlation. The genetic correlation rg of 24 diseases/trait were obtained from (Bulik-Sullivan et al., 2015a). Since Chao-Sorensen values are bounded from 0 to 1 and rg ranges from -1 to 1, we used the absolute value of rg here. Colored * indicates significant trait-pair for LDSC, iCPAGdb, or both at false discovery rate of 0.1.

(F) A model demonstrating how SNPs regulate uric acid levels to impact the development of kidney stones and gout.

(G) Riverplot of gout cross-phenotype associations generated from iCPAGdb output shows causal connections, comorbid outcomes, and regulators of disease. Mapped genes for SNPs associated with gout are shown on the left and connected to other NHGRI-EBI GWAS phenotypes grouped by EFO on the right.

Beyond this large expansion in traits and associations, we improved on the original CPAG algorithm by clumping GWAS data from each study (Fig. S1), creating a database of LD values based on 1000 Genomes (Genomes Project et al., 2015), allowing selection of either European, African, or Asian LD structure, and efficiently capturing cross-phenotype associations that are driven by LD proxy (Fig. 1B). For each trait pair, iCPAGdb first selects the lead SNPs from all associated loci at a selected p-value threshold (p < 5 x 10⁻⁸ was used for analysis of the NHGRI-EBI GWAS catalog; Table S1; Fig. S2). These lead SNPs are compared across the trait pair to count directly shared SNPs. For SNPs that are not directly shared, iCPAGdb then checks an LD database for overlap by LD proxy. For all directly or indirectly shared SNPs, iCPAGdb further forms them into bigger SNP blocks by recursively merging them until each SNP block has no LD proxy with R² >= 0.4 against all others. iCPAGdb improves memory efficiency with built-in functions connecting to SQL GWAS and LD proxy databases and improves computational efficiency and speed by utilizing multiple CPUs. For the NHGRI-EBI GWAS Catalog, the growth of GWAS findings and improvements of iCPAGdb over the previous version of CPAG led to a 27.7-fold increase in direct cross-phenotype associations and a 47.7-fold increase in indirect cross-phenotype associations, many of which would have been missed by the original CPAG algorithm (Fig. 1C, D). Indeed, analyzing the 2013 NHGRI-EBI GWAS catalog with iCPAGdb had little effect on direct associations but increased indirect associations by 76% (Fig. S3).

Results of iCPAGdb are consistent with results from the orthogonal approach of genetic correlation by LD score regression (Bulik-Sullivan et al., 2015b). Comparing the absolute values for genetic correlation of 24 phenotypes from (Bulik-Sullivan et al., 2015a) against a similarity index quantifying the degree of shared SNPs in iCPAGdb revealed that the two are significantly correlated (p=3.52 × 10⁻⁸; R² = 0.14) (Fig. 1E). Nearly all phenotypes (64 of 70) that showed significant correlation by LD score regression also demonstrated a significant excess of shared SNPs in iCPAGdb. The output of iCPAGdb provides the SNPs driving the similarity between the two phenotypes, facilitating follow-up studies. Interestingly, 61% of pairwise comparisons that had significant overlap based on iCPAGdb did not have significant genetic correlation based on LD-score regression. For example, LD-score regression did not detect significant genetic correlation between LDL and HDL cholesterol measurements, but iCPAGdb detected 92 shared SNPs, including 31 by direct overlap where the two phenotypes have the same lead SNPs (p=7.55×10⁻¹⁹⁵by Fisher’s exact test; p=1.49×10⁻¹⁹⁰ after Benjamini-Hochberg procedure. P-values from iCPAGdb in the remainder of the paper are FDR-corrected for all pairwise comparisons using Benjamini-Hochberg procedure).

GWAS of varying phenotypic scales reveals shared genetic architecture connecting molecular and cellular traits with human disease

In a previous study (Wang et al., 2015), we defined 4 categories of cross-phenotype associations: 1) SNP similarity between an intermediate trait/risk factor and disease, 2) SNP similarity between a disease and a consequence of disease, 3) SNP similarity between two traits affected by the same gene/pathway, and 4) SNP similarity between two traits affected by the same gene having effects in different tissues or on different pathways. Of these categories, perhaps the most clinically useful is the first category—shared SNPs that connect an intermediate trait to a disease may reveal how molecular or cellular phenotypes mediate some aspect of the pathophysiology of disease. While the NHGRI-EBI GWAS catalog is comprised primarily of case-control GWAS of disease, we detected numerous known shared associations linking a human disease with levels of a metabolite. Metabolites are the substrates, intermediates, and products of cellular metabolism and are routinely already used as biomarkers, such as measuring glucose in diabetes management.

Cross-phenotype associations involving the metabolite uric acid and gout, an inflammatory arthritis driven by excess levels of uric acid (Bodofsky et al., 2020), are illustrative of iCPAGdb’s usefulness. GWAS studies have been conducted on risk of gout (Chen et al., 2018; Lai et al., 2012; Lee et al., 2019; Li et al., 2015; Matsuo et al., 2016; Nakayama et al., 2017; Nakayama et al., 2020; Sulem et al., 2011) as well as uric acid or urate levels (Boocock et al., 2020; Dehghan et al., 2008; Doring et al., 2008; Kamatani et al., 2010; Kottgen et al., 2013; Li et al., 2007; Tin et al., 2019; Tin et al., 2011). Notably, of 31 GWAS loci for gout and 123 GWAS loci for serum uric acid levels at p<5×10⁻⁸, 13 loci overlap, including 9 loci identified only by LD proxy (nearly 6000-fold enrichment; p=5.9×10⁻⁴³). These loci are spread across 7 chromosomes and include several solute carrier (SLC) and ATP-binding class (ABC) transporters that control urate absorption and secretion. Some of the loci are in close proximity but are counted separately by iCPAGdb, as could occur if different GWAS studies locate nearby peaks that fall below our R²>0.4 threshold or if multiple causal signals are located in the same region. These data provide genetic evidence for the well-known causal role of excess uric acid in gout and further reveal multiple genes that may serve as therapeutic targets. Inhibitors of renal uric acid reabsorption through URAT1 (SLC22A12) are commonly used in treating gout (Dong et al., 2019), but additional transporters implicated through human genetics may also prove to be useful drug targets. Beyond uric acid levels, GWAS of kidney stones (Howles et al., 2019; Oddsson et al., 2015; Thorleifsson et al., 2009), a second manifestation of elevated uric acid levels, also share associated SNPs with gout (3 shared loci, all identified by proxy on chromosomes 2, 4, and 17; p=5.2×10⁻⁹). Finally, gout shares 2 loci (out of 5 from (Setoh et al., 2015; Suhre et al., 2017)) with levels of serum alpha-1-antitrypsin, an anti-inflammatory endogenous protease inhibitor (p=9.3×10⁻⁷), providing a human genetic rationale for the use of alpha-1-antitrypsin-based therapeutics in acute gouty flares (as has been demonstrated to be efficacious in mice (Joosten et al., 2016)). Thus, examining the gout cross-phenotype associations revealed causal connections, comorbid conditions with shared etiology, and factors that modulate inflammation in the disease (Fig. 1F, G).

Shared genetic associations reveal other well-known molecular and cellular disease relationships such as LDL cholesterol levels with cardiovascular disease (1.24×10⁻⁸¹) and Alzheimer’s disease (p=4.8×10⁻¹⁷) as well as glucose with type II diabetes mellitus (p=1.5×10⁻⁴⁰). Other cross-phenotype associations highlight genetic variation that can extend our knowledge. For example, cross-phenotype associations were found between malaria (Band et al., 2013; Jallow et al., 2009; Malaria Genomic Epidemiology, 2019; Malaria Genomic Epidemiology et al., 2015; Ravenhall et al., 2018; Timmann et al., 2012) and red blood cell distribution width (Astle et al., 2016; Chen et al., 2020; Chen et al., 2013; Fatumo et al., 2019; Kichaev et al., 2019) (p=1.3×10⁻⁹). This overlap is driven by well-known genetic variation in the beta-hemoglobin gene (HBB) and ABO blood type affecting malaria risk but also by genetic variation in ATP2B4 which encodes a calcium transporter. To the best of our knowledge, whether size of red blood cells impacts susceptibility to malaria parasites has not been examined. These cross-phenotype associations demonstrate the promise of this approach for revealing novel relationships that can be mined through iCPAGdb.

Expansion of iCPAGdb to additional datasets of molecular and cellular traits

The above examples of clinically relevant cross-phenotype associations involving metabolite and cellular phenotypes motivated expansion of iCPAGdb to additional datasets. We used three datasets to provide molecular and cellular traits to our analysis: 491 metabolites and xenobiotics in blood (Shin et al., 2014) and 55 metabolites in urine (Raffler et al., 2015), both from the Metabolomics GWAS Server (http://metabolomics.helmholtz-muenchen.de/gwas/index.php), and 79 cellular host-pathogen interaction traits from our dataset of cellular host-pathogen interaction GWAS, H2P2 (Wang et al., 2018). iCPAGdb revealed many connections between these molecular/cellular datasets and the NHGRI-EBI GWAS catalog (Fig. 2A; Table S2).

• Download figure
• Open in new tab

Figure 2. iCPAGdb integrates GWAS of different scales to reveal biological insight.

(A) Multi-dataset network of cross-phenotype associations detected by iCPAGdb. Phenotypes that demonstrated significant overlap (FDR ≤ 0.1) are color-coded in the indicated colors.

(B) Riverplot of macular telangiectasia type 2 (MacTel type 2) cross-phenotype associations generated from iCPAGdb shows causal connections, comorbid outcomes, and regulators of disease.

(C) Cross-phenotype associations connecting MacTel type 2 and serine. One locus demonstrated direct SNP overlap (rs715). A second locus demonstrated indirect overlap based on 4 SNPs in LD as visualized in the heatmap color-coded by LD.

(D) A model for how SNPs regulate serine levels to impact pathogenesis of MacTel type 2 based on iCPAGdb and prior work described in the text.

(E) Regional Miami colocalization plot demonstrates a genetic locus that impacts both CXCL10 level in lymphoblastoid cell lines following Chlamydia trachomatis infection and CXCL9 (MIG) levels in whole blood.

(F) Comparison of -log10(p value) for GWAS of CXCL10 following Chlamydia trachomatis infection and levels of CXCL9 (MIG) in whole blood. The lead SNP in the region for each phenotype is marked.

(G) Scatter plot demonstrates a highly positive correlation of the effect coefficients of cellular CXCL10 after Chlamydia trachomatis infection and of SNPs associated with blood CXCL9 levels. Each dot represents a SNP which has p value < 0.01 for both phenotypes. A total of 413 SNPs from a 4-mb window surrounding the leading SNP rs2869462 was selected. The blue vertical or red horizontal bar shows the standard error of the beta value for each SNP.

Cross-phenotype associations with macular telangiectasia (MacTel) type 2, a disease characterized by loss of central vision due to alterations in blood vessels in the macula of the eye, confirmed the importance of the amino acid serine (Fig. 2B). A GWAS of MacTel type 2 uncovered 3 genome-wide significant loci and the authors noted that two of these loci were involved in serine/glycine metabolism, with the alleles associated with low glycine and serine conferring increased risk of MacTel type 2 (Scerri et al., 2017). The authors speculated that the low serine levels could lead to high levels of ammonia and glutamate causing neurotoxicity and stress-induced angiogenesis (Scerri et al., 2017). Gantner et al. have since provided evidence that low serine levels result in elevated levels of deoxysphingolipids to trigger cell death in photoreceptors (Gantner et al., 2019). iCPAGdb rediscovered the connection of two loci being shared between serine in serum (measured by (Shin et al., 2014)) and risk of MacTel (Fig. 2C, D; p=4.0×10⁻⁷; 99,010-fold enrichment). iCPAGdb also revealed 7 other serum metabolites including glycine that shared an association with rs715 but not with the second MacTel locus. While serine was not part of the urine metabolomics dataset, iCPAGdb did detect overlap of glycine in urine and MacTel type 2 (p=0.01).

We also included host-pathogen traits from H2P2, a cellular GWAS we previously carried out using 528 lymphoblastoid cell lines (LCLs) exposed to 7 different pathogens (Wang et al., 2018). Notably, unlike the metabolomics datasets, H2P2 identified SNPs associated with traits at baseline and in response to stimuli. Further, as pathogens have likely been drivers of human evolution (Fumagalli et al., 2011; Pittman et al., 2016), comparing H2P2 to human disease GWAS may reveal unintended consequences of past pandemics on the human genome. Previously, we reported colocalization of a locus regulating CXCL10 levels following Chlamydia trachomatis infection (rs2869462) and risk of inflammatory bowel disease (Wang et al., 2018). iCPAGdb revealed shared genetic variants for this H2P2 phenotype and blood levels of CXCL9 (MIG) (Ahola-Olli et al., 2017) (Fig. 2E; p=0.04). P-values for the two associations are strongly correlated (Fig. 2F), and the effect size for SNPs associated with both chemokines are significantly positive correlated (Fig. 2G). We utilized COLOC, which uses a Bayesian framework to determine whether GWAS signals in the same region are likely due to the same causal variant (Giambartolomei et al., 2014). The posterior probability that both CXCL10 protein levels from cells and CXCL9 levels in blood share the same causal variant is 0.90 (Table 2), with rs2869462 identified as the most likely causal SNP (Table S3). The genes encoding these two chemokines are adjacent to each other on chromosome 4, and this result points to variants regulating expression of both genes that will make it challenging to disentangle their effects in disease.

Application of iCPAGdb to COVID-19 reveals susceptibility due to ABO may occur through regulation of CD209

We applied iCPAGdb to a recently published GWAS of severe COVID-19 with respiratory failure (Ellinghaus et al., 2020). While this study focused on two genome-wide significant associations at the ABO locus and in a cluster of chemokine receptors and other genes on Chromosome 3, we relaxed the p-value threshold for iCPAGdb to 1 × 10⁻⁵, resulting in 24 suggestive loci after LD clumping. Not surprisingly, iCPAGdb revealed that the genome-wide significant association near the blood type locus ABO is in LD with multiple other SNPs in this region associated with other human diseases and traits (Fig. 3A; Table S4). This included the classic association with malaria resistance (Timmann et al., 2012), but also less well known associations with duodenal ulcer (Tanikawa et al., 2012), pancreatic cancer (Amundadottir et al., 2009), and heart failure (Shah et al., 2020). Multiple studies have now reported the association of the ABO locus with risk of COVID-19 (Ellinghaus et al., 2020; Zhao et al., 2020). The causal effect on COVID-19 may involve A and B antigens on blood cells, antibodies against A and B antigens, the enzymatic activity of the ABO glycosyltransferase on possibly other glycoproteins, or even other genes in the region. Insight into these possible mechanisms was revealed by iCPAGdb, which identified association of this locus with levels of 8 individual proteins in the NHGRI-EBI GWAS catalog. These proteins, all encoded on different chromosomes than ABO, include IL-6, TNF-α, CD209 (DC-SIGN), Tie-1, mannose-binding protein C, FGF23, and clotting factors (factor VIII and vWF). In each of these cases, the association of the locus to both molecular trait and disease provides a plausible causal chain from SNP to cis-effect on ABO to trans effect on a protein to severe COVID-19 disease. For example, association with VWF and Factor VIII may indicate ABO affects COVID-19 through regulation of thrombosis, as patients with severe COVID-19 can have thromboembolic complications as part of a hyper-inflammatory state (Wool and Miller, 2020). In fact, both VWF and factor VIII are targets of glycosylation by ABO (Canis et al., 2018; Matsui et al., 1992; Sodetz et al., 1979) and levels of these proteins are reported to be regulated by ABO (Albanez et al., 2016; Gallinaro et al., 2008; Murray et al., 2020; Shima et al., 1995; Song et al., 2015). Further, regulation of levels of IL-6 and TNF-α suggest possible regulation of inflammation, as “cytokine storm” plays an important role during severe COVID-19 (Mangalmurti and Hunter, 2020). Most interestingly, the ABO locus is associated with both COVID-19 and CD209 (p=0.008). A preprint recently confirmed this association across populations, and these authors speculated that ABO may affect CD209 levels to regulate SARS-CoV-2 entry (Katz et al., 2020). Indeed, there has since been evidence from two preprints that CD209 can bind to SARS-CoV-2 and can act as a receptor for entry into immune cells (Amraie et al., 2020; Chen et al., 2021).

• Download figure
• Open in new tab

Figure 3. Cross-phenotype association of ABO reveals a possible role for CD209 in severe COVID-19.

(A) A network of genetic associations involving severe COVID-19. Each node represents either a disease/trait (filled circles) or a gene (dark blue diamond). The ABO locus was associated with multiple other diseases and levels of specific proteins, while DPP9 connects COVID-19 only with IPF and interstitial lung disease (idiopathic interstitial pneumonia).

(B) Regional Miami colocalization plot demonstrates the ABO locus impacts both CD209 protein levels and risk of severe COVID-19.

(C) A significant positive correlation for effect size of SNPs in the ABO locus on CD209 protein levels and risk of severe COVID-19.

(D) Model of how ABO may affect CD209 and severe COVID-19.

The “A” allele of rs657152 associated with increased risk of COVID-19 with respiratory failure is also associated with increased levels of CD209 (Fig. 3B). We performed colocalization analysis of the GWAS signals for COVID-19 (Ellinghaus et al., 2020) and CD209 protein levels (Suhre et al., 2017). This analysis indicated the two are likely driven by the same causal variants (Fig. 3C; COLOC posterior probability PP4 = 0.98 with the lead causal SNP of rs505922; Table S5). Thus, iCPAGdb and subsequent colocalization analysis support a model where ABO regulates CD209 protein levels to impact COVID-19 risk, though much future experimental and clinical studies will be required to fully test this hypothesis (Fig. 3D). The pleiotropic effects of ABO on levels of multiple proteins will make defining the mechanism challenging.

Application of iCPAGdb to COVID-19 reveals a role for DPP9 in regulation of both COVID-19 and idiopathic pulmonary fibrosis

Beyond ABO, a locus in the dipeptidyl peptidase 9 (DPP9) gene associated at p<1×10⁻⁵ with severe COVID-19 was identified as being shared with a GWAS of fibrotic idiopathic interstitial pneumonia (Fingerlin et al., 2013) and a recent GWAS of the most severe form of that group of diseases, idiopathic pulmonary fibrosis (IPF) (Allen et al., 2020). rs12610495 was the lead variant for each of these GWAS studies as well as the suggestive peak for severe COVID-19 (p=5.2×10⁻⁶; (Ellinghaus et al., 2020)). Much evidence has already accumulated that pulmonary fibrosis is a hallmark of severe COVID-19 (Ojo et al., 2020; Shi et al., 2020). While the association of rs12610495 with COVID-19 did not reach genome-wide significance in Ellinghaus et al. 2020, this SNP is in LD with the lead variant from a recent GWAS of critically ill COVID-19 patients that does surpass genome-wide significance (p=3.98×10⁻¹²; (Pairo-Castineira et al., 2020); R²=0.95 in 1000 Genomes European populations). Thus, iCPAGdb alerted us to the importance of a suggestive COVID-19 susceptibility locus that has since been validated in an independent cohort.

We determined that rs12610495 is an eQTL in lung tissue for the gene for DPP9 (and no other genes in the region) in GTEx (p=4.5×10⁻⁹; (Bao et al., 2015)), with the “G” allele being associated with lower expression (Fig. 4A). Interestingly, DPP9 is a protease in the same family as DPP4, the receptor for MERS-coronavirus (Raj et al., 2013). Additionally, DPP9 is an inhibitor of inflammasome activation by NLRP3 (Okondo et al., 2017; Okondo et al., 2018; Zhong et al., 2018). Colocalization analysis confirmed the signals from severe COVID-19 and IPF are likely driven by the same causal variant (Fig. 4B; COLOC posterior probability PP4 = 0.994, lead SNP rs12610495; Table S6). Based on these data and the known biology, we developed alternative hypotheses for how this SNP might be regulating risk of severe COVID-19: DPP9 may be acting as a previously unrecognized receptor for SARS-CoV-2 or it may be inhibiting inflammation during COVID-19 infection. Based on the directionality of effect of rs12610495 on DPP9 gene expression, the “G” allele should lead to lower DPP9 expression and less entry if the receptor model is correct. However, the “G” allele is instead associated with increased risk of severe COVID-19 (Fig. 4C). Alternatively, the “G” allele could lead to lower DPP9 to increase inflammasome activation in lung tissue, a model consistent with “G” increasing risk of severe COVID-19 and this allele also increasing risk of idiopathic pulmonary fibrosis (Fig. 4D).

• Download figure
• Open in new tab

Figure 4. Cross-phenotype analysis and COVID-19 patient transcriptomics reveals a role for DPP9 in severe COVID-19.

(A) Lung eQTL data from GTEx shows rs12610495 “G” allele is associated with reduced expression of DPP9.

(B) Regional Miami colocalization plot demonstrates the DPP9 locus impacts both idiopathic pulmonary fibrosis and risk of severe COVID-19.

(C) A significant positive correlation for effect size of SNPs in the DPP9 locus on idiopathic pulmonary fibrosis and risk of severe COVID-19.

(D) Model of how DPP9 may affect idiopathic pulmonary fibrosis and risk of severe COVID-19.

(E) DPP9 expression in peripheral blood is significantly higher in COVID-19 patients compared to healthy and bacteria-infected patients. The p values were calculated using the Wilcoxon rank-sum test.

(F) COVID-19 patients demonstrate significantly higher DPP9 expression compared to healthy controls during early (days 1-10), middle (days 11-20) and late (21+ days) stages of SARS-CoV-2 infection. The p values were calculated using the Wilcoxon rank-sum test.

(G) DPP9 demonstrates increased expression during recovery from COVID-19. A total of 11 patients were measured sequentially at enrollment (day 0), day 7, and day 14. The colored dash line connects measurements from the same patient across time points. P value was calculated using Friedman test.

(H) Decreased symptom severity scores of COVID-19 patients over time. The eleven subjects in G were assessed for symptom severity at day 0, 7 and 14. The colored dash line connects measurements from the same patient across time points. P value was calculated using Friedman test.

To further examine the role of DPP9 in COVID-19, we analyzed transcriptomics of peripheral blood from COVID-19 patients (McClain et al., 2020). Levels of DPP9 expression across 46 COVID-19 patients were compared to individuals with seasonal coronavirus, influenza, bacterial pneumonia, and healthy controls. DPP9 levels were significantly increased in COVID-19 patients compared to the other groups (fold-change = 1.15, p = 0.003 adjusted by Benjamini-Hochberg method). Comparing COVID-19 data vs. each comparator individually revealed that DPP9 levels were elevated vs. healthy controls (p = 0.0016) and bacterial infection (p = 0.0078) but not influenza or other coronavirus infection (Fig. 4E). This data supports a role for DPP9 in the host response to viral infections. In examining all samples in the cohort, increased DPP9 was observed both early and late in COVID-19 infection (Fig. 4F). However, eleven subjects that did not require hospitalization had repeated measurements at day 0 (initial enrollment into the study), day 7, and day 14 that revealed changes in DPP9 expression as infection resolved. While DPP9 expression increased from day 0 compared to 7 days and 14 days (Fig. 4G; p = 0.0089), symptom severity dramatically improved over this period (Fig. 4H; p = 0.00006). We speculate that DPP9 may be induced to effectively turn off the inflammatory response to SARS-CoV-2 to minimize tissue damage and fibrosis. Combined with our human genetic data, these findings suggest that insufficient induction of DPP9 expression could predispose to severe COVID-19.

Searching the iCPAGdb web server with user-provided GWAS summary statistics

As the above examples demonstrate, iCPAGdb analysis can rapidly generate hypotheses connecting molecular and cellular traits to human disease. The website allows quick access to the pre-calculated cross-phenotype associations results described in this manuscript. Users can also upload their own GWAS summary statistics for comparing against all 4414 GWAS traits in the iCPAGdb website, facilitating the discovery of new cross-phenotype relationships. Total time for uploading, clumping of summary statistics, and calculation of cross-phenotype associations is typically <2 minutes.

Collection of GWAS summary statistics

Publicly available GWAS summary statistics were downloaded from the following sources: 3793 traits from NHGRI-EBI GWAS Catalog (version 1.02, downloaded on 2020/08/05), 79 traits from H2P2 cellular GWAS (Wang et al., 2018), 587 traits from human blood circulating metabolites and urine metabolites GWAS (Raffler et al., 2015; Shin et al., 2014). NHGRI-EBI GWAS catalog traits included annotation by Experimental Factor Ontology (EFO). All GWAS data were harmonized to genome coordinates of HG19. In total, we collected 4,225 GWAS traits, and 104,247 Trait-SNPs pairs at a p value threshold of 5 × 10⁻⁸. A detailed list of trait-SNP pairs at varying p-value threshold can be found in Table 1.

Severe COVID-19 GWAS summary statistics were downloaded from the GRASP website (https://grasp.nhlbi.nih.gov/Covid19GWASResults.aspx) (download date 2020/07/15). Genome coordinates were converted from GRCh38 to HG19 using UCSC liftOver. GWAS summary statistics of IPF were kindly provided by Allen et al. 2020 after requesting access https://github.com/genomicsITER/PFgenetics.

LD clumping

GWAS summary statistics were individually pre-processed by LD clumping using PLINK v1.9 (Chang et al., 2015), based on genotypes from European populations from the 1000 Genome project. The general PLINK command was “--clump-p1 1e-5 --clump-p2 1 --clump-r2 0.4 --clump-kb 1000”. For NHGRI/EBI GWAS catalog, the index SNPs were selected using the genome-wide significant p value threshold of 5 × 10⁻⁸ (--clump-p1 5e-8). For molecular and cellular GWAS, we used a varying p-value cutoff from 1 × 10⁻³ to 1 × 10⁻⁵ for --clump-p1 parameter to choose the index SNPs.

For uploaded GWAS data, iCPAGdb calls on PLINK automatically to perform LD clumping. Users can define the p value for --clump-p1 to select the index SNPs and choose proper LD structure (European, African, or Asian) based on the ancestry of the GWAS.

LD proxy calculation

To maximize phenotypic associations due to indirect associations, pairwise LD R² values were computed for each leading SNP against its surrounding SNPs using the genotypes from the 1000 Genome project (Phase 3 genotypes). Prior to calculation, all SNPs with minor allele frequency less than 0.01 and missingness > 0.1 were removed. R² of pairwise SNPs within 10,000 bp windows were then calculated, and only LD proxies with R² > 0.4 were retained in further analysis. The PLINK parameters for calculating LD was “--ld-window-kb 1000 –ld-window 10000 –keep-allele-order –r2 in-phase with-freqs gz”.

Since GWAS may be performed on diverse populations from different ancestry or continents, we calculated ancestry-specific LD proxies for European, African, and Asian populations separately. European population included 503 samples from 5 populations (CEU, TSI, FIN, GBR, IBS), African included 661 samples from 7 populations (YRI, LWK, GWD, MSL, ESN, ASW, ACB), and Asian population included 504 samples from 5 populations (CHB, JPT, CHS, CDX, and KHV). We filtered genotypes for each ancestry population by minor allele frequency more than 0.01 and retained only biallelic SNPs. SNPs which have same genome coordinates were merged using “—merge-equal-pos”. For duplicated SNPs with same variant rsID, we kept only the first variant by using “--rm-dup force-first” using PLINK 2.0,

Cross-phenotype SNP analysis

Cross-phenotype SNPs were used to quantify the similarity of different traits. Cross-phenotype loci were identified as leading SNPs and/or their LD proxies having statistically significant associations with more than one trait/disease. If two traits shared a common leading SNP, we termed this “direct association”. If a leading SNP was associated with one trait, while its LD proxy SNPs were associated with another trait, we called this “indirect association”. If any shared SNP was in LD with another SNP with R² > 0.4, these SNPs were merged into a SNP block until no further LD was found across shared SNP/LD pairs.

The significant association among each trait pair were using the hypergeometric distribution.

Where N_e is the effective number of independent SNPs in the selected population, the n₁ and n₂ are the number of independent SNPs associated with trait 1 and trait2, and k is the number of independent SNPs blocks. The effective number of independent SNPs for European, African and Asian population were obtained from Table 4 from (Li et al., 2012).

The significance of associations for all trait pairs was further corrected for all possible pairwise comparisons using the Benjamini-Hochberg and Bonferroni methods for multiple test correction. A false discovery rate of 0.1 was chosen to identify significantly correlated trait pairs.

Comparison to LDSC

Bulik-Sullivan et al. (Bulik-Sullivan et al., 2015b) developed an innovative and unbiased method, LDSC, to estimate genetic correlation using GWAS summary statistic for all measured SNPs. Their model calculated the LD scores for a variant against all other variants in a 1 centimorgan window and hypothesized that SNPs with higher LD scores are tagged to a risk-conferring variant, and the genetic correlation among traits can be calculated by normalizing genetic covariance of SNP heritability. With this method, they estimated 276 genetic correlations for 24 diseases/traits based on full GWAS summary statistic (Bulik-Sullivan et al., 2015a). To evaluate the power of iCPAGdb, we calculated the genetic associations on the same 24 GWAS traits. For each trait pair, only SNPs associated with each trait passing genome wide significant threshold (5 x 10⁻⁸) were used by iCPAGdb. We quantified the strength of cross-phenotype similarity for each trait pair using Chao-Sorensen similarity index. Since the p values from (Bulik-Sullivan et al., 2015a) were not corrected by multiple test correction, we calculated the p values for r_g using R “p.adjust” function with a total number of 276 comparisons.

Colocalization analysis

To evaluate whether the associations of GWAS trait pairs identified by iCPAGdb were due to sharing the same causal variants, we performed colocalization analysis using R COLOC packages (Giambartolomei et al., 2014). COLOC uses a Bayesian framework to estimate the posterior probability that two GWAS traits share two independent causal signals (PP3) or shares a single casual variant (PP4) in the selected genome region. For each trait pair evaluated by COLOC, SNPs within 200 kb window from the lead SNP were included. Since COLOC requires minor allele frequency (MAF) for each SNP in both GWAS studies, when MAF was not available, we calculated the MAF using European populations from the 1000 Genome Project. We ran COLOC “coloc.abf” function using the default prior parameters, p1 = 1×10⁻⁴, p2 = 1×10⁻⁴, and p12 = 1×10⁻⁵. We also ran built-in “sensitivity” function to evaluate the robustness of predefined priors, and all tests suggested that default prior parameters are robust, therefore, we ran all colocalization analyses with default priors values.

COVID-19 transcriptomic analysis

As described in (McClain et al., 2020), samples were collected as part of the Molecular and Epidemiological Study of Suspected Infection (MESSI) which was conducted at Duke University Health System (DUHS) and the Durham Veterans Affairs Health Care System (DVAHCS). The study was approved by each institution’s IRB. Informed consent was obtained from all subjects or their legally authorized representatives, and informed consent were collected for all subjects. SARS-CoV-2 RT-PCR testing was used to confirm infection status. A total of 46 subjects were analyzed, 14 of which were assayed at more than 1 timepoint. In total, 77 samples were assayed. Subjects were divided into early (≤10 days), middle (11-21 days), and late (>21 days) stage based on duration of symptoms. Participant self-reported symptoms were recorded at each timepoint for 39 symptom categories. Each symptom was scored on a scale of 0–4, with 0 indicating not present, 1 mild, 2 moderate, 3 severe, and 4 very severe symptoms. Daily symptom severity (sum of symptom scores for all symptoms) was determined for each timepoint. At enrollment (Day 0), date of symptom onset was determined, and an initial symptom survey recorded maximum score for each symptom category between symptom onset and study enrollment. Total RNA was extracted from peripheral whole blood, and cDNA libraries prepared using NuGEN Universal Plus mRNA-seq with AnyDeplete Globin reduction were sequenced on the Illumina NovaSeq 6000, as described (McClain et al., 2020). In brief, STAR v 2.7.1 (Dobin et al., 2013) was used to align the short reads and generate the count matrix. The count matrix was further normalized using TMM method (Robinson and Oshlack, 2010) and log2 transformed. Associations were performed with generalized linear models (LIMMA, (Ritchie et al., 2015)) and corrected for multiple testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). Analysis of DPP9 was carried out in R, and p values were calculated using the Wilcoxon rank-sum test.

iCPAGdb software and website implementation

iCPAGdb is comprised of two core parts, the back-end computation and the front-end web browser. The back-end was written in python v3.6 with utilization of SQLite. SQLite tables were constructed for harmonized GWAS datasets and LD tables for different populations and are accessed using python sqlite3 package. The GWAS table stores clumped GWAS summary statistic, including trait name, trait sources, SNPs’ rsID, beta values, standard error/standard deviation of beta, effective allele, and p values. The ancestry-specific LD proxy tables contain pairwise SNPs’ rsID and R² values (R² >= 0.4) for different populations. All SQLite tables were indexed on unique combinations of SNP and trait or SNP pairs for LD proxy tables, which greatly reduces the searching time. To further increase calculation speed, the core cross-phenotype analysis part of iCPAGdb is parallelized by utilizing multiple threads.

Primary software components for the web portion of iCPAGdb are the R statistical programming language (Team, 2020), the R package Shiny (v1.5.0) for interaction of web pages with R scripts (Cheng et al., 2020), Shiny Server as a 24/7 multi-user platform to make Shiny apps publicly accessible (RStudio, 2020), the database environment SQLite for efficient cient querying of GWAS and CPAG results (Hipp, 2020), and the R package RSQLite to execute SQL queries from within R scripts (Muller et al., 2020). The results of a CPAG execution are read by the R script, processed, and presented to the viewer in various tables and graphs on a web page. The iCPAGdb website is currently loaded with associations across more than 4400 public GWAS datasets that can be browsed and searched in “Review” mode. The user requests an existing CPAG result set from which a corresponding table and heatmap are generated and displayed. Various filtering and graph construction controls are available for iterative sub-setting of data and selection of significance measure and number of top signicant phenotype pairs to plot. The “Download” button enables the researcher to make a local copy of records appearing in the currently displayed results table. Important packages used in this mode are DT for construction of and interaction with tables and ggplot2, plotly, and heatmaply for basic plotting, interactive plotting (hover labels), and heatmap generation, respectively. The web browser also allows users to upload their own GWAS summary data, and iCPAGdb will automatically perform LD clumping based on selected population and generate an atlas of connections for the user’s GWAS against > 4400 GWAS traits in the database. In this “Upload” mode, the user browses filles on a local computer, selects a properly formatted GWAS result file of interest (containing, for a single phenotype, SNP rsIDs and GWAS p-values), specifies format and column configuration, then uploads the file. Next, CPAG computation parameter values, including iCPAGdb GWAS set to be crossed with, significance thresholds for filtering, and linkage disequilibrium (LD) population are specified. When “Compute CPAG” is pressed, the R script composes a system level command to execute the CPAG (Python) function. The future() function of the R future package (Bengtsson, 2020) combined with a delaying pipe from the promises package execute CPAG operations asynchronously, waiting on completion before resuming R script execution. Typical run time for a single uploaded GWAS that is already clumped to lead variants is <30 seconds. For GWAS summary statistics including all SNPs in a study, run time is typically < 2 minutes. The results are available with downloadable tables and figures. Additional information on webapp is in Supplemental Note.

Web resources:

iCPAGdb: http://cpag.oit.duke.edu

NHGRI GWAS Catalog: https://www.ebi.ac.uk/gwas/

H2P2 cellular GWAS: http://h2p2.oit.duke.edu

Human metabolite GWAS summary statistics: http://metabolomics.helmholtz-muenchen.de/gwas/index.php?task=download

COVID-19 GWAS summary statistics from Ellinghaus et al. (2020): https://grasp.nhlbi.nih.gov/Covid19GWASResults.aspx

IPF GWAS: download link was obtained by applying for access following the collaborative protocol from https://github.com/genomicsITER/PFgenetics

Tools for visualization:

R packages:

ggplot2: https://cran.r-project.org/web/packages/ggplot2/

gggene: https://cran.r-project.org/web/packages/gggenes/index.html

tidygraph: https://cran.r-project.org/web/packages/tidygraph/

ggnetwork: https://cran.r-project.org/web/packages/ggnetwork/

circlize: https://cran.r-project.org/web/packages/circlize/

ggpubr: https://cran.r-project.org/web/packages/ggpubr/

DT: https://cran.r-project.org/web/packages/DT

plotly: https://cran.r-project.org/web/packages/plotly/

heatmaply: https://cran.r-project.org/web/packages/heatmaply/

promises: https://CRAN.R-project.org/package=promises

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Dennis C. Ko (dennis.ko{at}duke.edu). All iCPAGdb output described in this manuscript are available for browsing from http://cpag.oit.duke.edu. Supplemental files also contain iCPAGdb output and COLOC analysis results. Code is available at GitHub https://github.com/tbalmat/iCPAGdb.