Assuming conditional independence of haplotypes and exposure, we develop a test for the omnibus null hypothesis that any risk-related multi-locus diplotypes act independently of the exposure, in the sense that the joint effects of genotype and exposure on disease risk are multiplicative. Operationally this null means that, even though transmissions of specific haplotypes to affected offspring might violate Mendelian proportions, those transmissions should nonetheless be statistically independent of the exposure.
Under the no-interaction null, one might expect the covariance between a transmission score
at locus j
, defined by Dij
, and the exposure Xi
to have expected value 0 for every j
. We need, however, to account for the possibility that the exposure participates in bias-inducing population stratification, which can produce an overall correlation between transmissions to affected offspring and exposure, even under this null [1
]. Just as when studying a single SNP, protection against bias requires a strategy that groups together families with similar parental genotypes. Our approach calculates covariances within such groups of families.
For the moment imagine that we know phase and can stratify on parental diplotypes. Let g
index strata defined by particular pairs of parental diplotypes. Under the no-interaction null, the covariance between Dij
should be 0 within each stratum regardless of exposure-related genetic population stratification; under alternatives, transmission of some SNPs to affected offspring will instead be correlated with exposure. With that background as motivation, we designed our test procedure to estimate the average conditional-on-g
(within-stratum) covariance at each component SNP. We make use of the following decomposition (suppressing subscript i
) for locus j:
denotes the expected value taken over g
, and Cov
) denotes the covariance between two specified arguments conditional on parental diplotypes. The last term on the right in (1) can be nonzero under exposure-involved population stratification but does not reflect within-stratum correlation between Dij
. This decomposition indicates how to distinguish an interaction signal from population-stratification-induced bias and serves to motivate the development of our testing approach.
The mean of the cross product terms DijXi
across all triads, namely:
estimates the first term on the right in equation (1). We can generate a corresponding permutation distribution for locus j
by repeatedly permuting exposure values among families within each g
stratum and computing a cj
each time. Let
denote the average cj
across permutations. This permutation process recapitulates the null hypothesis by allowing transmissions to behave independently of exposure within each stratum. For the permutation distribution, then, the left side of (1) is zero by construction, and the two terms on the right must equal each other. Thus,
, which estimates E
) for permutations, also estimates Eg
)] for permutations. Consequently, we propose to regard permutation-based
as an estimate of Eg
)] under the null and to use it to adjust for bias in cj
due to exposure-involved population stratification.
To construct a multi-locus test of no haplotype-by-exposure interaction, we consider a vector with component cj
for each locus j
= 1, 2, 3, …, k
). Because the bias-related second term on the right in (1) will differ for different loci, from each component cj
, its permutation-based mean under the null, to form
. We then form a test statistic
as the sum of squared vj
, that is,
We use the within-stratum permutation procedure described above to generate the permutation distribution for
A modification of the statistic may enhance power under certain circumstances. The k
SNPs may contribute unequally to the genetic risks, hence to the interaction, especially when a disease-susceptibility mutation is one of the typed SNPs and that SNP not only confers risk among the unexposed but also interacts with the exposure. In that circumstance one could potentially increase the power by combining the vj
in a weighted fashion, using the estimated marginal genetic main effects of marker SNPs as weights. Accordingly, in a second procedure we used zj
, the j
-th component of the Z
statistic from TRIMM [11
] (the standardized mean of the j
-th column of D
) as a weight for SNP j
and formed an alternate test statistic
The use of Tmax–ZV
should also tend to favor finding interactions when the genotype main effect and a genotype-by-exposure interaction both involve the same haplotype.
In addition to the two test statistics
, we also considered a combined test, Tcombined
, that takes advantage of the complementary strengths of both. To do this, we used the sum_log(P) approach of Shi et al. [11
]. If the interaction effect is due to a single causative SNP, a SNP that was fortuitously included in our measured set of k
markers, then we expect Tmax–ZV
to do better than
. If instead a susceptibility haplotype is related to risk, or the SNPs studied behave collectively as markers for an untyped causative SNP or haplotype, then the advantage should go to
. Consequently, the combination of the two might perform well under either scenario. We refer to the test based on Tcombined
as GEI-TRIMM (gene-environment interaction triad multi-marker method).
With unphased genotype data, parental diplotypes are typically unknown, so instead of strata defined by parental diplotype pairs, we used strata defined by the vector of sums of parental genotypes (Mi,1 + Fi,1, Mi,2 + Fi,1, L, Mi,k + Fi,k). Strata defined this way represent a coarsening of strata based on parental diplotypes. Consequently, tests based on those strata are less than fully robust but remain resistant to bias from population structure. Other than this re-definition of strata, one generates the permutation distribution as described above. The empirical p value is the proportion of permutated-data-based test statistics that are larger than or equal to the observed-data-based test statistic.
We used simulations to evaluate the performance of our proposed tests in a realistic context. We simulated haplotypes based on HapMap-phased genotype data from the sample with European ancestry, using haplotypes and their frequencies for a 100-kb region around the replication factor C1 gene (RFC1). This haplotype set is defined by 12 SNPs, which specify 17 different haplotypes (Appendix). We then created a thirteenth SNP, representing a causative mutation, whose variant allele resided on haplotype 1. We also constructed a second haplotype set with fewer SNPs by using 5 LD tagging SNPS for RFC1. We again introduced a causative SNP, residing on the same haplotype. This second set contained 12 haplotypes (Appendix). In both these haplotype sets, only one haplotype was associated with risk and its frequency was typically 0.28, although we also examined the influence of risk-haplotype frequency on the performance of our tests. Our simulations always involved only a single risk-associated haplotype; in contrast, our test procedures neither presuppose a particular number of risk haplotypes nor designate any haplotypes as candidates.
We simulated both dichotomous and continuous exposures. Dichotomous exposures were generated via Bernoulli trials with outcome probability equal to the designated exposure prevalence. We also simulated continuous exposures from a zero-enriched truncated lognormal distribution. A proportion of individuals had no exposure (again generated via Bernoulli trials). For the remaining individuals, we generated exposure levels by exponentiating a standard normal random number and rejecting values above 6.
Besides specifying distributions of haplotypes and exposure for our simulations, we need a risk model to generate case-parent triads. Our simulations allowed the risk associated with haplotypes to depend on exposure level; but the risk model that we used, though convenient for simulations, is neither required for the validity of nor imposed by our tests. Because the case-parent triad design precludes estimation of effects of exposures, our simulations assumed for simplicity that exposure was unrelated to risk among offspring with no copies of the risk-associated haplotype. Suppose X0
is a designated reference level of the exposure and let H
denote the number of offspring-inherited copies of the risk-associated haplotype. In our simulations, X0
was always set to zero. We specify the genetic main effects parameters Rh
= 1 or 2 as:
Thus, haplotype relative risks at the reference exposure level depend on the inherited number of copies, H
, conditional on parental genotypes. For δ = X
, we define interaction parameters for the risk-associated haplotype as Ih
(δ) for h
= 1 or 2 as:
(δ) could be any function of δ, but for simulations we used a log-linear function, so that the log Ih
(δ) of was proportional to δ, a function we can characterize by specifying choices of Ih
(1) (simply denoted Ih
) for h
= 1 or 2. Consequently, risk models in our simulations, whether exposure is dichotomous or continuous, can be represented by the risk vector (R1
) for the designated risk-associated haplotype. Implicitly, all other simulated haplotypes have the risk vector (1, 1, 1, 1).
Extending the notation established above from a single risk haplotype to include separate functions Ih(δ) for each haplotype under study allows a more precise statement of the null hypothesis of interest, namely, H0: I1(δ) = 1 = I2(δ) for all δ and for all haplotypes studied. This null asserts that, conditional on parental diplotypes, the relative risks associated with having inherited copies of any haplotype do not vary with the level of the exposure. This null is equivalent to one where the effect on risk of a change in exposure level does not depend on how many copies of any haplotype were inherited. Alternatives would include situations where inheritance of one specific haplotype influences risk in the unexposed while a different haplotype may interact with the exposure. For simplicity, however, we simulate only scenarios where the same haplotype that contributes to risk in the unexposed might also interact with the exposure.
For each simulated scenario, we generated 1,000 data sets, each with 1,000 triads. We generated triads by sampling parental pairs (with population stratification, both parents came from the same subpopulation), randomly generating the offspring haplotypes based on Mendelian inheritance and an offspring exposure from the scenario's exposure model, and then generating disease status randomly from diplotype and exposure through the scenario's presumed risk model. Cases were retained until 1,000 case-parent triads were generated. For each risk scenario, we considered four situations: the mutation SNP was either typed or not and genotypes had either 0 or 20% of loci randomly missing. We applied our three proposed tests to each simulated data set. This tactic enabled us to examine whether Tcombined
(GEI-TRIMM) had power superior to Tmax–ZV
when the causative SNP was genotyped and power superior to
when the causative SNP was not genotyped. We used 1,000 permutations with each simulated data set to estimate p
To examine validity of our tests under population stratification, we simulated a no-interaction null scenario with a dichotomous exposure and two equal-sized subpopulations, each subpopulation having haplotypes in Hardy-Weinberg equilibrium (HWE). We used only the second haplotype set with 6 SNPs (Appendix). The disease ratio between subpopulations was 1:1.5 (ratio of risks for unexposed individuals with no copies of the risk haplotype). Risk-haplotype frequency and exposure prevalence were both 0.2 in one subpopulation and both 0.4 in the other. Thus this no-interaction scenario had both genetic-and exposure-related population stratification. We set (R1, R2, I1, I2) = (1.73, 3, 1, 1) in each subpopulation.
To further examine the Type I error behavior under population stratification, we used RFC1 with 6 SNPs again; but we reduced the risk haplotype frequencies and varied the proportional representation of the remaining haplotypes between the two subpopulations. We set the risk haplotype frequencies at 0.1 and 0.05 in the two subpopulations, respectively. In the first sub-population, all non-risk haplotypes retained their original HapMap frequencies rescaled to total 1–0.1; in the second sub-population, the non-risk haplotype frequencies were rescaled to total 1–0.05 and then re-assigned at random so that the frequency ranking of haplotypes differed between the two subpopulations. The exposure prevalences, disease ratio, and relative risk parameters were as in the previous paragraph. We constructed 12 distinct simulation scenarios by letting each haplotype in turn be the risk haplotype, and we simulated 1,000 data sets for each scenario.
We examined power for detecting interaction with a dichotomous exposure where either 13 or 6 SNPs (Appendix) were genotyped. Samples came from a homogeneous population, where the exposure prevalence was 0.2. For simplicity, we simulated genotypes in HWE (not required for test validity). We simulated data with four (R1, R2, I1, I2) settings: (1, 1, 2, 4), (1, 1, 1, 4) (1, 2, 2, 2), or (1.414, 2, 1.414, 2).
We also studied a quantitative exposure where the risk haplotype and exposure had a pure interaction, that is, haplotype was unrelated to risk among the unexposed. We drew triads from a population consisting of a mixture of two equal-sized subpopulations, each with haplotypes in HWE and with disease ratio 1.65. The proportions of unexposed individuals were 0.2 and 0.4, and the haplotype frequencies were 0.5 and 0.2 for the two subpopulations, respectively. We simulated scenarios with four (R1, R2, I1, I2) settings: a null scenario (1, 1, 1, 1) and three alterative scenarios (1, 1, 1.2, 1.44), (1, 1, 2.4, 2.4) and (1, 1, 1, 1.36).
To study the influence of risk-haplotype frequency on power, we generated triads either from a homogeneous population or from a stratified population formed from two equal-sized subpopulations each in HWE and with disease ratio 1.65. We used a quantitative exposure with unexposed prevalence of 0.3 for the homogeneous scenarios and of 0.2 and 0.4, respectively, for each subpopulation in the stratified scenarios. In these simulations, we set (R1
) at (1.2, 1.5, 1.2, 1.44). For the homogeneous scenarios, risk-haplotype frequencies ranged from 0.05 to 0.95. For the stratified scenarios, risk-haplotype frequencies (f1
, respectively) were set so that the odds ratio:
and the overall haplotype frequency, namely, (f1
)/2, ranged from 0.05 to 0.95. For example, for an overall haplotype frequency of 0.3, the frequencies in the two subpopulations are 0.5 and 0.2 respectively, which meets the odds ratio specification with (f1
)/2 = 0.3.
We compared the performance of our method on complete genotype data with the pseudocontrol approach [16
]. Because the current implementations of UNPHASED [8
] and PCPH [7
] do not provide for saturating the genetic main effects (i.e. for fitting a separate parameter for each diplotype), we excluded them from our comparison. We evaluated Type I error rates under a null scenario with population stratification in our second haplotype set with 5 SNPs (causative SNP not typed). This scenario used a quantitative exposure with unexposed prevalences 0.1 and 0.9 and risk haplotype frequencies 0.9 and 0.1 for the two equal-sized subpopulations, respectively. The disease ratio was set at 1.65 and (R1
) at (1.73, 3, 1, 1). We also simulated populations with three haplotypes, each with two SNPs. We simulated a no-interaction null scenario with population stratification: a quantitative exposure with unexposed prevalences 0.2 and 0.4 and risk haplotype frequencies 0.5 and 0.2 for the two equal-sized subpopulations, respectively, with disease ratio 1.65. (R1
) was set at (1, 3, 1, 1). We also simulated a non-null scenario with a homogeneous population: a quantitative exposure with unexposed prevalence 0.3, risk haplotype frequency 0.3, and (R1
) set at (1, 1, 1.2, 1.44). For the pseudocontrol approach [16
], we used ‘genotype’ coding, as described in the program's R-package documentation, to saturate main effects by fitting a parameter for each diplotype and tested for log-additive interaction with a likelihood ratio test.