MATS: a novel multi-ancestry transcriptome-wide association study to account for heterogeneity in the effects of cis-regulated gene expression on complex traits

Abstract

The Transcriptome-Wide Association Study (TWAS) is a widely used approach which integrates gene expression and Genome Wide Association Study (GWAS) data to study the role of cis-regulated gene expression (GEx) in complex traits. However, the genetic architecture of GEx varies across populations, and recent findings point to possible ancestral heterogeneity in the effects of GEx on complex traits, which may be amplified in TWAS by modeling GEx as a function of cis-eQTLs. Here, we present a novel extension to TWAS to account for heterogeneity in the effects of cis-regulated GEx which are correlated with ancestry. Our proposed Multi-Ancestry TwaS (MATS) framework jointly analyzes samples from multiple populations and distinguishes between shared, ancestry-specific and/or subject-specific expression-trait associations. As such, MATS amplifies power to detect shared GEx associations over ancestry-stratified TWAS through increased sample sizes, and facilitates the detection of genes with subgroup-specific associations which may be masked by standard TWAS. Our simulations highlight the improved Type-I error conservation and power of MATS compared with competing approaches. Our real data applications to Alzheimer’s disease (AD) case–control genotypes from the Alzheimer’s Disease Sequencing Project (ADSP) and continuous phenotypes from the UK Biobank (UKBB) identify a number of unique gene-trait associations which were not discovered through standard and/or ancestry-stratified TWAS. Ultimately, these findings promote MATS as a powerful method for detecting and estimating significant gene expression effects on complex traits within multi-ancestry cohorts and corroborates the mounting evidence for inter-population heterogeneity in gene–trait associations.

Introduction

The Transcriptome-Wide Association Study (TWAS) is a common approach to identify genes associated with complex traits and diseases. TWAS integrates expression and GWAS data to study the genetically regulated component of gene expression (also referred to as imputed gene expression), computed by taking a weighted sum of cis-eQTL genotypes, with weights derived from each Single Nucleotide Polymorphism's (SNP) estimated effect on gene expression (GEx). These weights can be obtained from GEx prediction models trained on independent data from a target case–control study, making TWAS applicable to large-scale genetic studies which do not have expression data. Unsurprisingly, imputed GEx is closely tied to (i) the underlying genetic architecture of gene expression (i.e. causal eQTLs and their effect sizes) and (ii) the accuracy and appropriateness of the prediction model for obtaining eQTL weights. Both of these factors can be heavily influenced by ancestry.

Research suggests that at least 25% of gene’s expression levels are significantly correlated with self-reported or local genetic ancestry, findings which are complemented by subsequent studies that identify variation in the underlying genetic architectures of GEx between populations (1–9). Population variation in cis-eQTLs can arise for a number of reasons, including variation in LD patterns and minor allele frequencies. A 2018 study by Mogil et al. revealed that the level of concordance between population eQTL effect sizes is proportional to their degree of shared ancestry, and that gene expression prediction and subsequent TWAS tests are optimized when the ancestry and MAF patterns of the prediction training sample and the target GWAS are similar (9).

The downstream effects of population genetic variation in GEx on complex traits are not fully understood to date, but it is increasingly clear that the impacts may be widespread throughout the genome and phenome (10,11). Bisogno et al. (12) showed that ancestral variation in gene expression can explain interpersonal differences in cell reprogramming, cell fate transitions and biological processes. Recent work by Blue et al. (13) revealed ancestry-linked heterogeneity in the effects of the protective (e2) and risk (e4) alleles for the APOE gene on Alzheimer’s disease (AD), suggesting that the effects of well-established loci involved in AD may be modified by ancestry.

Heterogeneity in GEx–trait interactions has also been directly observed through TWAS discoveries. Notably, TWAS may amplify inter-population heterogeneity by directly modeling GEx as a function of cis-eQTLs, which vary by ancestry. A race-stratified TWAS in breast cancer revealed heterogeneity in expression effects between African American and White women, including some ancestry-specific GEx associations which were not uncovered by unstratified standard TWAS (14). However, power may be lost in stratified analyses due to decreased sample sizes of each ancestry, particularly for genes with shared effects across populations. Therefore, while it is prudent to account for ancestry-specific GEx effects in TWAS, stratified TWAS may have limited power to detect homogeneous GEx–trait associations. The recent METRO model proposed by Li et al. incorporates the ancestry-specific expression prediction models to improve gene expression imputation (via a data determined weighted average of these prediction models). METRO estimates a single gene-trait effect for a GWAS dataset of a single ancestry and, therefore cannot estimate ancestry-specific gene-trait effects (15). As highlighted in our simulations, when applied to multi-ancestry GWAS data, METRO has low power to detect any gene expression effect in the presence of strong but oppositely signed effects between ancestries. Motivated by (16), we also explicitly account for individual-level heterogeneity of gene-trait associations, which is ignored in METRO.

In this study, we present the Multi-Ancestry TwaS (MATS), a novel approach to account for heterogeneity in the effects of imputed gene expression which are correlated with ancestry. By jointly analyzing samples from multiple populations, MATS is expected to improve power for detecting genes with shared expression–trait associations across populations through increased sample sizes, as compared with existing stratified TWAS approaches. Our proposed model of ancestry correlated GEx effects builds off of two existing and closely related approaches for modeling population substructure in genetic association testing, the generalized linear mixed model (GLMM) and principal component regression (PCR) (17). Specifically, we model population structure using a genetic relatedness matrix and its corresponding axes of genetic variation (i.e., Principal Components (PCs)), respectively, for GLMM- and PCR-based approaches (18–20). Using these models, we present a novel TWAS analysis pipeline that distinguishes between expression associations which are (i) shared across populations, (ii) population-specific and/or (iii) subject-specific. This pipeline relies on a sequence of tests to identify population-level and intra-population heterogeneity. We show the equivalence of these tests under both GLMM- and PCR-based modeling frameworks. Our test of intra-population heterogeneity is based on the empirical Bayes test for high-dimensional data of Goeman et al. (21) and has direct parallels to Kernel Machine Regression and related applications in genetic association studies (e.g. SKAT (22), C-alpha (23), SSU (24), MR-MEGA (16)).

In our primary analysis, we apply MATS to 2757 African American and 6710 Caucasian samples from the Alzheimer’s Disease Sequencing Project (ADSP), using ancestry-matched prediction models trained separately on African and European samples from the Multi-Ethnic Study of Atherosclerosis (MESA) study. We discover a number of genes with significant expression effects in AD that were undetected by standard TWAS, including SNX7, UGDH, HYAL3, CD68 and SUPT4H1. In a supplemental application, we apply MATS to a cohort of finescale population substructure, namely 37432 samples from the UK Biobank (UKBB) with neuroimaging phenotypes, which span four genetically determined subpopulations (African, Asian, European and Middle Eastern ancestry). Specifically, we study expression associations in 20 AD-associated imaging-derived phenotypes (IDPs) and 4 hematological traits. In the absence of ancestry-specific expression prediction models for each of the four UKBB subpopulations, we impute gene expression using the European expression predictions model in all samples. Nevertheless, we discover a handful of genes with suggestive population- and individual-level heterogeneity, highlighting the flexibility and power of MATS even when ancestry-specific expression models are unavailable.

These findings highlight the utility of our MATS framework for identifying expression-trait associations in multi-ancestry cohorts, even when the expression effects differ by subpopulation. Furthermore, our discovery of ancestrally heterogeneous expression–trait associations discourages against generalizing TWAS findings from predominantly European genetic studies to other populations and highlights the importance of expanding non-European cohorts in richly phenotyped genomic studies.

Results

Simulations

Here, we highlight key takeaways from our simulation study, which are outlined in greater detail in the Supplemental Material. Our first set of simulations assess the performance of MATS compared with METRO and standard TWAS in the absence of individual-level and population-level heterogeneity. As shown in Figure 1, the Type I errors for MATS are well controlled for the test of overall expression and population-level heterogeneity and are conservative for the test of individual-level heterogeneity. As is expected in the setting without heterogeneity, the TWAS Type-I errors are also well controlled, while METRO yields a conservative Type-I error rate.

Figure 1

Simulation Rejection Rates with varied total expression effect, in the absence of population level and individual level heterogeneity.

As the shared expression effect increases (α₁), we find that TWAS and MATS have greater power than METRO across all settings for α₁ (Fig. 1C). Unsurprisingly, TWAS is slightly more powerful than MATS when there is no heterogeneity, owing to fewer degrees of freedom and smaller standard errors (Supplementary Material, Fig. S1). However, when we introduce population-level heterogeneity (Fig. 2C; α₂ = 0.003, α₃ = −0.003), TWAS and METRO are unable to detect expression–trait associations when present within only two populations. This finding is expected since TWAS tests for a non-zero global mean expression effect, which is close to zero when α₂ = −α₃ and the shared expression effect (α₁) is small. In contrast, the MATS test of expression effects and population-level heterogeneity retains high power and unbiasedness across all magnitudes of α₁, even when only two populations have true non-zero expression effects (Fig. 2, Supplementary Material, Fig. S2). In fact, the MATS test of population-level heterogeneity and overall expression effects remains powerful across most reasonable magnitudes of α₂ and α₃ (Fig. 3A and B, α₁ = 0, vary α₂ = −α₃), even when only one population has a substantial expression effect (Fig. 3C and D, α₁ = 0, α₃ = −0.001, vary α₂). The same findings hold when one population’s sample size is significantly larger than the rest; we leave these simulation figures to the Supplemental Material.

Figure 2

Simulation Rejection Rates with varied total expression effect, population level heterogeneity and no individual-level heterogeneity.

Figure 3

Simulation Rejection Rates with varied population-level heterogeneity and no individual-level heterogeneity.

Figure 4 illustrates the rejection rates when perturbing individual-level heterogeneity. Specifically, we vary κ² when α₁ = α₂ = α₃ = 0. Notably, while the power for the test of individual-level heterogeneity increases as κ² increases, the rejection rates for population-level heterogeneity (H_0,4) and overall expression (H_0,1) also increase. While the true values for the main effects are zero (α₁ = α₂ = α₃ = 0), the rejection rates for these tests are not strictly Type-I error rates, given the assumption under both hypotheses that κ² = 0. Intuitively, the increased rejection rates for both H_0,1 and H_0,4 with growing κ² can be explained by spurious expression and inter-population signal introduced when the GEx–trait association is highly variable between subjects.

Figure 4

Simulation Rejection Rates with individual-level heterogeneity.

We further study the effects of extreme misspecification of K when α₁ = α₂ = α₃ = 0, specifically setting K = I, an n × n identity matrix (i.e. misspecifying that all subjects are entirely unrelated). As illustrated in Supplementary Material, Figure S4, we observe mild inflation in Type I errors for the test of population-level heterogeneity for high levels of individual-level heterogeneity, but the Type-I errors remain well conserved across moderate levels of inter-subject heterogeneity. Power for the test of individual-level heterogeneity remained high, in spite of the misspecification.

Importantly, the simulations illustrated in Figures 2 and 3 reflect settings where the non-negative weight constraint in METRO is violated, thus influencing its highly conservative performance. Nevertheless, we feel that the settings in which one or more populations may have opposite signed expression effects can be reasonably expected in real data, which we confirm in our real data applications and highlighted in a number of motivating examples throughout Introduction section. Furthermore, the simulations illustrated in Figures 1 and 4 reflect settings in which the METRO conditions are satisfied, since all population-specific effects are in the same direction.

To parse apart whether the association signal is coming from population-level vs. individual-level effects, we recommend fitting the PC-interaction model outlined in the Methods to follow when H_0,3 is rejected to adjust for individual-level heterogeneity by incorporating top PC-GEx interactions. We expect that the Type I error rates for H^*_0,4: α₂ = α₃ = 0 and H^*_0,1: α₁ = α₂ = α₃ = 0 under this full model will be better conserved when the leading PCs explain a large proportion of the variability. However, if a substantial proportion of variability is explained by the remaining PCs, we expect that there may still be some Type I error inflation. In Supplementary Material, Figure S3, we illustrate the rejection rates for LRTs of H^*_0,1 and H^*_0,4 under the full PC-interaction model using the top 10 PCs, compared with standard TWAS. As expected, the Type-I errors for population-level heterogeneity under the full model are improved over the model that does not include PC-interaction effects. However, we still observe notable Type-I error inflation as κ² increases, likely owing to missing signal when including only the top 10 PC-interactions, which explain 45% of the genetic variation in the simulated data.

Application to the Alzheimer’s disease sequencing project

Out of 1358 total prediction models from MESA, 1157 (African American (AFA)) and 1311 (Caucasian (CAU)) genes included at least one variant which was not in the target ADSP Whole Genome Sequencing (WGS) data (Fig. 5A). We applied our proposed method for integrating prediction and target model data using matched ancestry samples from 1000 genomes and successfully recovered complete information for all SNPs in 1160/1358 ≈ 85% and 1012/1358 ≈ 75% genes for the AFA and CAU prediction models, respectively. Figure 5B illustrates the number of recovered SNPs across all genes, highlighting the success of our proposed approach for salvaging information on missing prediction model variants.

Figure 5

Integrating GEx prediction and target data.

Figure 6A gives the PC plots for the top 6 PCs based on the realized genetic relatedness matrix of all ADSP samples, revealing some clustering between self-reported ancestries based on the top two PCs. The proportion of variance explained by the top 3 and 10 PCs is 10.4 and 16.4%, respectively. Figure 6B illustrates PC plots computed on a 9467 (n) × 1358 matrix of imputed gene expression, where each column represents an imputed gene. These plots show notable separation between the AFA and CAU samples based on imputed expression, which can be largely explained by the different prediction models used to imputed GEx expression in AFA and CAU populations, differences of which are driven by variability in the genetic architectures of GEx between populations, impacting downstream interactions between expression and the complex trait.

Figure 6

PC plots computed on A) genotypes and B) imputed IDPs in the ADSP cohort.

MATS identified 20 genes with significant expression effects on AD at a Bonferonni-corrected threshold of P ≤ 0.05/1358, and 363 genes with false discovery rate q-values ≤ 0.01. Notably, none of the 20 genes reaching the Bonferroni-corrected significance level were identified by standard TWAS, and 17 were also significant for the MATS test of individual-level heterogeneity. The remaining three genes included UGDH (p_MATS = 5.643 * 10⁻⁶, q_MATS = 0.00016), SNX7 (p_MATS = 1.618 * 10⁻⁶, q_MATS = 0.00016) and CAMK1D (p_MATS = 3.017 * 10⁻⁵, q_MATS = 0.00021). As to be discussed, both SNX7 and UGDH yielded a significant test for population-level heterogeneity at an unadjusted significance level of 0.05, suggesting that their discovery was likely driven by population-level effect differences. For CAMK1D, the MATS tests of population-level and individual-level heterogeneity yielded p-values of p_pop.het = 0.1048 and p_indiv.het = 7.758e−05, respectively, with log odds ratio (OR) estimates for the AFA and CAU populations of 0.0837 and −0.186, respectively. These findings suggest that the significant overall association test was driven by the strong (but not independently significant) signal from the inter- and intra-population expression effects and may warrant further investigation of CAMK1D as a putative gene involved in AD.

The Venn diagram in Figure 7 illustrates the number of genes identified by standard TWAS, CAU and AFA specific TWAS, and the MATS test of overall expression at an unadjusted P-value threshold of 0.05. Notably, we are limited to the comparison of genes identified at an unadjusted significance level given that no genes were detected by standard or stratified TWAS with p≤ 0.05/1358 or q ≤ 0.1, aside from SNX7 under the CAU-specific TWAS model (p_CAU-TWAS = 3.129*10⁻⁷, q_CAU-TWAS = 0.00035). MATS identified 344 genes which were not identified under the stratified or standard TWAS models (p < 0.05), 342 of which had suggestive evidence of individual-level heterogeneity and 8 of which had suggestive evidence of population-level heterogeneity. Of these eight genes, five remained significant for population-level effect differences (p < 0.05) after controlling for intra-population heterogeneity using three leading PC-interactions, namely PQLC3, ALS2, CHURC1, C17orf75 and DNTTIP1. Furthermore, 170 of the 344 genes identified uniquely by MATS had same-signed effect estimates, highlighting the benefit of MATS over standard TWAS to identify genes with expression-AD signal, even when the overall mean effect is not attenuated by opposite signed population-level effects. Many of the genes exclusively detected by MATS have existing links to AD, including (but not limited to) WDFY1, CCT5, CDKN2B (25), SIGLEC7 (26), and SLC22A15 and SLC22A16 (27).

Figure 7

Comparison of the number of significant genes identified by competing methods for ADSP.

MATS was also able to detect a large proportion of genes identified by a population-stratified TWAS model, specifically 29/59 and 57/86 genes discovered uniquely by the AFA- and CAU-specific TWAS tests (p< 0.05), respectively. Notably, a number of genes were uniquely identified by the AFA-specific (23 genes) and CAU-specific (12 genes) models which were not discovered by MATS. In spite of non-significant tests under the MATS framework (due to its larger degrees of freedom), the population-level effect estimates and 95% confidence intervals for these genes were highly concordant between MATS and stratified-TWAS for most genes (Fig. 8). On the other hand, we also note the power gains of MATS in detecting many more genes missed by the stratified (and standard) TWAS models, and yielding significant genes after correcting for multiple testing.

Figure 8

Estimated expression effects from standard TWAS vs. MATS for ADSP.

Forty genes were significant at a Bonferonni-corrected significance level (p < 0.05/1358) for individual-level heterogeneity in AD. The five genes which yielded the smallest p-values included ISG15 (p_indiv = 6.891 * 10⁻⁷), WDFY1 (p_indiv = 1.308 * 10⁻⁶), IFT88 (p_indiv = 1.308 * 10⁻⁶), PRELID1 (p_indiv = 2.369 * 10⁻⁶) and CCT5 (p_indiv = 2.402 * 10⁻⁶), none of which were significant for the MATS test of population-level heterogeneity standard TWAS, or stratified TWAS. A comparison of the population-level effect sizes for the full (with random effect for individual-level effect) and reduced model (without random effect) is illustrated in Supplementary Material, Figure S10. Interestingly, all of these genes were previously identified as susceptibility genes for AD or related dementias in the literature (ISG15 (28,29), WDFY1 (30), IFT88 (31), PRELID1 (32) and CCT5 (33)). Our findings suggest that the mechanisms underlying these expression-AD associations may be patient- or subgroup-specific, or may be influenced by environmental or batch effects within this sample. Ultimately, these genes can be promoted as candidate targets for personalized diagnosis and treatment which would likely go undetected in standard approaches focused on study or population-level mean effects.

At a Bonferroni-corrected significance level, no genes were identified by MATS for population-specific expression effects on AD, and only SNX7 was identified under the CAU-specific TWAS model. The lack of significant findings can be attributed to the small sample size of the ADSP data for detecting small population-level effect differences. At an unadjusted significance level of 0.05, 440 genes were identified by the MATS expression-association test, 401 of which were also significant (p < 0.05) for the MATS test of individual-level heterogeneity, and 43/440 were identified for population-level heterogeneity. The top 5 genes with the smallest p-value for population-level heterogeneity were SNX7 (p_pop = 0.00028), CD68 (p_pop = 0.0014), UGDH (p_pop = 0.00218), SUPT4H1 (p_pop = 0.00237) and HYAL3 (p_pop = 0.00247). None of these genes were significant for the MATS test of individual-level heterogeneity at an adjusted p-value threshold, and only UGDH yielded an individual-level heterogeneity test p-value below 0.05, suggesting that these findings are driven by population-level expression effect differences. Effect estimates for these genes without adjusting for PC-interactions are given in Table 1, along with the findings from standard and stratified TWAS. We discover opposite signed log odds ratios between the AFA and CAU populations for all five genes, suggesting that the current biological models describing these gene’s role in AD may not be generalized across populations. For example, the estimated log odds ratios (ORs) for AFA and CAU population-specific SNX7 expression from the MATS model are |${\hat{\alpha}}_1=0.0229$| and |${\hat{\alpha}}_1+{\hat{\alpha}}_2=-0.1861$|⁠, respectively. Overexpression of SNX7 has been linked to the reduction of the Aβ peptides whose accumulation is a key driver of AD (34,35), and our findings suggest that such protective mechanisms may not be generalized to individuals of AFA ancestry.

Table 1

Population-specific effects for the top 5 genes with significant population-level effects

Gene	MATS population   heterogeneity p-value	MATS AFA Effect   log OR/SE/p-value	MATS CAU Effect   log OR/SE/p-value	Standard TWAS   log OR/SE/p-value	AFA TWAS Effect   log OR/SE/p-value	CAU TWAS Effect   log OR/SE/p-value
SNX7	0.000283	0.023/0.0456/0.614	−0.186/0.0352/1.28e−07	−0.109/0.0278/9.17e−05	0.0213/0.0467/0.649	−0.179/0.0351/3.13e−07
CD68	0.001460	0.227/0.104/0.0288	−0.181/0.0751/0.016	−0.0416/0.061/0.496	0.224/0.106/0.0355	−0.177/0.0747/0.0175
UGDH	0.002180	0.325/0.138/0.0183	−0.119/0.0448/0.00801	−0.0763/0.0425/0.0728	0.328/0.141/0.02	−0.114/0.0445/0.0105
SUPT4H1	0.002380	−0.0781/0.082/0.341	0.252/0.0711/0.000398	0.112/0.0544/0.0392	−0.0956/0.0837/0.253	0.241/0.0707/0.000642
HYAL3	0.002480	−0.127/0.0542/0.0188	0.1/0.0521/0.0546	−0.0134/0.0363/0.712	−0.141/0.0563/0.0124	0.102/0.0518/0.0498

Given that UGDH reached suggestive significance (p < 0.05) for both population- and individual-level heterogeneity, we fit a full model with three PC-interactions to estimate the population-specific expression effects after accounting for intra-population heterogeneity. Under this model, the estimated effects (as log ORs) of UGDH in the AFA versus CAU population are |${\hat{\alpha}}_1=0.44482$| (SE = 0.14054, p = 0.0015) and |${\hat{\alpha}}_1+{\hat{\alpha}}_2=-0.05745$| (SE = 0.04931, p = 0.2440), respectively, where the standard error (SE) for |${\hat{\alpha}}_1+{\hat{\alpha}}_2$| was computed as|$\sqrt{\mathrm{Var}\big({\hat{\alpha}}_1\big)+\mathrm{Var}\big({\hat{\alpha}}_2\big)+2\mathrm{Cov}\big({\hat{\alpha}}_1,{\hat{\alpha}}_2\ \big)}$| and the p-value corresponds to a subsequent Wald test. The LRT testing of all PC-interactions was significant (p= 1.808* 10⁻⁶), and the p-value for α₂ (the difference between the AFA and CAU expression effect) is p= 0.000412. These findings suggest that the signal detected for UGDH was driven by both inter- and intra-population heterogeneity and that UGDH overexpression may be a risk factor for AD within the AFA population but incurs protective effects within CAU populations.

Application to imaging and hematological traits in UK Biobank

While we leave the details of our UKBB application to the Supplemental Material, we highlight some key findings here. In total, 63 GEx-IDP associations and 216 GEx-Hematological trait pairs were detected by MATS after Bonferroni correction. Only 24/63 and 10/216 of these pairs were also detected by TWAS, respectively, highlighting the improved power of MATS over standard TWAS to detect expression-trait associations. Among the 63 GEx-IDP associations, one pair was significant for individual-level heterogeneity and 12 were significant for population-level heterogeneity, with no overlap. Of particular note, we identified opposite signed effects in the expression of UTS2 on the insular cortex volume between the African and Asian populations (positive effects) versus the European and Middle Eastern populations (negative effects). Among the 216 GEx-Hematological trait pairs, 40 were significant for individual-level heterogeneity. One gene-trait pair, LRMDA on hemoglobin concentration, was significant for both individual- and population-level heterogeneity and remained significant for population-level heterogeneity after adjusting for the top 10 genetic PC-expression interactions, suggesting the presence of both inter- and intra-population heterogeneity.

Discussion

In this work, we presented the Multi-Ancestry TwaS, a novel extension to TWAS which models ancestral heterogeneity in the effects of cis-regulated gene expression on a complex trait. Within the MATS framework, we jointly model samples from multiple subpopulations from trans-ancestry studies, thereby boosting power over stratified TWAS to detect genes with shared expression-trait effects across ancestries. Furthermore, we proposed tests of both population- and subject-level heterogeneity to detect inter- and intra-population variation in GEx-trait effects, enabling MATS to detect association signal that would go undiscovered using standard TWAS, particularly in the most extreme cases where subpopulations have opposite signed expression effects. Our simulations highlight the superior power of MATS over existing methods to detect expression-trait association signals when the expression effects vary between and/or within subpopulations, and further reveal comparable power to TWAS in the absence of heterogeneity. This finding was validated through our applications to both the Alzheimer’s Disease Sequencing Project and UK Biobank. MATS identified many expression-trait associations which were undiscovered by stratified or standard TWAS, including genes with evidence of population- or subject-specific expression-trait associations. Of particular note, MATS confirmed previously supported protective effects of SNX7 overexpression in AD, but discovered possible risk conferring effects within African American populations. Furthermore, our application to the UK Biobank highlights the flexibility of MATS to study expression-trait associations within cohorts of finescale population substructure in which ancestry-matched expression prediction models are unavailable. Given the generalizability of MATS to both continuous and binary complex traits, we aim to raise the attention to more applications of MATS to both existing and future multi-ancestry genetic studies, which will be an important future direction.

This work, namely our discovery of genes with opposite signed expression effects between subpopulations, highlights the dangerous implications of generalizing findings from ancestrally homogeneous genetic studies to other populations. While our novel framework can detect such heterogeneity in the context of expression-trait association testing, the utility of MATS relies on the availability of trans-ancestry genetic studies with sufficiently sized subpopulations. In spite of the recent efforts of major sequencing studies to increase non-European populations (including ADSP (36)), the representation of non-EUR subjects in most major genetic biobanks remains disproportionate, as seen in our application to the UK Biobank (37). Furthermore, the reliability of MATS and standard TWAS are tied to the accuracy of the gene expression prediction models, which to date are trained on only a few hundred samples with available expression data in non-CAU samples. As such, our work emphasizes the need to prioritize ancestral diversity with well-powered subpopulations in large-scale biobanks with sequencing and expression data.

Materials and Methods

The transcriptome wide association study

The TWAS infers associations between genetically regulated gene expression (GEx) and a complex trait. TWAS is equivalent to the method of 2 Stage Least Squares (2SLS), so-called because it is performed in two stages. In the first stage, a prediction model for gene expression (x) is trained on a set of cis-eQTLs, commonly defined as SNPs associated with expression levels that lie within 1 Mb of the gene’s transcription start site. Stage 1 model coefficients |$(\hat{w})$| are then used to obtain fitted values for x (Eq. 1), denoted as |$\hat{x}$| and commonly interpreted as the genetically regulated component of gene expression. We use g_i (p x 1) to denote the vector of p cis-eQTL genotypes for subject i, i ∈ 1, …, n, with each SNP additively coded as 0, 1 or 2 effect alleles.

In Stage 2, the trait (Y) is regressed on imputed gene expression, as given in Eq. (2), where C represents a set of covariates (including the intercept) and h(·) represents some link function. The association between genetically regulated GEx and the trait is tested via divergence from the null hypothesis, H₀: β = 0, in the following model:

$$ \begin{equation} \hat{x}_i = g_{i}^{T} \hat{w},\end{equation} $$

(1)

$$ \begin{equation} h\left\{E\left[{Y}_i|{C}_i,{\hat{x}}_i\right]\right\}={C_i}^T{\beta}_C+{\hat{x}}_i\beta . \end{equation} $$

(2)

Conveniently, many studies have publicly available pre-trained GEx prediction models. This allows for imputation of GEx in cohorts that have available genotype and phenotype data, but not expression data, making it possible to apply TWAS to many large genetic case–control studies. Furthermore, TWAS benefits from improved biological interpretation over standard GEx-trait association tests by modeling GEx in terms of the additive effects from eQTLs (i.e. directly modeling the full pathway from SNPs to GEx to a complex trait). As described by Knutson et al. (2020), TWAS can also overcome potential biases introduced by confounding on the trait-GEx relationship and can thereby take on causal interpretations when certain conditions are met (38).

Modeling heterogeneous gene expression in TWAS

Here, we expand Stage 2 in TWAS to model possible heterogeneous effects of gene expression at both the intra- and inter-population levels. The proposed Stage 2 model is given in Eq. (3) for subject i in population j

$$ \begin{equation} h\{ E[Y_{i,j}|g_{i,j}, C_{i,j}] \} = C_{i,j}^T \beta_C + \hat{x}_{i,j}\beta_{i,j},\end{equation} $$

(3)

$$ \begin{equation} \beta_{i,j} = \alpha_1 + \alpha_j * D_j + \eta_{i,j},\end{equation} $$

(4)

$$ \begin{equation} h\left\{E\left[Y|g,C\right]\right\}=C{\beta}_C+\hat{x}{\alpha}_1+{\sum}_{j=2}^k\ \hat{x}\ast{\alpha}_j{D}_j+\hat{X}\eta, \end{equation} $$

(5)

where |${\hat{x}}_{i,j}$| represents the imputed gene expression for subject i in population j with j ∈ 1, …, k and i ∈ 1,... n_j such that n₁ + · · ·  + n_k = n. C represents a set of measured confounders, including an intercept and indicator variables for population. For now, we assume C also contains some leading genetic PCs to account for population stratification, though we consider this treatment in greater detail in the Supplemental Material.

A key feature of Eq. (3) is that the expression effect β_i,j is indexed at the population and individual levels. Specifically, we decompose β_i,j into three components, as given in Eq. (4): 1) the fixed expression effect for population 1, α₁, treated as the reference group, 2) the difference in the expression effect between the j-th ancestral group compared to the reference population, denoted as α₂, …, α_k, and 3) the subject-specific GEx effect, η_i,j. The final mixed-effects Stage 2 model in Eq. (5) follows from substituting Eq. (4) into Eq. (3), where |$\hat{X}=\mathit{{diag}}({\hat{x}}_1,\dots, {\hat{x}}_n)$| and D_j is an indicator for population j. We refer to the model in Eq. (5) as the Multi-Ancestry TwaS model, or MATS for short.

Under the MATS model, we treat α = [α₁, …, α_k]^T as fixed, and treat η = [η₁, …, η_n]^T as random such that E[η] = 0 and Var[η] = κ²K. Here, K represents an n x n similarity matrix and serves to account for genetic relatedness between subjects. There are a number of choices for K in practice, including a realized genetic relationship matrix (GRM), estimated as the sample covariance matrix computed on a matrix of genome-wide minor allele counts. Misspecification of K has been studied in detail in the related application of mixed modeling to adjust for population stratification in GWAS: if the extent of the misspecification is severe, power may be reduced and type I errors may be inflated. We explore a specific case of ‘extreme’ misspecification in our simulations, which is well known to inflate Type I error: namely, misspecifying K = I, where I is the n × n identity matrix. Notably, the use of K limits this model to individual-level data, though future work may extend the MATS model for summary data applications, similar to the work of (16).

The MATS model in Eq. (5) provides a framework to test for expression-trait associations while accounting for possible ancestral heterogeneity. As to be described in the following section, this MATS expression association test can detect association signals originating from 1) a shared expression effect across populations and subjects (as in standard TWAS), 2) population-specific effects and/or 3) subject-specific or intra-population expression-trait associations. In the sections to follow, we outline a pipeline of tests and corresponding models to identify the source of the association signal and appropriately estimate expression effects for the given setting.

Importantly, a key distinction between the MATS framework and existing models for Multi-Ancestry TWAS is that MATS models a distinct gene expression effect for each ancestral group, providing a way to estimate the effects of gene expression that differ at the population level (as has been observed in real data analysis). In contrast, the existing METRO method for multi-ancestry TWAS accounts for ancestry by weighting the relative contributions of each population’s gene expression model from Stage 1, but assumes a homogeneous gene expression effect across ancestral groups.

Given that MATS builds off of TWAS (or analogously, 2 Stage Least Squares), associations detected under MATS may be causal under certain conditions, particularly when instrumental variable conditions are met. However, the introduction of the random component and a non-linear link function complicates direct causal interpretations. Nevertheless, we highlight a number of MATS discoveries in our real data application that are closely related to the biology of disease, thereby validating the genes that are uniquely discovered by MATS (as compared to TWAS).

A novel test of expression-trait associations in multi-ancestry cohorts

Testing for expression-trait associations under the model in Eq. (5) corresponds to the null hypothesis H_0,1: α₁ = α₂ = · · · = α_k = 0 and κ² = 0, which will be rejected if any of the population-level or ancestry-adjusted subject-specific expression effects differ from zero. To test this hypothesis, we combine the two asymptotically independent tests of H_0,2: α₁ = α₂ = · · · = α_k = 0 given κ² = 0, and H_0,3: κ² = 0. A detailed outline of the score tests for H_0,2 and H_0,3 are given in the Supplemental Material. Briefly, the score statistic under H_0,2 is given as

$$ \begin{equation} {U}_2={S_2}^{\mathrm{T}}\ {\left({\Delta }^{\ast}\right)}^{-1}{S}_2\sim{\chi_k}^2, \end{equation} $$

(6)

with |${S}_2={\big[\hat{x},\hat{x}\odot D\big]}^{\mathrm{T}}\big(Y-{\mu}^{\ast}\big)\;\mathrm{and}\;{\varDelta}^{\ast }= Var\big[{S}_2\big]$|⁠, where |$\odot$| represents element-wise multiplication, D represents the n x k matrix of population indicators (i.e. D = [D₁, …, D_j, …, D_k]) and μ^* = (μ₁^*, … μ_n^*)^T represents the mean estimates under the null model for H_0,2.

The score test for H_0,3, i.e. testing for individual-level heterogeneity, builds off of the work of Goeman et al. (21) and Zhang and Lin (39) and corresponds to the score test statistic

$$ \begin{equation} {U}_3=\frac{1}{n}{\left(Y-\hat{\mu}\right)}^{\mathrm{T}}\hat{X}K{\hat{X}}^{\mathrm{T}}\left(Y-\hat{\mu}\right), \end{equation} $$

(7)

where |${\hat{\mu}}$| represents the mean estimates under H_0,3. We show that U₃ asymptotically follows a mixture of independent χ²₁ distributions, |${\sum}_{i=1}{\lambda}_i{\chi}_1^2$|⁠, where λ_i are eigenvalues of |$\frac{1}{n}{W}^{1/2}\hat{X}\;K\;{\hat{X}}^{\mathrm{T}}\;{W}^{1/2}$| (39). Several methods have been developed to approximate the mixture chi-square distribution to obtain p-values, including Davies (40), Liu et al. (41), Wang (42) and Chen et al. (43). However, these approaches require computing all eigenvalues λ, which can be very computationally expensive even in moderate sample sizes. As such, we approximate p-values using the highly accurate and computationally inexpensive fastSKAT algorithm proposed by Lumley et al. (44).

As outlined by Sun et al., U₂ and U₃ (the score test statistic for H_0,2 and H_0,3, respectively) are asymptotically independent, availing a number of methods for combining independent test statistics to test H_0,1 (45,46). Many of these approaches take the general form of a weighted sum of the independent tests, given in Eq. (8), where p_α and p_κ are the p-values resulting from the tests on the fixed and random components, respectively.

$$ \begin{equation} U_1 = \rho_\alpha[-2log(p_\alpha)] + \rho_\kappa[-2log(p_\kappa)] \end{equation} $$

(8)

The most straightforward approach is Fisher’s method. The test statistic of Fisher’s method is equivalent to an unweighted sum of the –2log(p-value) for the individual components, such that ρ_α = ρ_κ = 1 in Eq. (8). Fisher’s test statistic follows a χ²₄ under H_0,1.

Fisher’s method achieves high power when both |$\alpha \ne 0$| and κ > 0. However, when only one of the two alternatives is true, Fisher’s method loses power. Su et al. (47) propose a data adaptive approach for determining weights ρ_α, ρ_κ. However, their simulations show that Fisher’s approach generally yields power comparable to the adaptively weighted test across different scenarios. As such, we default to the standard Fisher’s method in practice.

Distinguishing between shared, population-specific and subject-specific expression effects

A significant test of H_0,1 suggests that gene expression has some effect on the complex trait, either at the global (shared) level, the population level and/or the individual level. A number of subsequent tests under the MATS framework can help identify the source of the association and inform estimation of either shared or population-specific effects, adjusted for individual-level heterogeneity as needed.

To start, the presence of individual-level heterogeneity can be assessed based on the significance of the score test of H_0,3. If H_0,3 is not rejected, we can operate from a reduced model without the subject-specific random effects, namely

$$ \begin{equation} h\left\{E\left[Y|g,C\right]\right\}=C{\beta}_C+\hat{x}{\alpha}_1+{\sum}_{j=2}^k\hat{x}\ast{\alpha}_j{D}_j. \end{equation} $$

(9)

Under this reduced model, we can test for the presence of heterogeneity in the population-level effects, corresponding to any standard test (e.g. Score, LRT or Wald test) of H_0,4: α₂ = · · · = α_k = 0. We outline the score test for H_0,4 in the Supplemental Material. Divergence from H_0,4 implies at least one population expression effect differs from the reference population. Notably, this null model is equivalent to the standard TWAS approach, and thus rejection of H_0,4 implies improvement of the population-specific MATS model over the standard TWAS approach with ‘pooled’ populations. If H_0,4 is not rejected, the standard TWAS model can be used to estimate the shared expression effect. Notably, due to the increased sample size, we expect this approach to have improved power over stratified TWAS to detect risk genes with shared effects across subpopulations, but may lose power compared to standard TWAS when no population heterogeneity is present due to the greater degrees of freedom of the MATS expression association test.

Fitting a full model in the presence of individual-level heterogeneity

When individual-level heterogeneity is detected (i.e. H_0,3 is rejected), it is appropriate to fit a full model to estimate and test for population-level expression effects which adjusts for individual-level heterogeneity. However, the full model in Eq. (5) cannot be readily estimated due to the difficulty in accurately finding the maximum likelihood estimate of κ², estimation which is avoided under the score test for H_0,3. To address this challenge, we propose a fixed effects approximation to the random individual-level effect in Eq. (5) which can be used to estimate an approximate full model in the presence of individual-level heterogeneity. A detailed outline of this approach is given in the Supplemental Material.

Briefly, our proposed full model replaces the random individual-level effect in Eq. (5) with interactions between leading genetic PCs and imputed genetic expression:

$$ \begin{align} h\left\{E\left[Y|g,C\right]\right\}= \ & C{\beta}_C+\hat{x}{\alpha}_1+{\sum}_{j=2}^k\hat{x}\ast{\alpha}_j{D}_j+\hat{x}\circ \left(\sqrt{\omega_1}{t}_1\right){\gamma}_1\nonumber\\&+\cdots +\hat{x}\circ \left(\sqrt{\omega_{\mathrm{L}}}{t}_L\right){\gamma}_L \end{align} $$

(10)

Here, t_ℓ (n x 1) represents the ℓ-th genetic PC with corresponding eigenvalue ω_ℓ, ℓ ∈ 1, …, L, where L is some arbitrarily selected number of top PCs. The last L terms in Eq. (10) correspond to the interactions between imputed gene expression and each weighted genetic PC, with each interaction effect denoted as γ_ℓ, ℓ ∈ 1, …, L. These interaction terms represent the effect of gene expression on the trait which is correlated with ancestry.

In the Supplemental Material, we show that if interaction effects for all PCs are included (i.e. L = rank(K)) and we treat all interaction effects as random, the corresponding test of heterogeneity under model in Eq. (10) is directly equivalent to the score test for H_0,2. It follows from this equivalence that a model which incorporates a few leading PC-expression interactions as fixed effects is a reasonable approximation to the GLMM MATS model with a random individual-level expression effect. While the ideal number of top PCs is unknown and this approximation may vary in strength depending on the degree of genetic differentiation captured in the top PCs, an LRT of H₀: γ₁ = · · · = γ_L^* = 0 can be used to confirm whether the intra-population variability is captured within the leading L^* PC-interactions. That is, rejection of H₀ suggests that non-ignorable intra-population signal is captured within the leading L^* terms. Less formally, additional PCs can be added and their importance assessed based on whether the results significantly change.

This model facilitates testing of global and population-level effects after adjusting for individual-level heterogeneity using a standard LRT of H^*_0,1: α₁ = · · · = α_k = 0 and H^*_0,4: α₂ = · · · = α_k = 0, respectively. This approach is suitable for fitting a heterogeneity-adjusted full model after rejecting the test for individual-level heterogeneity under the GLMM, which provides more accurate estimates of the population-level effects compared to a model which disregards the presence of intra-population heterogeneity altogether.

Data

GEx prediction models from MESA

We downloaded pre-trained gene expression models published by Mogil et al. (https://github. com/WheelerLab/DivPop) for 4674 genes (9). These models were trained using genetic and purified human monocytes expression data from The Multi-Ethnic Study of Atherosclerosis (MESA), which consisted of 6814 individuals, including 233 AFA, 352 Hispanic (HIS), and 578 CAU subjects. Ancestry data was based on self-report. Models were trained separately on data from each of the three ancestral groups. As detailed in Mogil et al., the expression models between AFA and CAU significantly varied at some genes. As such, we match the prediction model to the target GWAS ancestry (AFA and CAU only) when imputing gene expression in TWAS. Notably, a number of key risk genes in AD are not included in the MESA models (e.g. APOE, CLU, PICALM, TREM2, SORL1), which may have implications in our application to ADSP.

For each gene, an elastic net was trained to predict gene expression from cis-eQTL genotypes, specifically including common variants (MAF > 0.01) within 1 Mb of the transcription start site. Tuning parameters were selected using predictive R² and assessed via nested cross-validation for each ancestry model, where three specific mixing parameters were considered (α = 0.05, 0.5, 1). We limit our analysis to testing a set of 1358 genes with a prediction performance R² > 0.1 in both AFA and CAU populations.

Given that gene expression is cell-type specific, findings from MATS may differ depending on the tissue type of the gene expression model used. Investigators should carefully prioritize gene expression models of relevant cell-types for the trait in consideration. While the current model considers only expression for a single tissue, future extensions of MATS may integrate data from multiple tissues, which has been shown to boost power to detect associated genes for some traits (48).

Notably, the competing METRO method requires individual-level gene expression data or OLS estimates for eQTLs in the multi-ancestry GEx model. Neither of these data is presently available for our real data applications (ADSP or UKB) or for many other sources of multi-ancestry genetic data, thus highlighting another key advantage of MATS over METRO.

ADSP WGS data

Using the MESA GEx prediction models, we imputed gene expression using the individual-level data from the ADSP to test for heterogeneity-adjusted expression-AD associations. ADSP comprises three sub-studies (ADNI, ADSP discovery, and ADSP Extension), including both multiplex family and case–control design. ADSP includes subjects who self-reported as Black or African American (AFA), Non-Hispanic white (CAU), and Hispanic (HIS). However, due to limited sample sizes within the HIS subpopulation, we limited our analysis to the AFA and CAU samples from the case–control studies. Our sample selection procedure is illustrated in Figure 9. In total, we analyzed 4828 AD cases and 4558 controls (n = 9386 combined), each with available Whole Genome Sequencing data measured on the HisSeq X Ten platform. AD cases were identified as having possible, probable or definite AD based on cognitive assessment criteria from the National Institute of Neurological Diseases–Alzheimer’s (NINCDS-ADRDA). We excluded all subjects with uncertain AD statuses, or missing covariates for sex, age at onset (cases) or age at the last exam (controls). A breakdown of cases and controls by ancestry group is given in Table 2.

Figure 9

Sample selection procedure for real data applications.

Table 2

Sample sizes for ADSP case–control data by self reported ancestry

Ancestry	n cases	n controls	n total
Black/African American	1083	1666	2749
Caucasian	3745	2892	6637

To ensure concordance between all reference panels and the MESA prediction models, we performed liftover on the ADSP genotypes from hg38 to hg19. We computed a similarity matrix for all samples after pruning markers using an r² cutoff of 0.2, a window size of 1000 kb and a step size of 50 variants. We computed a genetic relatedness matrix using plink’s –make-rel function and performed PCA on the pruned genotypes using plink’s –pca option. For all MATS and TWAS models, we adjust for biological sex, age at baseline visit and APOE status.

UK Biobank GWAS

The UK Biobank (UKBB) is a rich resource comprising ∼500000 subjects with GWAS and phenotypic measurements. A subsample of 37438 subjects also provided Magnetic Resonance Images, from which the summarized brain imaging-derived phenotypes (IDPs) were computed, including many volumetric regions with previous evidence of causal associations in AD (38,49). While the UK Biobank is composed of primarily European individuals, population substructure and ancestry-linked gene-phenotype associations have been illustrated within UKBB, findings which motivate our application of MATS (50–52). To analyze population level heterogeneity, we derive four subpopulations based on genetically derived ancestry using the TRACE software (53), which included African, Asian, European and Middle Eastern populations. We apply MATS to study the expression associations in two sets of outcome phenotypes: (i) 20 IDPs with putative causal AD associations, analyses that complement our primary ADSP application and (ii) four blood cell traits, including platelet count, hemoglobin, hematocrit and white blood cell count. The latter analyses mirror the work of Wen et al., who studied ancestry correlated gene expression in these hematological traits using stratified TWAS between African, Hispanic and European ancestries (54). The details of these data and analyses are given in the Supplemental Material.

Integrating prediction models with target GWAS

When imputing GEx using the MESA prediction models, we found that many cis-eQTLs in the prediction models were not measured in the GWAS or WGS data, a common challenge in two-sample TWAS analyses. There are a number of common approaches to circumvent this issue, including (i) performing genotype imputation on the target sample to predict missing alleles, (ii) using a variant in high LD with the missing eQTL (i.e. an LD proxy) as a surrogate marker and (iii) excluding missing eQTLs when imputing gene expression. However, when genetic imputation of the target sample yields poorly predicted genotypes at missing markers or no variants in the target sample are in high LD with the missing variant, the first two approaches are not feasible. The latter approach may severely bias imputed GEx and the resulting TWAS when missing variants are strongly predictive of expression or many moderately predictive variants are missing.

Here, we propose an alternative approach for imputing gene expression using all prediction weights when some prediction model eQTLs are missing in the target data. As previously described, imputed gene expression is given as the fitted values from a regression of expression on its corresponding eQTLs, given as |$\hat{x}=E\big[x|G\big]=E\big[{g}^{(1)}{w}_1+\cdots +{g}^{(j)}{w}_j+{g}^{\big(j+1\big)}{w}_{j+1}+\cdots +{g}^{(p)}{w}_p|g\big]$|⁠, where g^(j) represents n × 1 vector of genotypes for the jth eQTL, j ∈ 1, …, p, w_j represents the estimated effect of the j-th eQTL and G represents the n x p matrix of genotypes, G = (g₁, …, g_p). Suppose SNPs 1, …, j are all missing from the target sample. Define the n × j matrix of missing genotypes as g_miss = (g₁, …, g_j) with weight vector |${\hat{w}}_{miss}$| and the n × (p–j) genotypes measured in the GWAS data as g_gwas = (g_j + 1, …, g_p) with weight vector |${\hat{w}}_{gwas}$|⁠, such that G = (g_miss,  g_gwas). Then, we can rewrite the stage 1 fitted values as

$$\begin{eqnarray}& \hat{x} = E[x|G] = E[g_{miss} \hat{w}_{miss}|G] + E[g_{gwas} \hat{w}_{gwas}|G] \nonumber\\ & \qquad \quad = E[g_{miss}|G]\hat{w}_{miss} + g_{gwas}\hat{w}_{gwas}. \end{eqnarray}$$

(11)

The unknown quantity in Eq. (11), E[g_miss|G], can be estimated by performing a multivariate regression of the genotypes g^* on G using a reference panel of similar ancestry, where we henceforth use the superscript * to denote genotypes from an external reference. That is, we obtain fitted values from the following regression:

$$ \begin{equation*} {\mathbf{g}}^{\ast}_{miss}={\mathbf{g}}^{\ast}_{gwas}\ \boldsymbol{\eta} +\epsilon^{\ast} .\end{equation*} $$

Reference panels such as the 1000 Genomes publicly report whole genome sequencing data on multiple populations, thereby making it likely that all prediction model variants are present in the reference data, and estimation of E[g_miss|G] can be performed separately using reference panel samples matched to the self-reported ancestry of target sample subjects.

Simulation setup

To evaluate the performance of our proposed test, we perform simulations under a variety of settings based on real data. We simulate genotype data using the sim1000G in R (https://github.com/adimitromanolakis/sim1000G), a genotype simulator that has been shown to successfully preserve ancestry-specific MAFs and LD patterns from a reference sample in simulated data (55). We simulate genotypes based on 50 cis-SNPs from a randomly selected gene (TOMM22 on chromosome 22) using the Phase III 1000 Genomes reference panel. For our primary simulations, we separately generate data on 2000 subjects for each of the EUR, AFR and AMR populations (total n = 6000).

We specify true Stage 1 GEx weights (i.e. marginal SNP effects on GEx) by sampling from a w_j ∼ N (0, 0.01). We allow 80% of these weights to be shared across all three subpopulations, and randomly sample the remaining 20% separately for each ancestry group to reflect the expected variability in gene expression prediction models across subpopulations (9). We generate gene expression from x = Gw + β_U U + ϵ_X, where ϵ_X,i is independently sampled from a N(0, 0.01) and U represents ‘unmeasured’ confounders randomly sampled from N(0, 0.1) with effect size β_U,X = 0.01. We then generate a continuous phenotype as |$Y={\sum}_{r=1}^{10}{c}_rP{C}_r+{\alpha}_1x+{\alpha}_2x\ast{D}_2+{\alpha}_3x\ast{D}_3+{\beta}_{U,Y}U+\mathbf{X}\eta +{\in}_Y$|⁠, where D_j represents an indicator for the jth population, ϵ_Y ∼ N (0, σ_ϵ) represents the random error and β_U,Y represents the effect of unmeasured confounders on Y, which we set to 0.01. Here, PC_r represents the rth genetic PC with effect c_r, sampled randomly from c_r ∼ N (0, 0.05). We sample the expression variance parameter, η, from MVN(0, κ²K), where K represents the genetic relatedness matrix based on the simulated genotypes.

To assess the performance of our proposed MATS model against standard TWAS, we first estimate weights and impute gene expression for each population separately, using the full samples for estimation (i.e. one sample TWAS). We assess the rejection rates for H_0,4 (population-level heterogeneity), H_0,3 (intra-population heterogeneity), H_0,1 (overall gene expression effect) against PC-adjusted TWAS. We vary the parameters α₁, α₂ and α₃ across the values [0, 0.001, 0.005, 0.01, 0.02, 0.03, 0.05], and vary κ² across [1e − 10, 5e − 04, 0.001, 0.00125, 0.0015, 0.00175, 0.002, 0.0025, 0.0035, 0.005, 0.01].

As it currently stands, most multi-ancestry GWAS or sequencing studies have imbalanced sample sizes between self-reported ancestries. To reflect this situation, we consider a supplementary simulation setting with n₁ = 4000, n₂ = 500 and n₃ = 500. Furthermore, while the Stage 1 sample sizes for our real data applications are small (n_AFA = 233, n_CAU = 578), recent work suggests that the cis-heritability of gene expression and Stage 2 sample size have substantially greater influences over the power of the Stage 2 TWAS test than does the Stage 1 sample size (56). This suggests our simulations reasonably reflect the expected power of our proposed Stage 2 MATS model in real data, and that the small Stage 1 sample size may not overly limit the power of our applied results (which are limited to heritable genes).

Code availability

MATS can be applied in R using the package available at https://github.com/kathalexknuts/MATS. This package integrates R and C++, along with a number of existing packages to improve computational efficiency, including bigQF (https://github.com/tslumley/bigQF) for p-value approximation in the MATS test of individual-level heterogeneity. For a total sample size of 6000, the suite of MATS tests took only 2.919 seconds to run.

Acknowledgements

Access to the ADSP and UKBB data was approved through the NIAGADS DAR ID #10141 and UKB Application #35107, respectively.

Conflict of Interest statement

The authors have no conflicts of interest to declare.

Funding

National Institutes of Health grants R01AG065636, R01AG069895, RF1 AG067924, U01AG073079, R01HL116720 and R01GM126002; Minnesota Supercomputing Institute at the University of Minnesota.