Multitrait genetic association test

The first step of our study consisted of ensuring the validity of the proposed statistical tests by studying potential bias, and assessing their statistical power under several simulation scenarios. Let us denote the vectors of single nucleotide polymorphism (SNP) Z-scores z = (z₁,…,z_K,), where K is the number of phenotypes (i.e., the number of GWASs analyzed jointly). The first model we used, which we refer to as sumZ, assumes that genetic effects across the phenotype analyzed follow a direction specified by a vector of weights w to form a weighted sum of Z-scores. Here we considered four weighting schemes: i) uniform weighting (sumZ₁); ii) weighting according to the first principal component of the phenotypic correlation matrix (sumZ_r); iii) weighting according to the first principal component of the genetic correlation matrix (sumZ_g); and iv) weighting according to the independent component analysis of the Z-scores matrix (sumZ_ica). The second approach, which we refer to as omnibus, does not rely on a prior specification on the direction or magnitude of the SNP effect across traits. In brief, it compares, for one SNP, the vector of genetic effects z with the expected multivariate normal distribution under the null. It is a standard omnibus test based on summary statistics that allows for one degree of freedom per outcome (here per phenotype). Both approaches (sumZ and omnibus) rely on a valid estimation of Σ_r, the variance-covariance matrix between z₁,…,z_K, under the null hypothesis of no association, and share similarities with previous approaches (see S1 Text).

We first performed an in-depth validation of each approach, starting with a series of simulations under an ideal situation, when there were no missing data and the true Z-score covariance matrix, Σ_r is known (S1–S3 Figs). We further show using both simulated data (S4 and S5 Figs and S1 Text) and real data from the UK Biobank cohort that in the specific case of complete sample overlap between GWASs, the omnibus test is asymptotically similar to a multivariate analysis of variance (MANOVA) applied to individual level data (S6 Fig and S1 Text). The only major potential source of bias we identified is the misspecification of Σ_r which can lead to severe type I error inflation (S7 Fig). Misspecification can affect all variants, if Σ_r is estimated naively using the z-score data as proposed in previous studies (S8 Fig). Comparing various approaches, we found that Σ_r can be accurately estimated using LDscore regression[9] (S9 Fig), and the approach was therefore used to estimate Σ_r along the genome-wide genetic correlation (Σ_g) for the 36 phenotypes analyzed (S2 and S3 Tables). Nevertheless, misspecification can also be SNP-specific when sample size varies across the SNPs analyzed. Per-SNP sample size heterogeneity can induce different proportions of sample overlap and potentially invalidate Σ_r for those variants. We illustrate this potential bias by applying omnibus tests for the joint analysis of the four GWASs from the GC consortium[33] (S10 Fig). To address this issue, we implemented additional tools to estimate sample size per SNP when missing and subsequently filtered the variants with small sample sizes (S11 andS12 Figs and S1 Text). Finally, out of 10 million variants reported for some GWAS, fewer than 1,000 had complete summary statistics for all 36 phenotypes analyzed. While methods exist to impute missing GWAS statistics, we found them to be inaccurate for multitrait analyses and we implemented the RAISS[27] approach we recently developed to ensure valid imputation for our context (S13 Fig). All preprocessing steps were recently incorporated into a publicly available toolset[28]. After applying our preprocessing pipeline to all 36 GWAS analyzed, there remained 6,978,319 SNPs with a missing data rate of 45% (59% before imputation).

Next, to illustrate the relative detection ability of each approach, we conducted a series of simulation studies across various scenarios. S14 Fig represents the null hypothesis rejection boundaries of each test using a simple example with two traits, when the genetic effect follows a single bivariate normal distribution and when it follows a mixture of two bivariate normal distributions (to reflect our working hypothesis). In this setting, the omnibus test displays the largest detection rate. Among the sumZ tests, sumZ_g and sumZ_ica show better performances than sumZ₁ and sumZ_r, especially when the structure of the genetic signal is heterogeneous. We then conducted more extensive simulations focusing on the omnibus and SumZ_g tests. Fig 1 lays out the results of simulations using 10 traits generated under contrasting scenarios. Again, the omnibus test shows the best performances, especially when all traits have a high heritability and if the genetic signal is not structured along a specific direction. The statistical power of the omnibus test also tends to increase with the sample overlap especially when the environmental correlation was not aligned with the genetic correlation. Interestingly this good performance was also observed for genetically uncorrelated traits when they have a high heritability. In the case of a highly structured genetic effect sampled along a specific direction SumZ_g and the omnibus test performed similarly. However, when only one of the 10 traits had a high heritability, the omnibus test underperformed compared to the SumZ_g, reflecting the cost of additional degrees of freedom in the omnibus test.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Statistical Power of the omnibus and sumZ_g tests under several simulation scenarios.

Power of the omnibus and sumZ_g tests with respect to sample overlap. Color of the line represents the test. Each panel correspond to a simulation scenario. Point shape indicates if the residual covariance was generated to be partially aligned with the genetic covariance or to be unconstrained (random). A 10 trait Z-score vector was generated as the sum of a genetic effect and a residual effect: Z = Z_g + Z_res. Z_res was sampled from a normal multivariate distribution with covariance matrix terms equal to where n₁ is the sample size of the first study, n₂ is the sample size of the second study and ρ is the phenotypic covariance among the n_s overlapping samples. Z_g varied depending on the simulation scenario: (eff ~ Random) Z_g were sampled from a uniform distribution with boundary [-6; 6], (eff ~ corG) Z_g was sampled from a normal multivariate distribution with a random covariance matrix, (eff ~ wG) Zg was sampled along a straight line blurred with a normal noise, (eff ~ high H) Z_g was sampled from a normal multivariate distribution simulating genetically uncorrelated traits with high heritability, (eff ~ het H) Z_g was samples from a normal multivariate distribution simulating genetically uncorrelated trait with only the first having a high heritability.

https://doi.org/10.1371/journal.pgen.1009713.g001

Comparison of multitrait association in real data

We analyzed the 36 GWASs of European ancestry (S1–S3 Tables) using the aforementioned multitrait approaches applied to seven phenotype sets: five medical-based sets (Immunity, Anthropometry, Metabolism, Cardiovascular and Psychiatric), a BMI related set including anthropometry traits and lipids (referred further as the Composite set), and finally all 36 phenotypes jointly (Fig 2). Note that we included bone mineral density traits in the Immunity set because an enrichment of BMD genome wide significant loci in immune pathways and immune cell regulatory regions has been previously reported[34, 35]. We derived the overlap of significant loci of the multitrait tests per phenotype set (S15–S21 Figs), and after merging all analyses (Fig 3A). We applied a Bonferroni correction to the joint tests and used a p-value threshold of 10⁻⁸ instead of the standard 5x10^-8. Univariate phenotype associations were included in the comparison using the minimum of univariate p-value across all outcomes (noted P_univ). Across all phenotype sets, 391 associations were identified by the multitrait tests only, 392 were identified by univariate association tests only, and 1557 were significant for both univariate and multitrait tests (see Fig 3A). The largest number of new associations was detected by the omnibus test. The performances of the sumZ tests varied substantially depending on the phenotype set. For example, the weighting scheme based on phenotypic correlation (SumZ_r), detects slightly more signals than other weights for the Immunity set (S19 Fig) but fewer associations in other phenotype sets (Fig 3A). While the Omnibus detected the largest number of new associations, the substantial share of signals found by other models suggests that applying several multivariate tests, especially the combination omnibus, sumZ_ica, and sumZ_g, could be an interesting solution to maximize detection. Finally, we checked the 392 associations identified by the multitrait test only in these data against previously reported associations from the GWAS catalog[1] for the same phenotypes. Altogether, we report a total of 322 new associations (S4–S10 Tables).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Analysis overview.

The diagram presents the overall analysis pipeline. A total of 36 GWAS were included covering several common diseases and quantitative traits. All GWAS summary statistics went through extensive pre-processing and quality control filtering, and missing single SNP statistics were imputed when possible. Multitrait approaches were then applied to all clean GWAS data and on each clinically based set (All, Immunity, Metabolism, Brain, Cardiovascular, Anthropometry, and Composite). After combining univariate and multivariate results, and merging SNPs within locus, a total of 6,767 associations were identified. After a comparison of results per approach, a clustering analysis was performed for variants within each set. Finally, we performed in-silico functional analysis of the clusters derived in the Metabolism set to assess their biological relevance.

https://doi.org/10.1371/journal.pgen.1009713.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Multitrait approach comparison.

Panel (A) shows independent variants detected across the six approaches: univariate test (univ), omnibus test (omni), weighted sum of Z-score with uniform weight (sumZ₁), weight defined as the loading of the first principal component of the phenotypic correlation (sumZ_r), the genetic correlation (sumZ_g), or defined using the loadings of an independent component analysis (sumZ_ica). Each line corresponds to a test and each column to a set of significant variants. For each set, the test for which variants are significant are represented with a black dot on the test line. The barplot at the left represents the total number of significant independent signals detected by each approach. The stack bar at the top represents the cardinality of the sets. The next panels show the link between strengths of univariate association signal and the relative performance (i.e. larger power) of the four most tests: univ, omni, sumZ_g, and sumZ_ica, for each phenotype set: anthropometry (B), cardiovascular (C), immunity (D), metabolism (E), brain (F), composite (G), and all phenotypes (H). Within each phenotype set, we split the top associated SNPs per region based on the most significant test, and derived the median chi-squared for each test. The radar plots show the derived median per test and illustrate the strong heterogeneity in patterns identified. For example, out of the 1605 SNPs from the anthropometry set, 1235 had stronger signal with univ as compared with other tests. The median chi-squares in that group were 49.1, 1.1, 2.0, 1.0, and 0.7 for height, body mass index (BMI), hip circumference (Hip), waist circumference (WC), and waist to hip ratio (WHR). Comparatively, the 267 SNPs harboring a stronger signal with omnibus, had median of 6.8, 20.1, 15.9, 11.2, and 7.2 for the same phenotypes.

https://doi.org/10.1371/journal.pgen.1009713.g003

To further understand the relative performances of those three tests (omnibus, sumZ_ica, sumZ_g) along the univariate test, we explored which multitrait signal was associated with the largest increase in detection per test. To this end, we listed all loci found associated with at least one of the four approaches, and assigned each locus to a test based on the lowest p-value. We then derived the median squared z-score by phenotype across the loci assigned to each test. As shown in Fig 3B–3H, the median pattern varied substantially across tests and phenotype sets. Higher power for the univariate test was, as expected, observed for strong association signals for a single phenotype, and mostly reflected a very large sample size for that phenotype and/or a strong heritability (e.g. height in the anthropometry set, Fig 3B, or atrial fibrillation in the cardiovascular set, Fig 3C). A strong association signal for the omnibus test was linked to the inclusion of correlated phenotypes and sample overlap, resulting in a high residual covariance (Σ_r, S2 Table). For example, the median squared z-score was elevated for any stroke (AS), any ischemic stroke (AIS) and cardioembolic stroke (CES) in the Cardiovascular set. The patterns preferentially detected by the sumZ_g test are harder to interpret. However, we noticed that sumZ_g displays a strong signal for SNPs associated with physiologically related traits (e.g., T2D and fasting glucose in the metabolism set, Fig 3E, or bone mineral density of neck and spine in the immunity set, Fig 3D).

To confirm the relevance of the associations detected by multivariate tests, we also conducted a tissue enrichment analysis to significant variants identified by the multitrait approaches and by the univariate analyses separately (S11 and S12 Tables). Overall, univariate variants and multitrait variants harbored a very similar functional enrichment landscape (S22 Fig). Most enriched tissues are already known to be involved in the phenotype in question, including liver, fat and pancreas for the Metabolism set, immune cell types and thymus for the Immunity set, and heart for the Cardiovascular set. Our enrichment study also confirmed less obvious observations, which have nevertheless been noted before: the involvement of immunity in brain-related traits (e.g. autisms and schizophrenia)[36, 37] and the overrepresentation of brain tissues in the Metabolism set[38, 39].

Distinct genetic association profiles correspond to distinct genetic correlations

Our comparison of approaches highlights that associated genetic variants display a broad range of multitrait association profiles. We investigated how these profiles can be broken down into groups of homogeneous multivariate genetic effects. This is directly related to the principle of genetic correlation, which quantifies the concordance of genetic effects across traits (e.g. [9]). The difference here, is that genetic correlation captures only the average over the whole genome, and as discussed in recent studies, more localized genetic structures likely exist for many pairs of traits[10]. To detect such a structure, we implemented a multivariate Gaussian mixture model (MGMM)[40] for the identification of clusters among SNPs found associated with at least one approach. We applied MGMM assuming between 2 and 10 clusters and used the BIC and silhouette criteria to determine the most relevant number of clusters. We further bootstrapped the computation of the clustering criteria to ensure the robustness of the selection (S1 Text). The best suited number of clusters is 6, 8, 8, 9, 3, 2 and 5 for the Metabolism, Immunity, Cardiovascular, Anthropometry, Psychiatric, Composite, and All sets, respectively (S23 Fig). As illustrated for the Metabolism set in S24 Fig, adding significant SNPs from the multitrait tests to those identified by the univariate tests enabled us to detect more clusters.

We assessed the uncertainty in cluster assignation by deriving the entropy for each variant (see S1 Text). We observed some heterogeneity in the distribution of entropy values across phenotype sets and clusters (Fig 4C and S16 and S17 Tables). The difficulty to attribute a cluster for certain variants might be due to the lack of representative cluster (i.e. sub-structure not captured in our analysis) or to shared functionalities between clusters that are not modelled in the GMM framework. Because differentiating those two possibilities using the available data would be very challenging, we decided to remove outlier variants (N = 28 across all sets) with entropy above 0.75 from further analyses.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Multitrait genetic association clusters for the Metabolism set.

The panels summarize the clustering of the 392 independent SNPs selected from the Metabolism set analysis. The set includes 10 phenotypes: triglyceride (TG), total cholesterol (TC), type 2 diabetes (T2D), low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C), glycated hemoglobin (HbA1c), Homeostasis model assessment of β-cell function (HOMA-B), homeostasis model assessment of insulin resistance (HOMA-IR), fasting insulin, and fasting glucose. The alluvial plot in panel A) represents the decomposition of univariate genetic association and its rewiring to the six inferred clusters. The flow widths represent the proportion of phenotype’s variance explained by the subset of SNPs assigned to each specific cluster, relative to the total genetic variance explained by all 392 SNPs. For example, SNPs from cluster 6 capture approximatively 41.7% and 54.6% of that genetic variance for TC and LDL, respectively. For clarity, flows explaining less than 0.1% of the variance are not represented. Panel B) shows the heatmap of normalized beta coefficients per phenotype within each cluster. Each column is a SNP, with blue and red colors indicating negative and positive beta, respectively. Coded alleles have been defined according to the per cluster first principal component. The boxplots in panel C) shows the distribution per cluster of SNP’s entropy, an indicator of the fitness of the SNP-cluster assignment. SNPs perfectly assigned are expected to have entropy close to zero.

https://doi.org/10.1371/journal.pgen.1009713.g004

The resulting clustering is presented in Fig 4 for the Metabolism set and in S25-S31 Figs for the other sets. Each figure presents a heatmap of Z-scores along with an alluvial plot displaying both the shared explained variance between phenotypes and the proportion of explained variance by clusters for each phenotype. The complete list of SNPs for the Metabolism set per cluster is presented in S16 Table. The multivariate effects vary substantially from one cluster to another. For instance, in Metabolism clusters, SNPs from the cluster 3 display increased HDL-C and decreased triglycerides, while SNPs from cluster 5 are more specific to triglycerides.

The alluvial figures and heatmaps provide an overview of the magnitude of genetic effects from one cluster to another. To further characterize the concordance or discordance of genetic contributions across phenotypes, we computed the pairwise SNP-based genetic correlations for each cluster (see S1 Text). Fig 5 presents those estimates for a subset of phenotypes within the Metabolism and Immunity phenotype sets. In the Immunity set, the correlations between rheumatoid arthritis (RA), ulcerative colitis (UC) and Crohn’s disease (CD) provide a striking illustration of how the genome-wide genetic correlation can be composed of smaller structures. The genome-wide genetic correlations between UC and CD is strong (0.41), but near 0 and not significant for RA (see S3 Table). In Fig 5B, we noticed a fairly large negative correlation in clusters 2 and 3 between RA both CD and UC whereas, the cluster 5 captures a group of variants displaying a strong positive correlation across the three traits. Similar negligible genome-wide correlations along with opposite genetic correlations across clusters were observed in the Metabolism set. For example, variants from cluster 1 display a strong concordant effect between LDL and T2D, but variants from cluster 6 harbor an equally strong negative correlation. Fig 5 also highlights that significant genome-wide genetic correlations across highly related phenotypes such as UC-CD and LDL-TG are not distributed evenly across variants.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Heterogeneity of genetic correlation across clusters for the Metabolism and Immunity sets.

We derived the genome-wide genetic correlation between phenotypes using LDscore regression and using Pearson correlation from all SNP Z-scores (top panels), and for SNPs within the identified clusters. Results for the Metabolism set are presented in panel (A) using only the four key traits, LDL, HDL, Triglyceride (TG) and type 2 diabetes (T2D). Results for the Immunity set are presented in the panel (B). For clarity only significant correlation are represented. The boldness of the line is proportional to the strength of the genetic correlation. Positive correlations are represented in blue and negative correlations in red. The values of the genetic correlation are indicated by the number next to the trait. Solid lines represent significant correlation (after Bonferroni correction) whereas dashed lines represent correlation significant only before Bonferroni correction. Note that because the clusters are inferred from the multivariate associations, the absolute value of the significance of the correlations is of limited interest. Nevertheless, it provides a useful descriptive statistic to identify the key structures within each cluster.

https://doi.org/10.1371/journal.pgen.1009713.g005

Biological meaning of genetic clusters

We conducted series of in silico functional analyses with the objective of mapping clusters to candidate biological functions. For each phenotype set, we evaluated two types of enrichment: tissue-specific chromatin mark enrichment per cluster (S13 Table) and a pathway enrichment framework (S14 and S15 Tables) which integrates multiple databases such as Gene Ontology (GO) and KEGG. Here, we focused on the Immunity and Metabolism sets as a case study. Given a phenotype set, while using the same set of traits for all clusters, we observed large differences in pathways, tissues and cell type enrichments between clusters.

For the Immunity set, clusters 1 and 4 predominantly capture genetic effects of bone-mineral density; clusters 2, 3 and 5 affect inflammatory bowel disorder (IBD); and clusters 6, 7 and 8 capture variants with pleiotropic effects on rheumatoid arthritis and IBD (S27 Fig). Both enrichment analyses pointed toward an overrepresentation of the immune system with all clusters–even those affecting primarily bone-mineral density–being enriched for at least one immunologic pathway or one immunological tissue. We highlight the top enriched tissues and top pathways in Table 1. Concerning pathway enrichment, immune related pathways regulating the shape of the immune response, such as cytokines and the JAK-STAT signaling pathway were recurrent. Interestingly, variants from those clusters map to a distinct set of cytokines and cluster of differentiation genes (e.g., IL4, IL13, IL33 for cluster 1 and IL3, IL5, IL10, IL19, IL20, IL21, IL27 for cluster 5), which suggests that they may impact different components of the immune system. Concerning tissue-specific active chromatin mark enrichment, clusters 2 and 3 contained multiple SNPs enriched primarily in transcriptionally active regions of “primary natural killer cells from peripheral blood” whereas cluster 7 and 8 are enriched for “primary T helper cells.” We also observed enrichment in the tissue where immune damage occurred for cluster 5 (colonic mucosa), which highlights the complex interaction between the immune system and the inflamed tissue.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Top tissue associations and Immune related Genes by Clusters for the Immunity set.

https://doi.org/10.1371/journal.pgen.1009713.t001

The Metabolism set includes several molecular phenotypes, which we expect to be closer to biological mechanisms than some of the macrophenotypes from other sets. Overall, cluster 1 is mostly associated with an increase in fasting glucose and impaired β-cell function; cluster 2 is highly pleiotropic and notably increased the risk of T2D; and clusters 3 to 6 are mostly associated with lipids (Fig 4). Accounting for the direction of effects, we also noted that the genetic associations in cluster 5 match the known phenotypic correlation with the inverse relationship between circulating levels of HDL-C with those of LDL-C and, more especially, TGs observed in epidemiological studies[41]. At the tissue level, we observed modest enrichment for adipocytes in clusters 1 and 2 (FDR p-value 0.028 and 0.01 respectively, S13 Table) and cluster 3 SNPs are up-regulated in the liver (FDR p-value 0.005).

As shown in S14 Table, each cluster was significantly enriched for a large number of GO terms. We report some specific and illustrative examples: cluster 1 is enriched for the carbohydrate homeostasis set (q-value = 2.5 x 10⁻³), cluster 3 is enriched for the reverse cholesterol transport set (q-value = 2.8 x 10⁻¹³), cluster 4 is enriched for the plasma lipoprotein clearance set (q-value = 1.7 x 10⁻⁵), cluster 5 is enriched for the protein lipid complex assembly set (q-value = 1.08x10^-9) and cluster 6 is enriched for the low density lipoprotein particle remodeling set (q-value = 1.07x10^-2). Cluster 4 also exhibits active chromatin tissue enrichment in immune T cells (q-value = 2.3x10^-3), highlighting the link between cholesterol and immunity. Indeed, cholesterol and modified forms of cholesterol, such as oxidized cholesterol and cholesterol crystals, promote inflammatory and immune responses through multiple pathways including the activation of the Toll-like receptor (TLR) signaling, the NLRP3 inflammasome and myelopoiesis[42, 43]. While the promotion of inflammation and immunity is carried by LDL particles, HDL particles were proposed to counteract this effect in part through reverse cholesterol transport[44]. However, cluster 3, which is enriched for reverse cholesterol transport did not exhibit such tissue enrichment in immune T cells, indicating that the link between HDL and immunity may be more complex, as recently pointed out by Madsen et al[45].

Metabolism pathways and diseases

To provide a perspective on the specificity of genetic variants across clusters and their potential contribution to human diseases, we investigated the lipids from the Metabolism set. We first projected each cluster gene onto KEGG pathways. Here, we used only maps corresponding to enriched GO gene sets identified or to tissue identified in the enrichment analysis at the previous stages (S14 and S15 Tables): fat digestion and absorption, cholesterol metabolism, and PPAR signaling pathways. We constructed a synthesis of these observations on the metabolic map presented in Fig 6A and 6B. Genes associated with clusters (S16 Table) had functions in agreement with their effects on blood lipid levels: cluster 3 is enriched in genes involved in HDL-C biogenesis and metabolism (LCAT, ABCA1, SR-B1, CETP, PLTP, LIPG, APOAx and APOCx), clusters 4 and 6 with genes related to LDL-C clearance (SORT1, PCSK9, LDLR, LDLRAP1, APOB and APOE), and cluster 5 to genes related to triglycerides and chylomicron transport (LPL, APOAx and APOCx).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Mapping clusters to pathways.

We projected cluster’s genes from the Metabolism phenotype set onto KEGG pathways and reconstructed a synthetic metabolic map. Panel A) presents the results for the lipoprotein component and panel B) for the lipid component. Gene names are highlighted by the colors of their associated clusters. When a gene is associated to several SNPs belonging to different clusters it is represented with several colors. To improve interpretation, we also present in panel C) a proxy for the relative contribution of each phenotype per cluster, defined as the loadings of the first principal component derived from the matrix of Z-score for the subset of SNPs in that cluster. Finally, panel D) shows the distribution of standardized beta for association between SNPs from each cluster and three diseases: any stroke (AS), coronary artery disease (CAD), and obesity (using body mass index as a proxy).

https://doi.org/10.1371/journal.pgen.1009713.g006

We then assessed the effect of variants from each cluster with three diseases known to be associated with serum lipids: coronary artery diseases (CAD), stroke, and obesity (defined as a BMI > 30) (S19 Table). Within each cluster, we aligned the SNP alleles with the main trend of the corresponding cluster so that all coded alleles fit the multitrait pattern defined in Fig 6C (see S1 Text). For example, all SNPs from cluster 5 were recoded to be associated with an increase in TGs, TC and LDL-C and a decrease in HDL-C. We plotted in Figs 6D and S31 the genetic effect of each SNP on the three diseases (using the effect on BMI as a proxy for obesity) after the aforementioned alignment, and we performed a sign test to assess the significance of the observed trend (S19 Table). SNPs from several clusters displayed a significant increase in the risk of CAD: cluster 2 (P = 6.6x10^-5), cluster 4 (P = 2.9x10^-2), cluster 5 (P = 3.9x10^-3) and cluster 6 (P = 2.8x10^-4). SNPs from cluster 2 also displayed a nominally significant increase in the risk of stroke (P-value = 1.6x10^-2). Finally, a large fraction of SNPs from cluster 3 had a negative effect on BMI (P = 6.4x10^-4). Interestingly, several SNPs from this cluster show association with CAD, but with heterogeneous effects–some associated with an increased risk and other associated with a decreased risk–so the absence of a global trend. The associations of clusters 4 and 6 with CAD add to the evidence of a causal effect of LDL-C on CAD[46], which has been established by prospective epidemiological studies[47], Mendelian randomization[48] and randomized clinical trials evaluating the effect of LDL-C reducing therapies[49]. Moreover, the association of cluster 5 with CAD risks corroborates a potential causal role of TGs[5] and remnant cholesterol[50, 51] in CAD. The role of TGs in CAD has also been confirmed by epidemiological studies[52], genome-wide association studies[5], Mendelian randomization studies[53] and randomized controlled trials aiming the lower of TGs[54]. Cluster 3, which is associated with increases in HDL-C, does not have a protective effect on CAD, which is again in agreement with Mendelian randomization analyses reporting no link between HDL and CAD[48, 55]. Finally, the association of cluster 2 with CAD and strokes further supports the potential causal effect of type 2 diabetes on CAD and stroke[56].

As a final exploratory analysis, we reported the cluster and multitrait genetic effect of genes targeted to mitigate hyperlipemia to prevent CAD (Table 2). Drug targets corresponding to cluster 3 (ABCA1, CETP, NR1H3) did not lead to successful clinical trials, whereas targets (PCSK9, NPC1L1, APOC3, HMGCR) in clusters 4, 5 and 6 were mostly successful or promising. The example of the CETP gene, that is classified in cluster 3, a cluster not associated with CAD, is of particular interest. CETP has been the target of failed clinical trials that attempted to prevent CAD by inhibiting CETP and consequently increasing circulating HDL-C[57–59]. Cholesteryl ester transfer protein (CETP) promotes the heteroexchange of cholesteryl esters and TGs between HDL-C and APOB-containing lipoproteins connecting HDL-C and TG metabolism[57]. Pharmacological inhibition of CETP was motivated by GWAS[60] and prospective cohorts[61] that indicated that CETP variants were associated with higher circulating HDL-C levels and lower LDL-C, TGs and CVD risk. However, although all CETP inhibitors achieved an effective increase in HDL-C, only anacetrapib led to a significantly lower incidence of major coronary events[62] in patients who were receiving statin therapy, an effect that might account for the reduction in ApoB (non-HDL-C) rather than the elevation of HDL, as suggested by Mendelian randomization analyses[63]. In addition to these well-known drugs, we provide a systematic listing of potential drug targets by cluster (S20 Table) based on the druggable genome database[64].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Drug target genes and associated SNPs in the metabolism set

https://doi.org/10.1371/journal.pgen.1009713.t002

Altogether, these results suggest that drug development might be more effective by accounting for the gene context, i.e., by selecting candidate genes not from their individual features but based on the disease association trend of genes displaying similar multitrait association profiles. Under this working hypothesis, the proposed inference of genetic functional groups can provide a means to identify those genes and therefore to select potential candidates.

Multivariate association test

Consider a vector z of K Z-scores statistics for a single nucleotide polymorphism (SNP) obtained from standard univariate genome-wide association screenings of K phenotypes. Under the null hypothesis, z = (z₁,…,z_K,) follows a normal distribution N(0,Σ_r), where Σ_r is the residual phenotypic covariance matrix (S1 Text), while under the alternative, z is expected to display additional covariance due to shared genetics (defined by a genetic correlation matrix Σ_g). We first considered an Omnibus test of the vector of Z-scores, which can be performed using the multivariate Wald statistics: where T_omni follows a chi-square with K degree of freedom (df) under the null hypothesis of no phenotype-genotype association. We also considered a classic weight-based test defined as: where w is a vector of K weights applied to the Z-score. Under the null, T_sumZ follows a chi-squared distribution with 1 degree of freedom. Note that this approach shares similarities with both standard fixed effect meta-analysis[14] and with dimensionality reduction methods (e.g. principal component analysis[25]). One can also note that the Omnibus statistics can be expressed as a combination of the sumZ statistics over all eigenvectors of Σ_r(S1 Text). We note v_i the i^th eigen vector of Σ_r:

We considered four weighting schemes for the sumZ tests: (i) in the SumZ_1, w is equal to the unit vector so all traits have the same weight; (ii) in the SumZ_r, w is equal to the first eigen vector of Σ_r so its direction represents phenotypic correlation between traits, (iii) in the SumZ_g, w is equal to the first eigen vector to Σ_g so its direction represents genetic correlation between traits, (iii) in the SumZ_ica w is computed by applying an Independent component analysis (ICA) to the complete matrix of Z-score. To compute the weight vector w of the SumZ_ica, for a given phenotype set, the genome wide Z-score matrix was extracted and an independent component analysis was performed with the scikit-learn python package. The component yielding the most novel association was selected as loadings. We verified that this selection procedure did not lead to an inflation under the null hypothesis by simulation (see S2 Fig).

Performing the omnibus test requires inverting the Z-score covariance matrix Σ_r. When this matrix does not have a full rank, we use a pseudo inverse of the matrix based on the singular value decomposition (S1 Text). Briefly, as Σ_r is a variance-covariance matrix, it can be written PDP^t where D = diag((λ_k)_{k = 1…K}), (λ_k)_{k = 1…K}) are the eigenvalues of Σ_r and P is the orthogonal matrix whose columns correspond to the eigenvectors of Σ_r. If it is not invertible, only K’ eigenvalues are different from 0 (where K’ denotes the rank of Σ_r) and an inverse of the matrix can be computed as , where and P_K′denotes the K×K′ matrix whose columns are the K′ eigenvectors corresponding to the eigenvalues different from 0. Note that the Omnibus statistics computed with follows a χ² distribution with K′ degree of freedom.

Robust estimation of Z-score covariance

The validity of the proposed multivariate tests mostly relies on the accurate estimation of Σ_r. In practice, the covariance between Z-scores from null SNPs from two GWAS will deviate from 0 when there is both sample overlap and correlation among the traits analyzed. When combining results from two independent studies, or when the trait analyzed has negligible correlation, Σ_r will be a diagonal matrix, so that the Omnibus test can be performed by summing chi-squared statistics for each SNP to form a K degree of freedom chi-square, and the sumZ test becomes a standard weighted meta-analysis of fixed effect. Yet, in the large-scale GWAS era, this situation is unlikely as most of the large GWAS are conducted in the consortium setting, where samples likely overlap across multiple GWAS. It follows that Σ_r can contain non-zero off-diagonal terms. Under the complete null model, the expected Z-score covariance for null SNPs between two traits equals where n₁ is the sample size of the first study, n₂ is the sample size of the second study and ρ is the phenotypic covariance among the n_s overlapping samples (see S1 Text and e.g. [3, 9]). In some specific cases, one can obtain these parameters directly from the data (e.g. when analyzing multivariate omics data). Conversely, obtaining all four parameters (ρ, n_s, n₁, n₂) from consortium GWAS based on dozen or even hundreds of cohorts can be a practically daunting and risky task. Moreover, accurate phenotypic covariance estimation would be particularly challenging when study-specific and trait–specific covariates adjustment has been performed. Recent studies proposed to estimate Σ_r using available SNPs from the GWAS in question using all available single SNPs Z-score[75] or using a random subset of pruned variants[14], though some discussed removing GWAS hits[15], focusing on a subset of SNPs in regions less likely to contain causal variants[76], or using tetrachoric estimator[16]. The validity of these estimators mostly relies on the assumption that the vast majority of the SNP effects in the genome are distributed under the null hypothesis. While this is likely to be true in some cases, associated variants can potentially lead to either upward or downward pairwise covariance between Z-scores. Instead, we leverage recent work by Bulik-Sullivan et al[3, 9] that allows for estimation of this covariance (and the diagonal variance terms) under a polygenic model and assuming multivariate normality of effect sizes across traits (see S1 Text). The estimation of Σ_r was performed on Z-scores before the imputation step described in the next section. For a few traits the estimated variance is markedly inferior to 1. As indicated in the LDSC regression method, this phenomenon happens when the original GWAS was corrected with a genomic control factor.

Simulation studies on test statistical power

To assess the statistical power of the Omnibus test and the SumZ tests, we designed the following simulation scenarios. A 10 trait z-score vector was generated for each of the 5000 causal SNPs as the sum of a genetic effect and a residual effect: Z = Z_g + Z_res. In the following, random covariance matrices were generated using the randcorr R package[77]. Z_res was sampled from a normal multivariate distribution with covariance matrix (Σ_r) terms equal to where n₁ is the sample size of the first study, n₂ is the sample size of the second study and ρ is the phenotypic correlation among the n_s overlapping samples. As demonstrated in[9], ρ = ρ_g+ρ_e where ρ_g is the genetic correlation and ρ_e the environmental residual correlation. We generated the residual covariance matrix Σ_r according to two scenarios: (random) ρ_r was randomly sampled with the only constraint that Σ_r remains a semi definite positive matrix, (aligned) ρ_r was constructed as with sampled at random. Z_g varied depending on the simulation scenario: (eff ~ Random) Z_g were sampled from a uniform distribution with boundary [-6; 6], (eff ~ corG) Z_g was sampled from a normal multivariate distribution with a random covariance matrix, (eff ~ wG) Zg was sampled along a straight line blurred with a normal noise, (eff ~ high H), Z_g was sampled from a normal multivariate distribution simulating genetically uncorrelated traits with high heritability, (eff ~ het H), Z_g was samples from a normal multivariate distribution simulating genetically uncorrelated trait with only the first having a high heritability.

Data preprocessing: an overview

The analysis of the 36 GWAS required substantial preprocessing, including the inference of several parameters. First, for many publicly available GWAS, sample size per SNP was not readily available and retrospectively collecting this information can be very challenging as it implies requesting this information from each individual cohort. For such a situation, we propose inferring a proxy for missing sample size as , where is the variance of , the estimated SNP effect, and the variance of the SNP, derived from the coded allele frequency which is either provided with the GWAS or extracted from a reference panel (see S1 Text). For linear regression this approximate , where N is the true sample size and is a residual variance of the outcome in the regression model. For logistic regression our estimator is a proxy for the term Np(1−p), where p is the in-sample proportion of cases, and it therefore assumes that the proportion of cases is relatively stable across SNPs with different sample size.

Another challenging issue was the merging of multiple GWAS with different set of assayed SNPs. Indeed, out of 10 million variants reported for some GWAS, fewer than 1,000 had complete summary statistics for all 36 phenotypes analyzed. We performed an imputation of missing Z-scores in each study using the RAISS[27] method we recently developed. The approach uses correlation between SNPs to predict Z-score at missing SNPs using available ones and achieves a level of imputation accuracy suitable for multitrait analysis (S1 Text). Here we used the European panels from the 1,000 Genomes project[78] as a reference for the estimation of the correlation between SNPs. Overall, imputation did not lead to any observable inflation of the omnibus statistic (S13 Fig). Nevertheless, as a supplementary quality control (QC), we excluded significant SNPs that were not surrounded by SNPs in linkage disequilibrium with significant or near significant p-values (P < 10⁻⁶).

These two parameter inferences were integrated along other preprocessing operations into a pipeline that is fully described here[28]. Given a reference panel with no ambiguous strand, it consists in the following steps (i) Extract, the coded and alternative alleles, signed statistics (regression coefficient or odds ratio), standard error, p-value, and sample size; (ii) Remove all SNPs that are not in the reference panel; (iii) Derive Z-score for each SNP from signed statistics and p-value; (iv) Infer sample size when not available; (v) Remove all SNPs whose sample size is less than 70% of the maximum sample size; and (vi) Infer missing Z-scores statistics based on the 1K genome reference panel. After applying our preprocessing pipeline to all 36 GWAS analyzed, there remained 6,978,319 SNPs with a missing rate of 45% (59% before imputation).

Characterization of new loci

To determine new and existing trait-associated loci we used genome regions formed by linkage disequilibrium (LD) blocks as defined in Berisa et al[79] using a reference panel of individuals of European ancestry. It included a total of 1,704 independent regions ranging from 10 kb to 26 Mb in length, with an average size of 1.6 Mb. For each independent LD region, we extracted the minimum p-value over all SNPs contained in the region, and a single univariate analysis p-value defined as the minimum across all single phenotype GWAS and all SNPs in the region. We consider that a region is newly detected by a multitrait test if the joint analysis p-value is genome-wide significant while its univariate p-value is not (joint analysis p-value < 1x10^-8 and univariate p-value > 1x 10⁻⁸). We determined SNPs carrying the signal inside significant region with the plink “clump” function using the following parameters:—clump-p1, 10⁻⁸;—clump-r2, 0.2. We kept the lead SNP by clump for further analysis (gene mapping and clustering).

To report associations exclusively detected in the current report (S4–S10 Tables), we filtered out association present in the GWAS catalogue[1] at the date of the 14^th of September 2020 (univariate p-value > 5x 10⁻⁸) for traits corresponding to our phenotype set. The following trait labels were used to retrieve associations: (Metabolism set) ‘Fasting blood glucose’, ‘Triglycerides’, ‘LDL cholesterol’, ‘LDL cholesterol levels’, ‘HDL cholesterol’, ‘HDL cholesterol levels’, ‘Total cholesterol levels’, ‘HOMA-B’, ‘HOMA-IR’, ‘Hemoglobin A1c levels’, ‘Type 2 diabetes’; (Psychiatric set) ‘Schizophrenia’, ‘Bipolar disorder’, ‘Major depressive disorder’, ‘Alzheimer’s disease’, ‘Educational attainment’; (Anthropometry set) ‘Height’, ‘Waist circumference’, ‘Waist-hip ratio’, ‘Body mass index’, ‘Hip circumference’; (Immunity set) ‘Bone mineral density’, ‘Rheumatoid arthritis’, ‘Ulcerative colitis’, ‘Inflammatory bowel disease’, ‘Crohn’s disease’, ‘Asthma’; (Cardiovascular set) ‘Coronary artery disease’, ‘Ischemic stroke’, ‘Large artery stroke’, ‘Stroke’, ‘Atrial fibrillation’, ‘Heart rate’, ‘Heart rate variability traits’; (Composite set) ‘Body mass index’, ‘Waist-hip ratio’, ‘Triglycerides’, ‘LDL cholesterol’, ‘LDL cholesterol levels’, ‘HDL cholesterol’, ‘HDL cholesterol levels’, ‘Total cholesterol levels’.

FUN-LDA tissue enrichment

We computed enrichment for SNPs belonging to regions of open chromatin (more likely to contain expressed genes[80, 81]) in specific tissues in three cases: i) when comparing results across phenotype sets, ii) when comparing univariate results, and iii) when comparing results across clusters. For all analyses we used functional annotations on 127 Roadmap tissues and cell lines defined by integrating activating histone marks (H3K4me1, H3K4me3, H3K9ac, and H3K27ac) with a latent Dirichlet allocation model as implemented in FUN-LDA[82]. The enrichment score for a tissue is based on the number significant SNPs compared with the total number of SNPs in open chromatin region (see S1 Text). Enrichment results are reported in S11–S13 Tables.

Multitrait genetic association clustering and selection of the optimal number of clusters

We performed a clustering of top associated SNPs for each phenotype set using a Gaussian Mixture model (GMM). Briefly, the assumed generative model is as follows. Consider n independent vectors Z_i = (Z_i1,…,Z_ij,…,Z_id)^t in ℝ^d, each arising from one of k distinct clusters. Each vector represents one SNPs. The value Z_ij is the Z-score of SNP i on trait j. Let I_ij = 1 if the i^th observation belongs to cluster j, and define the indicator vector I_i = (I_i1,…,I_ij,…,I_ik)^t. Conditional on membership to the j^th cluster, Z_i follows a multivariate normal distribution, with cluster specific mean μ_j and covariance Σ_j. Let denote π^j the marginal probability of membership to the j^th cluster. The observations can be viewed as arising from the following hierarchical model:

One major difficulty in applying the GMM was to deal with incomplete data. Indeed, even after imputation of some missing statistics, our datasets still contained some missing values. To solve the clustering, we implemented the statistical framework described by Ghahramani et al[83] which we recently implemented in a R package MGMM[40]. This method relies on EM optimization techniques enabling the inference of unobserved variables from observed variables and an assumed Gaussian mixture model. In classical GMM, the only variable inferred is the posterior probability cluster membership. In the MGMM algorithm, missing Z-scores are also inferred taking into account the observed Z-scores and the inferred probability to belong to a cluster.

The model gives for each SNP the posterior probabilities to belong to each cluster, and was therefore assigned to its most likely cluster, as long as its entropy was larger than 0.75. For a given variant SNP_i, the entropy was derived as follow: where k is the total number of clusters and P(SNP_i∈cluster_j) is the posterior probability of SNP_i to belong to cluster j. The higher the entropy the more the SNP attribution to one cluster is ambiguous. SNPs with an entropy higher than 0.75 were filtered out of the clustering results.

Clustering was performed on all independent significant SNPs. For each SNP, we defined three p-values on the phenotypic group traits: the minimum univariate p-value (P_univ), the SumZ_ica p-value and the omnibus p-value. All SNPs with at least one of the three p-value under 10⁻⁸ were selected for further analysis. For the Metabolism univariate clustering, we only considered the univariate p-value to perform the selection. We then applied the plink[84] clump function to retrieve practically independent associations using the 0.2 as clump-r2 parameter and 10⁻⁸ as clump-p1 parameter. For each clump we selected a representative SNPs as the one with the smallest p-value across the three tests and having more than 60% of its values observed. Note that for a negligible number of occurrences, the representative SNPs has a p-value above 10⁻⁸ (S15 and S16 Tables). We applied MGMM within each phenotype set and varied the pre-specified number of clusters between 2 and 10. To select the optimal number of clusters k, we performed the clustering 100 times on a random subset of 80% of the SNPs for each k. For each resulting clustering we computed the Bayesian Information Criteria and the Silhouette[85] (see S23 and S24 Figs). Except for the Metabolism set, the silhouette appears conservative and the BIC criterion anticonservative, i.e. the latter criteria tends to select a larger number of clusters. We decided to use the following ad hoc compounded criterion:

1. If the optimal number of clusters determined by the BIC criteria is higher than the one determined by the silhouette criteria, starting from the silhouette optimal, increase the number of clusters until one of these two conditions is met: 1) adding one cluster significantly decrease the silhouette criterion, 2) the BIC optimal number is reached.
2. In other cases, set the optimal number of clusters to the one determined by silhouette.

Cluster genetic correlation

We defined pairwise genetic covariance per cluster for as where β₁ and β₂ are the vector of genetic effects for the pair of phenotypes considered and M is the number of SNPs in the cluster. To estimate properly this quantity from the observed , we accounted for the bias introduced by sample overlap and phenotypic correlation using the following estimator (see S1 Text): where ρ_Y is the phenotypic covariance, and n_s, n₁ and n₂ are respectively the sample size shared between the two traits, for the trait 1, and for the trait 2. To assess whether the estimated genetic covariances are significantly different from zero, we performed for each pair of phenotypes within each cluster, a t-test on the vector of random variables (X₁, X₂,…,X_M), were is the contribution of SNP j to the covariance. Note that we used only independent SNPs selected using LD-clumping with squared-correlation parameter equals 0.2.

Functional enrichment of metabolism clusters

We used FUMA[86] SNP2GENE function to associate SNPs with genes based on two criteria, the physical position (in 30kb radius of a protein coding gene) and eQTLs (all significant cis-eQTL from GTEx up to a distance of 1Mb). Note that we restrained the eQTLs to the one that were found in relevant tissue for the Immunity and Metabolism set: immune cells for Immunity and adipose, intestine, liver and brain tissues for Metabolism (see S1 Data for complete parameters). After chaining genes to clusters based on SNPs, we performed a functional enrichment for pathways defined in KEGG[87] and GO[88] databases and derived report p-values using FUMA GENE2FUNC function. Here, cluster’s gene were compared against a background of protein coding genes. Finally, we used the R package pathview[89] to project genes onto KEGG pathways maps.

Disease-clusters association

For the metabolism phenotype set, to provide an indicator of the relative contribution of genetic variants to phenotypes in each cluster from the Metabolism set, we performed a principal component analysis (PCA) of the SNP-by-phenotype association matrix within each cluster. For this analysis, we used scaled beta coefficients, i.e. Z-scores divided by the square root of the phenotype GWAS sample size. To avoid bias due to the arbitrary choice of the coded allele, we randomly shuffled 20 times the coded allele, and repeated the PCA after each shuffling. We report in Fig 5, the average of the loadings of the first PC over all shuffling. Note that the first PC only provides the multidimensional direction explaining the largest variance and therefore do not fully capture the distribution of genetic effect within each cluster. Nevertheless, those first PCs explained a substantial amount of the total variance, equal to 75%, 38%, 53%, 64%, 80% and 93% of the variance in betas for cluster 1 to 6, respectively.

Then, we assessed the association between SNPs within the inferred cluster and three traits (none of which being included in the Metabolism set): cardiovascular diseases, any stokes and BMI. SNP alleles were aligned according to the first principal by clusters determined in the last section. We applied a sign test to assess the concordance of the sign of the projection on PC1 and the sign of Z-score for on the three tested additional traits. For this analysis we used more stringent criteria to ensure the SNPs independence. We selected the subset of Metabolism SNPs for which linkage disequilibrium does not exceed 0.2 (clump-r2 set to 0.05), which diminishes the number of SNPs considered from 391 to 285. Concerning the association of SNPs to drug target, we associated drug target to a representative SNPs by selecting the SNP with the lowest entropy and having a positive silhouette.

URL Resources

We developed the following R and python packages to perform our study:

JASS_Preprocessing: https://gitlab.pasteur.fr/statistical-genetics/jass_preprocessing

JASS: https://gitlab.pasteur.fr/statistical-genetics/jass

RAISS: https://gitlab.pasteur.fr/statistical-genetics/raiss

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

Figs 3–5 have been generated using the following R packages:

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

https://cran.r-project.org/web/packages/radarchart/README.html

https://cran.r-project.org/web/packages/ggalluvial/vignettes/ggalluvial.html

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

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