Cross-Population Joint Analysis of eQTLs: Fine Mapping and Functional Annotation
[['Expression quantitative trait loci , or eQTLs , are genetic variants that are associated with gene expression levels .', 'Mapping eQTLs can help in dissecting the molecular mechanisms by which genetic variants impact organismal phenotypes .', 'Recent studies [1–3] have revealed that there are substantial overlaps between eQTLs and genetic variants identified from genome-wide association studies ( GWAS ) of disease phenotypes .', 'In addition , eQTL mapping provides a powerful tool for investigating the regulatory machinery in different tissues [4 , 5] or cellular environments [6–8] .', 'In this paper , we jointly address three outstanding issues in eQTL mapping .', 'First , due to the high experimental cost , most available eQTL data sets typically have limited sample sizes .', 'To improve power of eQTL discovery , it becomes necessary to aggregate evidence across multiple data sets .', 'Second , because a gene is typically regulated by many regulatory elements , it is highly likely that there exist multiple independent eQTLs in its proximity ( i . e . , cis region ) .', 'In this scenario , a multi-SNP analysis is required to uncover all relevant cis acting genetic factors involved in the gene regulation process [9] .', 'Third , the availability of extensive functional annotations [10–12] now enables integration of functional genomic information into eQTL analysis , which can be useful to dissect the functional basis of eQTLsf .', 'Linking genomic annotations to eQTLs goes beyond genetic association analysis , and helps gain a better understanding of the underlying biological processes .', 'Individually , some of these three issues have been discussed by previous works .', 'For example , [3 , 9 , 13–16] discussed single SNP analysis of eQTLs jointly from different studies , populations or tissues .', 'But these methods do not naturally extend to multi-SNP analysis .', '[17–20] examined the enrichment of selected genomic features in cis-eQTLs , mostly based on single SNP association results .', 'To the best of our knowledge , there is no existing approach that jointly addresses all three issues in a systematic way .', 'In this paper , we demonstrate an integrative analysis framework to perform multi-SNP fine mapping analysis of eQTLs using cross-population samples .', 'Our statistical methods stem from an established Bayesian framework proposed by [14 , 15 , 21] , which has been successfully applied in mapping eQTLs from multiple tissues .', 'We apply our statistical framework to analyze the data from the GEUVADIS project [20] , where the expression-genotype data are collected from five population groups .', 'In GWAS , trans-ethnic meta-analysis of genetic association data from diverse populations has been shown to be a powerful tool in detecting novel complex trait loci and improving resolution of fine mapping of causal variants by leveraging population heterogeneity in local patterns of linkage disequilibrium ( LD ) and allele frequencies [22 , 23] .', 'This approach , to the extent of our knowledge , has not been applied to eQTL analysis .', 'Utilizing cross-population expression-genotype data , we are interested in identifying eQTL signals that show consistent effects in all populations .', 'Furthermore , we aim to examine whether we have sufficient statistical power to identify multiple independent cis-eQTL signals with the available aggregated sample size .', 'Finally , we set out to investigate whether the genetic variants that disrupt transcription factor ( TF ) binding are enriched in eQTLs , and we integrate such annotations to further improve fine mapping analysis of eQTLs .', 'Our three main aims are also inherently inter-related and reinforce each other .', 'With higher power through sample aggregation in mapping eQTLs , our method can identify at a high resolution potentially multiple genomic regions that harbor casual eQTLs .', 'Consequently , these efforts improve the statistical power and precision of localization for our functional analysis .'], ['We start with a brief description of our statistical framework and general strategy for multi-SNP fine mapping analysis in a meta-analytic setting across multiple populations .', 'We performed a series of simulation studies to evaluate the proposed Bayesian multi-SNP fine mapping procedures .', 'Specifically , we simulated under the setting where eQTL data are collected from multiple population groups .', 'This section briefly summarizes the main simulation results , and we leave the relevant details in the Materials and Methods section .', 'First , we examined the power of the proposed Bayesian multi-SNP analysis approach in identifying multiple independent association signals using cross-population samples .', 'In particular , we compared its performance with the commonly applied single SNP meta-analysis procedure and a ( step-wise ) conditional meta-analysis procedure [25 , 26] , treating each population group as a participating study in a meta-analysis .', 'To simulate the genotype-expression data in multiple population groups , we used the real genotypes of 2 , 500 SNPs from 100 randomly selected and distinct genomic regions ( i . e . , 25 consecutive SNPs per region ) in the GEUVADIS data .', 'As a consequence , there are generally modest to high levels of LD within each region , but unremarkable levels of LD between regions .', '( In other words , we assembled 100 relatively independent LD blocks for each gene . )', 'It is worth noting that within each region , the variation in patterns of LD between populations is maintained .', 'We randomly assigned one to four eQTL SNPs for each gene and simulated their expression levels according to a system of linear models ( see Material and Methods section and Section S . 5 of S1 Text ) .', 'For evaluation , we assessed the ability of each procedure in correctly identifying the 25-SNP regions ( i . e . , the LD blocks ) that harbor the true eQTLs .', 'The simulation results indicate that the proposed Bayesian multi-SNP analysis approach is the most powerful in identifying independent eQTL signals among the three methods compared ( Fig 1 ) .', 'Further investigation reveals that all three methods achieve similar power when genes contain only a single eQTL; however when a gene harbors multiple independent cis-eQTLs , the proposed Bayesian fine mapping approach becomes highly advantageous ( S1 Fig ) .', 'In addition , the simulation results show that the reported PIPs from the proposed Bayesian approach are well-calibrated ( S2 Fig ) .', 'We also performed additional simulations to evaluate the proposed computational procedures for enrichment testing and fine mapping analysis integrating genomic features .', 'In brief , we found that the proposed procedures have greatly improved power in testing enrichment of genomic features when compared with standard approaches based on single SNP testing results ( S2 Table ) .', 'In addition , we observed that the proposed computational procedure provides accurate estimates of the enrichment parameters ( S3 Fig ) .', 'In this paper , we focused on analyzing the expression and genotype data collected from the GEUVADIS project [20] .', 'More specifically , the data set consists of RNA-seq data on lymphoblastoid cell line ( LCL ) samples from five populations: the Yoruba ( YRI ) , CEPH ( CEU ) , Toscani ( TSI ) , British ( GBR ) and Finns ( FIN ) .', 'In our analysis , we selected 420 samples which were densely genotyped in the 1000 Genomes Phase I data release [27] and 11 , 838 protein coding genes and lincRNAs that are deemed expressed in all five population groups .', 'Throughout , our analysis focused on the SNPs that locate within a 200kb genomic region centered at the transcription start site ( TSS ) of each gene .', 'In contrast to the original eQTL mapping approach discussed in [20] , we treated each population as a single group and performed cis-eQTL analysis jointly across all five groups .', 'We have made the full analysis results available on the website http://www-personal . umich . edu/~xwen/geuvadis/ .'], ['We have presented an integrative analysis framework to perform fine mapping analysis jointly from cross-population samples while incorporating functional annotations .', 'Our core statistical methods are built on and naturally extend the established Bayesian framework for association analysis of genetic data from heterogeneous groups [14 , 15 , 21] .', 'The most notable features of our statistical framework are its ability of multi-SNP analysis and integration of SNP-level genomic features .', 'Through simulation studies , we have demonstrated the power and efficiency of the proposed approaches .', 'For the first time , we have applied this framework for multiple cis-eQTL fine mapping analysis in a cross-population meta-analytic setting using the GEUVADIS data .', 'Our analysis reveals that with an aggregated sample size of around 400 , multiple independent cis-eQTL signals can be confidently identified in many genes , which makes evident the necessity of considering multiple association signals in eQTL studies .', 'The commonly applied strategy for mapping additional association signals is to conduct conditional analysis , which can be viewed as a step-wise variable selection algorithm .', 'Besides lack of power , one additional disadvantage of conditional analysis is that it only reports a single “best” model in the end , and completely ignores its uncertainty .', 'From many of our examples shown in this paper , it is clear that in most cases there is typically a great deal of uncertainty in determining the truly causal SNPs due to LD , and the posterior probabilities of the “best” models are often unimpressive .', 'Consequently , it is generally inappropriate to solely rely on the information from the “best” model in the downstream analysis ( as we have illustrated in the power simulations of enrichment testing ) .', 'In comparison , in addition to gaining power in identifying multiple independent cis-eQTL signals , our Bayesian approach provides much more comprehensive information that fully conveys the uncertainty of the inference result , and the quantified uncertainty information is naturally propagated in our integrative analysis of genomic features .', 'In this paper , we have employed a two-step procedure that first screens eGenes by performing gene-level hypothesis testing , and then carries out multi-SNP analysis for the identified eGenes .', 'This procedure is analogous to the fine mapping procedure that is commonly used in GWAS , where interesting loci are ranked and selected by single SNP association testing before an in-depth analysis focusing on each flagged high priority locus .', 'We find this procedure not only yields considerable computational savings , but also provides a sound argument to specify the prior inclusion probability Pr ( γj ) .', 'It should be noted that our proposed analysis procedures are completely applicable for fine mapping of general QTLs in either a meta-analytic or single study setting .', 'They can be further extended to the applications of fine mapping analysis where subgroups of eQTL data are formed by different tissues [4 , 14] or cellular conditions [6 , 8] .', 'Comparing to the meta-analytic setting considered in this paper , the parameter space of {γj} in those applications is more complicated ( which includes 2s potential values , where s is the number of subgroups/tissues ) .', 'Nevertheless , [14] has provided a principled way to “learn” the priors on possible values that γj can take by pooling information across genes through a hierarchical model .', 'Our fine mapping results of eQTLs also demonstrate the benefit of utilizing cross-population samples in genetic association studies .', 'Most importantly , the population heterogeneity of local LD patterns serves as an efficient filter that narrows down the regions harboring casual eQTLs .', 'Nevertheless , varying LD patterns can cause some SNPs to display large degree of heterogeneity across populations in their estimated effect sizes from single SNP analysis: in the extreme cases , a SNP may appear to possess strong “population specific” effects .', 'As we acknowledge that genuine population specific eQTLs are certainly interesting phenomena and very much likely exist , we suggest interpreting highly heterogeneous eQTL signals from single SNP analysis with caution .', 'In our view , it may be necessary to carry out multi-SNP analysis , as we have demonstrated in this paper , to simply rule out the possibility that the seemingly population specific effects are artifacts due to varying LD patterns .', 'Our Bayesian inference framework naturally incorporates functional annotations in fine mapping eQTLs across population groups .', 'This feature allows us to quantitatively assess the enrichment of certain functional features in eQTLs , and in turn to use the quantified enrichment information to prioritize annotated SNPs for fine mapping analysis .', 'Overall , our model for integrative eQTL mapping analysis is similar to those presented in [17 , 18] .', 'However , these previous approaches make simplifying assumptions to restrict at most one cis-eQTL per gene , such that single SNP association results can be directly used .', 'Our method relaxes this assumption and is fully integrated into our multi-SNP analysis procedure .', 'In addition , our use of CENTIPEDE annotation to examine the relative enrichment of binding variants and footprint SNPs is also novel .', 'Although it is largely expected that binding variants are enriched in eQTLs , it is important to note that the level of enrichment for footprint SNPs is much lower than that for binding variants .', 'Interestingly , this finding seems concurring with the results reported by [30] where the relative enrichment of binding variants vs . footprint SNPs in other cellular and organismal phenotype QTLs is examined .', 'Last but not least , our multi-SNP fine mapping analysis of the GEUVADIS data has created a comprehensive resource for the community to gain better understanding of the genetic basis of gene regulation .', 'With proper uncertainty assessments , our results enable follow-up experimental validations and functional studies of causal genetic variants that alter gene regulation .', 'They also provide a unique and powerful resource to study population genetics of expression traits .', 'For example , we found that the distribution of Fst of eQTL signals consistent across populations has a shorter tail than the distributions of random SNPs ( Kolmogorov-Smirnov test p-value = 0 . 001 ) .', 'Nevertheless , when examining the Fst distributions between primary and secondary cis-eQTL signals ( stratified by the association strength in single SNP analysis with respect to their target genes ) , we did not find statistical evidence to differentiate the two ( Kolmogorov-Smirnov test p-value = 0 . 77 ) .', 'There are many interesting aspects to be explored with our analysis results , and we will leave the more in-depth investigations for our future work .'], ['The genotype and RNA-seq data were downloaded from the GEUVADIS project website .', 'We selected 420 samples whose genotypes are directly measured in the 1000 Genomes project .', 'The samples are evenly distributed in the five population groups , and the detailed breakdown of samples by population is shown in Table 1 .', 'For the RNA-seq data , we used a slightly more stringent threshold than the original analysis [20] to select genes that are expressed in all five populations .', 'Specifically , for each selected gene , we required > 90% individuals in each population group to have RPKM ≥ 0 . 1 .', 'From the 17 , 361 Ensembl genes that passed this filter , following the original analysis , we selected a subset of 11 , 838 genes consisting of annotated protein-coding genes and lincRNAs according to GENCODE [31] release 17 .', 'We log transformed the RPKM values , and used the pipeline employed in [29 , 32] to remove the effect of GC content on expression measurements .', 'We then followed the same strategy as described in [20] to remove latent confounding factors using the software package PEER [33] .', 'However , unlike the original analysis , we ran PEER for each population group separately .', 'In the end , we removed 15 , 13 , 15 , 20 and 20 PEER factors for samples from YRI , CEU , TSI , GBR and FIN , respectively .', 'Finally , the expression levels of each gene were quantile normalized across individuals separately in each population group .', 'For genotype data , we filtered out SNPs whose sample allele frequencies < 0 . 03 in the overall samples across population groups .', 'Note , we did not apply the allele frequency filter in each population group .', 'The SNPs passing this filter must have sample allele frequencies ≥ 0 . 03 in at least one population group .', 'In general , as discussed in the Results section , rare SNPs do not impose any statistical or computational problems for our analysis .', 'Nevertheless , removing SNPs that are not likely informative in any population group helps improve computational efficiency .', 'Following [17 , 32] , we defined a 200kb cis region for each gene centered at its TSS .', 'In total , the final data set contains 6 . 7 million gene-SNP pairs .', 'For each gene , we tested the null hypothesis that asserts no cis-eQTLs .', 'Specifically , we adopted the Bayesian hypothesis testing procedure discussed in [14] .', 'Essentially , [14] assumes a Bayesian model that is mostly similar to ( Eq 1 ) , except for an additional simplifying assumption , “at most one cis-eQTL per gene” .', 'Given a gene with p cis-SNPs , this additional assumption reduces the possible alternative scenarios into p single SNP association models , for which a gene-level Bayes factor can be easily computed analytically by Bayesian model averaging .', 'In the context of eQTL mapping in multiple tissues , [14] considered all possible configurations , i . e . , γj values , for each assumed associated SNP , whereas in our analysis of the GEUVADIS data we only allowed γj ∈ {0 , 1} to focus on discovering population consistent cis-eQTLs .', 'Briefly , for the alternative model where the j-th SNP is assumed the lone eQTL , the linear model ( Eq 1 ) is simplified to', 'y i = μ i 1 + β i , j g i , j + e i , e i ∼ N ( 0 , σ i 2 I ) , i = 1 , … , s .', '( 7 )', 'We modeled the correlation of genetic effects , βi , j’s across population groups through the following prior specification ,', 'β i , j ∼ N ( β ¯ j , ϕ 2 ) , β ¯ j ∼ N ( 0 , ω 2 ) .', '( 8 )', 'Equivalently , the above prior can be represented by a multivariate normal distribution by integrating out the average effect parameter β ‾ j , i . e . ,', 'β j = ( β 1 , j ⋮ β s , j ) ∼ N ( 0 , W ) , ( 9 )', 'where the s × s variance-covariance matrix W is given by', 'W = ( ϕ 2 + ω 2 ω 2 ⋯ ω 2 ω 2 ϕ 2 + ω 2 ⋯ ω 2 ⋮ ⋮ ⋱ ⋮ ω 2 ω 2 ⋯ ϕ 2 + ω 2 ) .', 'The parameter ϕ2 + ω2 characterizes the overall genetic effects of SNP j , and ϕ2/ ( ϕ2 + ω2 ) represents the degree of heterogeneity across population groups .', 'Following [14 , 15] we considered the values of ϕ2 + ω2 from a set E = {ϕ2 + ω2 : 0 . 12 , 0 . 22 , 0 . 42 , 0 . 82 , 1 . 62} which covers a wide range of plausible magnitude of genetic effects .', 'We allowed limited degree of heterogeneity by taking ϕ2/ ( ϕ2 + ω2 ) values from the set H = {ϕ2/ ( ϕ2 + ω2 ) : 0 , 0 . 1} which reflects our prior belief that effects of genuine eQTL signals should be highly consistent across population groups .', 'Overall , we considered a combination of |E| × |H| grid for ( ϕ2 , ω2 ) values for each alternative model .', 'Given this model , a single SNP Bayes factor , BFj , can be analytically evaluated following [14 , 15] , and the corresponding gene-level Bayes factor is obtained by Bayesian model averaging as', 'BF gene = 1 p ∑ j = 1 p BF j .', '( 10 )', 'Upon obtaining the gene-level Bayes factors , we used the methods implemented in the software package eQTLBMA [14] to select eGenes at FDR 5% levels .', 'eQTLBMA implements two types of FDR control procedures: one is a permutation based procedure which converts a gene-level Bayes factor to a corresponding p-value and control FDR using Storey’s procedure; the other procedure is based on the EBF procedure discussed in [34] which directly works with gene-level Bayes factors and avoids any permutations .', 'We found that the latter approach is much more computationally efficient , however slightly conservative .', 'The results reported in the Results section are based on the EBF procedure .', 'The power gain of the Bayesian gene-level testing procedure in identifying consistent association signals across multiple populations/subgroups and its comparison with commonly applied frequentist approaches have been fully demonstrated in [14] and [21] through simulations and real data analysis .', 'Here , we examined a few known factors that may potentially affect the power of the gene-level analysis in the GUEVADIS data .', 'Using the hierarchical model proposed in [14] , we performed the analysis described in [21] to estimate the effect size heterogeneity of eQTL signals across populations: our choice of the heterogeneity parameters is mostly informed by this analysis .', 'More specially , we used a full spectrum of grid values for the heterogeneity parameter ϕ2/ ( ϕ2 + ω2 ) from a comprehensive set H = {ϕ2/ ( ϕ2 + ω2 ) : 0 . 0 , 0 . 1 , 0 . 2 , … , 0 . 9 , 1 . 0} .', 'We then estimated the probability weight on each grid value by pooling information across all genes using the EM algorithm implemented in eQTLBMA .', 'The estimates show that the majority of the probability mass is concentrated in very low heterogeneity levels ( S1 Table ) , and the mean heterogeneity E[ϕ2/ ( ϕ2 + ω2 ) |Y , G] = 0 . 060 , i . e . , on average , the eQTL effects are very consistent across populations .', 'This result is also highly concordant to the findings in [21] where eQTL data from microarray experiments across European , African and Asian populations were examined , and further justifies our selection of heterogeneity parameter values in the analysis .', 'Among other factors impacting eGene discoveries , we noted that the number of cis-SNPs is negatively correlated with the gene-level local false discovery rate ( lfdr ) .', 'That is , on average , genes harboring more cis-SNPs tend to exhibit stronger gene-level association signals .', 'The Spearman’s rank correlation between the two quantities is modest ( -0 . 20 ) , nonetheless statistically highly significant ( p-value < 2 . 2 × 10−16 ) .', 'In our multi-SNP analysis , we no longer assumed “one cis-eQTL per gene” , and considered the full range of alternative scenarios described by model ( Eq 1 ) .', 'To make joint inference with respect to Γ = {γi , … , γp} , we further specified effect size distribution ,', 'β j | γ j = 1 , ϕ 2 , ω 2 ∼ N ( 0 , W ) , ( 11 )', 'where W is constructed in the same way as in the gene-level analysis .', 'In particular , we used the same set of grid values of ( ϕ2 , ω2 ) in the gene-level analysis .', 'Unconditional on γj , the prior on βj is a type of “spike-and-slab” , where the “slab” is represented by a mixture of multivariate normal distribution , and the vast majority of the prior probability mass , 1 − 1 p , is assigned to the “spike” ( i . e . , a point mass at 0 ) .', 'The described linear model system including the prior specification ( i . e . , SSLR ) , is a special case of the general system considered by [15] .', 'Given a specified Γ value , a Bayes factor , BF ( Γ ) , contrasting to the trivial null model , Γ ≡ 0 , can be analytically approximated by applying the result discussed therein .', 'With this result , the posterior probability , Pr ( Γ|Y , G ) , can be computed up to an unknown normalizing constant , i . e . , Pr ( Γ|Y , G ) ∝ Pr ( Γ ) BF ( Γ ) .', 'We implemented a Metropolis-Hastings algorithm , similar to what is discussed in [15] for multivariate linear regression model ( MVLR ) , to efficiently traverse the space of 2p possible Γ values .', 'In particular , we designed a novel proposal distribution that utilizes marginal and conditional analysis results to prioritize SNPs with strong marginal or conditional association signals .', 'In practice , we observed that the resulting MCMC algorithm achieves fast mixing .', 'The details of the algorithm is provided in S1 Text ( section S . 1 ) .', 'In the end , the MCMC algorithm yields a set of Γ samples from the posterior distribution , from which we computed the PIP for each SNP by marginalization .', 'The posterior expected number of independent cis-eQTLs and its variance for each gene are obtained by computing the sample mean and variance of the number of non-zero γj’s in each posterior model .', 'Equivalently , the posterior expected number of cis-eQTLs can be computed by the sum of PIPs , i . e . ,', 'E [ ∑ j = 1 p 1 ( γ j ≠ 0 ) | Y , G ] = ∑ j Pr ( γ j = 1 | Y , G ) .', '( 12 )', 'For the fine mapping analysis of the GEUVADIS data , we applied the MCMC algorithm for each identified eGene individually .', 'We carried out 25 , 000 burnin steps and 50 , 000 repeats for each MCMC run .', 'Taking advantage of parallel processing , we performed multiple-SNP analysis for multiple genes simultaneously in a distributed computing environment , which greatly reduced the overall computational time .', 'The MCMC run for the eGene with the most cis-SNPs , HLA-DRB1 ( Ensembl ID: ENSG00000196126 , with 11 , 400 cis-SNPs ) , took 30 minutes on a computer with a single Intel Xeon 2 . 13GHz 8-core CPU .', 'For average eGenes with ∼ 2 , 000 cis-SNPs the running time is approximately 3 to 4 minutes on the same machine .', 'We ran our fine mapping analysis for GEUVADIS data with 8 parallel threads , the overall computation took about 30 hours .', 'Although the MCMC output conveys the full information of the fine mapping results , parsing the information manually to identify independent eQTL signal clusters can still be non-trivial .', 'To tackle this problem , we designed and implemented a simple hierarchical clustering based algorithm that automatically parses the MCMC output and aides identifying independent eQTL signal clusters .', 'The details of the algorithm are described in the S1 Text ( section S . 4 ) .', 'To test the enrichment of a particular genomic annotation in eQTL signals , we examine the associations between SNP level annotations of interest and the PIPs from our multi-SNP fine mapping approach computed under the null hypothesis of no enrichment .', 'Without loss of generality , we consider a single annotation and use δ j g to denote the annotated value of SNP j in gene g ( the additional super-script for gene emphasizes the annotation is specific to each gene-SNP pair ) .', 'Our testing procedure starts by computing PIPs for all gene-SNP pairs using the proposed multi-SNP cis-eQTL mapping method with the exchangeable prior specification ( Eq 4 ) .', 'Under the null hypothesis , the resulting PIPs should be independent of the annotation information .', 'We then fit a logistic regression model', 'logit [ Pr ( γ j g | Y , G ) ] = α 0 + α 1 δ j g , ( 13 )', 'and perform the Wald test with respect to α1 .', 'The commonly applied enrichment testing approach typically classifies the latent association status of each SNP ( γ j g ) into 0 or 1 depending on the single SNP association test results ( usually p-values ) , and then test their associations with the annotations of interest .', 'In comparison , the use of PIPs in our proposed method presents at least two obvious advantages: first , it naturally accounts for potentially multiple association signals within a single gene; and second , it fully accounts for the uncertainty in determining γ j g which is ignored by the binary classification approach .', 'These advantages are reflected by the gain of powers in enrichment testing .', 'To incorporate genomic annotations into the fine mapping analysis of cis-eQTLs , we specify the prior distribution of Pr ( γj ) by the parametric function ( Eq 6 ) , whereas the other parts of the Bayesian multi-SNP analysis model remain intact .', 'For the g-th gene and its j-th SNP , we re-write ( Eq 6 ) using an equivalent vector form as follows', 'logit [ Pr ( γ j g = 1 ) ] = α ′ δ j g , ( 14 )', 'where the enrichment parameter α is assumed to be shared across all gene-SNP pairs .', 'Let Dg denote the collection of the annotation data for gene g .', 'For a total of q genes , it follows from the Bayes rule that', 'Pr ( α , Γ 1 , … , Γ q | { Y 1 , G 1 , D 1 } … { Y q , G q , D q } ) ∝ P ( α ) ∏ g = 1 q [ Pr ( Γ g | D g , α ) BF ( Γ g ) ] .', '( 15 )', 'Given a prior specification P ( α ) ( e . g . a flat prior ) , the quantities on the right hand side are individually straightforward to compute analytically .', 'It is conceptually easy to modify our MCMC algorithm for the multi-SNP analysis to jointly sample ( α , Γ1 , … , Γq ) .', 'Nevertheless , due to the extremely high dimensionality of the target space ( which is approximately the number of gene-SNP pairs , in the case of the GEUVADIS data the number is ∼ 6 . 7 million ) , the convergence of the MCMC within a reasonable time frame may be in doubt .', 'Furthermore , the computational resources , especially the memory usage , demanded by the MCMC algorithm may be too high to afford in a practical setting .', 'Alternatively , we consider an empirical Bayes approach which employs an EM algorithm to find the MLE of α by treating ( Γ1 , … , Γq ) as missing data , and perform a final round of multi-SNP analysis ( possibly only for eGenes ) conditional on the resulting MLE α ^ .', 'The derivation of the EM algorithm is mostly straightforward , we give the relevant details in S1 Text ( section S . 2 ) .', 'Briefly at ( t+1 ) -th iteration , in the Expectation step ( E-step ) , we evaluate E [ γ j g | Y g , G g , α ( t ) ] for each gene-SNP pair given the current estimate α ( t ) .', 'It should be noted that this conditional expectation is exactly the PIP of each gene-SNP pair which can be obtained by running the multi-SNP analysis for each gene separately given the hyperparameter α ( t ) .', 'In the Maximization step ( M-step ) , we find α ( t+1 ) by fitting a logistic regression model relating PIP of each gene-SNP pair to the genomic annotations of the corresponding SNP .', 'Overall , we describe the complete algorithm as an “MCMC-within-EM” algorithm , which is initiated at some arbitrary α ( 0 ) value , and iteratively performs multiple cis-eQTL mapping using the MCMC algorithm and maximization by fitting logistic regression models until convergence .', 'The main computational benefit of the “MCMC-within-EM” algorithm is that the E-steps involving MCMC runs can be processed in parallel on a distributed computing system , hence the memory requirement is much relaxed .', 'To illustrate the procedure for estimating α , we simulated a set of eQTL data using the same scheme in the enrichment testing .', 'In particular , we used a single binary annotation and set α1 = 1 . 50 .', 'We ran the MCMC-within-EM algorithm on this simulated data set with the starting point α1 = 0 , and plotted the estimated α1 values in each iteration in S3 Fig . The figure indicates that the estimates of α1 quickly converge to the close neighborhood of the true value after 6 to 7 iterations .', 'To reduce computation , in practice we analyzed the GEUVADIS data using only a single iteration of EM where all enrichment parameters are initialized at 0 .', 'As a consequence , our reported point estimates of the enrichment parameters are likely to be lower bounds of the true values .', 'Nevertheless , as we have shown in the analysis of the GEUVADIS data , this simplified version of the algorithm still improves the resolutions of the fine mapping analysis by up-weighting the SNPs with relevant annotations .']]
['Introduction', 'Results', 'Discussion', 'Materials and Methods']
['Mapping expression quantitative trait loci ( eQTLs ) has been shown as a powerful tool to uncover the genetic underpinnings of many complex traits at molecular level .', 'In this paper , we present an integrative analysis approach that leverages eQTL data collected from multiple population groups .', 'In particular , our approach effectively identifies multiple independent cis-eQTL signals that are consistent across populations , accounting for population heterogeneity in allele frequencies and linkage disequilibrium patterns .', 'Furthermore , by integrating genomic annotations , our analysis framework enables high-resolution functional analysis of eQTLs .', 'We applied our statistical approach to analyze the GEUVADIS data consisting of samples from five population groups .', 'From this analysis , we concluded that', 'i ) jointly analysis across population groups greatly improves the power of eQTL discovery and the resolution of fine mapping of causal eQTL', 'ii ) many genes harbor multiple independent eQTLs in their cis regions', 'iii ) genetic variants that disrupt transcription factor binding are significantly enriched in eQTLs ( p-value = 4 . 93 × 10-22 ) .']
['Expression quantitative trait loci ( eQTLs ) are genetic variants associated with gene expression phenotypes .', 'Mapping eQTLs enables us to study the genetic basis of gene expression variation across individuals .', 'In this study , we introduce a statistical framework for analyzing genotype-expression data collected from multiple population groups .', 'We show that our approach is particularly effective in identifying multiple independent eQTL signals that are consistently presented across populations in the proximity of a gene .', 'In addition , our analysis framework allows effective integration of genomic annotations into eQTL analysis , which is helpful in dissecting the functional basis of eQTLs .']
[]