## Abstract

High-throughput genomics allows genome-wide quantification of gene expression levels in tissues and cell types and, when combined with sequence variation data, permits the identification of genetic control points of expression (expression QTL or eQTL). Clusters of eQTL influenced by single genetic polymorphisms can inform on hotspots of regulation of pathways and networks, although very few hotspots have been robustly detected, replicated, or experimentally verified. Here we present a novel modeling strategy to estimate the propensity of a genetic marker to influence several expression traits at the same time, based on a hierarchical formulation of related regressions. We implement this hierarchical regression model in a Bayesian framework using a stochastic search algorithm, HESS, that efficiently probes sparse subsets of genetic markers in a high-dimensional data matrix to identify hotspots and to pinpoint the individual genetic effects (eQTL). Simulating complex regulatory scenarios, we demonstrate that our method outperforms current state-of-the-art approaches, in particular when the number of transcripts is large. We also illustrate the applicability of HESS to diverse real-case data sets, in mouse and human genetic settings, and show that it provides new insights into regulatory hotspots that were not detected by conventional methods. The results suggest that the combination of our modeling strategy and algorithmic implementation provides significant advantages for the identification of functional eQTL hotspots, revealing key regulators underlying pathways.

THE current focus of biological research has turned to high-throughput genomics, which encompasses large-scale data generation and a variety of integrated approaches that combine two or more -omics of data sets. An important example of integrative genomics analysis is the investigation of the genetic regulation of transcription, also called expression quantitative trait locus (eQTL) or “genetical genomics” studies (Cookson *et al.* 2009; Majewski and Pastinen 2011). A typical eQTL analysis follows a natural structure of parallel regressions between the large set of *q* responses (*i.e.*, expression phenotypes), and that of *p* explanatory variables (*i.e.*, genetic markers, often single nucleotide polymorphism, SNPs), where *p* is typically much larger than the number of observations *n*.

From a statistical point of view, the size and the complex multidimensional structure of eQTL data sets pose a significant challenge. Not only does one wish to detect the set of important genetic control points for each response (expression phenotype), including *cis*- and *trans*-acting control for the same transcript, but, ideally, one would wish to exploit the dependence between multiple expression phenotypes. This will facilitate the discovery of key regulatory markers, so-called hotspots (Breitling *et al.* 2008), *i.e.*, genetic loci or single polymorphisms that influence a large number of transcripts. Identification of hotspots can inform on network and pathways, which are likely to be controlled by major regulators or transcription factors (Yvert *et al.* 2003; Wu *et al.* 2008). Most importantly, there is mounting evidence that common diseases may be caused (or modulated) by changes at a few regulatory control points of the system (*i.e.*, hotspots), which can cause perturbations with large phenotypic effects (Chen *et al.* 2008; Schadt 2009).

In this article, we set out to perform hotspot and eQTL detection in an efficient manner, which exploits fully the multidimensional dependencies within the genome-wide gene expression and genetic data sets. We build upon our previous work (Bottolo and Richardson 2010), where we implemented Bayesian sparse multivariate regression for continuous response to search over the possible subsets of predictors in the large 2* ^{p}* model space. For each expression phenotype, this corresponds to carrying out multipoint mapping within an inference framework, Bayesian variable selection (BVS), where model uncertainty is fully integrated. Here, we propose a novel structure for linking parallel multivariate regressions that borrows information in a hierarchical manner between the phenotypes to highlight the hotspots. To be precise, we propose a new multiplicative decomposition of the joint matrix of selection probabilities ω

*that link marker*

_{kj}*j*to phenotype

*k*and demonstrate in a simulation study that this hierarchical structure and its Bayesian implementation (hierarchical evolutionary stochastic search or HESS algorithm) possess good characteristics in terms of sensitivity and specificity, outperforming current methods for hotspot and eQTL detection. Finally, we show the applicability of our method in two real-case eQTL studies, including animal models and human data. Our approach is broadly applicable and extendable to other high-dimensional genomic data sets and represents a first step toward a more reliable identification of functional eQTL hotspots and the underlying causal regulators.

Analysis models for eQTL data are linked to two strands of work: (i) methods for multiple mapping of QTL, where a single continuous response, referred to as a “trait,” is linked to DNA variation at multiple genetic loci by using a multivariate regression approach, and (ii) models that exploit the pattern of dependence between the sets of responses associated with a predictor (*i.e.*, genetic marker). There is a vast literature on multi-mapping QTL (see the review by Yi and Shriner 2008); some of the models have been extended to the analysis of a small number of traits simultaneously (Banerjee *et al.* 2008; Xu *et al.* 2008). Several styles of approaches have been adopted ranging from adaptive shrinkage (Yi and Xu 2008; Sun *et al.* 2010) to variable selection within a composite model space framework that sets an upper bound on the number of effects (Yi *et al.* 2007). Most of the implemented algorithms sample the regression coefficients via Gibbs sampling. However, these have not been used with a substantial set of markers in the “large *p* small *n*” paradigm, but mostly in case of candidate genes or in small experimental cross-animal studies. To face the challenges typical of larger eQTL studies, we have chosen to build our multi-mapping model using a recently developed Bayesian sparse regression approach (Bottolo and Richardson 2010). In this approach, subset selection is implemented in an efficient way for vast (potentially multi-modes) model space by integrating out the regression coefficients and by using a purposely designed MCMC variable selection algorithm that enhances the model search with ideas and moves inspired by evolutionary Monte Carlo algorithms.

The first eQTL modeling approach that explicitly set out to borrow information from all the transcripts was proposed by Kendziorski *et al.* (2006). In the mixture over markers (MOM) method, each response (expression phenotype) **y*** _{k}*, 1 ≤

*k*≤

*q*, is potentially linked to the marker

*j*with probability

*p*or not linked to any marker with probability

_{j}*p*

_{0}. All responses linked to the marker

*j*are then assumed to follow a common distribution

*f*(⋅), borrowing strength from each other, while those of nonmapping transcripts have distribution

_{j}*f*

_{0}. Inspired by models that have been successful for finding differential expression, the marginal distribution of the data for each response

**y**

*is thus given by a mixture model: . A basic assumption of the MOM model is that a response is associated with at most one predictor. For good identifiability of the mixture, MOM requires a sufficient number of transcripts to be associated with the markers. The authors use the EM algorithm to fit the mixture model and estimate the posterior probability of mapping nowhere or to any of the*

_{k}*p*locations. By combining information across the responses, MOM is more powerful and can achieve a better control of false discovery rates (FDR) by thresholding the posterior probabilities than pure univariate differential expression methods testing each transcript-marker pair. But as originally developed, it is not fully multivariate as it does not account for multiple effects of several markers on each expression trait.

To improve on identification of eQTL effects, Jia and Xu (2007) formulate a unifying *q* × *p* hierarchical model in which each transcript **y*** _{k}*, 1 ≤

*k*≤

*q*, is potentially linked to the complete set of

*p*markers

**X**through a full linear model with regression coefficients, β

*= (β*

_{k}

_{k}_{1},..., β

*,..., β*

_{kj}*)*

_{kp}^{T}. Inspired by Bayesian shrinkage approaches already used in conventional QTL mapping, they propose using a mixture prior on each of the β

*, also known as “spike and slab,” , (1)with a fixed very small δ for the spike and an independent prior for the variance of the slab in the*

_{kj}*j*th marker. They then link the

*q*responses through a hierarchical model of the Bernoulli indicators γ

*, establishing what we refer to as a hierarchical regression set-up. They assume that γ*

_{kj}*∼ Bernoulli(ω*

_{kj}*), 1 ≤*

_{j}*k*, ≤

*q*, and give ω

*a Dirichlet(1, 1) prior. In this model, to improve detection of transcript-marker associations, strength is borrowed across all the transcripts via the common latent probability ω*

_{j}*. Jia and Xu (2007) implement their hierarchical model in a fully Bayesian framework using an MCMC algorithm called BAYES, based solely on Gibbs sampling.*

_{j}The high dimensionality of both gene expression and marker space has been alternatively addressed through the use of data reduction methods. In particular, Chun and Keleş (2009) have proposed implementing sparse partial least-squares regression (M-SPLS eQTL) on preclustered group of transcripts. M-SPLS selects markers associated with each transcript cluster by evaluating the loadings on a set of latent components. As the dimension of each cluster is moderate, SPLS implements a multivariate formulation that takes into account the correlation between the transcripts in the same cluster. Sparsity of the latent direction vectors is achieved by imposing a combination of *L*_{1} and *L*_{2} penalties, similar to the elastic net. The tuning parameters *K* and η controlling the number of latent components and the convexity of the penalized likelihood are tuned together by cross-validation. The output of this method is the set of the regression coefficients of the markers belonging to the latent vectors that are significantly associated with a subset of transcripts, selected by bootstrap confidence interval.

Jia and Xu’s linked regression set-up and fully Bayesian formulation is a natural starting point for eQTL detection, which shares common features with our approach. Here, we present a novel model structure and state of the art implementation based upon evolutionary Monte Carlo. We report the results of a simulation study comparing our HESS method to BAYES (Jia and Xu 2007), as well as to two alternative approaches: MOM (Kendziorski *et al.* 2006) and M-SPLS (Chun and Keleş 2009). Finally, we show the application of our method to two diverse genomic experiments in mouse and human genetic contexts.

## Theory and Methods

### Hierarchical related sparse regression

Let the *n* × *q* matrix of responses, with **y*** _{k}* = (

*y*

_{k}_{1},…,

*y*)

_{kn}^{T}the sequence of

*n*observations of the

*k*th response, and let

**X**be the

*n*×

*p*design matrix with

*i*th row

**x**

*= (*

_{i}*x*

_{i}_{1},…,

*x*)

_{ip}^{T}. We assume throughout that

**x**

*is quantitative. It encompasses the case of continuous biomarkers, inbred genotypes {0, 1} for recombinant inbred (RI) strains and {0, 1, 2} genotype coding for F*

_{i}_{2}animal crosses or human data. A linear model for the

*k*th response can be described by the equationwhere α is an unknown constant,

**1**

*is a column vector of ones,*

_{n}**β**

*= (β*

_{k}

_{k}_{1},…, β

*)*

_{kp}^{T}is the

*p*× 1 vector of regression coefficients, and ε

*is the error term with , where*

_{k}**I**

*is the diagonal matrix of dimension*

_{n}*n*. BVS is performed by placing a constant prior density on α

*and a prior on*

_{k}**β**

*, which depends on a latent binary vector γ*

_{k}*= (γ*

_{k}

_{k}_{1},…, γ

_{k}_{j},…, γ

*)*

_{kp}^{T}, where γ

_{k}_{j}= 1 if β

*≠ 0 and γ*

_{kj}

_{k}_{j}= 0if β

*= 0,*

_{kj}*j*= 1,…,

*p*. Conditionally on the latent binary vector, the linear model becomeswhere

**β**

_{γ}

*is the nonzero vector of coefficients extracted from*

_{k}**β**

*,*

_{k}**X**

_{γ}

*is the design matrix of dimension*

_{k}*n*×

*p*

_{γ}

*, with columns corresponding to γ*

_{k}*= 1, and the number of selected covariates for the*

_{kj}*k*response. For every regression

*k*, we assume that, apart from the intercept α

*,*

_{k}**X**contains no variables that would be included in every possible model and that the columns of the design matrix have all been centered in 0.

Assuming independence of the *q* regression equations conditionally on the selected predictors modeled in the *q* × *p* latent binary matrix , the likelihood becomes (2)where φ* _{n}*(⋅) is the

*n*-variate normal density function.

The description of the joint likelihood as the product of *q* regression equations is similar to the one proposed by Jia and Xu (2007). However, one important difference is the assignment in (2) of a regression specific error variance , allowing for transcript-related residual heterogeneity and making our formulation more flexible. A more general model, seemingly unrelated regressions (SUR) introduces additional dependence between the responses **Y** through the noise ε* _{k}*, modeling the correlation between the residuals of different responses (Banerjee

*et al.*2008). However, it becomes computational unfeasible when the size of

*q*is large, which is typical in eQTL experiments.

### Prior set-up

For a given *k*, we follow the same prior set-up for the regression coefficients and error variance as described in Bottolo and Richardson (2010). First, we treat the intercept α* _{k}* separately, assigning it a constant prior,

*p*(α

*) ∝ 1. Second, conditionally on γ*

_{k}*, we assign a*

_{k}*g*-prior structure on the regression coefficients and an inverse-gamma (Inv Ga) density to the residual variance (3), (4)with

*a*

_{σ},

*b*

_{σ}> 0, and . This conjugate prior set-up has many advantages. The most important is that, for a given

*k*, the marginal likelihood

*p*(

**y**

*|*

_{k}**X**, γ

*, τ) can be written in a closed form that is particularly simple to compute once (3) and (4) are integrated out. Furthermore, it allows for more efficient MCMC exploration with correlated predictors than the nonconjugate case (*

_{k}*i.e.*, when the variance component in (3) is different from the error variance) and it provides more accurate identification of the high-probability models among those visited during the MCMC (George and McCulloch 1997). Finally it leads to a simple and interpretable expression, with the ordinary least-squares solution, of the level of shrinkage.

The hierarchical structure on the regression coefficients is completed by specifying a hyper-prior on the scaling coefficient τ, *p*(τ). We adopt the Zellner–Siow priors structure for the regression coefficients that can be thought as a scale mixture of *g*-priors and an inverse-gamma prior on τ, *p*(τ) = Inv Ga(1/2, *n*/2) with heavier tails than the normal distribution, . In general it has been observed (Bottolo and Richardson 2010) that data adaptivity of the degree of shrinkage conforms better to different variable selection scenarios than assuming standard fixed values (which can be easily implemented by using a point mass prior for τ). Since the level of shrinkage can influence the results of the variable selection procedure, in our model we force all the *q* regression equations to share the same common τ, linking the regression equations hierarchically through the variance of their nonzero coefficients.

The prior specification is concluded by assigning a Bernoulli prior on the latent binary value γ* _{kj}*,

*p*(γ

*|ω*

_{kj}*) = Bernoulli (ω*

_{kj}*). The prior chosen for ω*

_{kj}*is of paramount importance in BVS since it controls the level of sparsity,*

_{kj}*i.e.*, the association with a parsimonious set of important predictors. For a given response this task can be accomplished by specifying a common small-selection probability for all

*p*predictors, ω

*= ω*

_{kj}*and giving*

_{k}*p*(ω

*) = Beta(*

_{k}*a*,

_{k}*b*) (Bottolo and Richardson 2010). Inducing sparsity when all the responses are jointly considered is harder because further constraints are desirable. eQTL surveys (Cookson

_{k}*et al.*2009) suggest that only a fraction of expression traits are under genetic regulation and the number of their control points is usually small. This can be modeled by assigning a different probability for each marker ω

*= ω*

_{kj}*with an hyper-prior on ω*

_{j}*. This solution, first proposed by Jia and Xu (2007) with the conjugate prior*

_{j}*p*(ω

*) = Dirichlet(*

_{j}*d*

_{1}

*,*

_{j}*d*

_{2}

*), assumes that this selection probability is the same for all the responses. However, whatever the sensible choice of the hyperparameters*

_{j}*d*

_{1}

*and*

_{j}*d*

_{2}

_{j}_{,}

*d*

_{1}

*,*

_{j}*d*

_{2}

*= 1 or*

_{j}*d*

_{1}

*,*

_{j}*d*

_{2}

*= 0.5, the posterior density greatly depends on the ratio between the number of transcripts associated to the marker*

_{j}*j*,

*q*

_{j}_{,}and the total number the transcripts in the eQTL experiment,

*q*, since

*E*(ω

*|*

_{j}**) = (**

*Y**q*+

_{j}*d*

_{1}

*)/(*

_{j}*q*+

*d*

_{1}

*+*

_{j}*d*

_{2}

*), where*

_{j}*q*= # {

_{j}*j*: γ

*= 1}. In such formulation, the results are thus clearly influenced by the number of responses analyzed and sparsity of each*

_{kj}*k*th regression cannot be controlled in the prior specification adopted for ω

*of γ*

_{kj}

_{kj}_{.}

In this article we propose a novel way of specifying the selection probability ω* _{kj}* to synthetize the benefits of both approaches, Bottolo and Richardson (2010) and Jia and Xu (2007). We propose decomposing this probability into its marginal effects(5)with ω

*and ρ*

_{k}*the “row” and “column” effect, respectively, and 0 ≤ ω*

_{j}*≤ 1 and ρ*

_{k}*≥ 0, but constrained so that 0 ≤ ω*

_{j}_{j}

*≤ 1. The idea behind this decomposition is to control the level of sparsity for each row*

_{k}*k*through a suitable choice of the hyperparameters

*a*,

_{k}*b*of

_{k}*p*(ω

*) = Beta(*

_{k}*a*,

_{k}*b*), while the parameter ρ

_{k}*captures the “relative propensity” of predictor*

_{j}*j*to influence several responses at the same time. Large values of ρ

*indicate that predictor*

_{j}*j*has a marked influence on ω

*and thus likely to be a hotspot. The adopted multiplicative formulation has some similarity to the disease mapping paradigm where the relative risk level acts in a multiplicative fashion on an expected number of cases in a binomial or Poisson disease risk model. A gamma density on the*

_{jk}*j*th latent hotspot effect,

*p*(ρ

*) = Ga(*

_{j}*c*,

*d*), with

*E*(ρ

*) =*

_{j}*c*/

*d*, complete the hierarchical structure for the decomposition (5).

We conclude this section by describing the choice of the hyperparameters for ω* _{k}* and ρ

*. Since by construction ω*

_{j}*⊥ ρ*

_{k}*,*

_{j}*E*(ω

*) =*

_{jk}*E*(ω

*)*

_{k}*E*(ρ

*). If we assume*

_{j}*c*=

*d*, the hotspot propensity does not change the

*a priori*row marginal expectation,

*E*(ω

*) =*

_{jk}*E*(ω

*). However, it inflates the*

_{k}*a priori*row marginal variance Var(ω

*) > Var(ω*

_{jk}*), with Var(ω*

_{k}*) = Var(ω*

_{jk}*)(1 +*

_{k}*d*

^{−1}) +

*d*

^{−1}

*E*

^{2}(ω

*). For the specification of the hyperparameters*

_{k}*a*and

_{k}*b*, we use the Beta-binomial approach illustrated in Kohn

_{k}*et al.*(2001), after marginalizing over the column effect in (5). The two hyperparameters can be worked out once

*E*(

*p*

_{γ}

*) and Var(*

_{k}*p*

_{γ}

*), the expected number and the variance of the number of genetic control points for each response, are specified.*

_{k}### Posterior inference

After integrating out the intercepts, the regression coefficients and the error variances, the joint density can be factorized as (6)where , , and . Posterior inference is carried out on the *q* × *p* latent binary matrix **Γ** = {γ* _{kj}*, 1 ≤

*k*≤

*q*, 1 ≤

*j*≤

*p*}, on the

*q*×

*p*selection probability matrix

**Ω**= {ω

*, 1 ≤*

_{kj}*k*≤

*q*, 1 ≤

*j*≤

*p*}, and on the scaling coefficient τ, if not fixed.

Sampling **Γ** is extremely challenging since complex dependence structures in the **X** create well-known problems of multimodality of the model space even for a single regression equation. Here the computational challenge is higher since we are aiming to explore a huge model space of dimension (2* ^{p}*)

*. For this reason vanilla MCMC (MC*

^{q}^{3}, Gibbs sampler, simple dimension changing moves) cannot guarantee a reliable exploration of the model space in a limited number of iterations. In this article we use a sampling scheme introduced by Bottolo and Richardson (2010), evolutionary stochastic search (ESS) as a building block for our new algorithm HESS. For each response, HESS relies on running multiple chains with different “temperature” in parallel, which exchange information about the set of covariates that are selected in each chain. Since chains with higher temperatures flatten the posterior density, global moves (between chains) allow the algorithm to jump from one local mode to another. Local moves (within chains) permit the fine exploration of alternative models, resulting in a combined algorithm that ensures that the chains mix efficiently and do not become trapped in local modes. Specific modifications of ESS were introduced to comply with the structure of ω

*, which are sampled with the decomposition (5) rather than integrating them out as in ESS. This requires some modifications in the local moves (details in Supporting Information, File S1, Section S.2.2.)*

_{kj}For sampling the selection probability matrix **Ω**, we implemented a Metropolis-within-Gibbs algorithm for each element of the row effect **ω** = (ω_{1}, …, ω* _{q}*)

^{T}and column effects

**ρ**= (ρ

_{1}, …, ρ

*)*

_{q}^{T}, rejecting proposed values outside the range [0, 1]. However, since the dimension of

**ω**and

**ρ**is very large, tuning the proposal for each element of the two vectors is prohibitive. To make HESS fully automatic, we use the adaptive MCMC scheme proposed by Roberts and Rosenthal (2009), where the variance of the proposal density is tuned during the MCMC to reach a specified acceptance rate. To satisfy the asymptotic convergence of the adaptive MCMC scheme, mild conditions are imposed (details in File S1, Section S.2.3).

If not fixed, the scaling coefficient τ, which is common for all the *q* regression equations and all the *L* chains, is sampled using a Metropolis-within-Gibbs algorithm with random walk proposal and fixed proposal variance (details in File S1, Section S.2.4).

Finally, we describe a complete sweep of our algorithm. We assume that the design matrix is fully known. If missing values are present, these can be imputed in a preprocessing imputation step (for instance using the fill.geno function from the *qtl* R package for genetic crosses (Broman and Sen 2009) or *FastPhase* (Scheet and Stephens 2006) for human data). Without loss of generality, we assume that the responses and the design matrix have both been centered. The same notation is used when τ is fixed or given a prior distribution. For simplicity of notation we do not index variables by the chain index, but we emphasize that the description below applies to each chain:

Given

**Ω**and τ we update γ, according to the ESS procedure, using global and local moves. During the burn-in, we sample the latent binary vector γ_{k}for each_{k}*k*to tune the regression specific temperature ladder (details in File S1, Section S.2.5). After the burn-in, at each sweep, we select at random without replacement a fraction φ of the regressions where to update γ._{k}Given

**Γ**and τ, we sample**ω**and**ρ**with a random walk Metropolis with adaptive proposals.Given

**Γ**and**Ω**, we sample τ with a random walk Metropolis with a fixed proposal. To balance the number of updates of the latent binary values γwith those of the scaling coefficient, at each sweep, the number of times we sample τ is proportional to_{kj}*q*×*p*×*L*.

The Matlab implementation of the HESS algorithm is available upon request from the authors.

### Postprocessing analysis

In this section we present some of the postprocessing operations required to extract useful information from the rich output of our model. We stress that, while here for simplicity we are not using the output of the heated chains, following Gramacy *et al.* (2010), posterior inference could also be carried out using the information contained in all the chains.

The primary quantity of interest is the posterior propensity of each predictor to be a hotspot. In the spirit of cluster detection rules in disease mapping (Richardson *et al.* 2004), we use tail posterior probabilities of the propensities ρ* _{j}*,

*i.e.*, declare the

*j*th predictor to be a hotspot if, (7)where

*t*is a chosen threshold. We have found by empirical exploration and simulations that choosing a posterior threshold of

*t*= 0.8 gives good performance across different scenarios with varying dimensions (data not shown).

The next quantity of interest is the posterior probability of inclusion for the pair (*k*, *j*). Following Petretto *et al.* (2010), the *marginal probability of inclusion* is (8)where is the latent binary vector sampled at iteration *s*, is the model posterior probability obtained through inexpensive numerical integration in the full output (see File S1, Section S.3) and is the constant of normalization. The Bayes factor (BF) for the pair (*k*, *j*) is derived from (8) as the ratio between posterior odds and prior odds (9)where *E*(*p _{γk}*) is the

*a priori*expected number of genetic control points for the

*k*th transcript.

Similarly to (8), if of interest, we can further evaluate the joint posterior probability of the set of predictors declared as hotspots as (10)with *C ^{k}* as before and ℋ the set of markers identified as hotspots.

Finally, the *best model visited* is defined as (11)Note that the configuration posterior probability *p*(**Γ**|**Y**) (see File S1, section S.3) can be used as an alternative weight in (8) and (10) or to derive the max *a posteriori (MAP) configuration visited*

## Results

### Simulation studies

We carried out a simulation study to compare our algorithm with recently proposed multiple response models: MOM (Kendziorski *et al.* 2006), BAYES (Jia and Xu 2007), and M-SPLS (Chun and Keleş 2009).

To create more realistic examples, we decided not to simulate the **X** matrix, but to use real human-phased genotype data spanning 500 kb, region ENm014 (chromosome 7: 126,368,183–126,865,324 bp), from the Yoruba population used in the HapMap project (Altshuler *et al.* 2005) as the design matrix. After removing redundant variables, the set of SNPs is reduced to *P* = 498, with *n* = 120, giving a 120 × 498 **X** matrix. As noted by Chun and Keleş (2009), high correlations between markers might affect the performance of variable selection procedures that do not explicitly consider such a grouping structure. The benefit of using real human data are to test competing algorithms when the pattern of correlation, *i.e.*, linkage disequilibrium (LD), is complex and blocks of LD are not artificial, but they derive naturally from genetic forces, with a slow decay of the level of correlation between SNPs (see Figure S1).

In the simulated examples, we carefully select the SNPs that represent the hotspots (Figure S1): (i) all hotspots are inside blocks of correlated variables; (ii) the first four SNPs are weakly dependent (*r*^{2} < 0.1); and (iii) the remaining two SNPs are correlated with each other (*r*^{2} = 0.46) and also linked to one of the previous simulated hotspots (*r*^{2} = 0.52 and *r*^{2} = 0.44, respectively), creating a potential masking effect difficult to detect. Apart from the hotspots, no other SNPs are used to simulate transcript–SNP associations. We simulated four cases:

SIM1: In this example we simulated

*q*= 100 transcripts from the selected six hotspots, with some transcripts predicted by multiple correlated markers (polygenic control): for instance transcripts 17–20 are regulated by three SNPs at the same time (see Figure S2). Altogether we simulated 94 transcript–SNP associations in 50 distinct transcripts. The effects were simulated from a normal density with smaller variance than in Jia and Xu (2007), β∼_{kj}*N*(0,0.3^{2})with**ε**_{k}∼*N*(**0**, 0.1^{2}**I**) to mimic the smaller signal-to-noise ratio expected in genetically heterogeneous human data._{n}SIM2: As in the previous example, we simulated 100 responses, but there are only three hotspots with the same simulated association as before, leading to 64 transcript–SNP associations in 30 distinct transcripts. Moreover we created potential false-positive associations by simulating transcripts 81–90 and 91–100 using a linear transformation of transcript 20 with a mild negative correlation (in the interval [−0.5, −0.4]) and of transcript 80 with a strong positive correlation (in the interval [0.8, 0.9]), respectively. Since we create false-positive associations, the scenario will inform on how different algorithms behave when correlations among some transcripts are not due to SNPs.

SIM3: This simulation set-up is identical to the first scenario for the first 100 responses, but we increase the number of simulated responses to

*q*= 1,000, simulating the further 900 transcripts from the noise.SIM4: This is the same as the second simulated data set for the first 100 responses, with additional 900 responses simulated from the noise, giving altogether

*q*= 1000 responses.

We discuss here the hyperparameters set-up. Since *a priori*, in addition to a large effect of a SNP that is located close to the transcript (*cis*-eQTL), we expect only a few additional control points associated with the variation of gene expression (typically *trans*-eQTL); in HESS we set *E*(*p*_{γ}* _{k}*) = 2 and Var(

*p*

_{γk}) = 2, meaning the prior model size for each transcript response is likely to range from 0 to 6 (Petretto

*et al.*2010). Following Kohn

*et al.*(2001), we fixed

*a*

_{σ}= 10

^{−10}and

*b*

_{σ}= 10

^{−3}, giving rise to a noninformative prior on the error variance. We run the HESS algorithm for 6000 sweeps with 1000 as burn-in with three chains and ϕ = 1/4. Computational time is similar for the first two simulated examples, 6 hr, and 10 times greater for the last two simulated scenarios on a Intel Xeon CPU at 3.33 GHz with 24 Gb RAM.

We run BAYES for 15,000 sweeps with 5000 as burn-in, recording sampled values every 5 sweeps. The variance δ of the spike component 1 is set 10^{−4}, which is 100 times lower than the noise variance. Since the code available from the authors was written in SAS/IML, we recoded their Gibbs sampler in Matlab. We used the default parameters for MOM, while in M-SPLS the two tuning parameters are obtained through cross-validation selected in the interval *K* = 1,…, 10 and η = 0.01, …, 0.99. Each simulated example was replicated 25 times and we run the four algorithms on each replicate.

#### Power to detect hotspots:

The identification of the hotspots is of primary interest for all the algorithms we are comparing. In HESS using the tail posterior probability Pr(ρ* _{j}* >1|

**Y**), we can rank the predictors according to their propensity to be a hotspot, while in BAYES the posterior mean of the common latent probability ω

_{j}_{,}

*E*(ω

*|*

_{j}**Y**) is utilized to prioritize important markers. In MOM the strength for a predictor to be a hotspot is not directly available but, as suggested by the authors, given a marker, it can be obtained by taking a suitable quantile of the transcript–marker associations distribution across responses. We use their R function get.hotspots recording the average of the distribution for each predictor. M-SPLS, after cross-validation, provides a list of latent components that predicts most of the variability of

**Y**. While this group cannot be interpreted directly as the list of hotspots, we use it to test the existing overlap between the simulated hotspots and the latent components. Finally, differently from the analyses presented in Jia and Xu (2007) and Chun and Keleş (2009), in the HESS power calculation, we simply rank the evidence for being a hotspot provided by each algorithm across the 25 replicates. Therefore we are not using any method-specific procedure to call a hotspot, based for instance on FDR considerations, that could influence the comparison results.

Figure 1 shows the ROC curves for the four algorithms considered. HESS (blue lines) outperforms all the other methods with sizeable power on the simulated examples. It is not significantly affected by the dimension of the eQTL experiment (top, *q* = 100; bottom, *q* = 1000). This is somehow expected since the hotspot propensity does not depend directly on the number of transcripts analyzed (see File S1, section S.2.3). Spurious correlations among transcripts not due to SNPs (right) have a negligible effect on the HESS power, showing robust properties of our algorithm in detecting hotspots under different scenarios.

The other methods show good properties when *q* = 100, but their power degrades sensibly when *q* = 1000. This is expected for BAYES (green lines) since *E*(ω* _{j}*|

**Y**) is affected by

*q*, while the performance of MOM (red lines) is more stable. MOM (red dashed lines) shows good power in the simulated examples even when the number of markers is larger than the number of traits, a situation that MOM is not designed for (top). Finally M-SPLS has greater power than BAYES, but it is outperformed by both HESS and MOM in the more sparse scenarios (bottom). Looking more closely at the list of latent vectors identified by M-SPLS (data not shown), we noted that the simulated hotspots at SNP 362 and 466 that are linked to SNP 239 were rarely selected (false negative) in both SIM1 and SIM3. On the contrary, SNP 75 is very often included (false positive) in the list of latent vectors in all the scenarios. This might reflect the high correlation between SNP 362 and 466 with SNP 239, as well as the strong dependency between SNP 75 and 30 where we simulated a hotspot. This suggests that M-SPLS has limited efficiency in the presence of complex correlation patterns in the design matrix. In Figure S3, Figure S4, Figure S5, and Figure S6, interested readers can find the visual representation of the signals detected by each algorithm and averaged across the 25 replicates.

#### Power to detect transcript-marker associations:

Figure 2 shows the ROC curves of the transcript–marker marginal associations detected by each method. Also in this case, to perform the power calculation, we are not using any method-specific way to declare a significant association, since we simply record the output from each algorithm and rank it across the 25 replicates. In particular we use the marginal probability of inclusion *p*(γ* _{kj}* = 1|

**y**

*) (8) for HESS; the posterior frequency for BAYES, where is the value recorded at iteration*

_{k}*s*; the transcript–marker association provided the MOM object momObj; and finally the associations selected by bootstrap confidence interval at different type I error levels (α = 10

^{−4}, 10

^{−3}, 10

^{−2}, 0.05) for M-SPLS.

For transcript–marker association detection, we find that HESS has higher power than that of the other methods in all the simulated scenarios. As expected when more responses are included, the power decreases slightly (bottom), while spurious associations due to the correlation between transcripts do not seem to affect the ability of HESS to distinguish between true and false signals (right). MOM is quite stable across scenarios, but it reaches only half of the power of HESS. BAYES and M-SPLS have similar behavior and their performance degrades when *q* = 1000. BAYES, in particular, has very low power since the shrinkage to the null effect, caused by common latent probabilityω_{j}_{,} is particularly strong in SIM3 and SIM4.

The different power of the methods considered can be better understood by looking at Figure S3, Figure S4, Figure S5, and Figure S6, where, for each simulated example, we averaged the evidence of transcript–marker association across replicates. HESS is able to identify the correct simulated pattern, with very few false positives. When false-positive associations are simulated, for instance, in SIM3 and SIM4, HESS assigns on average lower posterior probability of inclusion than for the true positive ones (Figure S4 and Figure S6, top left). While MOM is able to identify the simulated hotspots, it finds it difficult to separate the true transcript–marker associations from the spurious ones (Figure S3, Figure S4, Figure S5, and Figure S6, bottom left). The main limitation of M-SPLS is the correct identification of the latent vectors when highly correlated predictors are considered (Figure S3, Figure S4, Figure S5, and Figure S6, top right). Finally BAYES is able to identify the simulated pattern when *q* = 100 (Figure S3 and Figure S4, bottom right), but it seems to be too conservative when the number of responses is large, *q* = 1000 (Figure S5 and Figure S6, bottom right). The higher false-negative rate in BAYES may depend on the poor efficiency of the MCMC sampler (which is based exclusively on the Gibbs sampling that is not able to jump between distant competing models) and on the spike and slab prior that is not integrated out. The latter influences the sampling of γ* _{kj}* since the latent binary vector depends on the regression coefficients (see Figure S7 for an illustration).

### Real case studies

Here we present two applications of HESS to: (i) mouse gene expression data published in Lan *et al.* (2006) that is commonly used as a benchmark data set for detection of eQTL (Chun and Keleş 2009) and eQTL hotspots (Kendziorski *et al.* 2006; Jia and Xu 2007) and (ii) human monocytes expression data set recently analyzed for disease susceptibility by Zeller *et al.* (2010).

#### Mouse data set:

The mouse data set has been previously described in detail (Lan *et al.* 2006), and it consists of 45,265 probe sets the expression of which has been measured in the liver of 60 mice. Mice were collected from the F_{2}*-ob/ob* cross (B6 × BTBR) and genotype data were available for 145 microsatellite markers from 19 autosomal chromosomes. To make our analysis comparable with previously reported studies (Jia and Xu 2007; Chun and Keleş 2009), we focused on 1573 probe sets showing sizeable variation in gene expression in the mouse population (sample variance >0.12). Running HESS for 12,000 sweeps with 2000 as burn-in and the same choice of the hyperparameters described in the simulation studies, among the 145 markers 16 were identified with posterior tail probability >0.8, regulating a significant number of probe sets (Table S1). We report the genome location of the identified hotspots in Figure 3 and show transcript–marker associations in Figure 4. Since large hotspot propensity reveals that multiple traits are controlled by the same marker, we focused on biologically meaningful transcript–marker associations by using marginal probability of association >0.95 (corresponding to local FDR 5%, Ghosh *et al.* 2006). Six markers were found to control more than 5% of all analyzed probe sets as shown in Figure 3. While marker D15Mit63 was previously detected by BAYES and M-SPLS, three other major regulatory points were identified solely by our method: D13Mit91, D18Mit9, and D18Mit202, controlling 14.1, 10.6, and 9.7% of all analyzed probe sets, respectively (Table S2).

The regulatory hotspot at marker D13Mit91 is located within the *Kif13a* (kinesin family member 13A) gene, which is involved in intracellular protein transport and microtubule motor activity via direct interaction with the AP-1 adaptor complex (Nakagawa *et al.* 2000). This hotspot is associated with 222 probe sets, representing 190 distinct well-annotated genes, that are enriched for specific gene ontology (GO) terms, including “protein localization” (*P* = 4.2 × 10^{−6}), “protein transport” (*P* = 5.7 × 10^{−6}), and “establishment of protein localization” (*P* = 6.4 × 10^{−6}). Hence, given its molecular function *Kif13a* is likely to be a candidate master regulator of the genes implicated with protein transport, and whose expression is associated with marker D13Mit91.

The other two newly identified markers, D18Mit9 and D18Mit202, are located on mouse chromosome 18. D18Mit9 resides within a known QTL (*Hdlq30*) involved in HDL cholesterol levels (Korstanje *et al.* 2004) whereas D18Mit202 resides within a known diabetes susceptibility/resistance locus (Idd21, insulin-dependent diabetes susceptibility 21) (Hall *et al.* 2003).

#### Human data set:

The human data set included 648 probe sets, representing 516 unique and well-annotated genes (Ensemble GRCh37), that were found to be coexpressed in monocytes, delineating a network driven by the *IRF7* transcription factor in 1490 individuals from the Gutenberg Heart Study (GHS) (for details on the network analysis, see Heinig *et al.* (2010)). This *IRF7*-driven inflammatory network (IDIN) was also reconstructed in a distinct population cohort: 758 subjects from the Cardiogenics Study showing significant overlap with the network in GHS. The “core” of the network consisted of a small gene set (*q* = 17), including *IRF7* and coregulated target genes, the expression of which was found to be *trans*-regulated by a locus on human chromosome 13q32 using MANOVA in Cardiogenics (Heinig *et al.* 2010). However, this *trans*-regulation was not found in the GHS study, using similar MANOVA analysis.

Here we take a new look and use HESS to analyze the larger IDIN with 648 probe sets in the GHS population (*n* = 1490 individuals) and the SNP set (*P* = 209) spanning 1 Mb on chromosome 13q32 (data available upon request from Stefan Blankenberg under the framework of a formalized collaboration via a Memorandum Transfer Agreement). While MOM and BAYES fail to detect any signal at this locus, using HESS we found two SNPs, rs9557207 and rs11616269, which are detected as hotspots for the IDIN expression with tail posterior probability 0.83 and 0.91, respectively (Figure 5). These SNPs are located 45.3 kb (rs9557207) and 25.1 kb (rs11616269) from SNP rs9585056, which previously showed significant *trans*-effect on the core gene set of the network in Cardiogenics (*P* = 5.0 × 10^{−3}). This region was also associated with *EBI2* expression (*P* = 2.2 × 10^{−8}), the candidate gene at this locus, and with type I diabetes (T1D) (*P* = 7.0 × 10^{−10}) (Heinig *et al.* 2010). For the two identified hotspots, we looked in detail at each transcript–marker association and compute their BF as given in (9). We observe that 26 and 13 transcripts show clear evidence of associations (BF > 10, Kass and Raftery 2007) in the two hotspots identified (Table S3) delineating the extent of regulatory effects. To further calibrate this evidence, we investigated BF for marker–transcript associations in a comparable simulated set-up, that of SIM3. Using the threshold BF > 10 would lead to declaring <5% false positive marker–transcript associations in the identified hotspots (data not shown). Note that most of these transcripts (80%) are found only in the network inferred in GHS and not with the Cardiogenics network, suggesting a complex pattern of regulatory effects around locus rs9585056 which is highlighted in a specific manner in each population. These population-specific regulatory effects could reflect differences in monocytes selection protocols between GHS and Cardiogenics (see Heinig *et al.* (2010) for details). However, the identification of hotspots at the 13q32 locus by HESS in GHS represents a significant replication of the findings previously reported, which reflects the increased power of HESS over alternative methods.

## Discussion

We have presented a new hierarchical model and algorithm, HESS, for regression analysis of a large number of responses and predictors and have applied this to hotspot discovery in eQTL experiments. Simulating a variety of complex scenarios, we have demonstrated that our approach outperforms currently used algorithms. In particular, HESS shows increased power to detect hotspots when a large number of transcripts are jointly analyzed. This is due to the propensity measure ρ* _{j}* that we use, which quantifies the latent hotspot effect independently of the response dimensionality. One improvement of HESS over vanilla MCMC-based algorithms is in the search procedure that efficiently probes alternative models and assesses their importance, thus providing a reliable model space exploration (Bottolo and Richardson 2010). We have also illustrated the potential of HESS to discover regulatory hotspots in two eQTL studies that encompass diverse genetic contexts (animal model and human data). In contrast to other methods, using HESS, we were able to replicate an established regulatory control of a large inflammatory network in humans (Heinig

*et al.*2010). Moreover, in the mouse data set, we identified a new candidate (

*Kif13a*) for the regulation of a set of genes implicated in protein transport, which was not detected by other approaches.

Our model is embedded in the linear regression framework with additive effects. One distinct feature of our formulation is the multiplicative decomposition of the selection probabilities and its hierarchical set-up, which allows other structures and/or different types of prior information to be included. For example, specific weights π* _{kj}* (suitable normalized) could be added in (5), ω

*= ω*

_{jk}*× ρ*

_{k}*× π*

_{j}*, to provide additional prior information about the regulation of*

_{kj}*k*th transcript. This may include

*cis*-acting genetic control or auxiliary information on regulatory effects of the

*j*th SNP (

*i.e.*, evolutionary conservation, coding, noncoding, genomic location, etc.) (Lee

*et al.*2009). Likewise, additional structure on the responses (

*e.g.*, KEGG pathways membership, predicted targets of transcription factors, protein complexes, etc.) could be included using

*k*-indexed weights, ω

*, to favor detection of hotspots for similar responses.*

_{k}Another possible extension of our method is the inclusion of interactions in the linear model and their efficient detection. Recent advances in this direction have employed either a stepwise search for interactions between preselected main effects (Wang *et al.* 2011) or partition models with which to discover modules or clusters of transcript–marker responses (Zhang *et al.* 2010). Such approaches could be embedded in our variable selection algorithm.

The current Matlab version of HESS represents a first step toward a more efficient implementation in high-level coding languages (currently undergoing), taking advantage of the existing C++ version of ESS algorithm (Bottolo *et al.* 2011). The approach that we propose here is ideally suited after prioritizing candidate genomic regions or gene networks, as shown in the discussed human case study. The flexibility to incorporate prior biological knowledge makes our method suitable for a wide range of analyses beyond eQTL hotspots detection, including genetic regulation of miRNA targets and metabolic and epigenetic phenotypes.

## Acknowledgments

We thank two anonymous referees and the associate editor for their comments that improved the manuscript. L.B. received funding from the Wellcome Trust Value-in-People award. S.R. gratefully acknowledges support from ESRC National Centre for Research Methods (BIAS II node, grant RES-576–25–0015). L.B. and S.R. acknowledge support from the Medical Research Council (grant G1002319). E.P. and S.A.C. acknowledge support from the Medical Research Council and the British Heart Foundation (grant FS/11/25/28740).

## Footnotes

*Communicating editor: G. A. Churchill*

- Received July 5, 2011.
- Accepted August 23, 2011.

- Copyright © 2011 by the Genetics Society of America