GRAT: Speeding up Expression Quantitative Trail Loci (eQTL) Studies

Fig. 1. An example of applying GRAT in two hypothetical regions. First, the proxy SNP (rectangle) is tested and its statistics is compared to the threshold (dashed line). If the statistic is above the threshold, the remaining SNPs in the region are tested.

Fig. 1. An example of applying GRAT in two hypothetical regions. First, the proxy SNP (rectangle) is tested and its statistics is compared to the threshold (dashed line). If the statistic is above the threshold, the remaining SNPs in the region are tested.

Emrah Kostem in our group has recently published a new method for association analysis targeted toward eQTL analysis called GRAT(10.1007/978-3-642-37195-0_10).  GRAT is designed to speed up eQTL studies and the software is available at

Over the past few years, the genome wide association study (GWAS) approach has been applied to identify regions of the genome, which harbor genetic variation that affects gene expression levels. These regions are referred to as expression quantitative trait loci (eQTL)(10.1038/nrg2537).(10.1038/nrg1964). In a typical eQTL study, the GWAS approach is applied to tens of thousands of gene expression levels using millions of SNPs, resulting in billions of association statistics to be computed. This results in a tremendous computational burden, which is only increasing with sequencing technology collecting more genetic variations and high-throughput genomic data collecting more phenotypic data such as isoform expression(10.1016/j.tig.2010.10.006). This problem is compounded by the fact that some of the statistical techniques for analyzing eQTLs utilize mixed models and themselves are computationally expensive(10.1038/ng.548),(10.1038/nmeth.1681).(10.1038/ng.2310).

We recently published a paper on a method, GRAT, to perform association analysis in high-throughput phenotype datasets, such as the eQTL studies.
The key idea behind GRAT is that we first test a subset of the SNPs and only in regions where the statistic is above a threshold, we test the remaining regions. In contrast to testing all SNPs, our approach tests around 10% of the SNPs in two-stages and guarantees to identify all significant associations with a very high accuracy.

Here is a description of our method from the paper:

Genome-Wide Rapid Association Testing (GRAT)
In Figure 1, we consider two possible scenarios for a genomic region in a GWAS. In (a) the region contains no significant associations and in (b) the region con- tains a causal SNP. In (a) and (b), the statistics for each SNP are shown, denoting what could have been observed in each scenario had all the SNPs in the region been tested. Let m2 be the proxy SNP for this region to decide whether or not to test the rest of the SNPs. We refer to the SNPs other than the proxy SNP ( m1, m3, m4, m5, m6 and m7 ) as the “remainder SNPs”. If the observed statistic of the proxy SNP is stronger than a threshold value, which in this example is 3.0, the remainder SNPs are tested.
In the first-stage, only the proxy SNP is tested and its association statistic is observed. In (a), where the region contains no associations, the statistic of the proxy SNP is 0.7. The observed statistic of the proxy is less than the threshold value ( 0.7 < 3.0 ) and hence none of the remainder SNPs within the region are tested. In (b), the region contains associations and the proxy SNP captures this information. The observed statistic of the proxy SNP is stronger than the thresh- old value ( 5.0 > 3.0 ), which leads to testing each of the remainder SNPs in the region. This results in identifying all the significant SNPs ( m3, m4 and m5 ).

In the paper, we introduce a novel approach for choosing the proxy SNPs and the threshold values, which provide guarantees that all statistically significant associations will be discovered while computing the least amount of association tests. Due to the complexity of linkage disequilibrium (LD) across the genome, we use a separate threshold value for each remainder SNP rather than using a common threshold value for all the remainders SNPs in an LD region. This is performed by pairing each remainder SNP with its most strongly correlated proxy SNP and a threshold value is used for the pair to decide whether or not to test the remainder SNP. We have precomputed the proxy SNPs for the 1000 Genomes Project and studies imputing to SNPs in this reference can benefit from our method. Even though the LD structure among the SNPs in the study and the reference dataset may be different, our method guarantees to discover all significant associations with high-probability. This is achieved by updating the threshold values using the LD structure observed in the study. We term our novel two-stage testing procedure as Genome-wide Rapid Association Testing (GRAT).

GRAT can be applied to a wide range of statistical models, such as case- control studies, quantitative traits and linear mixed models (LMM). In particu- lar, the LMM approach has recently become popular due to its effective control of population structure. Computing the LMM association statistic is compu- tationally expensive and recently its efficient computation has attracted great interest(10.1038/ng.548),(10.1038/nmeth.1681).(10.1038/ng.2310). The speed-up due to GRAT is cumulative with these efforts.

There are some interesting aspects to the computational method. For given proxy SNPs, our method solves an optimization problem that minimizes the number of SNPs tested that will result in a given recall rate. We prove that this problem is convex and show that it can be solved very efficiently.
Furthermore, we propose a greedy algorithm to search for the optimal proxy SNPs.

The full citation is:

Kostem, Emrah; Eskin, Eleazar

Efficiently Identifying Significant Associations in Genome-Wide Association Studies Conference

Research in Computational Molecular Biology, University of California Springer Berlin Heidelberg, 2013.

Abstract | Links | BibTeX


Correcting Population Structure using Mixed Models Webcast

Golden Helix yesterday hosted a excellent webcast on correcting population structure in association studies using mixed models and they highlighted our EMMA(10.1534/genetics.107.080101) and EMMAX(10.1038/ng.548) algorithms.  The webcast was given by Greta Peterson and is available at  The webcast is a great overview of mixed models applied to population structure in general as well as specifically how to use the Golden Helix software to use mixed models in association studies.

A interesting aspect of the story is that we found out about the webcast from an email advertising that they will cover the EMMAX algorithm.  It turns out that there were 863 people registered for the webcast which surpassed their previous record (for a webcast on NGS) by almost 100!  It is exciting to see how much interest there is in mixed models and in our EMMA paper which we published in 2008.

On our website, we have a bunch of resources for mixed models including the EMMA, EMMAX and ICE softwares.  We recently posted an overview of mixed models here.  Below is a list of our papers related to mixed models.


Sul, Jae Hoon; Martin, Lana S; Eskin, Eleazar

Population structure in genetic studies: Confounding factors and mixed models. Journal Article

In: PLoS Genet, 14 (12), pp. e1007309, 2018, ISSN: 1553-7404.

Abstract | Links | BibTeX


Sul, Jae Hoon; Bilow, Michael; Yang, Wen-Yun Y; Kostem, Emrah; Furlotte, Nick; He, Dan; Eskin, Eleazar

Accounting for Population Structure in Gene-by-Environment Interactions in Genome-Wide Association Studies Using Mixed Models. Journal Article

In: PLoS Genet, 12 (3), pp. e1005849, 2016, ISSN: 1553-7404.

Abstract | Links | BibTeX


Joo, Jong Wha J; Hormozdiari, Farhad; Han, Buhm; Eskin, Eleazar

Multiple testing correction in linear mixed models. Journal Article

In: Genome Biol, 17 (1), pp. 62, 2016, ISSN: 1474-760X.

Abstract | Links | BibTeX


Joo, Jong Wha J; Kang, Eun Yong; Org, Elin; Furlotte, Nick; Parks, Brian; Hormozdiari, Farhad; Lusis, Aldons J; Eskin, Eleazar

Efficient and Accurate Multiple-Phenotype Regression Method for High Dimensional Data Considering Population Structure. Journal Article

In: Genetics, 204 (4), pp. 1379-1390, 2016, ISSN: 1943-2631.

Abstract | Links | BibTeX


Schweiger, Regev; Kaufman, Shachar; Laaksonen, Reijo; Kleber, Marcus E; März, Winfried; Eskin, Eleazar; Rosset, Saharon; Halperin, Eran

Fast and Accurate Construction of Confidence Intervals for Heritability. Journal Article

In: Am J Hum Genet, 98 (6), pp. 1181-92, 2016, ISSN: 1537-6605.

Abstract | Links | BibTeX


Furlotte, Nicholas A; Eskin, Eleazar

Efficient Multiple Trait Association and Estimation of Genetic Correlation Using the Matrix-Variate Linear Mixed-Model. Journal Article

In: Genetics, 200 (1), pp. 59-68, 2015, ISSN: 1943-2631.

Abstract | Links | BibTeX


Kostem, Emrah; Eskin, Eleazar

Improving the accuracy and efficiency of partitioning heritability into the contributions of genomic regions. Journal Article

In: Am J Hum Genet, 92 (4), pp. 558-64, 2013, ISSN: 1537-6605.

Abstract | Links | BibTeX


Sul, Jae Hoon; Eskin, Eleazar

Mixed models can correct for population structure for genomic regions under selection. Journal Article

In: Nat Rev Genet, 2013, ISSN: 1471-0064.

Links | BibTeX


Sul, Jae Hoon; Han, Buhm ; Ye, Chun ; Choi, Ted ; Eskin, Eleazar

Effectively Identifying eQTLs from Multiple Tissues by Combining Mixed Model and Meta-analytic Approaches Journal Article

In: PLoS Genet, 9 (6), pp. e1003491, 2013, ISSN: 1553-7404.

Abstract | Links | BibTeX


Listgarten, Jennifer; Lippert, Christoph ; Kadie, Carl M; Davidson, Robert I; Eskin, Eleazar ; Heckerman, David

Improved linear mixed models for genome-wide association studies. Journal Article

In: Nat Methods, 9 (6), pp. 525-6, 2012, ISSN: 1548-7105.

Links | BibTeX


Kenny, Eimear E; Kim, Minseung ; Gusev, Alexander ; Lowe, Jennifer K; Salit, Jacqueline ; Smith, Gustav J; Kovvali, Sirisha ; Kang, Hyun Min ; Newton-Cheh, Christopher ; Daly, Mark J; Stoffel, Markus ; Altshuler, David M; Friedman, Jeffrey M; Eskin, Eleazar ; Breslow, Jan L; Pe'er, Itsik

Increased power of mixed models facilitates association mapping of 10 loci for metabolic traits in an isolated population. Journal Article

In: Hum Mol Genet, 20 (4), pp. 827-39, 2010, ISSN: 1460-2083.

Abstract | Links | BibTeX


Kang, Hyun Min; Sul, Jae Hoon ; Service, Susan K; Zaitlen, Noah A; Kong, Sit-Yee Y; Freimer, Nelson B; Sabatti, Chiara ; Eskin, Eleazar

Variance component model to account for sample structure in genome-wide association studies. Journal Article

In: Nat Genet, 42 (4), pp. 348-54, 2010, ISSN: 1546-1718.

Abstract | Links | BibTeX


Kang, Hyun Min; Ye, Chun ; Eskin, Eleazar

Accurate discovery of expression quantitative trait loci under confounding from spurious and genuine regulatory hotspots. Journal Article

In: Genetics, 180 (4), pp. 1909-25, 2008, ISSN: 0016-6731.

Abstract | Links | BibTeX


Mixed models can correct for population structure for genomic regions under selection

Genome-wide association studies (GWAS) collect people with a disease (called  “cases”) and people without a disease (called “controls”) and compare allele frequencies between cases and controls to identify genomic locations associated the disease. An underlying assumption of GWAS is that cases and controls are sampled from the same population. If they are not, then a phenomenon called “population structure” may cause spurious associations. Correcting for population structure in GWAS has been a very important problem in model organism such as mouse and in human genetics.

In 2010, our group proposed a method called “EMMAX” (10.1038/ng.548) that uses a linear mixed model to correct for population structure in human GWAS. EMMAX computes the relationship between every pair of individuals from SNP data (called “kinship matrix”) and uses this kinship matrix to control population structure. We showed that our method removes effects of population structure better than previous methods using two human GWAS datasets. However, Price et al. showed in this paper (10.1038/nrg2813) that EMMAX may be susceptible to spurious associations for genomic regions under selection; these are regions where two populations have significantly different allele frequencies.

We investigated this issue further and found that by using an appropriate kinship matrix (or matrices), EMMAX can correct for population structure for genomic regions under selection. We showed in the paper that by computing the kinship matrix only from SNPs whose allele frequencies are very different between two populations, we can successfully remove effects of population structure. We also proposed using two kinship matrices; one kinship computed from SNPs under selection and the other kinship from the rest of SNPs. This also correctly controls population structure. Lastly, we looked at whether SNPs under selection actually cause this problem in two human GWAS datasets, but did not identify the problem in both datasets.

Full Citation:

Sul, Jae Hoon, and Eleazar Eskin. 2013. Mixed models can correct for population structure for genomic regions under selection. Nature Reviews Genetics 14, no. 4 (February 26): 300–300.