Cross-GWAS coherence test at the gene and pathway level

2.1 Coherence test

Consider the index I = Σ_i w_iz_i, with w_i and z_i N independent samples of two random variables z, w ∼ 𝒩 (0, 1). The index can also be written as I = N 𝔼(wz), with 𝔼 denoting the expectation. Clearly, I is proportional to the standard empirical covariance of w and z.

For independent pairs of samples, the sampling distribution of I is simply a sum of independent product-normal distributions. Hence one can infer that for identically correlated random variables with correlation coefficient ϱ (cf., Methods, Product-Normal distribution) In particular, for ϱ = 0, we have that I ∼ V G(N, 1), with V G being the variancegamma distribution discussed in more detail in the supplement. One should note that the difference is in the distributional sense and, therefore, generally non-vanishing. A null hypothesis of zero correlation (or some other fixed value) can therefore be tested, as the cumulative distribution function (cdf) for I can be calculated explicitly and efficiently with Davies’ algorithm (cf., Methods, Linear combination of χ² distributions). The main advantage of the above significance test is that it is straightforward to relax the requirement of sample independence. That is, we can view the index I as a scalar product of random samples of w ∼ 𝒩 (0, Σ_w) and z ∼ 𝒩 (0, Σ_z), with 𝒩 denoting here the multivariate Gaussian distribution and Σ_w|z covariance matrices. For Σ_w = Σ_z = 1, the corresponding distribution of I is given by (2.1). In the general case, the inter-dependencies can be corrected via linear decorrelation, cf., Methods Coherence test decorrelation. The case of interest for this paper is Σ := Σ_w = Σ_z. In this case, one can show that under the null of w and z being independent with λ_i the ith eigenvalue of Σ. Hence, I is distributed according to a linear combination of distributions with positive and negative coefficients, and therefore the cdf and tail probability can be calculated with Davies’ algorithm. Note that the above discussion can be extended to non-standardized variables (for ϱ = 0) via the results given in supplementary section I.

As discussed in detail in [2], a GWAS (cf., Methods, GWAS) gene enrichment test can be performed via testing against I with w = z. This effectively tests against the expected variance of SNPs’ significances in the gene. Here we propose to use I = w^T z with w and z resulting from two different GWAS phenotypes to test for the co-significance of a gene for two GWAS. A significance test can be performed either against the right tail of the null distribution (2.2) (coherence) or against the left tail (anti-coherence). For a single SNP in the gene window, the test reduces to the product-normal test, whose properties will be discussed in more detail in the following section 2.2.

Note that we do not centralize w and z over the gene SNPs. Hence, we do not test for covariance but for a non-vanishing second cross-moment. After decorrelation, it is best to interpret this as testing each joint SNP within the gene independently for a coherent deviation from the null hypothesis. Therefore we refer to this test as testing for genetic coherence or simply as cross-scoring. For simplicity, we only consider a fixed effect size model and assume that the correlation matrix Σ obtained from an external reference panel is a good approximation for both GWAS populations.

The above coherence test assumes that there is no sample overlap between the populations of the pair of GWAS considered. However, the presence of sample overlap can be corrected for, as discussed in supplementary section 4.

2.2 Simulation study

It is useful to discuss the single SNP case in more detail. The distribution (2.2) simplifies for N = 1 to , which corresponds to the uncorrelated product normal distribution, cf., Methods eq. (4.4). The index I for N = 1 is a measure of coherence (or anti-coherence) between z and w. The significance thresh-old curve for a fixed desired p-value, say p_I = 10⁻⁷, is illustrated in figure 1. The curve corresponding to a given p_I is unbounded. For a given w, there is always a corresponding z such that the resulting product I is significant. This differs from Fisher’s exact method which combines two p-values p_w,z into a combined one (p_F) via . Since Fisher’s method combines significance by addition, the corresponding combined significance curve is bounded. Specifically, in the extreme case of one of the p-values being equal to one, the significance threshold is finite and fixed by the other p-value. In contrast, in the case of the product-normal, the divergence between the p-values is penalized. If one of the p-values is large, say p_w ≃ 1, the other one has to be extremely small to achieve a given combined significance value, i.e. p_z ≪ p_w (cf., figure 1). Therefore, the product-normal-based method becomes more and more conservative for increasingly diverging p-values. For very similar p-values, the method has a lower threshold than Fisher’s, p-values reaching the same minimum. Therefore, one should see Fisher’s method as additive in the evidence, while the product-normal-based method is multiplicative.

• Download figure
• Open in new tab

Figure 1:

Significance threshold curves in the one element case for the product-normal (dashed) and Fisher’s method (dotted) for various p-values. The diagonal is indicated by a gray dash-dotted line and corresponds to equal − log₁₀ p-values.

The importance of correcting for the inter-dependence between w and z elements in the index I is already visible for N = 2, as discussed in supplement section III. For higher N, this becomes even more pronounced. For example, consider a correlation matrix Σ of dimension one hundred with off-diagonal elements identically set to 0.2. We draw 1000 pairs of independent samples of 𝒩 (0, Σ) and calculate I for each pair. A p-value is then obtained for each index value for the linear combination of χ² distributions (2.2) (also referred to as weighted χ²), and for the variance-gamma distribution. Recall that the latter does not correct for the off-diagonal correlations. We repeat the experiment with the off-diagonal elements set to 0.8, resulting in a stronger element-wise correlation. The resulting qq-plots for both cases are shown in figure 2. We observe that the variance-gamma distribution (2.1) (with ϱ = 0) indeed becomes unsuitable for increasing element-wise correlation of the data sample elements. Not correcting for the inter-sample correlation leads to more and more false positives with increasing correlation strength. In contrast, the weighted χ² distribution (2.2) yields stable results in both the weakly and strongly correlated regime, as is evident in figure 2.

• Download figure
• Open in new tab

Figure 2:

QQ-plots of observed p-values resulting from the index I for 1000 pairs of samples of 𝒩 (0, Σ) against uniform p-values. Top: Off-diagonal elements of Σ set to 0.2. Bottom: Off-diagonal set to 0.8. The blue curve is obtained using the variance-gamma distribution to perform the statistical test, while the orange curve is obtained via the weighted χ² distribution. The latter corrects for the correlation and therefore is well calibrated.

2.3 Ratio test

Consider the normalized index with w and z as in the previous sections. In particular, being independent random variables. The cdf for R can be calculated to be given by With , Λ the matrix of eigenvalues of Σ, and the linear combination of cdf evaluated at the origin, see Methods, Ratio test derivation (Similar results can be obtained for w and z interchanged). Hence, the cdf of the ratio (2.3) can also be calculated with Davies’ algorithm. A consistency check follows from the case of one dimension, where R has to be Cauchy distributed. The corresponding cdf is given by F_C = 0.5 + arctan(r)/π. Evaluation for various r shows agreement with values calculated from (4.10) in the one-dimensional case. Note that similar expressions can be derived for non-standardized w_i and z_i, albeit in terms of the non-central χ² distribution, cf., the supplement.

Note that the ratio with and i.i.d. 𝒩 (0, 1), and λ_i the ith eigenvalue of Σ, can be interpreted as the weighted least squares solution in the case of heteroscedasticity for the regression coefficient of the linear regression between the de-correlated effect sizes of the two GWAS. Therefore, with the cdf for R derived above, we can test for a significant deviation from the null expectation of no relation. In general R is not invariant when swapping w and z and may be used under certain conditions to make inference about the causal direction, cf., multi-instrument Mendelian randomization, in particular [13]. Specifically, we consider the trait related to w as exposure and that of z as the outcome. Then the ratio test tries to detect genes that are strongly (anti)-coherent between exposure and outcome but at the same time only exhibit relatively low variance in the outcome, as is clear from the definition of R. Thus the ratio test effectively normalizes the alignment of the outcome to the exposure, with respect to the outcome’s own variation. The causal direction from exposure to outcome is implied if the exposure is confirmed to be associated with the gene via a significant gene enrichment p-value, p_V, obtainable from the usual χ² test of [2], and if confounding factors can be excluded. The overall scheme is illustrated in figure 3.

• Download figure
• Open in new tab

Figure 3:

Interpretation of the ratio test. The test detects potential gene-wise causal relations between a GWAS trait viewed as exposure (E) and a GWAS trait viewed as outcome (O). The association of the exposure has to be confirmed independently via a standard gene enrichment test. Potential confounders (C) have to be excluded by other means.

Similar to the previously introduced coherence test, we assumed above that the populations of the pair of GWAS considered present no overlap. However, sample overlap can be corrected for as well, cf., supplementary section 4.

2.4 COVID-19 GWAS application

To demonstrate the usefulness of our methods in the context of actual and topical GWAS data (cf., Methods, GWAS), we considered the recent meta-GWAS on very severe respiratory confirmed COVID-19 [14, 15] as the primary phenotype and co-analyzed it with several other related traits (see also Data availability statement). Specifically, we used the summary statistics resulting from European subjects, excluding those from UK Biobank (A2 ALL eur leave ukbb 23andme) in order to avoid overlap with the secondary trait GWAS. This GWAS shows significant gene enrichment on chromosomes 3 and 12, as can be inferred via testing for gene enrichment following [2]. The Manhattan plot for the resulting gene p-values is shown in supplementary figure 7.

We cross-scored this GWAS against a panel of GWAS on medication within the UK Biobank [16], which take prescription (self-reported intake) of 23 common types of medications as traits. Hence, in cross-scoring against the severe COVID-19 GWAS, we sought to uncover whether there is a shared genetic architecture of predisposition for severe COVID and being prescribed specific medications. The description of all medication class codes used can be found in supplementary table 2. Note that, as usual, these GWAS have been performed with age as one of the covariates. Therefore, in our pair-wise analysis age-dependent effects are already regressed out and do not play a significant role in our analysis.

All calculations were performed with the python package PascalX [17], which incorporates the methods detailed in the Methods section (see also Code availability statement). The association directions were extracted from the signs of the raw effect sizes (β). As a reference panel, we used the European sub-population of the high coverage release of the 1K Genome project [52]. Note that we took only SNPs with matching alleles between the GWAS pairs and the reference panel into account. We transformed the raw GWAS p-values via joint rank transformation, for the reason discussed in Methods under SNP normalisation. We tested all genes with at least one SNP present in both traits, and defined the gene region as the transcription side plus a window of 50kb on both sides in order to capture as well potential regulatory effects. We verified at hand of one drug class (M05B) that the choice of gene window does not lead to an inflation of signal, see figure 9 in the supplement.

The cross-scoring results for the test are shown in figure 4. Cross scoring revealed joint coherent signals between GWAS effects from the prescription of group M05B medications (drugs affecting bone structure and mineralization) and that for very severe COVID-19. We show the Manhattan plot resulting from the null model (2.2) for the medication group M05B in figure 5 (for the corresponding qq-plot, see supplementary figure 8). We observe that SNP-wise effects in genes in the well-known COVID-19 peak locus on chromosome 3 appear to be coherent with those from being prescribed M05B medication, with Bonferroni significance for the chemokine receptor genes CCR1, CCR3 and the gene LZTFL1 in the region chr3p21 (we Bonferroni corrected for number of genes and drug classes tested). We tested the orientation of the aggregated associations of these genes via the D-test, cf., Methods, Direction of association, over the COVID-19 GWAS and find a positive direction (right tail) with p_D ≃ 1.7 × 10⁻⁴, 8.6 × 10⁻³ and 0.03, respectively.

• Download figure
• Open in new tab

Figure 4:

Resulting p-values for cross-scoring 23 drug classes GWAS against very severe COVID-19 GWAS for anti-coherence and (top figure) coherence (bottom figure). Each data point corresponds to a gene. The dotted red line marks a Bonferroni significance threshold of 1.18 × 10⁻⁷ (0.05 divided by the 18453 genes tested and 23 drug classes). Left: Anti-coherence. Right: Coherence. Note that the drug class M05B shows the most significant enrichment in coherence with severe COVID-19.

• Download figure
• Open in new tab

Figure 5:

Manhattan plot for cross scoring very severe confirmed COVID-19 with medication class M05B for coherence. Data points correspond to genes. The dotted red line marks a Bonferroni significance threshold of 1.18 × 10⁻⁷ (0.05 divided by the 18453 genes tested and 23 drug classes). The inlay plot shows a zoom into the relevant locus chr3p21.31.

For illustration, we show the spectrum of SNPs considered in the CCR3 region and their SNP-SNP correlation matrix in figure 6. One can see a large block of SNPs in high LD, encompassing some of CCR3’s 5’UTR and most of its gene body, all having positive associations with both severe COVID-19 and M05B drug prescription. Furthermore, a somewhat weaker coherence signal can be seen for SNPs in its 3’UTR, some of which share negative associations. Using our coherence test described above, we calculate a p-value of p ≃ 8.2 × 10⁻⁹ for CCR1 and p ≃ 3.3 × 10⁻⁹ for CCR3. The SNP spectrum in the LZTFL1 region is shown in supplementary figure 10. The p-value for the significance of the coherence is p ≃ 2.7 × 10⁻⁸.

• Download figure
• Open in new tab

Figure 6:

Left: SNP p-values after rank transform for the CCR3 gene. The x-axis is numbered according to the ith SNP in the gene window ordered by increasing position. The dotted black lines indicate the transcription start and end positions (first and last SNP). Up bars correspond to M05B and down bars to severe COVID-19. Green indicates positive and violet negative association with the trait. Right: SNP-SNP correlation matrix inferred from the 1KG reference panel. The color scale is shown on the right with dark red (+1) corresponding to perfect SNP-SNP correlation and dark blue (−1) to perfect anti-correlation over the reference population.

The directions of associations inferred from the D-test and their coherence suggest that genetic predispositions leading to M05B prescription may carry a higher risk for severe COVID-19. We list the Bonferroni significant pathways for (anti)-coherent M05B and severe COVID-19 genes detected via a pathway enrichment test (Methods, Pathway enrichment) in table 1. Note that the enrichment test of pathways is performed with the (anti)-coherent gene enrichment p-values, and therefore tests for pathways that are relevant for both traits simultaneously (we Bonferroni corrected for number of pathways and number of traits we tested for pathway enrichment)

The cross-scoring test also detected the gene HLA-DQA1 to be Bonferroni significant under anti-coherence for the drug classes H03A (thyroid preparations), as well as R03A and R03BA (drugs for obstructive airway diseases). Further, we note that the drug classes C10AA (HMG CoA reductase inhibitors) and L04 (immunosuppressants) have gene hits with a p-value < 1 × 10⁻⁵. Interestingly, all these drug classes possess indications in autoimmune diseases and allergies. In particular, of the non-significant genes we detect the gene FUT1 with p ≃ 2.9 × 10⁻⁷ for C10AA to be closest to Bonferroni significance. The D-test shows that FUT1 is located in the left-tail under C10AA (p_D ≃ 9.8 × 10⁻⁴), and therefore this gene’s variations may be a possible risk factor for severe COVID-19. It is known that FUT1 is involved in creating a precursor of the H antigen, itself a precursor to each of the ABO blood group antigens. It has been hypothesized before in the literature that the ABO blood system correlates to COVID-19 severity [18].

We also find the observed leading TRIM genes for L04 of interest. For illustration of the anti-coherence case, the SNP spectrum for TRIM26 (p ≃ 1.1 × 10⁻⁶) under L04 is shown in supplement figure 12. The pathway enrichment test in the anti-coherent case shows that L04 cross-scored severe COVID-19 has Bonferroni significant enrichment in interferon gamma-related pathways, see table 1. Note that TRIM proteins are expressed in response to interferons [19], and therefore the detected pathways are consistent with the enrichment of TRIM genes.

2.5 Osteoporosis related GWAS

The main application of M05B medications is the treatment of osteoporosis. In order to investigate this potential link further, we cross-scored the COVID-19 GWAS against a selection of GWAS with phenotypes related to osteoporosis, namely, bone mineral density (BMD) estimated from quantitative heel ultrasounds and fractures [20], estrogen levels in men (estradiol and estrone) [21], calcium concentration [22], vitamin D (25OHD) concentration [23] and rheumatoid arthritis (RA) [24]. The inferred gene enrichments for coherence and anti-coherence with the COVID-19 GWAS are shown in figure 7. We found Bonferroni significant enrichment for RA and enrichment with gene p-values < 10⁻⁵ for calcium and vitamin D (we Bonferroni corrected for number of genes and traits tested). The Manhattan plot for the Bonferroni significant trait is shown in figure 8 (for the qq-plot, see supplementary figure 11). Note that ≈ 47% of the SNP alleles between the RA and COVID-19 GWAS were not matching and that we discarded all non-matching SNPs for the analysis. A list of the most significant detected genes under the test is given in table 2. The table also contains a brief description of the potential relation of the respective gene to COVID-19, if known.

• Download figure
• Open in new tab

Figure 7:

Resulting p-values for cross-scoring several GWAS related to osteoporosis against the severe COVID-19 GWAS. The top figure shows the anti-coherent case and the bottom figure the coherent case. The red dotted line marks the Bonferroni significance threshold of 3.9 × 10⁻⁷ (0.05 divided by the 18453 genes tested and 7 traits). We observe Bonferroni significant enrichment for RA plus enrichments in vitamin D and calcium with p-values < 10⁻⁵.

• Download figure
• Open in new tab

Figure 8:

Manhattan plot for cross-scoring severe confirmed COVID-19 with rheumatoid arthritis for anti-coherence. Data points correspond to genes. The dotted red line marks a Bonferroni significance threshold of 3.9 × 10⁻⁷ (0.05 divided by the 18453 genes tested and 7 traits). The labels denote the leading gene hit for the corresponding peak.

For RA, several TRIM genes are Bonferroni significant. TRIM proteins are associated with innate immunity, and are in particular involved in pathogen recognition and host defense [19, 25]. This is consistent with the (weaker) TRIM signals we detected for the immunosuppressants L04, which are indicated in the treatment of RA. The D-test shows that the leading TRIM genes are in the right tail under the RA GWAS, with p_D < 10⁻⁶. Therefore, the anti-coherence suggests a protective effect of variants related to a predisposition for RA.

For the calcium anti-coherent case, the top gene HGFAC with p ≃ 4.0 × 10⁻⁷ is only slightly below the Bonferroni significance threshold. This gene sits in the right tail (p_D ≃ 8.9 × 10⁻⁴) under the D-test (4.8) applied to the calcium GWAS. Therefore, the aggregated variants in this gene imply that a predisposition for high calcium concentration may implicate a reduced risk for severe COVID-19, i.e., a protective effect.

Let us also briefly discuss an application of the ratio-based causality test. We ratio tested the COVID-19 GWAS against M05B, RA, vitamin D, and calcium in the coherent and anti-coherent case, using eq. (2.4). The resulting gene hits of interest are summarized in table 3, with a brief description of the potential COVID-19 context of the gene hits. Unless not indicated otherwise in the table, confounding between calcium and vitamin D could be excluded. The Manhattan plot resulting from the ratio test for Vitamin D as exposure and COVID-19 as outcome is shown in figure 9 (the corresponding qq-plot is shown in supplementary figure 13). We found that Vitamin D carries two interesting hits suggesting causal pathways from Vitamin D concentration to the severity of COVID-19, namely, the genes KLC1 and ZFYVE21. The observed p_R-values suggest a causal flow from a genetic predisposition for vitamin D concentration to the severity of COVID-19, mediated via these genes.

• Download figure
• Open in new tab

Figure 9:

Manhattan plot for ratio scoring Vitamin D concentration against severe COVID-19 in the coherent case, with ratio denominator given by COVID-19. The Bonferroni significance threshold of p = 3.4 × 10⁻⁷ is indicated by the red dotted line (0.05 divided by the 18453 genes and 2 times 4 traits tested).

In the anti-coherent case, the gene HOXC4 is detected for M05B. Modulo potential confounders, a causal relation from a predisposition to take M05B to the severity of COVID-19 via HOXC4 is suggested. Interestingly, we also find an inverse causal relation, namely, a predisposition for severe COVID-19 implying a predisposition for calcium concentration, mediated by the gene DPP9. For a brief discussion of the potential role played by these genes in the COVID-19 context, we refer to table 3.

We also tested all cases against pathway enrichment (see Methods, Pathway enrichment). We only detected Bonferroni significant pathways for RA, see table 4. Note that while we did not detect genes of significance at the gene level under the ratio test for RA, the non-significant genes combine to Bonferroni significant pathways for the trait (we Bonferroni corrected for number of pathways and number of traits pathway enrichment has been tested for). We observed the detected pathways already in table 1 to be co-enriched for immune system-related drug classes and severe COVID-19. The ratio test implies that these pathways may play a causal role.

2.6 Multiplicative meta-analysis

We investigated how the gene-wise coherence test based on the product-normal statistic compares to results from the more simple procedure of first computing gene-wise significance for each trait (via the usual χ²-test, cf., [2]), and then testing the corresponding combined scores via the product-normal as if there was only a single genetic element affecting both traits (cf., fig. 1 for an illustration of the corresponding combined product normal p-values). However, in doing so, we loose the information on the direction of coherence. As before, we used jointly qq-normalized GWAS p-values to avoid the risk of unintended uplift and only kept SNPs with consistent alleles between the GWAS pairs.

In general, we expect such a simple meta-analysis approach to generate both more false positives and false negatives. It is clear that ignoring the coherence of the contributing SNP-wise effects in the gene window can create false positives. However, false negatives may also occur: In our proposed novel coherence test, introduced in the Results section, several independent weak signals may combine into a stronger signal if the directions of associations are consistent over the gene window. Such signals are invisible in the simple approach.

Figure 10 shows scatter plots for gene-wise log-transformed p-values for some of the GWAS we discussed above, namely, the severe COVID-19 GWAS against the M05B medication class and rheumatoid arthritis GWAS. Significance thresholds for the (simple) product-normal are also indicated. We note that the boundaries of the scattered p-value distributions appear to follow product-normal significance curves, justifying our choice to take the single-element product-normal to define significance thresholds for the combined χ² statistic.

• Download figure
• Open in new tab

Figure 10:

GWAS gene scores, obtained separately for two GWAS using joint qq-normalization and standard χ² based enrichment test, plotted against each other. Each point corresponds to a gene. The gray dotted line marks the diagonal corresponding to equal − log₁₀ p-values. The orange and red dashed curves mark the simple gene-wise product-normal threshold curves for significance of p < 10⁻⁶ and Bonferroni significance, respectively. The color of a gene (orange and red) indicates the full SNP-wise cross scored p-value (p < 10⁻⁶ and Bonferroni significant, respectively).

The left plot of figure 10 shows the co-analysis of COVID-19 and M05B GWAS signals. We observe that CCR1 and CCR3 are also Bonferroni significant under the simple meta-analysis test, however not LZTFL1. Also, while CCR9 and CXCR6 cross the p < 10⁻⁶ threshold under the coherence test, not so under the simple single-element product-normal test. This confirms that the coherence test can improve the signal strength and may lead to fewer false negatives. The right plot of figure 10 shows severe COVID-19 against RA. While, again, all Bonferroni significant genes are also Bonferroni significant under the single-element test, we observe several additional genes which appear to be Bonferroni significant under the simple test. For instance, it appears that the single element test cannot completely resolve the TRIM gene cluster and hence may lead to false positives.

4.1 Linear combination of χ² distributions

We denote the χ²-distribution with n degrees of freedom as . It is well known that the sum of N independent distributed random variables v_i is also χ² distributed, i.e., However, no closed analytic expression is known for the distribution Ξ of a general linear combination where a_i are real coefficients. Nevertheless, various numerical algorithms exist to compute the cdf of Ξ, denoted as F_Ξ, up to a desired precision. Perhaps most well known are Ruben’s algorithm [45, 46, 47] and Davies’ algorithm [11, 12]. The latter is the most relevant for this work, as it allows for negative a_i.

With F_Ξ at hand, for a given real x a right tail probability p (p-value) can be calculated as

4.2 Product-Normal distribution

The product-normal distribution, the distribution of the product of two normaldistributed random variables w and z, plays a central role in this work. The moment generating function for joint normal samples with correlation ϱ reads [48] Our key observation is that the above moment generating function factorizes into moment generating functions of the gamma distribution, , i.e., Therefore, For general parameters of the two gamma distributions, the corresponding difference distribution is known as the bilateral gamma distribution [49]. (Note that the subtraction in equation (4.3) is in the distributional sense, so, even for ϱ = 0, the corresponding distribution does not vanish.)

Due to the well known relation between the gamma and χ² distributions, one can also express the product-normal distribution in terms of the χ² distribution introduced in the previous section. In detail, The cdf of the product-normal can therefore be efficiently calculated using Davies’ algorithm, as the distribution (4.4) is simply a linear combination of distributions. A similar relation can be derived for the product distribution of non-standardized Gaussian variables, albeit in terms of the non-central χ² distribution, cf., supplementary section I. Note that the relation (4.3) allows for a simple analytic derivation of a closed-form solution for the product-normal pdf, but not for the cdf. For completeness, details can be found in the supplement, section II.

4.3 Coherence test decorrelation

We make use of the eigenvalue decompositions and , with Λ. the diagonal matrix of eigenvalues of Σ., to decorrelate the elements of each set.

The index I can then be written as with and . In components, I reads For Σ_w = Σ_z =: Σ the matrix K is the identity matrix, and the result given in equation (2.2) follows.

4.4 GWAS

As explained in the introduction, a prime example of very strong inter-element correlations are SNPs in LD. Recall that the univariate least squares estimates of the effect sizes β in genome-wide association studies (GWAS) reads with x_i the ith column of the genotype matrix X of dimension (n, p) and y the phenotype vector of dimension n. Both x and y are mean-centered and standardized. n is the number of samples and p the number of SNPs. The central limit theorem and standardization ensures that for n sufficiently large .

As a multi-variate model, we have with α the vector of p true effect sizes and ϵ the n-dimensional vector of residuals with components assumed to be ϵ_i ∼ 𝒩 (0, 1) and independent. Substituting the multivariate model into (4.5), yields As a fixed effect size model, under the null assumption that α = 0 (no effects), we infer that We can stack an arbitrary collection of such z_i to a vector z via stacking the x_i to a matrix x, such that with . Note that we made use of the affine transformation property of the multi-variate normal distribution.

It is important to be aware that z is only a component-wise univariate estimation, and hence the null model (4.6) is for a collection of SNPs with effect sizes estimated via independent regressions.

We can also take the effects to be random variables themselves. Let us assume that independently α_i ∼ 𝒩 (0, h²/p), with h referred to as heritability. We then have that such that with , and where we assumed that α and ϵ are independent. Two remarks are in order. As X runs over all SNPs, the calculation of L usually requires an approximation, for instance, via a cutoff. Furthermore, the null model (4.7) requires an estimate of the heritability. Such an estimate can be obtained for via LD score regression [9].

4.5 Direction of association

The direction of effect of the aggregated gene SNPs can be estimated via the index In detail, making use of the Cholesky decomposition Σ = CC^T and the affine transform property of the multi-variate Gaussian, we have that as null with |.|_F the Frobenius norm. Testing for deviations from D in the right or left tail indicates the direction of the aggregated effect size. Note that the (anti)-coherence test between pairs of GWAS introduced above is, alone, not sufficient to determine the direction for a GWAS pair but requires, in addition, testing at least one GWAS via (4.8) to determine the base direction of the aggregated gene effect. Furthermore, at least one GWAS needs to carry a sufficiently oriented signal in the gene such that (4.8) can succeed. We will also refer to testing for deviations from (4.9) as D-test. (Note that we regularize Σ via thresholding eigenvalues smaller than 10⁻⁸ for the D-test.)

4.6 Ratio test derivation

Clearly, for Σ_w = Σ_z = Σ we have that We define such that . Note that the component-wise correlation coefficient ϱ between and reads Hence, from eq. (4.4) and the supplementary section I, we deduce that We conclude that (2.4) holds.

4.7 Pathway enrichment

The gene scores resulting from the above coherence or ratio test can be utilized instead of the usual gene scores to perform a gene set (pathway) enrichment test. A pathway is thereby tested for enrichment in coherent or causal genes for a GWAS trait pair. We follow the pathway scoring methodology of [2]. Genes of a pathway close to each other are fused to so-called meta-genes, and coherence or ratio test-based gene scores are re-computed for the fused genes. The purpose of the fusion is to correct for dependencies between the gene scores due to LD. The resulting gene scores (p-values) are qq-normalized and inverse transformed to distributed random variables. This is followed by testing against the distribution, with n the total number of (meta)-genes in the pathway. In this work, we tested for pathway enrichment against MSigDB 7.4 [50, 51].

4.8 SNP normalisation

As discussed in the Results section, the product normal combines evidence for coherent association in a multiplicative manner. A potential challenge to the proposed method arises when one of the two GWAS has associations with very low p-values. Such highly significant associations are common for GWAS with very large sample sizes. Without moderating these p-values, such associations may appear nominally co-significant as soon as the other GWAS provides a mild significance level. We propose two possible strategies to mitigate this.

One strategy is introducing a hard cutoff for very small SNP-wise p-values. The precise cutoff depends on the desired co-significance to achieve, and the amount of possible uplift of large p-values one finds acceptable. The dynamic is clear from figure 1. If we target a co-significance of p = 10⁻⁸ and accept to consider SNPs with a p-value of 0.05 or less in one GWAS to be sufficiently significant, we have to cut off p-values around 10⁻¹⁶. While such a cutoff ensures that no SNPs with p-values above 0.05 in one GWAS can become co-significant due to very high significance (i.e., p-value below 10⁻¹⁶) in the other GWAS, applying such a hard cutoff point hampers distinguishing differences in co-significance.

An alternative strategy is to transform all p-values, rather than only the most significant ones. For example, the so-called qq-normalisation, re-assigns uniformly distributed p-values according to the rank r, i.e. p = (r + 1)/(N + 1), where N is the number of p-values. (This approach implies the strongest p-value moderation and therefore is very conservative.) The product-normal statistic in (2.2) is then computed with w and z according to the inverse χ² cdf of the respective p-values. Since for GWAS usually N ≃ 10⁶, the most significant transformed p-value is ∼ 10⁻⁶. According to figure 1 the other p-value has then to be smaller than 10⁻³ − 10⁻⁴ to achieve (genome-wide) co-significance.

Since the qq-normalisation allows for combining two GWAS with significantly different signal strengths without the need to introduce an adhoc cutoff, it is our approach of choice and will be used in this work. We also prefer to apply this transformation for the ratio test in order to compensate for different signal strengths between the GWAS.