Introduction to generalized linear models in genetics.

In this subsection, we introduce the theory of GLMs and their application to genetic data. We then describe how a so-called parameterization can be used to incorporate genetic assumptions about a single variant. Finally, we describe how a multi-variant parameterization can be constructed from several single variant parameterizations.

In this work, we use GLMs to describe the relation between predictor variables, which here typically are genotypes, and an outcome variable, that is, the phenotype. For each individual i, Y_i is a measured phenotype and x_i is a vector of predictor variables. The observed phenotype y_i is modeled by its expected value E[Y_i ∣ X_i] = μ_i along with a distribution from the exponential family that captures the stochastic variation around μ_i. We write y_i ∼ f(μ_i) and refer to f as the dispersion distribution. The expected value μ_i is in turn related to the linear predictor ψ(x_i)β by g(μ_i) = ψ(x_i)β, where g is the link function, ψ is an encoding of the predictor variables, and β is a vector of parameters. A combination of ψ and β is called a parameterization. The parameters are generally estimated according to the maximum likelihood principle by applying the iteratively reweighted least squares algorithm. The commonly used linear regression is a special case of GLMs, obtained by using the identity link function (g(μ_i) = μ_i) and the Normal dispersion distribution. Moreover, logistic regression is a GLM with the logit link function () and a Binomial dispersion distribution.

As our focus here is interactions, we restrict the predictor variables x_i to be genotypes. Because the set of possible genotypes for a set of variants are finite, and in practice often fewer than the number of individuals, many individuals will have identical x_i. It is therefore convenient to define a model that is indexed by an enumeration of the genotypes h instead of by individuals i. In this notation, p_h is the design vector of the genotype that corresponds to h, and, for an individual i with genotype h, μ_h, the expected value of y_i, is modeled by g(μ_h) = p_hβ. Here h can be the genotype of either a single variant or multiple variants. Furthermore, we can define a design matrix P, so that p_h are the rows of P; the complete relationship between the mean levels of the genotypes and the parameterization can then expressed g(μ) = Pβ.

In general, we will call any model where the number of parameters equals the number of unique design vectors a saturated model, otherwise it is called unsaturated. Additionally, we will call a model full if the number of unique design vectors equals the number of possible genotypes, otherwise it is called constrained. A design matrix that recodes genotypes into dominant and recessive states is an example of a constrained model.

Single-variant models. We first introduce a notation for the parameterization for a single variant and, in the next section, we extend this notation to multiple variants. We let β = (α, β₁, β₂)^T, where α is referred to as the reference level, and β₁ and β₂ are deviations from the reference (the ^T indicates the transpose, i.e., β is a column vector). For example, if p_h = (1, 2, 0) then g(μ_h) = α + 2β₁. For the single variant case, let A be the most frequent haploid genotype in the population considered, and a be the less frequent. We have three possible diploid genotypes {AA, Aa, aa} that are enumerated by h ∈ {0, 1, 2}. We now describe the parameterizations corresponding to the most common single variant genetic models.

• Saturated full parameterizations
  • G—Genotypic. Let β₁ be the mean difference in phenotype between the reference and the heterozygote, and β₂ be the mean difference in phenotype between the reference and the minor homozygote. Then the genotypic can be expressed as follows:
  • AD—Additive/deviation. Letting β₁ denote the additive component, and β₂ denote the deviation from additivity, this can be expressed as follows:
• Saturated constrained parameterizations
  • R—Recessive. This model assumes that the effect allele has a recessive effect:
  • D—Dominant. This model assumes that the effect allele has a dominant effect:
  • H—Heterozygote. This model assumes an heterozygote advantage effect:
• Unsaturated full parameterizations
  • A—Additive. This model assumes that only additive effects are present:

We note that G alternatively could be perceived as either a Heterozygote/deviation, HD, or a Recessive/deviation, RD.

From single to multi-variant models. We will now show how a multi-variant parameterization with relevant interaction parameters can be obtained from a set of single variant parameterizations. This provides a intuitive way to construct higher order interaction parameterizations. Thereafter, we provide two examples of the construction of bi-variant models (see Fig 1).

Download:

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

Fig 1. Two examples of the construction of bi-variant interaction models as the Kronecker product (⊗) of two uni-variate models.

(a) The construction of the P matrix of the G × G model, (b) The construction of the P matrix of the AD × AD model, (c) the G × G model in bi-variant notation, (d) the AD × AD model in bi-variant notation. For both models, the resulting parameter vector β = (α, β₁, β₂, γ₁, δ₁₁, δ₁₂, γ₂, δ₂₁, δ₂₂). In c) and d) a ∈ {0, 1, 2} is the genotype for the first variant and b ∈ {0, 1, 2} is the genotype for the second variant, implicitly β₀ = γ₀ = δ_0* = δ_*0 = 0, and, lastly, I(x) is an indicator function taking the value 1 if x is true and 0 otherwise.

https://doi.org/10.1371/journal.pcbi.1005556.g001

A multi-variant parameterization can be constructed from single variant parameterizations using the Kronecker-product (⊗) [17]. Given a design matrix, P⁽ⁱ⁾, for each variant i ∈ {1…k}, in which each row parameterizes the corresponding genotype, then the multi-variant design matrix, P, for these variants is (1) The parameters β for the multi-variant model can similarly be constructed by taking the Kronecker-product of the parameter vectors (2) where, in each Kronecker operation β^(k) ⊗ β^(k−1) (i) the intercept has been replaced by 1 in all parameter vectors β^(⋅) on the right hand side of Eq (2), (ii) the first element of the resulting vector β on the left hand side of Eq (2) is the new intercept, and (iii) any product between two different parameters is replaced by a new parameter that represents interaction. The number of factors in each element represents the order of interaction and factors of order 1 are the main effects. The matrix P and the parameters β are now related to the mean values by

We will below mainly focus on bi-variant models resulting from applying Eqs (1) and (2) to two uni-variant models, say M_i and M_j, of those enumerated above in Section Introduction to generalized linear models in genetics. We denote the resulting bi-variant model M_i × M_j. We detail the construction of G × G and AD × AD in Fig 1. Notice that the Kronecker-product of two saturated single-variant models yield a saturated bi-variant model.

This concludes our description of the GLM components used to model genotype-phenotype association. Using this uniform GLM framework, we formulate, in the next section, efficient tests for interaction applicable to very general families of GLMs, allowing any number of variants and a wide variety of phenotypes and link functions. Moreover, in Section Relating different GLMs, we use this framework in a comparative analysis of major GLM families with respect to differences in parameterization, link function and dispersion distribution.

Fast estimation and testing of interaction in generalized linear models.

In this section, we first derive closed-form estimates of the parameters of a GLM and the corresponding covariance matrix for any saturated parameterization. We then build on this to design an computationally efficient test for interaction. A complete discussion and proofs can be found in S1 Text. We generalize these results to unsaturated models using a two-step estimation procedure. Finally, we show how the Wald tests can be used in a meta-analysis.

A general Wald test for saturated models. We first describe how to estimate the parameters. Let N be a diagonal matrix in which the diagonal element N_hh is the number of individuals with the genotype corresponding to h. Let t be a vector in which element t_h contains the sum of the phenotypic values of individuals with the genotype corresponding to h. Maximum likelihood estimates of β is then obtained by (3) where the inverse of P exists because P is a full rank square matrix when the parameterization is saturated (see S1 Text). The time complexity of this estimation is, assuming that matrix inversion can be performed in cubic time, where T is the number of variants in the model and n is the number of samples. The first term is for computing the sufficient statistics, the second for the matrix inversion of P, and the third for computing the product between the transformed sufficient statistics and P⁻¹. Here, P⁻¹ only needs to be computed once, and T is typically small (T = 2 in this article), so that the dominant time complexity is linear in the number of samples .

The covariance matrix of the saturated full parameterization can be similarly derived, and has the following simple form where is the Fisher information matrix. The inverse of the Fisher information matrix is a diagonal matrix with elements on the form where v_h = Var(y_h ∣ p_h, ϕ = 1), which can be calculated from the specific dispersion distribution that is used, is the derivative of the link function evaluated at (the estimate of μ_h), n_h is the number of samples with genotype h, and ϕ is the dispersion parameter of the GLM (see S1 Text).

The general test for interaction in association analysis evaluates whether the relevant interaction parameters are significantly different from zero. There are three asymptotically equivalent tests that can be used to test the interaction parameters: the score test, the Wald test, and the explicit likelihood ratio test (LRT). The score test requires estimation only under the null hypothesis, the Wald test requires estimation under the alternative hypothesis, and the standard likelihood ratio test under both. For saturated GLMs, Eq (3) provides an efficient parameter estimation under the alternative hypothesis, and, consequently, we base our test for interaction on the multivariate Wald test. Let denote the vector of estimated interaction parameters, and denote the sub-matrix of the estimated covariance matrix that corresponds to the interaction parameters. The Wald test statistic is a quadratic form in constructed using the estimated covariance matrix, and the test is defined as This Wald test requires computing the inverse of a typically small covariance matrix C_δ. Consequently, the total time complexity of estimation and testing is . Again, as T is small, the dominant time complexity is linear in the number of samples n. Our Wald test is applicable to any saturated GLM, and contains, as a special case, the previously described Wald test for logistic regression [21].

The score test and the LRT additionally requires estimation of β under the null model of no interaction, for which no closed form expression exists. Hence, in this case, we need to resort to the slower iteratively reweighted least squares (IRLS) algorithm. The average time complexity for estimation and testing based on the IRLS algorithm is where A(n) is the expected number of iterations until convergence. Since the time complexity for computing the actual test statistics is similar to that of the Wald test, these two tests will have a theoretical complexity that is dominated by a factor that is linear in n. The practical difference is shown in S1 Fig, in which the Wald test derived in this paper gives a speedup of around 20-100 compared to the likelihood ratio test implemented using the ILRS algorithm. The speed of the Wald test also compares favorably with the less general interaction test in the highly optimized Plink 1.9 software, we note that this test is limited to an A × A model with linear or logistic link function.

Unsaturated models. For an unsaturated model the design matrix P is no longer full rank, and for parameter estimation one usually must resort to an iterative algorithm. We develop a two-step closed-form estimation, in which we first transform the model to a convenient saturated model (e.g., G × G), by extending P and β with corresponding columns and parameters, respectively. We estimate all parameters in this model, and then back-transform them to the unsaturated model.

Let Pβ be the parameterization of the unsaturated model, and P_Sβ_S be the parameterization of a corresponding saturated model. For convenience of notation let ; we then obtain the following relationship between the parameters of the saturated and unsaturated model: β_S = Sβ. Given an estimate of the covariance matrix C_S of the saturated parameters we can estimate β and the corresponding covariance matrix C by and The idea underlying these estimators is that we can view the estimated parameters of the saturated model as an outcome of a multivariate Normal distribution with covariance matrix C_S and mean Sβ. Maximum likelihood estimation of β in this model leads to the equations above (this can be equivalently viewed as a generalized least squares estimation). We show in S1 Text that this leads to a consistent estimator.

The time complexity involves estimating the saturated model, transforming these estimates into estimates for the unsaturated model, and computing the final Wald test statistic. This time complexity is bounded by the estimation of the saturated model, so the total time complexity is still . For the score and LRT tests the time complexity is, assuming one parameter per variant, .

Fixed effect meta-analysis of multivariate Wald tests. We now describe how to combine association results, based on the Wald test, from multiple studies to perform a meta-analysis. Let denote the vector of estimated coefficients from study i and the corresponding covariance matrix. The combined across studies is estimated with and the combined covariance matrix with The combined Wald statistic across studies is then which follows a χ²-distribution with four degrees of freedom. Of note, each covariance matrix is commonly small and the time for computing the combined covariance matrix is shorter than the time for estimating each individual covariance matrix.

A practical problem for large-scale meta-analysis is that the result files from the interaction analysis of each study will be very large, typically in the order of 1 Terabyte. We suggest one possible solution to this issue. The analysis can be split into two stages. In the first stage, each study reports all variant pairs below some p-value threshold. A reduced set of candidate variant pairs are then created by taking the union of the significant pairs over all studies. In the second stage, all studies perform a second analysis of the reduced set of variant pairs. Finally, meta-analysis is performed on the result from the second stage. This effectively limits the storage space required, but may miss variant pairs with intermediate effects in all studies.

Relating different GLMs.

In this section, we demonstrate some aspects of how interaction models can be related to each other within the GLM framework.

Linear reparameterizations of saturated models are equivalent. We first show that a certain class of parameterizations are equivalent with respect to the Wald test, providing some insight to when two parameterizations are identical in terms of the inference of epistasis.

We start by demonstrating an important property of the Wald test; if two parameterizations with interaction parameters δ and δ′, respectively, are linear transformations of each other, i.e. δ′ = Bδ, then the corresponding Wald tests are equivalent, We observe that an important corollary of this is that a joint test, i.e. testing all interaction parameters simultaneously, for the G × G and AD × AD parameterizations is equivalent. In general, any joint test of interaction for saturated reparameterizations, will have identical results. However, this is not true if parameterizations are non-linear transformations of each other, i.e. if they use different link functions.

Of note, the degrees of freedom can be reduced if some of the genotypes contains very few samples. We can then avoid estimating the interaction parameters for these cells and adjust the degrees of freedom accordingly. This could, however, cause the GLMs to be non-equivalent in terms of the joint likelihood, and an interaction in one parameterization may not be reflected in another.

Relating parameter estimates. We show above that two saturated full GLM reparameterizations are equivalent in terms of the joint Wald test. However, different saturated full parameterizations may emphasize different interaction components, which potentially reflect different biological mechanisms. We now investigate this further.

The relationship between two saturated full parameterizations Pβ and Qβ′ is Importantly, the interaction parameters of β are a function solely of the interaction parameters of β′. This is in contrast to the main effect parameters of β, which may depend on both main effect and interaction parameters of β′. The relationship between the interaction parameters of G × G and AD × AD is, and We observe that the two parameterizations distribute the total interaction effect differently over the interaction parameters. For example, a single double homozygote interaction is apparent in G × G: δ^G×G = {0, 0, 0, 1}, but not in AD × AD: δ^AD×AD = {0.25, 0.5, 0.5, 0.25}. Conversely, a single additive-additive interaction is apparent in AD × AD: δ^AD×AD = {1, 0, 0, 0}, but not in G × G: δ^G×G = {1, 2, 2, 4}. We note that, for certain questions, exploring multiple parameterizations would be beneficial to allow the formulation of most-parsimonious constrained parameterizations, or to provide hypotheses on possible biological mechanisms for inferred interactions.

When performing the corresponding investigation of parameter relations of a non-saturated full parameterization, Pβ, to another parameterization, Qβ′, we find that the interaction parameters of β are no longer guaranteed to be a function of the β′ interaction parameters alone, but can depend on its main effects as well. More generally, constrained/non-saturated models reduce the number of interaction parameters (and sometimes also the main effect parameters) of saturated full parameterizations, either by constraining some parameters to be zero or to be functions of other parameters; some examples of the relation between the G × G parameterization and selected constrained or unsaturated parameterizations are given in Table 1.

Download:

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

Table 1. Relation between G × G interaction parameters and those of selected constrained or unsaturated models.

https://doi.org/10.1371/journal.pcbi.1005556.t001

The link function determines how fast the phenotype mean, μ, changes with the genotype. The choice among major classes of link function is further discussed, below, in Section Relating different GLMs and we will here focus on understanding the effect of a small change in the link function has on the parameters β. If the link function g(μ) is perturbed to g(μ) + ϵ (while keeping the design matrix P constant), we have the following approximate relationship between the parameters in the two models, where ∇_x f(x) denotes the gradient of f(x) with respect to x. This means that the incorrect link function will introduce a bias in the parameter estimates, and the magnitude of that bias will depend on the rate of change of the inverse of the link function, as well as on the design matrix (some examples are given in S1 Table).

Finally, the dispersion distribution models the variation around the phenotype mean μ (given by the link function and the parameterization). However, for a given data set, the dispersion distribution is seldom varied because different dispersion distributions are often tightly connected to specific types of phenotypes (see examples in S2 Table). In particular, if, for a given dispersion distribution, f, with (shape) parameter θ, there exists a link function g(μ) such that θ = μ, then g(μ) is called the canonical link function of f (S2 Table).

Examining the effect of different link functions. In the previous sections, we observed that a small perturbation of the link function may bias the parameter estimates, and, for general changes in the link function the corresponding Wald tests may not be equivalent. This suggests that some interactions may not be consistent between link functions. In fact, it is well-known that, under some circumstances, interactions inferred with one link function might be absent when another link function [25, 26] is used. Loftus [25] referred to these interactions as “uninterpretable”, because although they represent an interaction we cannot be sure unless we know the true underlying link function. However, Loftus also defined a class of “interpretable” interactions whose inference is invariant of the link function.

To make a robust statement of the existence of a particular interaction, we must therefore determine whether the interaction seems to be invariant of the link function. One previously proposed link function test investigates systematic bias in the residuals [24]. The test is called a goodness-of-link test and is constructed by first fitting GLM with the canonical link, to obtain a first estimate of , and then fitting the same GLM with as an additional covariate to investigate non-linearity. However, for saturated GLMs, the goodness-of-link test becomes over-parameterized, and therefore meaningless, because there are no degrees of freedom left after fitting parameters to accommodate (it does work in the special case of unsaturated GLMs).

We will instead use another approach that tests interaction over multiple link functions in a family of link functions, parameterized by some parameter λ [22, 23]; we will refer to this test as the link family test. Many such families have been suggested, but the most important property is to be able to model functions increasing either faster or slower than the canonical link function. For simplicity we have selected the following families for the Normal and Binomial dispersion distributions, respectively:

The first one is a symmetric version of the Box-Cox family of power transforms that lies between the -link (λ = 0) and the exp-link (λ = 2). It also contains the identity link as a special case λ = 1. The second family, the Pregibon family, was suggested for studies of link functions for Binomial outcomes which contains the logit-link (λ = 0) as a special case [24]. These families are by no means a complete representation of all possible links, but serve as an attempt to broaden the family of null models.

Evaluation of link function invariance is a computer-intensive approach as it requires re-analyses while stepping through different values of λ to identify any critical link function. In practice, we therefore apply this test on significant pairs, only, as an a posteriori measure of link appropriateness.

In the following sections we investigate how the false positive rate and statistical power is affected by different assumptions on the genetic effects. In both cases, we evaluate the following parameterizations: the G × G joint, AD × AD joint, G × G separate, A × A, and R/D × R/D tests (described below). We then apply the G × G joint test to two biological data sets to demonstrate the efficacy and plausibility of our method. In the following sections, we will (unless it is clear from context) discriminate between a generative model referring to the model from which data is generated and a test model, which is used in the statistical test. For details on data generation, please see Materials and methods.

Lipoprotein(a).

We performed a classic discovery-replication analysis using quantitative data for circulating Lp(a). Lp(a) is an LDL-like lipoprotein particle also containing apolipoprotein [30], circulating in the blood. The concentration of circulating Lp(a) has been associated to coronary heart disease [31] and has attracted attention as a possible target for lipid-lowering therapies [32]. For the discovery we used the PROCARDIS cohort [33], comprising in total 8,112 individuals genotyped for 566,865 variants, while replication was done in the SCARF-SHEEP cohort [34, 35], comprising 2,345 individuals and 116,540 variants.

A full chip-wide analysis, testing all possible pairs from L, the set of variants in approximate linkage equilibrium (≈10¹⁰ pairs), resulted in no variant pair passing the very severe Bonferroni-corrected significance level (p < 10⁻¹²).

We therefore investigated a variance heterogeneity-based genome-wide association study (vGWAS) [36] approach as an a priori selection procedure for single variant interaction candidates The vGWAS analysis is claimed to be independent of, and not biasing ensuing interaction analysis; in fact a theoretical proof is provided [36]. In a vGWAS, each variant is first tested with a heterogeneity test (in this case the Brown-Forsythe test [37]). Variants that display sufficient variance heterogeneity comprise the set V. We used a liberal p-value threshold of 10⁻⁴. The final interaction analysis was performed on all pairs between V and L. A single variant pair (rs13209686-rs1884480) was significant from this analysis (adjusted p = 0.0306; Table 2).

Download:

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

Table 2. Association results for the discovery and fine-mapping analyses.

https://doi.org/10.1371/journal.pcbi.1005556.t002

Because the variant pairs in the first analysis were first pruned for LD, we fine-mapped the association signal. For each of the variants i ∈ (1, 2) of the significant variant pair from the discovery analysis, we defined the set F_i of all variants in close proximity of i, and then analyzed all pairs of variants from F₁ and F₂. The top three pairs in the region are shown in Table 2. Only one of these variant pairs, comprising rs3103353 and rs9458157, was available for replication in our replication cohort SCARF-SHEEP. The replication result (Table 3) suggests that the discovery signal indeed replicates (p = 4.7378e − 07); the value of λ is consistent between discovery and replication and indicates that link should increase faster than linearly. Lastly, we performed a fixed-effect meta-analysis of the PROCARDIS and SCARF-SHEEP results for the rs3103353-rs9458157 pair. This resulted in a substantially stronger signal (p = 1.06 ⋅ 10⁻¹²; Table 3), indicating consistent effects in the two cohorts. This is also strengthened by the consistency of the estimated models for the two cohorts (S2 Fig), showing markedly lower expected Lp(a) levels for the double minor homozygote.

Download:

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

Table 3. Replication analysis of LP(a).

https://doi.org/10.1371/journal.pcbi.1005556.t003

Both rs3103353 and rs9458157 are located in the proximity of the previously associated LPA locus [29], but are not in LD (r² = 0.0219008 and D′ = 0.186003). Functional annotation analysis using HaploReg [38](S3 Table) showed that there is previous evidence that rs3103352 is an eQTL for SLC22A2 [39], while rs9458157 is situated in an intron of AGPAT4. Both variants are situated in regions with histone or methylation patterns associated with enhancer activity—in adipose tissue and in blood or cardiovascular tissue, respectively. Gene Ontology analysis [40](S4 Table) indicated that SLC22A2 is a membrane transporter of organic cations, including histamine, choline and dopamine, while AGPAT4 is involved in metabolic processes related to triglyceride and phospholipid metabolism.

Because of the proximity to the LPA locus, it is of interest to investigate the relation of this interaction to other single SNP associations in LPA, in particular those identified in previous studies (e.g., rs3798220 and rs10455872 [29] and rs41272114 [41]). We performed an additional analysis in the replication cohort, where we first identified all single variants associated to Lp(a), and then included each of these variants in turn as a covariate in the interaction analysis of the rs3103353-rs9458157 pair. The result (Supplementary table S10 Table) showed that the interaction association is not diminished by the additional covariate variants, suggesting that the identified interaction and neither of the associated single variants are proxies for the same association. However, the previously published associated variants rs10455872 and rs41272114 were not genotyped in our cohorts. Thus, we cannot test the possibility that the discovered interaction reflects the same association as one of these two variants or, for that matter, of any unsampled variant associated to Lp(a) in our cohort (i.e., either the interaction being a proxy for the single-variant [42] or vice versa [43]). We therefore performed a simulation study where we, using the PROCARDIS genotype data for chromosome 6, generated a continuous phenotype from a randomly picked single variant and then applied the same interaction analysis as for the biological data, but with different combinations of LD-pruning and vGWAS selection of variants; this was repeated 200 times (for details see S2 Text). The results (S8 Fig) shows that, while both approaches involving no pruning or LD-pruning controls the family-wise error rate (FWER), vGWAS does indeed inflate the FWER substantially. This contradicts the claim by [36] that vGWAS does not affect ensuing interaction analysis, this is likely caused by LD between tested variants which is assumed to be absent in the theoretical proof of Pare et al. [36]. To conclude, we cannot exclude the possibility of the existence of unsampled variants that display the same association as the discovered interaction between rs3103353 and rs9458157.

Myocardial infarction.

MI is one of the common endpoints of coronary artery disease (CAD). Typically for complex diseases, such as CAD (and MI), very large sample sizes are needed for GWAS. Major collaborative efforts aimed at genetic analyses of CAD are therefore being performed, typically in large-scale meta-analyses driven by international consortia. As a proof-of-concept, we performed a two-cohort meta-analysis of gene-gene interactions associated to myocardial infarction (MI).

The first cohort was obtained by augmenting the PROCARDIS MI case-control data with control data from the WTCCC [44] (resulting in a total of 10,139 samples for analysis), the second cohort was the MIGEN cohort [45] (6,042 samples). This analysis only included variants that were genotyped in both cohorts (132,181 variants). We tested for interaction between all pairs among variants (≈ 8.74 ⋅ 10⁹ pairs) using a fixed effect meta-analysis, as described above in Section Fast estimation and testing of interaction in generalized linear models. No variant pair was declared significant using the Bonferroni correction (p < 5.72 ⋅ 10⁻¹²), the top five variant pairs (best p-value p = 1.91 ⋅ 10⁻¹⁰) are shown in S9 Table.