## Abstract

We devise a significance test for covariance of samples not drawn independently, but with known inter-sample covariance structure. We propose a test distribution which is a linear combination of *χ*^{2} distributions, with positive and negative coefficients. The corresponding cumulative distribution function can be efficiently calculated with Davies’ algorithm at high precision. As an application, we suggest a test for dependence between SNP-wise effect sizes of two genome-wide association studies at the level of genes. This test can be extended to detect gene-wise causal links. We illustrate this method by uncovering potential shared genetic links between severity of COVID-19 and (1) being prescribed class M05B medication (drugs affecting bone structure and mineralization), (2) vitamin D (25OHD) and (3) serum calcium concentrations. Our method detects a potential role played by chemokine receptor genes linked to *T*_{H}1 versus *T*_{H}2 immune reaction, a gene related to integrin beta-1 cell surface expression, and other genes potentially impacting severity of COVID-19.

## I. INTRODUCTION

Pearson’s sample correlation is one of the most used techniques in data analysis, independent of the specific field of science. It is defined as the covariance divided by the standard deviations of two random variables and gives a measure of linear dependence. The core underlying requirement is that the observed samples are drawn independently from a joint bivariate distribution. However, not all applications possess independently drawn samples.

One particular application with dependent samples occurs in genome-wide association studies (GWAS). Specifically, GWAS correlate genotypes, most commonly single nucleotide polymorphisms (SNPs), with a phenotype of interest, both measured in the same study population. For human studies usually between 1 and 10 million SNPs are considered, and in most GWAS each of these SNPs is tested independently for correlation with the phenotype. By now, thousands of such GWAS have been conducted and identified a plethora of statistically significant associations of SNPs with complex traits. For most traits, in particular complex diseases or their risk factors, that have been assessed in very large populations (100k-1M subjects), usually hundreds of SNPs turn out to be significant, even after stringent correction for multiple hypotheses testing. Individual SNP-wise effect sizes are often very small, but add up to sizable narrow sense heritability, pointing to a polygenic genetic architecture [1]. Mapping SNP-wise effects to genes and annotated gene-sets (*pathways*) [2, 3] can yield valuable insights into the genetic underpinning and potential pathomechanisms of complex diseases and aid drug discovery and re-purposing.

However, due to Linkage Disequilibrium (LD) (*cf*., [4]), many SNPs in close proximity are not independent, and this leads to dependencies between observed SNPs’ effect sizes. This is of particular relevance when aggregating SNP-wise effects to genes or pathways. Specifically, gene-wise effects are computed by adding up the (squared) effects of all SNPs within the transcript region of a gene of interest, as well as sizable up- and downstream-regions that may contain regulatory elements of this gene. Path-way effects are computed from the gene-wise effects [2]. Thus, if not corrected for LD, the resulting gene and pathway scores may reflect the level of importance for the phenotype inaccurately. For this reason, techniques and tools have been developed to correct for LD in the aggregation process, like for example *Pascal* [2] *and MAGMA* [3]. These tools mainly differ in how they map SNP effect sizes to genes, how the LD structure is accounted for, and details of the numerical procedure to estimate significance.

Some SNPs are significantly associated with more than one trait. The phenomenon of a single genetic variant affecting two or more traits is called *pleiotropy* [5]. In the case of disease traits such a shared genetic component hints at the same functional pathology contributing to several diseases, *cf*., [6]. At the gene level, usually a gene is considered to be relevant for two different traits if the respective gene effect sizes are significant in both traits. However, this may not be a good criterion under all circumstances. Protein coding genes often contain several independent LD blocks. Therefore, two traits may associate with genetic variation in two or more functionally different blocks of SNPs within the same gene, which may independently be significant. Hence, even though the two traits share the same significant gene, they may not share the same functional mechanism. To call a gene pleitropic, one should therefore move beyond comparing single variants, and take all SNPs in the gene region into account.

Several methods have been proposed to uncover shared genetic origin of two traits from GWAS summary statistics: One early method is a test of co-localisation between GWAS pairs based on Bayesian statistics [7]. This method assumes that at most one association is present for each trait in the region of interest. The extension to the general case of multiple associations (which is usually the case) appears, however, to be non-trivial. A more recent method is cross-trait *LD score regression* [8], *an extension of single-trait LD score regression* (LDSR) [9], which is a method to estimate heritability and confounding biases from GWAS summary statistics. Like single-trait LDSR, cross-trait LDSR considers the effect sizes as random variables and uses *LD-scores* (*i*.*e*., the sum of genetic correlations between a given SNP and all other SNPs) to estimate the genetic correlation between two traits. Yet this estimate is at the level of the entire genome, and it is difficult to use this approach to obtain estimates at the level of genes. The reason is that SNPs in the same genomic region are often in high LD, so the variables entering the linear correlation may be highly dependent, and this has to be corrected for. Similarly, the standard errors and *p*-values are estimated via resampling (jackknife). But this requires independent samples, which is not the case for strong LD. Another method to estimate local genetic covariance and correlation has been introduced in [10]. However, this method requires the computation of the inverse SNP-SNP correlation matrix, which is in general problematic as it often would need regularization.

Here, we propose to consider the second cross moment between the effects of SNPs for two traits within a gene region, corrected for LD, as a measure for pleiotropy. This gives us a simple, but systematic definition of cosignificance of a gene for two traits, *i*.*e*., if the SNP-wise cross-moment is significant. For a single SNP this reduces to testing against the product-normal distribution, and hence corresponds to a multiplicative meta-analysis, rather than an additive one like Fisher’s. We show that for multiple SNPs with a non-trivial covariance structure, one can express the underlying test distribution in terms of a linear combination of *χ*^{2} distributions, with a mixture of positive and negative coefficients. The corresponding cumulative distribution function can be efficiently calculated with Davies’ algorithm at high precision. Thus, our method takes into account not only isolated significant SNPs, but all common SNPs within the gene region to call a gene co-significant for two traits. The assumption we have to make is that there is an underlying joint normal distribution, and the dependency pattern of the SNPs can be jointly inferred from a reference population. Furthermore, we assume that the study populations of the two traits do not have overlapping samples.

A limitation of our approach is that common genetic effects to different traits cannot be disentangled from potential other joint residual contributions. Yet, using GWAS data from different populations such contributions are likely to be negligible. Using the notion of Mendelian randomization, our statistic can be extended to test for a possible causal relation between the two traits, mediated by the tested gene. However, possible confounders have to be excluded by other means.

The outline of this paper is as follows: In sections II and III we introduce the probability distributions of relevance for this work. We briefly review the *χ*^{2} distribution in section II, since it is essential for our work, as all the other distributions we consider can be expressed as a linear combination of *χ*^{2} distributions. This holds in particular for the product-normal distribution, as we show in section III, as well as for the Variance-Gamma distribution discussed in appendix B. In section IV we explain how to perform a significance test of independence under the conditions described above. In section V we derive how to calculate the cumulative distribution function of a particular ratio distribution of relevance for a causality test. Simulated examples are discussed in section VI, followed by an application to real GWAS data in section VII. We demonstrate the utility of our method by a timely co-analysis of GWAS summary statistics on the severity of COVID-19, being prescribed class M05B medication, and vitamin D and calcium concentrations. In particular, we detect a potential role for severity of COVID-19 played by chemokine receptor genes linked to *T*_{H}1 versus *T*_{H}2 immune reaction, and a gene related to integrin beta-1 cell surface expression. Several other genes related to COVID-19, discussed before elsewhere in the literature, are replicated. In addition we uncover hints for a potential protective and/or therapeutic pathway related to being prescribed specific immune related medications (H03A, L04 and M01A).

## II. LINEAR COMBINATION OF *χ*^{2} DISTRIBUTIONS

We denote the *χ*^{2}-distribution with *n* degrees of freedom as . It is well known that the sum of *N* independent distributed random variables *v*_{i} is again *χ*^{2} distributed, *i*.*e*.,

However, no closed analytic expression for the distribution Ξ of a general linear combination,
where *a*_{i} are real coefficients, is known. Nevertheless, a variety of numerical algorithms exist to compute the cumulative distribution function (cdf) of Ξ, denoted as *F*_{Ξ}, up to a desired precision. Perhaps most well known are Ruben’s algorithm [11–13] and Davies’ algorithm [14]. The latter is of most relevance 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

A variety of approximations to the distribution Ξ exists, *cf*., [15]. One of the more well known methods is the Satterthwaite-Welch approximation [16, 17], which matches the first two moments of Ξ to a Gamma distribution, denoted as Γ. In detail [18],
with

Here, *d*_{i} denotes the *i*th degree-of-freedom parameter of the *i*th *χ*^{2} in equation (1) and *i ∈ {*1, …, *N}*.

To our knowledge, the quality of second order approximations to Ξ at very high precision (below double precision) has not been evaluated in the literature, *cf*., [15]. However, it is expected that in general the approximation breaks down far in the tail.

## III. PRODUCT-NORMAL DISTRIBUTION

A central role in this work is played by the product-normal distribution, the distribution of the product of two normal-distributed random-variables *w* and *z*. The moment generating function for joint normal samples with correlation *ϱ* reads [19]

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 [20]. (Note that the subtraction in equation (4) 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 *χ*^{2} distributions, one can also express the product-normal distribution in terms of the *χ*^{2} 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 (5) 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 *χ*^{2} distribution, *cf*., appendix A. Note that the relation (4) 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 appendix B.

We conclude that the cdf of the product normal distribution can be calculated with Davies algorithm. In particular, note that for *ϱ* = 0 we can not make use of the Satterthwaite-Welch approximation to calculate the cdf, as the first cumulant vanishes and therefore (3) is not well defined.

## IV. COVARIANCE SIGNIFICANCE TEST

Consider the index
with *w*_{i} and *z*_{i} *N* independent samples of two random variables *z, w* ∼ 𝒩(0, 1). The index can also be written as
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 we infer from the previous sections that for identically correlated random variables

In particular, for *ϱ* = 0, we have that *I* ∼ *X*(*N*, 1), with *X* being the Variance-Gamma distribution discussed in more detail in appendix B. As mentioned already before, one should note that the difference is in the distributional sense and therefore generally non-vanishing. A null assumption of zero correlation (or some other fixed value) can therefore be tested via (2), as the cdf for *I* can be calculated explicitly and efficiently with Davies algorithm.

We have now all ingredients in place to discuss the problem we want to address with this paper. The main advantage of the above significance test is that it is straight-forward 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 Σ. covariance matrices. For Σ. = 1, the corresponding distribution of *I* is given by (7). In the general case, the inter-dependencies can be corrected for as follows.

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

The case of interest for this paper is Σ_{w} = Σ_{z} =: Σ such that *K* is the identity matrix. In this case, making use of the moment generating function of the product-normal, as in section III, we can show that under the null of *w* and *z* being independent
with *λ*_{i} the *i*th 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 making use of the result of appendix A.

## V. RATIO SIGNIFICANCE

Consider the normalized index
with *w* and *z* as in the previous section, in particular independent. The cdf for *R* can be calculated as follows (similarly for *w* and *z* interchanged). Clearly, for Σ_{w} = Σ_{z} = Σ we have that

We define such that . Note that the component-wise correlation coefficient *ϱ* between and reads

Hence, from (5) and (A1) we deduce that

We conclude that
with the linear combination of cdf evaluated at the origin. Hence, the cdf of the ratio (9) can be as well calculated with Davies algorithm. A consistency check can be performed as follows. In one-dimension, we have that *R* has to be Cauchy distributed. The corresponding cdf is given by . Evaluation for various *r* shows agreement with values calculated from (10) in the one-dimensional case. Note that for non-standardized *w*_{i} and *z*_{i} similar expressions can be derived, albeit in terms of the non-central *χ*^{2} distribution, *cf*., appendix A.

## VI. SIMULATIONS

In the following, several examples will be discussed, illustrating more specific aspects and applications of the theoretical material presented so far.

### A. Single element

It is useful to discuss the single element case in more detail. The distribution (8) simplifies for *N* = 1 to
which corresponds to the uncorrelated product normal distribution, *cf*., (5). The index *I* for *N* = 1 is a measure of coherence (or anti-coherence) between *z* and *w*. The significance threshold curve for a fixed desired *p*-value, say *p*_{I} = 10^{−7}, is illustrated in figure 1. The curve corresponding to a given *p*_{I} is unbounded. That is, 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 . 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 simply corresponds to that of the other *p*-value. In contrast, for the product-normal 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, one should see Fisher’s method as being additive in the evidence, while the product-normal based method as being multiplicative.

### B. Two elements

The importance of correcting for the inter-dependence between the elements of *w* and *z* can be seen easily in the *N* = 2 case. Consider the covariance matrices,
such that *w* ∼ 𝒩(0, Σ) and *z* ∼ 𝒩(0, Σ), and with *r* varying. In figure 2 we show various significance threshold curves of *I* for varying *r*, as calculated from the distribution (8). Clearly, for increasing inter-element correlation, the significance threshold level rises. The magnitude of the effect increases with the desired level of significance.

### C. Normal draws

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 for each pair *I* defined in (6). A *p*-value is then obtained for each index value for the linear combination of *χ*^{2} distributions (8) (also referred to as weighted *χ*^{2}), and for the gamma-variance distribution (B4). 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. Resulting qq-plots for both cases are shown in figure 3.

We observe that the gamma-variance distribution (7) (with *ϱ* = 0) indeed becomes unsuitable for increasing element-wise correlation of the data sample elements. In detail, not correcting for the inter-sample correlation leads to more and more false positives with increasing correlation strength. In contrast, the weighted *χ*^{2} distribution (8) yields stable results in both the weakly and strongly correlated regime, as is evident in figure 3.

## VII. GENETIC COHERENCE

### A. Generalities

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 *β* reads
with *x*_{i} the *i*th 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 multi-variate model into (12), yields

#### 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 (13) is for a collection of SNPs with effect sizes estimated via independent regressions.

#### b. Effect sizes as random variables

We can also take the effects to be random variables themselves. Let us assume that independently *α*_{i} ∼ 𝒩(0, *h*^{2}*/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, calculation of *L* usually requires an approximation, for instance, via a cutoff. Furthermore, the null model (14) requires an estimate of the heritability. Such an estimate can be obtained for via LD score regression [9].

### B. Product normal test

For this paper, it is only of importance that in both cases *a*. and *b*., the null model for *z* is a multi-variate Gaussian.

Therefore,
with *λ*_{i} the *i*th eigenvalue of the covariance matrix of the Gaussian. As discussed in detail in [2], a GWAS gene enrichment test can be performed via testing against (15). This effectively tests against the expected variance of SNPs significances in the gene.

What we propose here, is to use (8) of section IV, *i*.*e*.,
for *w* and *z* resulting from two different GWAS pheno-types, to test for co-significance of a gene for two GWAS. In more detail, a significance test can be performed either against the right tail of the null distribution (16) (*coherence*), or against the left tail (*anti-coherence*).

A clarifying remark is in order. 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 de-correlation, it is best to interpret this as testing independently each joint SNP in the gene for a coherent deviation from the null expectation. Therefore we refer to this test as testing for genetic coherence, or simply as cross-scoring. An example will be discussed below. For simplicity, in this work we only consider the fixed effect size model *a*. and assume that the correlation matrix Σ obtained from an external reference panel is a good approximation for both GWAS populations.

#### a. 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} (requires regularization of the estimated Σ), 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, gives an indication of 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 to test at least one GWAS via (17) to determine the base direction of the aggregated gene effect. Furthermore, at least one GWAS needs to carry sufficiently oriented signal in the gene such that (17) can succeed. We will also refer to testing for deviations from (18) as D-test.

### C. Ratio test

Note that the ratio
with and i.i.d. 𝒩(0, 1), and *λ*_{i} the *i*th 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 in section V, we can test for a significant deviation from the null expectation of no relation. Note that in general *R* is not invariant under interchange of response and explanatory variable, and may be used under certain conditions to make inference about the causal direction, *cf*., multi-instrument Mendelian randomization, in particular [21]. The main idea can be readily sketched.

Note first that the ratio test tries to detect genes that are maximal (anti)-coherent over the SNP effects between two GWAS, but at the same time carry minimal variance in one of the GWAS, as is clear from (19). The purpose of this, at first sight somewhat non-intuitive, statistical test becomes more clear if one thinks in terms of one GWAS trait as being the exposure and the other as outcome. If the ratio tested gene does carry minimal SNP variance only over the outcome, but on the same time SNP coherence between the exposure and outcome, the causal direction from exposure to outcome is implied, if the exposure is confirmed to be associated to the gene via *p*_{V} obtainable from (15), and confounding factors can be excluded. The setup is illustrated in figure 4.

### D. Pathway enrichment

The gene scores resulting from the above coherence or ratio test can be utilized instead of the usual gene scores resulting from (15) 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]. In detail, genes of a pathway which are in close proximity 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.

### E. SNP normalisation

As we have seen in section VI, the product normal combines evidence for coherent association in a multiplicative manner. A potential challenge to the method we propose 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 level of significance. We propose two possible strategies to mitigate this.

One strategy is to introduce 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 dynamics is clear from figure 1. If we target a co-significance of *p* = 10^{−8} and accept to consider SNPs with a *p*-value of 0.05 or less in one GWAS to be sufficiently significant, we have to cutoff *p*-values around 10^{−16}. 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^{−16}) 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 (16) is then computed with *w* and *z* according to the inverse *χ*^{2} cdf of the respective *p*-values. Since for GWAS *N* ≃ 10^{6} the most significant transformed *p*-value is ∼ 10^{−6}, according to figure 1 the other *p*-value has to be smaller than 10^{−3} − 10^{−4} 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 we will use it in the following example, and we also prefer to apply this transformation for the ratio test in order to compensate for different signal strengths between the GWAS.

## VIII. APPLICATION TO COVID-19 GWAS

### A. Coherence test

To demonstrate the usefulness of our methods in the context of real and topical GWAS data we consider the recent meta-GWAS on very severe respiratory confirmed COVID-19 [22, 23] as main phenotype and co-analyse it with several other related traits. Specifically, we use the summary statistics resulting from European subjects excluding those from UK Biobank [22] (A2 ALL eur leave ukbb 23andme, release 7. Jan. 2021), 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 against the null model (15) with fixed effect sizes. The Manhattan plot for the resulting gene *p*-values is shown in figure 5.

We cross-scored this GWAS against a panel of GWAS on medication within the UK Biobank [24]. These GWAS on medication take prescription (self-reported intake) of 23 common types of medications as traits. Hence, in cross-scoring against the severe COVID-19 GWAS we hope to uncover whether there is a shared genetic architecture between severity and pre-disposition to being prescribed specific medications.

All calculations have been performed with the python package *PascalX* [25], which incorporates the methods described in section VII. Default settings and the European subpopulation of the 1000 Genome Project [26] as reference panel to estimate the SNP-SNP correlations Σ were used. We only considered protein coding genes and SNPs within a gene window extending 50*kb* beyond the transcription start and end position for the cross-scoring. We transformed the raw GWAS *p*-values via *joint* rank transformation, for the reason discussed in section VII E. The directions of association are extracted from the signs of the raw effect sizes (*β*). All summary statistics data we made use of in this work have been made publicly available by the respective authors of the cited corresponding publications. Pathway enrichment has been tested against the gene sets included in MSigDB version 7.1 [27, 28].

#### a. Coherence

The cross-scoring results for the coherence test are shown in figure 6. Cross scoring of GWAS effects from prescription of medication group M05B (drugs affecting bone structure and mineralization) reveals joint signals with that for very severe COVID-19. We show the Manhattan plot resulting from the null model (13) and the index (16) for the medication group M05B in figure 7. 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 strong significance for the chemokine receptor genes *CCR1, CCR3* and the gene *LZTFL1* in the region chr3p21. Note that the two chemokine receptors are Bonferroni significant under the number of genes and drugs tested. However, *LZTFL1* slightly missed the threshold. We tested the orientation of the aggregated associations of these genes via (17) over the COVID-19 GWAS and find a positive direction with *p*_{D} ≃ 1.6 × 10^{−4}, 7.5 × 10^{−3} and 0.03, respectively.

For illustration, we show the spectrum of SNPs considered in the *CCR3* region and their SNP-SNP correlation in figure 8. One can see that there is a large block of SNPs in high LD, encompassing some of its 5’UTR and most of its gene body, all having positive associations with both COVID-19 and M05B drug prescription. Furthermore, a somewhat weaker signal of coherence can be seen for SNPs in the 3’UTR, some of which share negative associations. Using the method described in section IV, we calculate a *p*-value of *p* ≃ 1.50 × 10^{−8} for the significance of the coherence.

Note that the chemokine receptor of type 1 (*CCR1*) is involved in regulation of bone mineralization and immune/inflammatory response. In particular, chemokines and their receptors are critical for recruitment of effector immune cells to the location of inflammation. Mouse studies suggest that this gene plays a role in protection from inflammatory response and host defence [29, 30]. The chemokine receptor *CCR3* is of importance for regulation of eosinophils, a leukocyte involved in many inflammatory pathologies [31]. In particular, mouse models suggest a complex role of *CCR3* in allergic diseases [32]. In general, chemokines are suspected to be a direct cause of acute respiratory disease syndrome, a major cause of death in severe COVID-19. For a review of chemokines and their receptors in the COVID-19 context we refer to [33].

The SNP spectrum in the *LZTFL1* region is shown in figure 9. From the SNP-SNP correlation matrix on the right one can see that the region including the *LZTFL1* gene contains at least three sizable LD blocks. The *p*-value for the significance of the covariance is calculated to be given by *p* ≃ 1.19×10^{−7}. The observed co-significance of *LZTFL1* for the two traits is in line with the fact that the gene *LZTFL1* modulates T-cell activation and enhances IL-5 production [34]. In particular, mouse models suggest that expression of IL-5 alters bone metabolism [35].

Direction of association and coherence suggests that genetic predispositions leading to M05B prescription may carry a higher risk for severe COVID-19, with shared functional pathways related to the genes discussed above. We list the most significantly enriched pathways for (anti)-coherent M05B and severe COVID-19 genes detected via the pathway enrichment test described in section VII D in table I. Note that we detect several Bonferroni significant pathways, and that most of the leading coherent pathways appear to be of interest in the COVID-19 context, as immune system related. It is also noteworthy that in looking at common pathways between M05B and COVID-19 we increase the signal strength, as only the lead pathway (*Roeth tert targets dn*) is also detected to be Bonferroni significant under the COVID-19 GWAS alone, *cf*., table III in appendix C. The significant pathway of Roeth et al. [36] corresponds to genes that are significantly down regulated in T lymphocytes that overexpress human telomerase reverse transcriptase (hTERT). The also significant pathways of Kurozumi et al. [37] correspond to inflammatory cytokines and their receptors modulated in brain tumors after treatment with an oncocytic virus.

#### b. Anti-coherence

We perform a similar test between COVID-19 severity and drug classes as above for anti-coherence. The corresponding results for all medication classes are shown in figure 6. Despite the fact that no drug received a Bonferroni significant hit, we note that the drug classes H03A (Thyroid preparations), C10AA (HMG CoA reductase inhibitors), L04 (Immuno-suppressants) and M01A (Anti-inflammatory and anti-rheumatic products) have gene hits with a *p*-value < 1 × 10^{−5}. Interestingly, all of these drugs find applications in auto-immune diseases and allergies. The genes with a *p*-value below 1 × 10^{−5} for H03A are *HLA-DQB1*, for C10AA *SMARCA4* and *BCAT2*, for L04 *LZTFL1, TRIM10, TRIM15* and *TRIM26*, and for M01A *TRIM10* and *TRIM31*. Note that TRIM proteins are associated with innate immunity, and are in particular involved in pathogen-recognition and host defence [38, 39]. For illustration of the anti-coherence case, we plot the SNP spectrum for *TRIM10* under L04 in figure 10. We tested the above genes for direction of aggregated effect sizes under the COVID-19 GWAS via (17) and found that the detected *TRIM* genes tend to be localized in the left tail. In particular, *p*_{D} ≃ 0.04 for *TRIM10* and *TRIM15* under L04 (right tail). Hence, a protective effect is suggested for these genes, and therefore conditions leading to prescription of these medications may imply a lower genetic risk for severe COVID-19. The pathway enrichment test for the anti-coherent L04 cross-scored severe COVID-19 genes show Bonferroni significant enrichment in interferon gamma related pathways, see table I. Note that *TRIM* proteins are expressed in response to Interferons, *cf*., [38], and therefore the detected pathways are consistent with the above observed enrichment of TRIM genes. Remarkably, the individually non-Bonferroni significant genes aggregate to Bonferroni significant path-ways. While two of the detected pathways are also Bon-ferroni significant under L04 alone, one pathway (*GO interferon gamma mediated signaling pathway*) is only Bon-ferroni significant under anti-coherence with COVID-19. Note that we detected for L04 individually 35 Bonferroni significant pathways, *cf*., table III. Therefore, the coherence test singles out a small subset of potentially shared pathways with COVID-19.

#### c. M05B related traits

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 [40], estrogene levels in men (estradiol and estrone) [41], calcium concentration [42], vitamin D (25OHD) concentration [43] and rheumatoid arthritis [44]. The inferred gene enrichment for coherence and anti-coherence with the COVID-19 GWAS is shown in figure 11. We found enrichment with gene *p*-values < 10^{−5} for vitamin D and calcium. In particular, in the anti-coherent case the most significant genes for vitamin D are *OAS3, OAS2, OAS1, FYCO1, CXCR6* and *LZTFL1*. We show the corresponding Manhattan plot in figure 12. The D-test shows that all these genes are in the right tail (*p*_{D} < 3 × 10^{−2}), in particular *OAS3* with *p*_{D} ≃ 7.21 × 10^{−5}, under the COVID-19 GWAS, and therefore suggests an increase in the risk for severe COVID-19 in the presence of effect alleles in these genes.

The *OAS* family are essential proteins involved in the innate immune response to viral infection. They are involved in viral RNA degradation and the inhibition of viral replication [45]. It has been observed that vitamin D can increase the expression of the *OAS* genes [46].

FYCO1 plays a role in microtubule plus end-directed transport of autophagic vesicles [47]. It has been demonstrated that SARS-CoV-2 inhibits autophagy activity [48]. A regulatory role of vitamin D on autophagy at different steps, including induction, nucleation and degradation has been suggested, *cf*., [49].

The gene *CXCR6* (*p* ≃ 7.96 × 10^{−6}) is expressed by subsets of T_{H}1 cells, but not by T_{H}2 cells, and may be important in trafficking of effector T cells that mediate type-1 inflammation [50]. It has been shown that the vitamin D analog TX527 promotes the surface expression of *CXCR6* on T-cells and inhibiting effector T cell reactivity while inducing regulatory T cell characteristics, promoting migration to sites of inflammation [51]. It is therefore suggestive that vitamin D concentration influences differences in immune system reaction, *i*.*e*., between type-1 (inflammatory) or type-2 (anti-inflammatory), and thereby impacting the severity of COVID-19. Indeed, possible links between severity of COVID-19 and vitamin D are actively discussed in the current literature. See for instance [52–54] and references therein.

For calcium in the anti-coherent case, the top genes are *HGFAC* (*p* ≃ 4.15 × 10^{−7}) and *DOK7* (*p* ≃ 1.61 × 10^{−6}). Note that *HGFAC* is Bonferroni significant under the seven traits tested. Both genes sit in the right tail under the D-test (17) applied to the calcium GWAS (*p*_{D} ≃ 5.87 × 10^{−4}, 3.5 × 10^{−5}). Therefore, the aggregated variants in these genes imply that a predisposition for high calcium concentration may implicate a reduced risk for severe COVID-19, *i*.*e*., a protective effect. Note that in general it is known that viruses appropriate or interrupt *Ca*^{2+} signaling pathways and dependent processes, *cf*., [55].

The gene *HGFAC* plays a role in converting hepatocyte growth factor (HGF) to its active form. In particular, binding of HGF causes the up-regulation of *CXCR3*, which is primarily activated on T lymphocytes and NK cells. *CXCR3* is preferentially expressed on T_{H}1 cells, while *CCR3* on T_{H}2 cells [56]. In detail, *CXCR3* binds the chemokine receptor *CCR3* and prevents an activation of T_{H}2-lymphocytes. Thereby, a towards T_{H}1 biased in-flammation immune reaction is triggered [57]. Note that *CXCR3* is able to increase intracellular *Ca*^{2+} levels [58].

The gene *DOK7* is of importance for neuromuscular synaptogenesis [59]. It activates *MuSK*, which is essential for maintenance of the neuromuscular junction as it is involved in concentrating *AChR* in the muscle membrane at the neuromuscular junction. The latter protein is critical for signaling between nerve and muscle cells, a necessity for movement, and is influenced by intra-cellular calcium [60]. Hence, it may be possible that *DOK7* could be involved in a genetic explanation of the muscle weakness impacting some of the severe COVID-19 patients [61].

Note that for the vitamin D and calcium GWAS, we observe cross enrichment in, both, the coherent and anti-coherent case with the severe COVID-19 GWAS, *cf*., figure 11. For example, for vitamin D the leading gene hit is *AGAP2* (*p* ≃ 2.84 × 10^{−6}). *AGAP2* modulates the transforming growth factor beta-1 (*TGF-β1*), one of the most potent pro-fibrotic cytokine known to date, currently accepted as the principal mediator of the fibrotic response in liver, lung, and kidney [62]. The D-test under the COVID-19 GWAS is however inconclusive (left-tail *p*_{D} ≃ 0.15).

Another implied coherent gene locus of potential interest is *OS9* (*p* ≃ 4.2 × 10^{−6}). The corresponding protein binds to the hypoxia-inducible factor 1 (*HIF-1*). *HIF-1* is a key regulator of the hypoxic response [63, 64]. In particular, regulation of *HIF-1* interpolates between regeneration and scaring of injured tissue [65]. It is known that severe COVID-19 may lead to lung tissue fibrosis [66, 67]. The D-test applied to the COVID-19 GWAS shows that *OS9* is located in the left tail (*p*_{D} ≃ 0.03), and therefore a protective property of effect alleles in this gene are suggested.

### B. Can coherence make a difference in identifying co-significant genes?

We next wanted to explore whether gene-wise significance based on the product normal statistic that tests for consistently coherent or incoherent SNP-wise associations within the gene-window gives different results than when first computing gene-wise significance for each trait (via the usual *χ*^{2}-test, *cf*., eq. (15)), and then combining the corresponding scores, as if there was only a *single* genetic element affecting both traits (*cf*., fig. 1 for the corresponding combined product normal *p*-values). As before, we use jointly qq-normalized GWAS *p*-values to avoid the risk of unintended uplift.

In general, we would expect that the simple approach described in this section would lead to more false positives, but would also lead to false negatives. The danger for false positives is clear, as the simple approach is blind to incoherent directions of association between the SNPs of the two GWAS in the gene window under consideration. However, also false negatives may occur: In the proposed novel approach using relation (16), several independent weak signals may combine to a stronger signal, if the directions of association are consistent over the gene window. Such signals are invisible in the simple approach of comparing aggregated gene *p*-values.

Figure 13 (top plot) shows the scatter plot for genewise log-transformed *p*-values for the severe COVID-19 GWAS against the medication class M05B GWAS, and significance thresholds for the product-normal. We observe that in this example all of the genes which are significant (*p <* 10^{−5}) based on the SNP-wise product normal coherence test, are also significant under the simple single-element test. Yet, interestingly, two genes that pass the simple test, i.e. *CXCR6* and *CCR2*, are not significant under the coherence test. While these genes obtain significant association for both traits, their SNP-wise effects appear to be partly coherent and partly incoherent, explaining why the coherence test is not as significant, *cf*., figure 14, which compares the SNP spectrum of *CCR2* against *CCR5*.

Nevertheless the corresponding coherence *p*-values are just above the 10^{−5} threshold. Since these genes are also cytokine receptors, it seems likely that their joint association with COVID-19 and M05B points to a common mechanism relevant for these traits that is modulated by variants in these genes.

A similar picture emerges for the genes with significant cross-scores of severe COVID-19 against vitamin D concentration, see figure 13 (middle plot): Their simple single-element statistic is also significant and conversely only very few genes whose cross-score is insignificant, *i*.*e*., with *p <* 10^{−5}, have a significant single-element statistic, and if so only by a small margin.

A somewhat different pictures emerges from the co-analysis of COVID-19 and calcium GWAS signals: Here a number of genes with a significant simple single-element statistic do not achieve a significant cross-score, and there are several genes with the reverse signature of a significant cross-score, but insignificant single-element statistic. The latter include *LRPAP1* and *DOK7*, which both do not exhibit a strong association with COVID-19, yet the available annotation (*LRPAP1* codes for a LDL-receptorrelated protein and *DOK7* has been assigned to the GO category “lipid binding”) makes a role for these genes in severe COVID-19 plausible, as dislipidemia has been associated with the severity of COVID-19 [68].

### C. Ratio test

Let us also briefly discuss an application of the ratio based causality test introduced in sections V and VII C. We ratio tested the COVID-19 GWAS against the vitamin D and calcium concentration GWAS in the coherent case, making use of eq. (11) and (19), and found that Vitamin D carries two interesting hits which suggest causal pathways from genetic predisposition for Vitamin D concentration to severity of COVID-19. The ratio test Manhattan plot is shown in figure 15. The leading gene is *KLC1* (*p*_{R} ≃ 9.18 × 10^{−7}). It is known that Kinesin-1 uncoats viral DNA and compromises nuclear pore complex integrity, allowing viral genomes nuclear access to promote infection [69]. Note that this gene also features in the COVID-19 virus-host protein interactions, belonging to the functional group of viral trafficking, as recently discussed in [70]. The gene carries under the outcome (COVID-19) a *p*-value of *p*_{V} ≃ 0.95 and under the exposure (vitamin D) a *p*-value of *p*_{V} ≃ 3.23 × 10^{−9}. Under calcium we have *p*_{V} ≃ 0.25, hence calcium can be excluded as potential confounder. The gene is in the right tail (*p*_{D} ≃ 0.07) under the vitamin D GWAS. We conclude that the gene may mediate a causal relation.

Furthermore, besides the paralog gene *AL139300*.*1* of *KLC1*, the gene *ZFYVE21* (*p*_{R} ≃ 3.9 × 10^{−6}) stands out. Under the COVID-19 GWAS (outcome) the gene carries a *p*-value of *p*_{V} ≃ 0.99 while under the vitamin D GWAS (exposure) a *p*-value of *p*_{V} ≃ 6.5 × 10^{−4}. This implies a causal flow from genetic predisposition for vitamin D concentration to severity of COVID-19. We confirmed that calcium concentration is not a potential confounder for *ZFYVE21* (*p*_{V} ≃ 0.26). The D-test shows that the gene is in the right tail (*p*_{D} ≃ 0.03) and hence is associated with high vitamin D concentration. A possible explanation of this hit may go as follows. *ZFYVE21* regulates microtubule-induced PTK2/FAK1 dephosphorylation, which is important for integrin beta-1/ITGB1 cell surface expression. It has been discussed before in the literature that integrins in host cells may play the role of alternative receptors to ACE2 for SARS-CoV-2 [71, 72]. Hence, it is tempting to speculate that genetic predisposition for vitamin D concentration may also influence expression of the integrin receptors, and thereby the outcome of COVID-19 via an increased risk of cellular infection.

Another gene we observe to be of relevance in the coherent case is the Ring Finger Protein 217 (*RNF217*) with *p*_{R} ≃ 1.83 × 10^{−6} and calcium as exposure. We have *p*_{V} ≃ 2.15 × 10^{−5} under calcium and *p*_{V} ≃ 0.82 under COVID-19. Vitamin D has for this gene *p*_{V} ≃ 0.26, hence can be excluded as potential confounder. *RNF217* is a member of the E3 ubiquitin protein ligase family. Ubiquitination of proteins is a post-translational modification process with different cellular functions, including antiviral functions and virus replication [73]. A potential COVID-19 therapeutic pathway based on E3 ubiquitin ligases has been recently proposed in [74]. As the gene shows an illustrative coherence pattern (multiple LD blocks), we show the corresponding SNP correlation plot in figure 16.

We also tested the anti-coherent case. The gene *HOXC4* with *p*_{R} ≃ 4.5 × 10^{−6} is detected for M05B. Individual scored *p*-values for the gene read *p*_{V} ≃ 2.04 × 10^{−5} under M05B and *p*_{V} ≃ 0.75 under COVID-19. The D-test under M05B shows that the gene signal is located in the right tail, but not significantly. Modulo potential confounders, a causal relation from predisposition to take M05B to a protective property towards severity of COVID-19 via *HOXC4* is suggested. We have *p*_{V} ≃ 0.27 under calcium and *p*_{V} ≃ 0.06 under vitamin D, hence, a potential confounding role of vitamin D is possible. We note that *HOXC4* has been discussed before in the COVID-19 context [75]. The gene is related to an enhanced antibody response under the regulation of estrogens. This matches the observation of anti-coherence with severity of COVID-19.

Interestingly, also for the inverse causal relation, *i*.*e*., a predisposition for severe COVID-19 implying a predisposition for the need to take M05B medication, a relevant gene is singled out. Namely, the gene *DPP9* with *p*_{R} ≃ 3.88 × 10^{−5}, and *p*_{V} ≃ 0.74 under M05B, respectively, *p*_{V} ≃ 2.57 × 10^{−5} under COVID-19. For this gene *p*_{V} for calcium and vitamin D is not significant. This gene has been implicated before to be involved in the genetic mechanisms for critical illness due to COVID-19 [76]. Note that we also observe *DPP9* (*p*_{R} ≃ 5.82 × 10^{−7}) as top hit for calcium concentration as outcome, *cf*., the corresponding Manhattan plot shown in figure 17. Under both exposures the gene signal tends to be more in the right tail of the D-test, but is too weak for calling. Among the genes we detect in this ratio test, the gene *IFNAR2* (*p*_{R} ≃ 9.70 × 10^{−6}) appears to be of particular relevance, as it is a known possible drug target against COVID-19 [76].

For vitamin D as exposure, we find in the anti-coherent case that the genes *MED23* (*p*_{R} ≃ 5.22 × 10^{−5} and *GATA4* (*p*_{R} ≃ 7.88 × 10^{−5}) are the leading hits with a *p*-value < 10^{−4}. Both genes are not significant under calcium. The D-test shows that *MED23* is in the right tail under the exposure (*p*_{D} ≃ 1.22 × 10^{−3}), whereas the test is inconclusive for *GATA4*. The Mediator subunit gene *MED23* has emerged before in SARS-CoV-2 screens, with a critical role of this complex during infection and death suggested [77]. *GATA4* has been observed previously in the context of the top significant biological processes likely associated with respiratory failure in COVID-19 patients [78].

Finally, note that the genes receiving very small *p*-values for our anti-coherent ratio test with calcium as exposure, i.e. *HGFAC* (*p*_{R} ≃ 4.2 × 10^{−5}) and *DOK7* (*p*_{R} ≃ 1.27 × 10^{−5}), were already detected above via the coherence test. Therefore, it is suggested that the corresponding pathways could be causal. However, for both genes we have that *p*_{V} ≃ 0.01 under vitamin D, and therefore vitamin D could be a potential confounder. In addition, we observe the leading gene *PPP2R3A* (*p*_{R} ≃ 1.0 × 10^{−5}) with *p*_{V} ≃ 1.9 × 10^{−4} under calcium, but with a potential confounding role of vitamin D (*p*_{V} ≃ 0.03). Note that *PPP2R3A* interacts with *CDC6* [79], which is up-regulated at the early stages of human coronavirus 229E infection [80], and initiates assembly of a pre-replication complex [81].

We also tested all cases against pathway enrichment as outlined in section VII D. The resulting enriched path-ways with *p <* 10^{−4} are listed in table IV of appendix C. Note that most of the detected enriched pathways are immune system related, and may be of interest for further investigations in the COVID-19 context.

## IX. SUMMARY AND CONCLUSION

The work presented here relies on the novel mathematical insight that the null distribution of the normal product of two effect size vectors with a common known covariance structure, can be expressed in terms of a weighted *χ*^{2} difference distribution (*cf*., eq. (8)). Moreover, we showed that the ratio between this product and the square of one of the vectors, which allows for testing causal relationships, also relates to this distribution. In both cases the corresponding cdf and tail probabilities can be computed efficiently with Davies’ algorithm.

This insight has application for GWAS, because it enables testing for *coherent* SNP-wise effects within a gene window for two different traits, *cf*., section VII. Importantly, our test is different from an additive test looking merely for co-significance. Indeed testing whether the signs of the effects sizes tend to be the same (or opposite) within the window, has potentially more power to identify genes that modulate both traits. Our test is statistically rigorous, properly taking into account the SNPs’ correlation structure (i.e. LD), while still being able to be computed rapidly and accurately without any approximations.

We applied our method to identify genes with a role both in severe COVID-19 and medical conditions leading to the prescription of any of 23 medication classes. Our analysis revealed a particular strong signal of coherence between COVID-19 and M05B medications, with chemokine receptor genes *CCR1* and *CCR3* as lead hits. We then searched for coherence signals between COVID-19 and additional GWAS traits, such as vitamin D and calcium concentrations, known to be related to diseases treated with M05B drugs. This analysis implicated genes related to differentiation between type-1 and type-2 immune response. A possible explanation could be that many patients being prescribed class M05B medication suffer from osteoporosis. Vitamin D stimulates the absorption of calcium and therefore is often prescribed to these patients in order to increase their bone mineral density [82]. While data clearly support the function of vitamin D in bone growth and maintenance, evidence supporting the role of vitamin D in other health and disease processes, in particular in acute respiratory tract infections [83], often observed for patients with severe COVID-19, is less clear cut: Vitamin D is thought to reduce the risk of infection, mainly due to factors involving physical barriers, cellular natural immunity and adaptive immunity [84]. Furthermore, a lack of plasma vitamin D level has been associated with the risk of infection [85]. Yet, direct conclusive evidence for its proposed protective function specific to severe COVID-19 is still lacking [86]. In this work we found potential genetic hints that link the severity of COVID-19 to vitamin D concentration through specific genes operating in immune response pathways, including those coding for chemokine receptors. In addition, we detected a potential causal link from genetic predisposition for vitamin D deficiency to severity of COVID-19, mediated by the *ZFYVE21* gene, which is related to integrin beta-1/ITGB1 cell surface expression and thereby potentially impacting disease outcome by influencing alternative receptors to ACE2 for host cell entry.

Our causal analysis using the ratio test (*cf*., eq. (19)) suggests that the genes *HOXC4* and *DPP9*, which have already been associated with severe COVID-19 [75, 76], may contribute to the protective role of conditions leading to M05B medication prescription, and that the genes *MED23* and *GATA4* may play a role in the causal link between vitamin D deficiencies and severe COVID-19. Note that the novel aspect of our method is its capacity to identify candidate genes mediating the causal effects, but we cannot exclude potential confounders in this causal analysis.

We also found that genetic variants in genes related to both the adaptive (HLA) and innate immune system (TRIM genes) that have a higher frequency in subjects being prescribed medications indicated for specific auto-immune disorders tend to reduce the risk of severe COVID-19. Indeed, it is well known that auto-immune disorders are more common in females [87] who also have a smaller risk of severe COVID-19 in comparison to men [88]. While it is of course reasonable to expect that subjects with an increased risk for auto-immune disorders will tend to fight off infections more efficiently, the added value of our analysis is to pinpoint specific genes that are potentially involved in mediating this effect, and which may hint towards a protective pathway and/or therapeutic targets against severe COVID-19.

Our novel localized gene-centred approach was able to single out several plausible candidate causal genes, discussed in part already elsewhere in the literature in the COVID-19 context, and worthy of further investigations. Note that recent work using standard GWAS methods could not detect genetic evidence linking vitamin D to severity of COVID-19 [89]. This illustrates the power of the methods developed in this work as a discovery engine. We therefore believe that these methods are likely to be of utility for other studies trying to identify the genetic player mediating pleiotropic effects. An expected increase in the power of severe COVID GWAS and other relevant traits is likely to refine the picture starting to emerge in our analysis even further.

We consider it natural that the same genes, or even their domains, may exhibit effects with opposite signs for GWAS of different traits. Similarly pathways, or subsets thereof, may include sets of genes whose effects either coherently increase or reduce the effect on (or risk of) the trait. Coherence may occur at different levels of granularity (domain, gene, gene-pair, subunit or entire pathway), and these levels depend on which GWAS trait pair is co-analysed, specifically the complicated interplay between the common genetic components in each trait. Full genome based methods like genetic correlation or usual applications of Mendelian randomization can not resolve this fine structure, but rather yield only an aggregated trend for the traits. We therefore like to stress that statements made in this work regarding trait risk implications are not made for the traits in general, but rather for the respective gene under discussion. For a given pair of traits, the implication may differ between different genes.

In this work we assumed that there is no sample overlap between the GWAS under consideration and that the GWAS under consideration either have similar sample sizes, or that effects of sample sizes differences can be compensated by qq-normalization. Further work is needed to study the potential impact of sample overlap in more detail and to develop potential corrections along the lines of [90, 91]. Also our method would profit from better ways to deal with GWAS pairs whose traits exhibit very different effect size distributions or sample sizes.

Finally, we can envisage other, more general, applications of the methods discussed in this paper. For instance, our approach could be used to correct for the auto-correlation structures in estimating the significance of correlation between time-series. Hence, the technical results of this paper may be as well of interest for other domains.

## Data Availability

All data is publicly available from third-parties in terms of summary statistics and references given.

## Acknowledgments

We like to thank A. L. Button, A. Brümmer, Z. Kutalik and S. O. Vela for valuable comments on an earlier draft of the manuscript.

## Appendix A: General product-normal

Let us briefly consider as well the distribution for *xy* with and , *i*.*e*., the product-normal for non-standardized Gaussian random variables. The moment generating function reads [19]

Note that for, either, *µ*_{x} = 0 or *µ*_{y} = 0, the factorization of the main text still holds. Defining and , the exponent of the exponential in can be split for the general non-correlated case (*ϱ* = 0) into

We deduce that we can still factorize the moment generating function, such that
with the non-central *χ*^{2} distribution with one degree of freedom. Hence, the general product-normal can be expressed as a linear combination of non-central distributions, for *ϱ* = 0.

## Appendix B: Probability density functions

### a. Product-Normal

Making use of the relation (4), the pdf, denoted as *f*, of the product-normal distribution can be calculated analytically via convolution (*cf*., [92])

Completing the square and invoking hyperbolic substitution, we arrive at
with *K*_{0} the modified Bessel function of second kind at zero order. The result above for *f* is in agreement with the previous derivations of [93, 94]. Note that the analytic calculation of the corresponding cdf requires the solution of an integral of the type . We are not aware of a known closed form solution.

For illustration, we plot the pdf for *ϱ* = 1*/*2 together with the corresponding histogram sampled from (4) in figure 18. We also show the cdf obtained via numerical integration of (B1). We verified that the numerical integration matches the results obtained via Davies algorithm for the cdf calculation for various *ϱ, cf*., table II.

Note that since the pdf of the non-central *χ*^{2} distribution includes a Bessel function, analytic calculation of the pdf of *xy* in the more general case of section A is more complicated than in (B1), and will not be discussed here.

### b. Variance-Gamma

Consider a random variable *X* distributed according to

The corresponding pdf can be calculated similar as above via convolution. We infer

Completing the square and using as before hyperbolic substitution, we arrive at
with *K*_{n} a modified Bessel function of second kind at order *n*. Using the integral , we easily verify that . Hence, due to symmetry the pdf is well normalized. However, we are not aware of closed form solutions for Bessel function integrals of the type , which are needed to provide a closed form expression for the cdf.

The distribution given by (B3) occurred before in the finance domain as a special case of the *Variance-Gamma* distribution [95]. It can be traced back further to the distribution of the bivariate correlation. In detail, the gamma-variance corresponds to the off-diagonal marginal of a two-dimensional Wishart distribution, which models the covariance matrix [96]. However, what is new, to the best of our knowledge, is the expression in terms of the difference distribution in equation (B2).

For *h* = *n* and *g* = 1, we have

Hence, for *n* = 1 we obtain the product-normal distribution with *ϱ* = 0 as discussed in the previous section. For general *n* we can view *X*(*n*, 1) as the distribution of a sum of *n* independently distributed product-normal random variables. In particular, we can make use of Davies algorithm to calculate the cdf for *X*(*n*, 1) exactly at a desired precision.

## Appendix C: Enriched pathways

The enriched pathways (*p <* 10^{−5}) for the individual COVID-19, M05B and L04 GWAS, calculated as described in section VII D, are given in table III. Note that we detect for the COVID-19 trait only one Bonferroni significant pathway (*Roeth tert targets dn*).

We also tested for pathway enrichment of the genes under the ratio test for the combinations of traits investigated in section VIII C. We detected no Bonferroni significant hit. However, several of the lead pathways may be of potential interest in the COVID-19 context, as immune system related. Therefore, the top pathway scores (*p <* 10^{−4}) are listed in table IV.

## Footnotes

Pathway analysis added; Discussion extended; Minor modifications and corrections

## References

- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].↵
- [10].↵
- [11].↵
- [12].
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].↵
- [51].↵
- [52].↵
- [53].
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
- [59].↵
- [60].↵
- [61].↵
- [62].↵
- [63].↵
- [64].↵
- [65].↵
- [66].↵
- [67].↵
- [68].↵
- [69].↵
- [70].↵
- [71].↵
- [72].↵
- [73].↵
- [74].↵
- [75].↵
- [76].↵
- [77].↵
- [78].↵
- [79].↵
- [80].↵
- [81].↵
- [82].↵
- [83].↵
- [84].↵
- [85].↵
- [86].↵
- [87].↵
- [88].↵
- [89].↵
- [90].↵
- [91].↵
- [92].↵
- [93].↵
- [94].↵
- [95].↵
- [96].↵