The null distribution of stochastic search gene suggestion: a Bayesian approach to gene mapping

With the advancement of computers over the past few decades comes a rise in the application of Bayesian methods and Markov-chain Monte Carlo methods (MCMC) to genetic data (see [1–6] for examples). Many of these MCMC methods capitalize on Bayesian priors to model additional complexity in the problems facing complex disease mapping. Bayesian priors allow prior information about the loci in question to statistically enter the inference regarding the current data set. It is well known that with Bayesian methods, the more data that is collected, the less impact the prior has. Also, using uninformative priors when modeling large data sets generate posteriors that can be approximated by the same asymptotic normal distribution of the corresponding maximum likelihood estimates [7], leading to Bayesian and Frequentist agreement under these conditions.

But how do Bayesian methods behave when there is not a gene to find? Or, what impact does the prior have on inference if, in reality, there are no genes to find? We use simulated data of Problem 3 for the Genetic Analysis Workshop 15 to evaluate the performance of stochastic search gene suggestion (SSGS), a Bayesian method introduced by Swartz et al. [3], in a region of chromosome 6 with no association to disease.

SSGS combines hierarchical priors to model alleles within loci with the conditional logistic regression likelihood to model the probability of transmission to a diseased child in case-parent triad data. The full details are given in Swartz et al. [3], and we give a brief overview of the method here.

where g*|G _m , G _f represents all possible pairings of transmitted and non-transmitted alleles consistent with the parental genotypes and RR [(g _m , g _f )] is a relative risk function defined previously as a conditional logistic regression function [3].

Using generalized transmission-disequilibrium test (GTDT) coding as described by Schaid [8], we must omit a reference allele for model identifiability. Using calculations from Thomas et al. [4], we select the most prevalent allele as our reference allele. Then, using l = 1,..., L to index the loci, and a = 1,..., A _l -1 to index the alleles within a locus (and A _l denotes the maximum number of alleles at locus l), each β _la refers to an allele main effect for allele a at locus l, with the usual interpretation as the log relative risk for transmission to the case. Note that using GTDT coding assumes additive allelic effects [3].

Next, we assume a multivariate normal prior for the main effects, and hierarchical latent indicators to stochastically search through the main effects that impact the probability of transmission [3, 9]. We denote the first latent indicator as the vector λ= (λ₁,..., λ _L ). Each element of λhas a Bernoulli distribution with probability p _l , where p _l is the prior probability that locus l is associated with the disease. Similarly, we define a partitioned allele indicator vector, conditional on the loci selected as α|λ= (α₁,..., α_L) with each partition ${\alpha }_l=\left({\alpha }_{l1},\mathrm{...},{\alpha }_{l{A}_l-1}\right)$ indicating the alleles at locus l associated with disease, and each element in α₁ has a Bernoulli distribution with probability q _la as its prior, where q _la is the prior probability that locus l and allele a at locus l are associated with the disease. This dual hierarchical prior structure restricts the stochastic search to models that select alleles within loci.

Conditional on the indicators, we can define the prior for main effects. Let γ _la = λ _l α _la . Therefore, γindicates the allele main effects when both the locus and allele are indicated. We can define the conditional prior for the coefficients, β|γ, as

π (β|γ) = MVN (0, D _γ RD _γ ),

where c is large, τ is small, and R is either an identity matrix or a covariance matrix defined by genetic correlation. (For more details concerning c and τ and the flexibility of SSGS using these parameters, see Swartz et al. [3] and references therein.)

To define R, we start with a blocked matrix, L, whose diagonal blocks represent within loci covariance and off-diagonal blocks represent between loci covariance [3]. We model the within-locus covariance by modeling the probability of an allele's presence as a multinomial distribution. The covariance is then defined in the usual way for a multinomial distribution, using the normalized allele frequencies as cell probabilities for the distribution. To model between locus covariance, each element is simply the allele-wise linkage disequilibrium value: ${\delta }_{i_a{j}_b}=p_{i_a{j}_b}-p_{i_a}{p}_{j_b}$. Once the elements of the blocks are defined, we set R = L^-1. More details are given in Swartz et al. [3].

Our posterior distribution is intractable, and therefore, we use MCMC simulations to sample from the posterior distribution. The parameters λand αcan be calculated using a Gibbs sampler, while the βvalues are updated by locus using a multivariate Metropolis Hastings step [3]. Software to implement SSGS is available at http://www.epigenetic.org/Linkage/ssgs-public/.

We are mainly interested in the marginal posterior of the γ _la values given the data. We use the proportion of iterations for each γ _la = 1 to estimate the posterior probability of each allele being associated with transmission to the case. Once the posterior probabilities are calculated, we use the median model decision rule developed by Barbieri and Berger [10]: select genes with posterior probability of inclusion greater than or equal to 0.5.

Recall that conditional logistic regression models the probability of transmission to the affected offspring, and can include other covariates, such as environmental factors in the model. We compare SSGS with the TDT, implemented in TDTEX program in the Statistical Analysis for Genetic Epidemiology release 5.0 (SAGE 5.0) [11]. TDTEX performs McNemar's test of counts of transmitted versus non-transmitted alleles, which is a powerful χ² text of association that cannot include additional covariates.

The methods ranked differently when analyzing the dense microsatellite loci. SSGS without including a dependence structure indicated all alleles at dense microsatellite Locus 2, having posterior probability for all alleles greater than 99% across all prior values of p and q, yet did not indicate any alleles at dense microsatellite locus 1 or 3. The conditional logistic regression maximum-likelihood estimates also indicated all alleles at dense microsatellite 2 with p-values less than 0.001, without indicating any alleles for Loci 1 or 3. When incorporating covariance structure defined as a function of allele frequencies and linkage disequilibrium through R, SSGS only indicated one allele at dense microsatellite 2, (allele 14) with posterior probability of 64% when the prior values for p and q were set equal to 50%. For other smaller values of prior probabilities p and q, SSGS did not indicate any alleles at any loci. TDTEX indicated both dense microsatellite 2 (p-value < 3.3 × 10^-28) and dense microsatellite 3 (p-value < 0.024) to be associated with the simulated disease locus. Here, at some point, all four methods found the marker closest to the disease locus. However, SSGS excluding a dependence structure and the conditional logistic regression gave the strongest evidence for the location of the disease locus, defining the narrowest region; while SSGS including the dependence structure gave the weakest signal, and was also the least robust to different prior probabilities.

The results of this paper differ from the results of our previous work [3]. Previously, we concluded that including the dependence structure in the model was best. However, in these data sets, we see that excluding the dependence structure performed better than including it when looking for a signal. So why did the dependent prior do so well in our previous research? To answer this question we looked at the LD pattern between markers among the dense microsatellites. When looking at the D' values for the dense microsatellites (calculated by ldmax [14]), we found that the LD falls in the range of 0 to 0.24. Meanwhile, the D' values (again calculated by ldmax) for the data for our previous work [3] was much higher with values ranging from 0.25 to 0.52. Therefore, we see a greater improvement gained by including a covariance structure that is a function of allele frequencies and LD when the LD between markers is higher.

We have revealed important features for SSGS, a new method introduced in Swartz et al. [3]. We found that in null regions, SSGS can be sensitive to the magnitude of prior probability of inclusion, especially when not including dependence structure in the model. To calibrate the prior, this study suggests using a prior probability of inclusion of 0.25 to control for false positives, while not erroneously excluding true loci. In general, we showed that using smaller prior probability of inclusion for loci and alleles reduces false-positive values, particularly in null regions. Most importantly, we found that the usefulness of including genetic dependence between markers in the model depends on the amount of LD between the markers. Therefore, before including the dependence structure in the matrix, it is very important to examine the LD pattern between the marker loci. If the D' values for the markers is less than 0.25, then incorporating a covariance structure based on the LD seems to shrink the signal drastically toward the null and increase the sensitivity to the prior probabilities of inclusion. Therefore, it is better to exclude the dependence structure from the model when using SSGS to minimize false positives in null regions. However, if the LD between the markers is high (greater than 0.25), as in those analyzed in Swartz et al. [3], we see a clear advantage in including dependence structure based on LD into the model.