Mapping Quantitative Trait Loci From a Single-Tail Sample of the Phenotype Distribution Including Survival Data

Mikko J. Sillanpää; Fabian Hoti

doi:10.1534/genetics.107.081299

Abstract

A new effective Bayesian quantitative trait locus (QTL) mapping approach for the analysis of single-tail selected samples of the phenotype distribution is presented. The approach extends the affected-only tests to single-tail sampling with quantitative traits such as the log-normal survival time or censored/selected traits. A great benefit of the approach is that it enables the utilization of multiple-QTL models, is easy to incorporate into different data designs (experimental and outbred populations), and can potentially be extended to epistatic models. In inbred lines, the method exploits the fact that the parental mating type and the linkage phases (haplotypes) are known by definition. In outbred populations, two-generation data are needed, for example, selected offspring and one of the parents (the sires) in breeding material. The idea is to statistically (computationally) generate a fully complementary, maximally dissimilar, observation for each offspring in the sample. Bayesian data augmentation is then used to sample the space of possible trait values for the pseudoobservations. The benefits of the approach are illustrated using simulated data sets and a real data set on the survival of F₂ mice following infection with Listeria monocytogenes.

QUANTITATIVE trait locus (QTL) mapping methods often assume that the trait, conditionally on the effects of the QTL, follows a normal distribution. However, nonrandom missing data patterns resulting from single-tail sampling may violate this assumption. The target in single-tail sampling is to increase the expected genotype–phenotype correlation of a sample with respect to the original population parameters. By sampling (ascertaining) individuals from the right tail of the phenotype distribution, the genotype frequencies for QTL with positive phenotype effects are potentially enriched. Similarly, sampling individuals from the left tail of the phenotype distribution can increase our chances to find QTL with negative effects. Single-tail sampling may also arise from censoring or if a quantitative trait exhibits measurable values only for a portion of the individuals, i.e., there is a spike in the phenotype distribution (Broman 2003). However, due to single-tail sampling, the phenotypic variation of a sample may become too small for standard QTL mapping methods to work properly, i.e., the signal is totally masked by the error. Therefore current approaches to QTL mapping of data resulting from single-tail sampling of the phenotype distribution consider the deviation of the allele- (or genotype-) frequency distribution at the marker loci from their Mendelian expectation, use logistic regression-based analysis strategies, or combine both of these approaches (Henshall and Goddard 1999; Beasley et al. 2004; Tenesa et al. 2005). Alternatively one can apply nonparametric/semiparametric methods, rank-based statistical procedures, or a robust mixture model to analyze such data (Kruglyak and Lander 1995; Zou et al. 2002, 2003; Broman 2003; Feenstra and Skovgaard 2004). A disadvantage of these approaches is that a single-QTL model is implicitly assumed, since only a single chromosomal position is tested at a time.

As stated in Luo et al. (2005), the viability (survival) of an individual can be simply defined as a binary phenotype indicating whether an individual has survived (y = 1) or not (y = 0). For continuous survival (or failure) time data, such as time to tumor or time to death (measured in logarithmic scale), the single-tail sampling approach can be considered (Broman 2003). Alternatively, methods exist for survival phenotypes (Diao et al. 2004; Moreno et al. 2005). In controlled crosses, several methods have been designed specially to map viability loci, the gene positions that have an influence on the fitness or the survival of an individual (e.g., Vogl and Xu 2000; Luo and Xu 2003; Luo et al. 2005; Nixon 2006). In outbred populations, similar/related methods are adopted to locate the signatures of selection—the genomic regions having been under selective pressure (subject to natural or artificial selection). It is well known that (1) the variability (diversity) is reduced, (2) the linkage disequilibrium is enriched, and (3) the segregation ratios depart from their Mendelian expectations in the genomic regions at the immediate surroundings of the gene positions that influence survival. The size of the effect (selection intensity) the position has on the survival can be indirectly monitored via the extent of the above influences and their decay as a function of the genetic distance. Thus, the general rationale behind the mapping methods of such loci is in testing distorted segregation, testing linkage disequilibrium patterns, or comparing levels of genetic variability between the particular genomic position and other parts of the genome or between species. Again, as a drawback, a single-QTL model is usually implicitly assumed in these methods. Moreover, a common difficulty in applying these methods in outbred populations is that the demographic history (population growth or recent expansion) leaves kinds of local signs in the genome similar to those of selection (e.g., Schlötterer 2003).

For case–control and association studies of binary traits in human genetics, it is common that one samples affected individuals only (see Greenland 1999). There is the affecteds-only test for trios, where genotyped or haplotyped parents and their affected offspring are both collected (e.g., Falk and Rubinstein 1987; Terwilliger and Ott 1992; Lander and Schork 1994; Gauderman et al. 1999). Such a test is constructed between the cases and the controls where the individuals of the control sample, so-called pseudocontrols, “artificial controls,” or “antisibs,” are created from the genetic material that was not transmitted from the parents to the cases (their chromosomes are mirror images of the case chromosomes). The information needed to generate the chromosomes for the antisibs is obtained from the parental haplotypes by taking the complement of the genetic material of the parents that was transmitted to the affected offspring (see Figure 1). The ability to derive such a complement on the basis of the parental genotype data depends only on the parental mating type for the marker (marker informativeness); e.g., it is easy to derive complemental observations for mating type AB × CD. Some single-locus tests also utilize genotypic pseudoobservations in 1:3 proportions. By adopting the pseudocontrol approach, one can obtain well-matched controls and avoid spurious associations due to ethnic confounding, i.e., closer kinship (higher degree of background linkage disequilibrium) in the affected sample (Terwilliger and Weiss 1998).

Figure 1.—

A general representation on how chromosomes of the “antisib” (pseudoobservation) are created on the basis of parental information. The genetic material is numbered (and shaded) according to the four original grandparental sources and the recombination points are indicated with vertical lines. Note that the genetic material of the antisib is complementary (each allele originating from the other grandparent) to that of the real offspring. Also the recombinations occur in constant places and make chromosomes (inherited from the same parent) equally probable between real and pseudoobservations. Because the offspring and its antisib have different QTL genotypes ({1, 3} and {2, 4}) and they do not share any QTL alleles with each other, they also are expected (in the presence of the 1-QTL model) to have a different status in their binary phenotypes (ovals).

Here we bring this antisib idea to mapping QTL in experimental crosses (e.g., backcross and F₂) of inbred lines as well as provide a theoretical basis for applying this method to outbred populations using a multiple-QTL model. To illustrate the methodology, we extend the affected-only tests to single-tail sampling with quantitative traits. In the method, continuous-trait values are generated for all of the antisibs on the basis of Bayesian hierarchical modeling and data augmentation. For data augmentation, see Albert and Chib (1993), Rubin (1996), and Van Dyk and Meng (2001). Unlike many others who consider mapping QTL (for quantitative, viability, or survival data) from single-tail samples, we use a multiple-QTL model.

To map signatures of selection, viability, and other binary traits from case-only data (where only selected/survived individuals are in our sample and have phenotypic value one), it is straightforward to genotype individuals from the single-phenotypic group only (i.e., survived individuals). In the continuous-trait case, one can selectively sample backcross or F₂ progenies so that the (case) individuals from only one tail of the phenotype distribution are genotyped. Pseudocontrols, corresponding to observations from the other tail of the phenotype distribution, can then be created as mirror images for each case individual. The binary phenotypic value of the pseudocontrol individuals is zero and their continuous-scale observations (so-called liability values) can be predicted using data augmentation. (Note that the liability values of the case individuals are already observed.) The genotype data for these artificial observations can be created on the basis of the genotypes and linkage phases of the parents, which in inbred line-cross designs are known by definition. For example, in backcross, one can follow the principle that the mirror image of genotype AA is AB and vice versa. In F₂, the mirror images of the three genotypes AA, AB, BB are BB, AB, AA, respectively. The same principle applies for outbred crosses and populations where the parental genotypes and haplotypes are known or estimated.

MODEL

Notation:

Let us consider an inbred line-cross experiment (e.g., backcross, double haploids, or F₂) with N_gen possible genotypes. We assume that measurements of a quantitative trait and marker genotypes at N loci have been obtained from N_ind individuals sampled from a single tail of the phenotype distribution. To consider a pseudocontrol idea (Figure 1) for quantitative traits, it is easy to adopt a liability and threshold model framework (e.g., Albert and Chib 1993). Further, we assume that each sample has a mirror image (hidden observation) in the unsampled part of the phenotype distribution. Using case–control and threshold model terminology this assumption means that for each case individual (whose liability value is measured) we need one control individual (whose liability value is systematically lower/greater than that in case individuals).

Denote the phenotype and marker data of the observed offspring as (y^o, M^o) and the hidden phenotype and marker measurements of the offspring data as (y^h, M^h). From here on we refer to individual i and its unobserved counterpart as pair i. (Note that survival data, where genotypes M^h are observed for censored individuals, are a special case of this setting.) The observed and hidden phenotype vectors, $\text{[math]}$ and $\text{[math]}$ include the observed and hidden phenotypes of pair i, $\text{[math]}$ and $\text{[math]}$ respectively. Similarly, $\text{[math]}$ and $\text{[math]}$ are the observed and hidden marker matrices where the elements, $\text{[math]}$ and $\text{[math]}$ , are the coded genotypes from the set [1, … , N_gen] for pair i on marker j. To allow that some marker genotypes may be missing among the observed half of the individuals, the incomplete form of M^o is denoted by $\text{[math]}$ Note that if there are no missing marker genotypes, then M^o = M*.

To exploit the “mirroring idea” by applying data augmentation on the hidden observations (y^h, M^h), it is helpful to consider that the observed phenotypes (y^o) give rise to the discrete auxiliary variables (discrete phenotypes). Let $\text{[math]}$ denote a discrete phenotype vector where for individual i, $\text{[math]}$ Here T is a (known) discretization threshold and the binary phenotype $\text{[math]}$ obtains value 1 if the underlying (observed) continuous phenotype $\text{[math]}$ is higher than the threshold T and $\text{[math]}$ = 0 otherwise. Similarly, for the hidden observations y^h, we have a discrete vector $\text{[math]}$ where $\text{[math]}$ To conclude, in single-tail sampling, (1) the threshold T uniquely determines the proportion of individuals selected from the phenotype distribution, (2) all elements in vector z^o are either 0 or 1, and (3) all elements in vector z^h are either 0 or 1 and opposite to z^o. In the case that T is unknown, depending on which tail has been sampled we can define T as the smallest or the highest phenotype value.

Phenotype model:

We adopt the additive multiple-QTL model considered earlier by Xu (2003). This model is closely related to the model of Meuwissen et al. (2001) and can be viewed as a submodel of Hoti and Sillanpää (2006). Although it is straightforward to include also pairwise epistatic interaction terms in the design matrix of this model (see Zhang and Xu 2005; Xu 2007), we omit such extensions here. For considering other models and designs, see the discussion. In the model, let us assume that the putative QTL can be placed only at marker points. However, this is not a very restrictive assumption because in experimental designs, arbitrary map positions (putative QTL) can be included into the analysis as pseudomarkers (Sen and Churchill 2001). Given the overall mean α and the effect-specific coefficients β = (β_j,k) of marker j, the phenotypes $\text{[math]}$ s = o, h of the pair i can be expressed as $\text{[math]}$ (1)where the residuals (the phenotypes after correcting for QTL effects) are assumed to be normally distributed, $\text{[math]}$ with unknown variance $\text{[math]}$ Note that this same model is assumed for both the observed $\text{[math]}$ and the hidden $\text{[math]}$ phenotypes. The indicator variable $\text{[math]}$ if the marker observation $\text{[math]}$ equals genotype code k and is 0 otherwise. For each marker j we introduce the constraint β_j,1 = 0. Thus for a backcross or double haploids, where N_gen = 2, only a single coefficient β_j,2 at each marker is needed to capture the contrast between the two genotypes. Similar treatment for F₂, where N_gen = 3, leads to two coefficients, β_j,2 and β_j,3, that can be estimated for each marker. Here we use a random-variance model, where the genetic coefficients β_j,k, for k > 1, are assumed to be normally distributed $\text{[math]}$ with unknown variances $\text{[math]}$ In the following, we denote all unknown QTL parameters together as $\text{[math]}$ where σ² = $\text{[math]}$ is a vector of the effect-specific variances.

Key assumptions (which should be considered jointly because assumption 1 is a necessary condition for assumption 2):

Given haplotypes for parents, genotype data of the hidden observations is (generated to be) maximally dissimilar to the observed data (cf. Terwilliger and Ott 1992): When we create a mirror image (see Figure 1), we actually produce a maximally dissimilar pseudoindividual for each observation with respect to the marker data. This means that by doing so we maximize the information content of the sample (see O'Brien and Funk 2003). Note the related approaches that try to maximize the information in the sample by selecting individuals or pairs of individuals to be phenotyped on the basis of their genetic dissimilarity (Jin et al. 2004; Jannink 2005; Xu et al. 2005; Fu and Jansen 2006).
The phenotypes of the hidden observations are (generated using data augmentation) either on the left or the right side of the truncation point and the observed phenotypes: We briefly consider what is assumed at the genetic level, in the presence of the additive multiple-QTL model (1), when “mirroring” of the genotypes is performed. Now, we assume ordering of the phenotypes with respect to the threshold T: $\text{[math]}$ (2)

To see what this means at the genetic level, we substitute model (1) into both sides of Equation 2 so that α cancels out, and we obtain $\text{[math]}$ (3)where residuals $\text{[math]}$ and $\text{[math]}$ are both independently normally distributed with mean zero. Let us consider the following cases:

When the residuals are ordered $\text{[math]}$ in the same way as phenotypes of Equation 2: For such pairs i, Equation 3 imposes an ordering constraint for the sum of the QTL effects that the mirrored genotypes at the QTL loci need to fulfill. However, this ordering constraint is relaxed by the nonnegative factor $\text{[math]}$ This means that the model can cope with some proportion of phenocopies (i.e., such data individuals whose phenotype is not in agreement with the QTL model).
When the residuals are ordered $\text{[math]}$ in the opposite way as the phenotypes of Equation 2: For such pairs i, the ordering constraint is adjusted by the negative factor $\text{[math]}$ This means that for some of the phenotypes it is required that the ordering constraint is fulfilled in a tighter form. We now represent Equation 3 as $\text{[math]}$ (4)where $\text{[math]}$ is a relaxation/adjustment factor of the ordering constraint, whose sign and size depend on the rank and the difference of the two residuals, respectively. Further understanding of this question requires simulation studies that are not in the scope of this article.

Assumptions 1 and 2 together imply that each sample has a mirror image (hidden observation) in the unsampled part of the phenotype distribution.

Hierarchical model:

In Bayesian analysis, the aim is to obtain an estimate for the posterior distribution of the model parameters given the data, p(θ, y^h, M^h, M^o, z^o, z^h | y^o, M*). This can be achieved by using Markov chain Monte Carlo (MCMC) methods, exploiting the fact that the posterior is proportional to the joint distribution of the parameters and the data, p(θ, y^h, M^h, M^o, z^o, z^h, y^o, M*). By adopting suitable conditional independence assumptions (leading to the graphical model of Figure 2), the joint distribution of the parameters and the data can be presented as $\text{[math]}$ where the likelihood can be factorized as $\text{[math]}$ The functional forms of p( $\text{[math]}$ | θ, $\text{[math]}$ ) and p( $\text{[math]}$ ) are normal densities of the residuals $\text{[math]}$ and $\text{[math]}$ of model (1) with mean zero and variance $\text{[math]}$ (see, e.g., Sillanpää and Arjas 1998).

Figure 2.—

A graphical display of the hierarchical structure of the model. The data and the priors are presented using boxes and the unknown variables are shown as ovals. The hierarchical dependency and the conditional independence assumptions are visible in the graph. The dependency can be either deterministic (dashed line) or stochastic (solid line) and its direction is indicated with an arrow.

Constraining priors:

This model includes an exceptionally high number of constraining priors, which introduce restrictions into the (MCMC) sampling scheme and take care of the consistency between variables. Their specific forms are given below. The prior for the discretized hidden observations is $\text{[math]}$ where p( $\text{[math]}$ ). Here, $\text{[math]}$ and $\text{[math]}$ Moreover, the prior for the discretized “nonhidden” observations is $\text{[math]}$ where $\text{[math]}$ and $\text{[math]}$ The markers of the pseudoobservations are created according to the mirroring prior $\text{[math]}$ where $\text{[math]}$ and $\text{[math]}$ Also the indicator function prior $\text{[math]}$ is used to ensure that the complete marker observations are compatible with the observed data.

Other priors:

In inbred line-cross data, the prior used to handle missing observations can be presented as a Markov chain $\text{[math]}$ For the actual forms of the transition probabilities $\text{[math]}$ ) in various designs, see Jiang and Zeng (1997) and Sillanpää and Arjas (1998). The prior for the QTL parameters can be factorized as p(θ) = p(α)p(β | σ²)p(σ²)p( $\text{[math]}$ ), where p(α) ∝ 1, $\text{[math]}$ , and $\text{[math]}$ Further, p(β_j,k | $\text{[math]}$ ) is the density function of a normal distribution with mean zero and variance $\text{[math]}$ and $\text{[math]}$ is the Jeffreys scale invariant prior having most of the support (mass) in values near zero. Also, for the residual variance, we choose $\text{[math]}$ The use of effect-specific variance components together with Jeffreys' prior is well justified because the prior adaptively shrinks QTL variances at unlinked positions to zero—which then leads to the positioning of QTL with nonnegligible effects (see Xu 2003; Hoti and Sillanpää 2006). Note that Meuwissen et al. (2001) fitted a common variance for all coefficients at single locus, which, however, does not lead to an equally sparse solution. It is good to know that even if Jeffreys' prior in this context seems to work extremely nicely, theoretically the posterior is improper because Jeffreys' prior has an infinite amount of mass near zero (e.g., Hopert and Casella 1996; ter Braak et al. 2005). One way to avoid the theoretical problem is to specify a small positive number as a lower bound for the parameter in the prior, which we, however, did not apply here.

Parameter estimation:

We use a MCMC algorithm (e.g., Casella and George 1992; Chib and Greenberg 1995) to estimate the posterior distribution of the unknown model parameters. Here we assume that the truncation point T of the original population is known, equals the smallest (or the highest) phenotypic value in the data, or has been successfully estimated before the analysis. Also if the phenotypic mean $\text{[math]}$ of the original population is available, we can utilize it as a starting value for α; otherwise we initialize α to zero (i.e., we set $\text{[math]}$ ). We use nonzero starting values for the variances so that nonzero values are proposed for all effects—all positions initially explain the phenotype. In the following, we outline the MCMC sampling scheme used for continuous traits. For survival data, if genotypes M^h are observed/available for censored individuals, we can omit the generation of mirror images in steps 1 and 3. For binary traits, step 4 below is replaced by the step “the updating liabilities for a binary trait” found in earlier articles (Kilpikari and Sillanpää 2003; Hoti and Sillanpää 2006):

Specify initial values $\text{[math]}$ , (β_j,k = 0, $\text{[math]}$ = 0.5, j = 1, … , N, k = 2, … , N_gen), $\text{[math]}$ = 0.5; initialize the missing genotypes in M^o from their prior distribution; and generate mirror images M^h conditionally on M^o.
Update the QTL parameters (θ) needed in the phenotype model according to the Gibbs sampling scheme outlined elsewhere (Xu 2003; Hoti and Sillanpää 2006).
Update missing values in M^o and the corresponding mirror images in M^h using a separate Metropolis–Hastings step for each individual and each marker. Propose the genotypes from their prior distribution p(M^o). The acceptance ratio contains only the likelihood p(y^o, y^h | θ, M^o, M^h). Note, however, that each change in M^o also changes the value of p(y^h | θ, M^h) because M^h contains the mirror image of the proposed value.
Update y^h using Gibbs sampling. A new $\text{[math]}$ is sampled (for each individual separately) from p( $\text{[math]}$ | θ, M_i^h, $\text{[math]}$ > T) if the observations y^o have been collected from the right tail of the phenotype distribution and $\text{[math]}$ is sampled from p(y_i^h | θ, M_i^h, y_i^o ≤ T) if y^o are from the left tail. The fully conditional posterior distributions are $\text{[math]}$ and $\text{[math]}$ , where $\text{[math]}$ is the predictive mean and ϕ(·) is the cumulative distribution function of the standard normal distribution. Similarly, as in data augmentation algorithms for binary traits (Albert and Chib 1993; Hoti and Sillanpää 2006), this Gibbs sampling step requires sampling from a truncated normal distribution (for algorithms, see Devroye 1986). Note that the condition $\text{[math]}$ uniquely determines the values for $\text{[math]}$ and $\text{[math]}$ which again imply the constraint for the possible values of $\text{[math]}$ .
Repeat steps 2–4 until a prespecified number of rounds have been reached.

DATA ANALYSIS

In the following we present example analyses and comparisons of our method under different sampling schemes (random, single-tail, and two-tail sampling), using simulated backcross data in cases of unlinked and linked QTL. We consider both the average performance (assessed by analyzing 50 or 100 data replicates) and performance under a single realization of a data set (assessed by analyzing several single data sets with small heritability in each). For reasons why correction methods based on truncated normal distribution as “incomplete-data” likelihood are not used here, see the discussion. We used unrealistically large (QTL) heritabilities and small sample sizes in our example data sets to reduce computation time when analyzing data replicates. Albeit this treatment may appear to be unrealistic, the analyses presented here arguably correspond to the analyses with smaller heritabilities and larger samples. One can use existing power tables to find rough correspondence between the two cases (e.g., Van Ooijen 1992; Carbonell et al. 1993; Beavis 1998). For example, using a traditional approach and backcross data, the probability of success to find a QTL with heritability 0.05 in a sample of 400 is roughly comparable to that to find a QTL with heritability 0.16 in a sample of 100 (Lander and Botstein 1989). Additionally, we illustrate the performance of our method with survival data and censored observations, using previously analyzed real F₂ mice data that have some degree of randomly missing genotypes (Broman 2003).

SIMULATION ANALYSIS OF UNLINKED QTL

Simulated data:

The performance of the new approach was tested using simulated data, which were generated in two phases. First, linked marker data for a population of 250 backcross individuals were generated using the QTL Cartographer software (Basten et al. 1996). The produced offspring data consisted of 33 markers that span the area on three 100-cM-long chromosomes. Each chromosome had 11 equidistant markers, one every 10 cM. To generate phenotypes, we selected 3 markers (nos. 3, 17, and 30) of 33 as QTL with additive genetic effects as β₃ = 3, β₁₇ = −2, β₃₀ = 1, respectively. For each individual i, a quantitative phenotype y(i) was generated using an additive genetic model $\text{[math]}$ (5)where the indicator functions $\text{[math]}$ and $\text{[math]}$ take value 1 if individual i has genotype AB at positions 3, 17, and 30, respectively. The additive error e(i) was generated from the normal distribution with mean zero and variance 4. This resulted in a heritability ∼0.5. Note that there were no missing values in the marker data.

Sampled subdata:

In the following, to distinguish between the original simulated data set (250 individuals) and a smaller sample (∼40 individuals), which is obtained by sampling according to the sampling scheme, we call these two alternatives data and subdata, respectively.

Analyses:

To demonstrate the efficiency of the proposed pseudocontrol approach and to address the sampling variation around the estimates, six different sampling schemes were compared by analyzing simulated subdata replicates. Sampling of the individuals and analysis of the resulting subdata were repeated 100 times under each sampling scheme, using 100 different simulated data sets. All simulated data sets included the same genotype data of 250 offspring (see above) but a different set of phenotypes was generated each time by using the same generating model (Equation 5). (The method is expected to be more sensitive to sampling variation due to phenotypes than due to genotypes.) In each repetition, subdata were sampled from the same phenotype distribution of the 250 individuals according to different sampling schemes. Analyses using six different schemes of sampling from the phenotype distribution were considered: (A) an analysis using a random subdata sample, (B) an analysis using a sample from both the left and the right tails, (C) an analysis using the right tail sample without doing any correction with respect to truncation, (D) the pseudoobservation analysis using the right tail sample, (E) an analysis using the left tail sample only without doing any correction with respect to truncation, and (F) the pseudoobservation analysis using the left tail sample. Analyses D and F were carried out for all the subsamples, using the approach presented in the model section. Analyses A–C, and E were carried out using our implementation of the approach of Xu (2003); see details from Hoti and Sillanpää (2006). Analyses C, D and E, F correspond to the same analysis with and without generating pseudoobservations, respectively. The truncation threshold T, which determines the sampled individuals in the subdata, was defined as $\text{[math]}$ for the right-tail sampling analyses (C and D) and as $\text{[math]}$ for the left-tail sampling analyses (E and F). In the two-tail sampling analysis (B), both thresholds T_l and T_r were used. Due to the resampling of the phenotype, the sample sizes used in analyses C–F varied in each repetition. Therefore, to maintain comparability between the schemes (in each repetition), the size of the random sample (A) was chosen to equal the mean of the sample sizes of the left and right tail samples (rounded upward). In the two-tail sampling analysis (B), we randomly sampled half of the individuals (rounded upward) from both tails. The sample sizes varied in the range [35, 51] with median 41, which closely coincides with the theoretical expectation (16% of the samples in a normal distribution should be beyond one standard deviation, here corresponding to 0.16 × 250 ≈ 40 individuals).

Results:

We implemented the methods using Matlab software on a personal computer. The posterior estimation (of the effects) for each of the 100 repetitions was based on 10,000 Markov chain Monte Carlo cycles. In each MCMC run the first 1000 initial cycles were discarded from the chain as “burn-in” rounds and thinning of 10 was applied (by saving the values at every 10th cycle) to reduce autocorrelation between the samples. Due to the rather simple data generation model, the MCMC sampler converged rapidly in all 100 cases.

Instead of using the estimated effect size directly to summarize the results, we use a standardized form of the effect size, because then the selected QTL threshold 0.1 (giving definition for QTL as in Hoti and Sillanpää 2006) is directly comparable/applicable to other traits (sampling schemes) and marker data. For a backcross, at marker j, the standardized effect is $\text{[math]}$ where $\text{[math]}$ is the empirical standard deviation of the genotypes at marker j and $\text{[math]}$ is the empirical standard deviation of the phenotype (calculated from augmented data in D and F). Following Hoti and Sillanpää (2006), we define an indicator variable for the event, that the absolute value of the standardized effect size is larger than the given QTL threshold 0.1. This enables us to estimate the QTL occupancy probability (as a function of the posterior distribution of the QTL effects) in models such as those of Xu (2003), which do not originally include model selection indicators. Thus, we present the results in the form of the posterior probability of the QTL occupancy P(|θ_j| > 0.1 | data), which we calculate as the proportion of MCMC rounds where |θ_j| > 0.1. To summarize the posterior QTL occupancy over the 100 repetitions of data for each of the six sampling schemes, we calculated the mean value of the estimated posterior QTL occupancy probability at each locus j as $\text{[math]}$ by taking the average over the 100 subdata analyses in the different sampling schemes (Figure 3). In Figure 3, the corresponding means of the standardized effects (calculated over the MCMC samples in each subdata analysis where |θ_j| > 0.1) are shown using the curve. At each repetition, the calculation of the QTL occupancy and the posterior mean of the estimated (standardized) effect size was based on 900 effective MCMC samples. [The reason for not taking the average over repetitions with respect to the correctly identified QTL was that we wanted especially to monitor the magnitude of the signals (cf. Broman and Speed 2002).]

Figure 3.—

Unlinked QTL: 100 data sets. (A–F) The posterior QTL occupancy probability (indicated with a bar at each marker locus j) averaged over 100 subdata analyses in the different sampling schemes. The corresponding mean of the standardized effect (calculated only over such MCMC rounds where |θ_j| > 0.1) is shown using the curve. The true positions of the simulated QTL are marked with an asterisk. The right column indicates the sampling scheme used/part of the distribution sampled and arrows indicate the utilization of the pseudoobservations (mirroring idea) in the analysis. Shown are the analysis of the random subdata sample (A), the two-tail sample analysis (B), the direct and the mirror analysis of the sample from the right tail of the phenotype distribution (C and D), and the direct and the mirror analysis of the sample from the left tail of the phenotype distribution (E and F).

In Figure 3, A–F, the QTL occupancy probability (bars) and the mean standardized QTL effect (curve) summarize the QTL evidence at each locus over the 100 repetitions. Note that the same scales of the y-axis are used throughout. As expected, the analyses without doing any correction for truncated data performed very badly, showing practically no signals (Figure 3, C and E). It becomes evident from the graphs that the single-tail sampling analyses with pseudoobservations (Figure 3, D and F) have clearly more power than the analyses based on random subdata samples (Figure 3A) or the analyses that do not utilize the pseudoobservations (Figure 3, C and E). This is in the sense that analyses (Figure 3, D and F) on average show higher or elevated signals (QTL occupancy probabilities) around the true positions (3, 17, and 30). Surprisingly, in the case of the unlinked QTL, one can even conclude that the single-tail sampling analyses with pseudoobservations (Figure 3, D and F) showed power comparable to the analysis with two-tail sampling (Figure 3B). On the other hand, the mean of the estimated standardized effect size is practically at the same level in most of the analyses, which indicates that standardized effects seem to be (on average) comparable across the different analyses. In Figure 3, D and F, note the negligible bias of the QTL position (around locus 17) in the opposite direction from the direction of sampling. The simulated effect size at position 30 was apparently very small because the position (on average) stayed undetected in most of the cases.

Heritability estimation:

In schemes A–F in Table 1, the posterior mean heritability was estimated using the formula $\text{[math]}$ where $\text{[math]}$ is the empirical phenotypic variance (from augmented data in D and F), $\text{[math]}$ is the residual variance at round t, and r is the total number of MCMC rounds (after burn-in). The mean and the standard deviation were used to summarize the distribution of the estimated heritability over the 100 repetitions for each of the six sampling schemes. The random-sampling analysis (scheme A) seemed to give (on average) slightly small heritability estimates and the standard deviation (sampling variance) was very large. The heritability estimates were negligible (or negative) in almost all analyses based on the single-tail sampling without doing any correction with respect to truncation (schemes C and E). (The negative values arise because the residual variance was not restricted in the prior to be smaller than the phenotypic variance.) On the other hand, the heritability was overestimated in most of the analyses and the sampling variance was relatively small for analyses using single- or two-tail sampling (schemes B, D, and F) and pseudoobservations.

View this table:

TABLE 1

Heritability estimates

To enable comparison to more traditional methods (based on hypothesis testing) and to obtain an understanding of their potential performance on these data, Figure 4 shows the deviation of the marker genotype frequencies from their expected values in a typical realization of the data after right-tail sampling. One can clearly see the dependence of the frequencies over linked loci and the potential difficulty to control false positives by choosing the significance threshold. By looking at these data, it is easy to understand also the value of the multiple-QTL model.

Figure 4.—

Unlinked QTL: a typical realization of the data obtained by sampling only the subset belonging to the right tail of the phenotype distribution of 250 simulated backcross individuals. At each locus the bar indicates the deviation of the number of individuals with genotype AB from the Mendelian expectation. The true positions of the simulated QTL are marked with an asterisk.