## Abstract

In 2001, Sen and Churchill reported a general Bayesian framework for quantitative trait loci (QTL) mapping in inbred line crosses. The framework is a powerful one, as many QTL mapping methods can be represented as special cases and many important considerations are accommodated. These considerations include accounting for covariates, nonstandard crosses, missing genotypes, genotyping errors, multiple interacting QTL, and nonnormal as well as multivariate phenotypes. The dimension of a multivariate phenotype easily handled within the framework is bounded by the number of subjects, as a full-rank covariance matrix describing correlations across the phenotypes is required. We address this limitation and extend the Sen–Churchill framework to accommodate expression quantitative trait loci (eQTL) mapping studies, where high-dimensional gene-expression phenotypes are obtained via microarrays. Doing so allows for the precise comparison of existing eQTL mapping approaches and facilitates the development of an eQTL interval-mapping approach that shares information across transcripts and improves localization of eQTL. Evaluations are based on simulation studies and a study of diabetes in mice.

THE quantitative trait loci (QTL) mapping framework developed by Sen and Churchill (2001), referred to hereinafter as the Sen–Churchill framework, unifies many methods for QTL mapping in inbred line crosses. The seminal work of Lander and Botstein (1989) and subsequent methods including Haley–Knott regression (1992), composite interval mapping, and multiple QTL mapping (Jansen 1993; Zeng 1993, 1994; Jansen and Stam 1994), are all represented, at least approximately, as special cases of the framework. The framework also accounts for covariates, nonstandard cross designs, missing genotype data, genotyping errors, multiple interacting QTL, and nonnormal as well as multivariate phenotypes. As a result, it provides a powerful approach to localize the genetic basis of quantitative traits.

There has been much interest recently in identifying the genetic basis of thousands of gene- expression traits measured via microarrays (Brem *et al*. 2002; Schadt *et al*. 2003; Yvert *et al*. 2003; Cox 2004). The multi-trait version of the Sen–Churchill framework is based on the multivariate normal distribution. This approach becomes problematic when the number of traits is larger than the number of subjects, as the estimated covariance matrix will have less than full rank. To address this, we here extend the Sen–Churchill framework to accommodate expression phenotypes. We first highlight aspects of the Sen–Churchill framework important to our development, and then detail the extension. We show that the extended framework generalizes the currently available expression QTL (eQTL) mapping methods and facilitates the development of an approach that allows for both interval mapping of eQTL and information sharing across transcripts. Evaluations are based on simulation studies and a study of diabetes in mice. Generalizations of the framework are also discussed. Many of the technical details can be found in the appendixes.

## A FRAMEWORK FOR EXPRESSION QTL INFERENCE

#### The Sen–Churchill framework:

The Sen–Churchill framework supports a Bayesian approach to QTL mapping that accommodates a variety of phenotypes and data structures. Much of the flexibility of the approach is due to two main features. The first is marked separation of the genetic model, which relates phenotype to genotype, and the linkage model, which relates putative QTL genotype to the marker map. The second feature is that computation relies on an efficient Monte Carlo component instead of a more complex MCMC procedure as employed in a number of other Bayesian QTL methods (Satagopan *et al*. 1996; Yi and Xu 2000; Yi 2004). As we discuss in detail below, these two features allow for accommodation of microarray data as a phenotype within the framework. We here provide an overview of the framework, focusing on aspects important to our extension.

Suppose that quantitative traits are measured for *n* members of an inbred line cross. Denote the traits by and denote the corresponding marker data by the *n* × *M* matrix *m*, where *M* denotes the total number of markers. Marker location and genetic distances are assumed known, although in practice these are estimated. A genetic model *H* describes the way in which QTL genotypes determine a phenotype; it is prescribed by the number of QTL, their locations, and the way in which they act and interact to affect the phenotype. Assuming *p* QTL in a genetic model, let γ denote the *p*-dimensional vector of QTL locations and *g* denote the *n* × *p* matrix of QTL genotypes. The parameters of the genetic model are denoted by μ.

Of primary interest is the posterior distribution of QTL location, *p*(γ | *y*, *m*), given by(1)where modes of *p*(γ | *y*, *m*) estimate QTL position. An exact evaluation of Equation 1 is computationally prohibitive, but an approximation can be obtained by sampling multiple versions of the putative QTL genotypes *g* and averaging as follows:

Select a regularly spaced grid

*G*of pseudomarker locations, locations for which genotypes are not known, and create*q*realizations of the pseudomarkers by sampling from*p*(*g*|*m*). Assuming known genetic distances and no crossover interference, a Markov chain sampling scheme can be used. Each realization of pseudomarker genotypes is an*n*×*G*matrix.For the assumed genetic model

*H*, a*p*-dimensional vector of pseudomarker locations corresponding to the QTL, γ, is prescribed; and the^{H}*i*th realization of pseudomarker genotypes provides*g*(γ_{i}), an^{H}*n*×*p*matrix of pseudomarker genotypes at the QTL locations.For each realization, calculate a weight under the assumed genetic model

*H*. The weight for the*i*th realization isAn average over

*q*of these weights approximates (1), according to the principle of importance samplingfor some constant of proportionality*C*.

#### Extensions to eQTL mapping:

Consider for simplicity a backcross population genotyped as *aa*(0) or *Aa*(1) at *M* markers (this simplification to a backcross is not required and is relaxed in the *Applications to data from a study of diabetes*). For eQTL mapping, the observed phenotype data *y* are no longer a vector as above, but rather a *T* × *n* matrix of expression levels. Specifically, *y* = (*y*_{1}, *y*_{2},…, *y _{T}*)′, where vector

*y*= (

_{t}*y*

_{t}_{1},…,

*y*) denotes the (possibly transformed) expression levels for transcript

_{tn}*t*measured in

*n*animals. As in the univariate phenotype case,

*m*denotes an

*n*×

*M*matrix containing genotypes on

*M*markers.

Of most interest is the identification of significant linkages between transcripts and genome locations. To be precise, a transcript *t* is linked to location *l* if where denotes the latent mean level of expression for transcript *t* for the population of animals with genotype 0(1) at location *l*. Two *T* × *G* matrices, θ^{0} and θ^{1}, contain the latent mean levels of expression ; and, as above, *G* denotes the total number of locations considered. In the Sen–Churchill framework, of primary interest is the posterior evaluation of γ, a vector of QTL locations. In this context, γ is transcript specific. For example, for transcript *t*, γ_{t} would contain indexes *l*′ such that .

#### Single eQTL mapping methods:

Suppose that a transcript is affected by at most one genotype location *l* (this assumption can be relaxed as discussed later) and consider inference at location *l*. Of most interest is the posterior probability that transcript *t* is linked to location *l*. We show in appendix b that(2)where is the marginal density describing the data in the case of linkage to *l*.

Equation 2 is similar in form to (1), but there are some important differences. In Equation 2, conditioning is done on the full set of transcripts. An assumption of conditional independence across transcripts (see appendixes) yields a right-hand side (RHS) that is *evaluated* only at the transcript of interest *t*. The form of *f*_{P1} determines whether or not information from other transcripts affects the evaluation. For example, if *f*_{P1} is taken to be a univariate Gaussian (or other parametric) distribution, then the RHS is completely determined by the data at *t* since the parameters of *f*_{P1} do not depend on other transcripts. An application of the extended Sen–Churchill framework in this case would consist of a repeated application of a single-transcript analysis to each expression trait in isolation. This has been done in a number of eQTL studies to yield effective results. However, with this approach, there is no information shared across transcripts. As pointed out in a number of articles on microarray data analysis (Newton *et al*. 2001; Tusher *et al*. 2001; Kendziorski *et al*. 2003; Smyth 2004; Cui *et al*. 2005), information sharing is important to improve sensitivity and moderate test statistics that are otherwise prone to inflated error. Kendziorski *et al*. (2006) demonstrated this in the context of eQTL studies and proposed the mixture-over-markers (MOM) model to localize eQTL while allowing for information sharing across transcripts. The MOM model is further developed in Gelfond *et al*. (2007).

The MOM model is represented as a special case of the extended framework when *f*_{P1} is taken to be a certain predictive density. In short, assume measurements of transcript *t* for animal *r*, denoted *y _{tr}*, arise as conditionally independent random deviations from an observation distribution

*f*

_{obs}(· | , θ) with the 's as random effects described by a distribution π(μ). The model is assumed to be the same across locations and so dependence on

*l*is suppressed. In this model, an equivalently expressed transcript

*t*presents data

*y*according to the distribution(3)where and

_{t}*f*

_{P1}(y

_{t}) =

*f*

_{P0}()

*f*

_{P0}() describes the data for mapping transcripts, owing to the fact that different mean values, and , govern the different subsets and of samples and are considered independent draws from π(μ) (see appendixes). Here, and denote the collection of expression values from subjects with genotypes

*aa*and

*Aa*, respectively. As detailed in Kendziorski

*et al*. (2006), a Gaussian model is assumed for

*f*

_{obs}(·) and π(·). We also allowed for the possibility that different clusters of transcripts could present data with different variances.

Specification of the denominator *p*(*y _{t}* |

*m*) of Equation 2 is not required if closed forms for (or good approximations of) parameter estimates are available and estimation of the false discovery rate (FDR) is not of interest. When closed forms are not available and/or calculation of estimated FDR is of interest,

*p*(

*y*|

_{t}*m*) must be evaluated. Note that , where

*p*(γ

_{t}= 0) implies that the transcript does not map to any of the

*G*locations. We do not assume any specific priors on the mixing proportions. They will be estimated using the data. As detailed in appendix b, Equation 2 then becomes(4)Note that conditioning on genotype is dropped if the transcript is not linked to location

*l*as all measurements arise from a distribution with common mean and so genotype information, which prescribes groups in the case of a transcript mapping to

*l*, is not required.

When evaluated at markers only, where genotypes are known, Equation 4 is identical to the MOM model. Extensions of MOM to interval mapping have been difficult to date, as evaluation of Equation 4 can be prohibitive in between markers. Since the *l*th column of *g*, denoted *g ^{l}*, is a vector of length

*n*, there are 2

*possible genotypes (for a backcross); and as a result, the integral in Equation 4 is a very large mixture, when*

^{n}*n*is even moderately large. In practice, one could potentially restrict to fewer possibilities since many genotype vectors have very small probabilities. However, as the number of individuals in the study gets large (>200), this quickly becomes computationally infeasible even with the restriction. Fortunately, pseudomarkers can be used, as in the Sen–Churchill framework, to overcome this problem.

In the extended framework, multiple versions of pseudomarkers are sampled from *p*(*g ^{l}* |

*m*). Suppose for each location

*l*(

*l*= 1,…,

*G*),

*q*genotype vectors are sampled from the proposal distribution

*p*(

*g*|

*m*) to yield (, ,…, ). Then Equation 4 is approximated by(5)and modes of this distribution are used to estimate eQTL positions. One can apply this approach to grids of varying sizes (

*i.e*., varying

*G*) to localize eQTL at and in between markers. We refer to this approach as pseudomarker MOM (psMOM).

#### Simulations:

We conducted a small set of simulations to compare psMOM with traditional interval mapping (IM) applied to each transcript in isolation. The simulations are not designed to capture the many complexities of eQTL data, but rather they provide some preliminary information on operating characteristics in simple settings. Marker genotype data were simulated for four chromosomes, each of length 100 cM and having 11 equally spaced markers (10-cM spacing). We assumed that 15% of all transcripts map to at least one genomic location; 5% map to a single location on chromosome 1 (26 cM); 5% map to two locations on chromosome 2 (44 and 56 cM); the remaining 5% map to two locations on chromosome 3 (22 and 82 cM). No transcripts are affected by alleles on chromosome 4.

Backcross data were simulated for 200 animals and 4000 transcripts. Simulated intensities follow the approach described in Kendziorski *et al*. (2006). Briefly, we assume log intensities are normally distributed, which is consistent with the assumptions of both IM and psMOM. Transcript-specific means and variances are sampled from the empirical means and variances of the F_{2} cross described previously. The latent means of transcripts mapping to a single location satisfy . For the transcripts mapping to two locations *l* = (*l*_{1}, *l*_{2}), their latent means satisfy . Twenty simulated data sets were generated.

#### Implementation of IM:

For IM, we consider *f*_{P1} as univariate Gaussian with transcript-specific parameters obtained via the method of moments using only data at that transcript; 500 sets of pseudomarkers are generated every 2 cM and LOD scores are computed. To compare with the highest posterior density (HPD) regions of psMOM (see below), likelihood ratios (LRs) are derived from the LOD scores, normalized, and converted to quantities similar to posterior probabilities. For example, if *L*(*H*_{1}, *l*′)/*L*(*H*_{0}) denotes the likelihood ratio at location *l*′, we consider and as evidence of equivalent and differential expression at *l*′, respectively. We refer to these as LOD posterior probabilities. Transcripts with LOD posterior probability of differential expression exceeding some threshold are considered mapping transcripts. As shown in Tables 1–3⇓⇓, IM is evaluated for varying thresholds.

For some examples (noted in subsequent text), to compare with the HPD regions derived from LOD posterior probabilities, we also considered 1.5-LOD drop-support intervals around peak LOD scores (Mangin *et al*. 1994; Dupuis and Siegmund 1999). They are designed to target confidence regions of level 95%, but in general, these intervals are known to be biased in that they are too small (Visscher *et al*. 1996). On the other hand, confidence intervals that are slightly too small favor IM as eQTL appear to be better localized. To give IM the best results, we consider a 10-cM window around the true eQTL positions and define the respective LOD peaks as the highest LODs within the windows. The 1.5-LOD support intervals are then constructed. Of course, in practice, one does not have the luxury of knowing where to choose these peaks and perhaps only the largest peak would be identified. In this way, the results of this approach are further biased in favor of IM.

#### Implementation of psMOM:

Equation 4 is first evaluated at the genotyped markers and the *M*′ markers with posterior probabilities ≥0.9 define an HPD region. In particular, the posterior probabilities at the identified *M*′ markers are averaged across the mapping transcripts and, using this *M*′ vector, an HPD region is identified. Basically, the HPD region contains the minimum number of support points with corresponding posterior probabilities having a sum exceeding 1 – α. More precisely, the HPD region of level 1 − α is constructed by rank ordering posterior probabilities *p*_{(1)} ≤ *p*_{(2)} ≤ … ≤ *p*_{(n)} where and identifying the largest (*n*′) such that . The HPD region then consists of the support points corresponding to *p*_{(n′)}, *p*_{(n′+1)},…, *p*_{(n)}.

A pseudomarker grid is set up within the HPD region and multiple versions of pseudomarkers are generated. The model fit is carried out utilizing all genotyped markers *plus* pseudomarkers in the HPD regions. The procedure provides a matrix of posterior probabilities for every transcript at every test point. As in IM, 500 sets of pseudomarkers are generated every 2 cM; unlike IM, the pseudomarkers are only generated within the HPD region and thus the dimension of *G* here is reduced, thereby reducing the computational burden. The procedure gives a matrix of posterior probabilities for every transcript and every test point. Transcript *t* is defined to map to location *l* if the posterior probability in the second stage exceeds some threshold. As in IM, psMOM is evaluated for varying thresholds.

#### Choice of threshold:

A list of mapping transcripts with target FDR α can be constructed by taking those with posterior probability of equivalent expression less than α (Newton *et al*. 2004). This specifies transcripts that likely harbor at least one eQTL, but does not provide information on the total number of eQTL per transcript. For the latter, a linkage threshold must be set. The thresholds evaluated here for both IM and psMOM are varied from 0.1 to 0.4. For example, recall that the LOD posterior probability profiles from IM and the posterior probability profiles from psMOM each sum to 1. When a single eQTL is simulated, often the (LOD) posterior probability of linkage is quite large at the location of the eQTL (*e.g*., >0.95). However, with two eQTL, individual posterior probabilities are rarely that large since evidence is spread out across multiple locations. Thresholds could also be chosen on the basis of transcript-specific HPD regions (see supplemental material at http://www.genetics.org/supplemental/ for an example).

#### Power and false discovery rate calculation:

A call is said to be “correct” if the genome location identified is within 5 cM of the true eQTL (*i.e*., within the 10-cM window centered at the true eQTL location). In the case of two eQTL, at least one location has to be within the 10-cM window of a true eQTL for the identified eQTL to be deemed as correct. Power measures the ability to correctly identify mapping transcripts. It is calculated to be the ratio of the number of correct calls to the total number of eQTL. FDR is calculated as the ratio of the number of incorrect calls to the total number of calls. Incorrect calls consist here of nonmapping transcripts that are identified to map. Our calculation of FDR *does not* consider mapping transcripts that map outside the 10-cM window of the eQTL, since this led to greatly inflated FDR for IM. Specificity represents the proportion of nonmapping transcripts that are identified as nonmapping.

#### Tests for enrichment:

A number of efforts utilize information from multiple sources to annotate transcripts; and it is informative to identify sets of transcripts that are enriched for some annotation compared with a randomly sampled set of the same size. A hypergeometric calculation is often used to assess evidence of enrichment, but interpretation of resulting *P*-values is not straightforward due to the many dependent hypotheses tested. Furthermore, the hypergeometric calculation tends to result in small *P*-values when few transcripts are considered. For these reasons, it has been suggested that one consider only interesting small *P*-values obtained from a relatively large set of transcripts (>10) (Gentleman 2004). That is the practice we follow here, considering lists of at least size 10 and setting *P*-value thresholds for enrichment at 0.05. We considered biological functions annotated in the Gene Ontology (GO) database.

#### Software:

All calculations were carried out in R (http://www.r-project.org). The IM method was performed using the scanone function with the “imp” option in R/qtl (Broman *et al*. 2003).

## RESULTS

#### Simulation results:

The results from a single simulation are shown in Figure 1 (results are representative of those observed in the other 19 simulations). The top graphs show results from psMOM and the bottom graphs show those from IM. The left graphs demonstrate the average linkage evidence -posterior probabilities in psMOM and LOD posterior probabilities in IM; and the right graphs show transcript-specific HPD regions.

From the average linkage evidence shown in the left graphs, the two approaches considered are very similar: they each identify the regions for the single eQTL on chromosome 1 and the two unlinked eQTL on chromosome 3. Both MOM and marker regression miss the two linked eQTL, identifying a wide peak only in the middle region of chromosome 2. The interval-mapping approaches (psMOM and IM) further refined the eQTL underneath this wide peak.

Differences between the approaches are more pronounced when transcript-specific linkage evidence is considered. As shown in the right graphs of Figure 1, psMOM precisely identifies the eQTL correctly for most of the mapping transcripts. In contrast, the regions surrounding the IM identifications are relatively wide. This result could be due to the way the LR normalization was done, to the fact that information shared across transcripts is not accounted for, or to both. To test the former, instead of HPD regions constructed from LOD posterior probabilities, we considered confidence intervals constructed using 1.5-LOD drop intervals around the peak LOD score. As detailed in the *Implementation of IM* section, this procedure favors IM. Even so, the approach still provided very imprecise estimates of eQTL location, much worse than those shown in Figure 1, and we do not recommend this in practice.

The results for this single simulation hold across simulations as shown in Tables 1–3⇓⇓. Table 1 reports power at varying thresholds averaged over 20 simulated data sets for each eQTL and each chromosome. For low thresholds, power is similar for both approaches, with psMOM showing slightly higher power. Table 2 shows that FDR from psMOM is well controlled for the three chromosomes. The level stays the same under different thresholds. On the other hand, the FDR from IM is quite high at the 0.1 cutoff point; it decreases with increasing threshold, but with reduced power. Table 3 shows the specificity calculated over an average of 20 simulated data sets. The specificities are very similar for both approaches and they are satisfactory.

#### Applications to data from a study of diabetes:

The data set considered here is discussed in detail in Lan *et al*. (2006) and is available at GEO (Barrett *et al*. 2007), accession no. GSE3330. Briefly, 60 mice (29 males and 31 females) were selected from an F_{2} population segregating for phenotypes associated with diabetes and obesity (Stoehr *et al*. 2000). The population was derived from B6 male and BTBR female parents. Selection was based on the selective phenotyping algorithm developed in Jin *et al*. (2004), which can substantially improve sensitivity for QTL localization compared with random sampling of the same sample size (Jin *et al*. 2004). The marker map consists of 145 microsatellite markers spanning the 19 mouse autosomes, with an average intermarker distance of 13 cM. Over 90% of the animals are genotyped at any given marker.

Liver total RNA was extracted from frozen tissue samples with RNAzol reagent (Tel-Test). Crude RNA samples were purified with RNeasy mini columns (QIAGEN, Valencia, CA) before hybridization. The RNA samples were processed according to the Affymetrix Expression Analysis technical manual. Expression levels for 45,265 probe sets (referred to hereinafter as transcripts) were measured using the MOE430A and MOE430B chips for each of the 60 F_{2} mice. Preprocessing and normalization was done using robust multi-array average (RMA) (Irizarry *et al*. 2003) to obtain a single normalized summary score of expression for each gene in each animal.

Both IM and psMOM were applied to the data; psMOM accommodates F_{2} populations by increasing the number of expression patterns. For example, with three genotype groups (0, 1, and 2) there are three latent means of interest and four non-null expression patterns for each transcript *t* at each location .

Posterior probabilities from psMOM and LOD posterior probabilities from IM were averaged across transcripts to identify the genomic regions of most interest. As in the simulation study, the average results from IM and from psMOM largely agreed (Figure 2a). However, when looking at a finer scale, one does observe important differences (Figure 2b). There are two locations in particular (on the distal regions of chromosomes 2 and 5) where psMOM shows some evidence for linkage but IM does not. To test whether the regions identified by psMOM might be meaningful, we consider the biological functions of identified transcripts.

We tested for functional enrichment among the transcripts mapping to the two subpeaks (call these *trans1a* and *trans1b*) on the distal region of chromosome 2. Both sets show significant enrichment of lipid-metabolism and fatty-acid-metabolism genes (*P*-values are 0.0016 and 0.0044, respectively). The lipid-metabolism group on the distal region of chromosome 2 coincides with the lipid-metabolism cluster discovered in Lan *et al*. (2006). Several QTL for obesity and related traits have been mapped to this region (Stoehr *et al*. 2000). As shown in Figure 2b, there is some evidence of linkage (a single peak) on the distal region of chromosome 2 provided by MOM (the first pass of psMOM at markers only). The transcripts mapping to this region did not show significant linkages for lipid metabolism or fatty acid metabolism (*P*-values are 0.2465 and 0.2302, respectively) or for any categories that appeared to be related to our diabetes or obesity phenotypes of interest.

In addition to the chromosome 2 linkages, we considered two linkage peaks on the distal region of chromosome 5, near the marker *D5Mit240*, as these peaks are identified by psMOM alone. As on chromosome 2, tests here show enrichment for lipid-transport and fatty-acid-metabolism genes. In addition, the enrichment of genes responsible for positive regulation of metabolism is highly significant (*P*-value = 0.002). A closer look at the mapping list reveals some interesting members. They include PPARα and PPARγ, two major lipid-metabolism transcription factors (Attie and Kendziorski 2003). Other interesting genes include fatty-acid-synthase genes (Fasn, Elovl6, Elovl5, and Fads2), lipid-transport genes (Scp2, Pltp, and Apoa4), and two fatty-acid-metabolism genes (Gpam and CD36). Taken together, these results provide some support for the peaks uniquely identified by psMOM.

## DISCUSSION

We have extended the QTL mapping framework of Sen and Churchill (2001) to accommodate expression phenotypes. The Bayesian formulation prescribed by Sen and Churchill (2001) and the pseudomarker-sampling approach developed there is maintained. Our extension relies on specifying a more general form for the genetic model, the model relating phenotype to genotype. By fitting the full genetic model to all phenotypes simultaneously, information can be shared across transcripts through the estimated hyperparameters, which in many cases leads to improved inference.

The extended framework generalizes most eQTL mapping approaches and in doing so facilitates their understanding, evaluation, and precise comparison by revealing their specific characteristics in the context of a common notation, which in turn provides an improved environment for addressing open questions and developing ideas for future methods. As an example, we considered a deficiency of the MOM model, namely that no information is provided between markers. Viewing MOM as a special case of the extended framework clarified how to address this deficiency using pseudomarkers. We expect that other open questions in the area of eQTL mapping can be more readily addressed in the context of this unified framework.

One such question might be the choice of an appropriate threshold. In most eQTL studies to date, thresholds are varied and the one that yields a list with many transcripts while controlling some measure of false positives at a reasonable level is used. A number of false positive measures have been considered; and, clearly, investigators define “many” and “reasonable” quite differently in different contexts (Kendziorski and Wang 2006). The framework presented here can be used to investigate common approaches and perhaps to rigorously address this question.

Without exact knowledge of appropriate thresholds for either psMOM or IM, our evaluations were based on varying thresholds. There appears to be a slight advantage of psMOM over IM, likely due to the information shared across transcripts. For some thresholds, the advantage is negligible; while for others, it is much more pronounced. A bigger advantage of psMOM is the precision provided by the HPD regions. Analogous regions were constructed from the LOD profiles, but these provided much less precise localization than the standard HPD regions of psMOM. Indeed, as we noted in *Power and false discovery rate calculation* section, the FDRs shown in Table 2 *did not* consider differentially expressed transcripts that mapped outside the 10-cM window of the eQTL. Considering these transcripts greatly inflated the FDR for IM. Considering alternative methods for LOD profile normalization did not yield results better than those shown here.

Finally, we have focused on representation of approaches for single-eQTL mapping. The simulation results show that the approaches can work well, even for two eQTL settings, much like single-QTL models can provide information on multiple QTL. Of course, single-QTL models will not work well when multiple QTL are tightly linked; and we here note that extensions of psMOM as presented are possible. In the context of the extended framework, the extension is seen as changes in *f*_{Pk}, where now *k* > 2. In particular, if a transcript *t* is affected by two genotype locations *l*_{1} and *l*_{2}, then four latent means are of interest: , and . Here denotes the latent mean level of expression for transcript *t* for the populations of animals with genotype (*g*_{1}, *g*_{2}) at locations *l*_{1} and *l*_{2}. These latent means can be arranged into 15 possible expression patterns, all of which may be of interest (see supplemental materials at http://www.genetics.org/supplemental/ for pattern specification and further detail). As before, of primary interest is the posterior probability of particular expression patterns. These can be calculated for any pattern of interest.

The approach was applied to the simulated data described previously using all 15 expression patterns (*i.e*., *k* = 0, 1, … , 14). Results for one data set (representative of the other 19) are presented in Figure 3. Figure 3a shows the posterior probabilities of P6 (additive model with equal effects) calculated for each marker pair averaged over all the transcripts. Because of symmetry, only the bottom triangle was plotted. The posterior probabilities from single-eQTL psMOM are shown on the diagonal. The two eQTL on chromosomes 2 and 3 are located between markers 5 and 7 and markers 3 and 9, respectively; psMOM identified them with fairly strong evidence.

Figure 3b shows the LOD scores derived from a standard two-QTL IM approach. The top triangle contains LOD scores for epistatic interactions; the diagonal shows LOD scores from a single QTL model; the bottom triangle shows LOD scores for the additive model. Because the simulations did not include an interaction, the top triangle correctly shows very little linkage signal. In the bottom half, however, the entire path between markers 5 and 7 on chromosome 2 and markers 3 and 9 on chromosome 3 is highlighted with the highest LOD scores occurring at the marker pair regions between chromosomes 1 and 2 and 2 and 3. In contrast, the graph from the two eQTL psMOM model gives improved localization of the true eQTL. This is promising evidence for the utility of extending psMOM to multiple eQTL.

One of the main obstacles in the multiple-eQTL model extension is the computational burden. The number of components in the mixture model grows rapidly with the number of loci under investigation. Fitting the full model can be a daunting task. A Dirichlet process mixture model for which the number of components is no longer a bottleneck has been introduced (Chen 2006) and is currently under investigation.

## APPENDIX A

Assume transcript *t* is linked to location *l*. For a backcross, there are then two distinct latent expression means, one for each genotype group, denoted by and . Consider a conditional distribution of measurements for animals with genotype 0 given by , *r* = 1,…, *n* and a prior distribution on given by | θ ∼ π_{θ}(μ),*t* =

1,…, *T*. The notation for dependence on *l* has been suppressed. Under this model, the marginal distribution of measurements is given by(A1)The same form holds for . The marginal distribution of measurements *y _{t}* is then given by

*f*

_{P1}(

*y*) =

_{t}*f*

_{P0}( )

*f*

_{P0}( ), assuming conditional independence of and given the latent expression means and .

For calculations presented here, we evaluate expression measurements on the log scale and assume a Gaussian model for *f*_{obs}(), with variance σ^{2}; π() is also Gaussian with mean μ_{0} and variance . Hence, the hyperparameters shared by all transcripts are σ^{2}, μ_{0}, and . The joint predictive density, *f*_{P0}, is then also Gaussian with mean *n* vector (μ_{0}, μ_{0},…, μ_{0}) and exchangeable covariance matrixwhere **I _{n}** is an

*n*×

*n*identity matrix and

**M**is an

_{n}*n*×

*n*matrix of ones. Further detail can be found in Kendziorski

*et al*. (2003).

## APPENDIX B

The posterior distribution of the QTL location for transcript *t* is given by(B1)In detail,where *y*_{−t} denotes the matrix of expression phenotypes with transcript *t* omitted. We assume here (second equality) that the distribution of the marker genotypes is independent of the eQTL location and latent expression means (this is analogous to the assumption made in Sen and Churchill (2001) in their Appendix A in justifying their final equality with latent expression means here corresponding to their model parameters). We further assume that *y _{t}* is conditionally independent of

*y*

_{−t}given the latent expression mean for the

*t*th transcript, μ

_{t}, and that the latent expression means are independent across transcripts (third equality). A similar derivation gives

*p*(

*y*,

*m*) =

*p*(

*y*|

_{t}*m*)

*p*(

*y*

_{−t}|

*m*)

*p*(

*m*).

Substituting these quantities into (B1), we have(B2)

Note that *p*(*y _{t}* |

*m*, γ

_{t}=

*l*) in the numerator of (B2) can be further written asHere,

*g*represents the eQTL genotype of the

*t*th transcript. We once again assume that the distribution of the marker genotypes is independent of the eQTL location (second equality). The third equality follows from the assumption that the expression levels of the

*t*th transcript are independent of eQTL locations and marker genotypes given the eQTL genotype and latent expression means (this is similar to the second equality of Sen and Churchill 2001, their Appendix A, with latent expression means corresponding to their model parameters) and that the eQTL genotype is independent of the latent expression means for transcript

*t*given the eQTL locations and markers (similar to Sen and Churchill 2001, their Appendix A, third equality). The last equality is given by the definition of

*f*

_{P1}.

In summary, we have(B3)for .

The notation π implies that integration with respect to is a two-dimensional integral over the joint distribution of the latent means in the two genotype conditions.

## Acknowledgments

The authors thank Gary Churchill, Jessica Flowers, Michael Newton, Saunak Sen, and Ping Wang for useful discussions. This work was supported in part by National Institute of General Medical Sciences (R01GM076274-01) (G.M.) and National Institute of Diabetes and Digestive and Kidney Disease (R01DK066369-03) (C.K.).

## Footnotes

Communicating editor: R. W. Doerge

- Received January 31, 2007.
- Accepted July 21, 2007.

- Copyright © 2007 by the Genetics Society of America