Simulation framework

We created sets of simulated phenotypes based on the genotypes of the HDMP RI panel, which is an admixed population in which the B6 strain, on average, contributes 50% of each strain’s DNA. For each simulated phenotypes, we drew a SNP (75% > MAF > 25%) at random from the HMDP genotypes and created a phenotype based on β, the genetic effect size, ϕ the effect size of the interaction between global ancestry (θ) and the chosen SNP and a multivariate normal (mvn) derived from three variance terms: σg2, the proportion of variance attributable to genetic effects σθ2, the proportion of variance attributable to Gxθ effects and σε2, the residual proportional variance attributable to all combined sources of error and variation not considered in this study. Phenotypes were generated both with and without the Gxθ variance term to ascertain the necessity of incorporating a second GRM (K^A), which would correct for relatedness in ancestry, into the algorithm.

We simulated four distinct phenotypes for our analysis

Phenotypes generated by including a SNP Effect

Phenotypes generated by including a Gxθ Effect

In each phenotype, σg2 was set to 0.4. When incorporated, σθ2 was set to 0.2 and σε2 was set to the remainder of the variance (0.6 or 0.4). The power of our model and independence of our β and β_Gxθ terms were queried by varying either β or β_Gxθ from 0 to 1 (200 values set 0.005 apart) with 1,000 simulated phenotypes at each step (200,000 total simulations per phenotype).

Using the same panel and set of SNPs described above, we drew one SNP at random to be our test SNP and one to ten additional SNPs to be simulated epistatically interacting SNPs (10,000 simulations per set of interacting SNPs, 100,000 simulations in total). For each simulation we created a composite SNP in which only strains with the minor (non-B6) allele in every one of the tests and interacting SNPs had that allele in the composite SNP. We used this composite SNP to generate a phenotype as described above with a χ² test statistic for the 2–11 interacting SNPs set to a large value of 20 to ensure a consistent and observable effect that could be recovered using either regular GWAS with pyLMM or our approach.