 Proceedings
 Open access
 Published:
A Bayesian genome screening of maximum number of drinks as an alcoholism phenotype with the new HasemanElston method
BMC Genetics volume 6, Article number: S116 (2005)
Abstract
Common human disorders, such as alcoholism, may be the result of interactions of many genes as well as environmental risk factors. Therefore, it is important to incorporate gene × gene and gene × environment interactions in complex disease gene mapping. In this study, we applied a robust Bayesian genome screening method that can incorporate interaction effects to map genes underlying alcoholism through its application to the data of the Collaborative Studies on Genetics of Alcoholism provided by Genetic Analysis Workshop 14. Our Bayesian genome screening method uses the regressionbased stochastic variable selection, coupled with the new HasemanElston method to identify markers linked to phenotypes of interest. Compared to traditional linkage methods based on singlegene disease models, our method allows for multilocus disease models for simultaneous screening including both main and interaction (epistatic) effects. It is conceptually simple and computationally efficient through the use of Gibbs sampler. We conducted genomewide analysis and comparison between scans based on microsatellites and singlenucleotide polymorphisms. A total of 328 microsatellites and 11,560 singlenucleotide polymorphisms (by Affymetrix) on 22 autosomal chromosomes and sex chromosome were used.
Background
Alcohol dependence is a complex disorder that is influenced by many genetic and environmental factors. Identifying genes associated with alcohol dependence is critical to understand its etiology and to develop efficient methods for prevention and treatment. However, this effort has been hampered by the complexity underlying alcohol dependence: rather than there being one or a few major genes affecting alcohol dependence, it is likely that multiple genes interact with each other, together with environmental factors, to affect susceptibility to alcohol dependence. In this paper, we describe analyses of the Collaborative Study on the Genetics of Alcoholism (COGA) data (Problem 1), using selfreported "maximum number of drinks consumed in a 24hour period" (denoted by M) as a quantitative trait, to map genes underlying alcohol dependence. The measure M is closely related to alcoholism diagnosis and provides a quantitative measure for alcohol dependence. For genome screens for this trait, we use the modified HasemanElston regression method [1] along with the Bayesian variable selection methods developed by Oh [2] to locate susceptibility genes for alcohol dependence and assess their epistatic effects.
The HasemanElston method and its derivatives allow one to apply linear regression methods in linkage analysis. For each sibling pair, these methods use the number of alleles identical by descent (IBD) at each marker as the explanatory variable and a statistic measuring similarity of the quantitative traits in the sibling pair, squared difference, or crossproduct, as the response variable.
In practice, the number of markers and their possible epistatic effects are often larger than the number of observations (patients or sibpairs), where the design model is referred as being "supersaturated". As we often have hundreds of markers to consider, we must deal with the problem of multiple testing in this context. Besides, if one would like to take epistasis into account, the number of tests can easily exceed tens of thousands. Performing hypothesis tests for linkage for all of these possibilities without appropriate adjustment of multiple comparisons can lead to the identification of spurious genetic effects or the masking of real effects. The supersaturated nature of the design model also makes the conventional best subset model selection methods [3] practically infeasible. To overcome this infeasibility, efficient and robust Bayesian variable selection methods for multiple regressions were proposed by Oh [2] utilize the stochastic search variable selection (SSVS) methodology developed by George and McCulloch [4] in gene mapping. SSVS uses Gibbs sampling to sample from the posterior distribution of the model space and obtain the "best" subsets from the sampled posteriors. Therefore, it is capable of considering a large number of candidate variables without evaluating all possible models. Oh et al. [5] applied this approach to genetic linkage studies and successfully estimated main and epistatic effects. George and McCulloch's method [4], yet, is known to be sensitive to the choice of the priors. Oh [2] modified and extended SSVS, focusing only on statistics robust to the choice of the prior. These methods were shown to be computationally more efficient than the original SSVS [4] in handling many hundreds of thousands epistatic effects and gene × gender interactions through the simulation study by Oh [2].
In this study, we apply the method developed by Oh [2] to the COGA data of Genetic Analysis Workshop 14 (GAW14) to screen disease susceptibility genes related to the quantitative phenotype, "maximum number of drinks consumed in a 24hour period", and compare the results from singlenucleotide polymorphisms (SNPs) with those from microsatellites. Because SNPs offer denser and more automated genotyping than microsatellites, we expect that SNPbased linkage analysis may have better performance in the mapping of disease loci.
Methods
HasemanElston method
The original HasemanElston method [6] is a general modelfree method for testing linkage between candidate markers and quantitative trait loci on a sample of sibpairs. This method involves the regression of the squared trait difference D^{2} = (X_{1}  X_{2})^{2}) in pairs of siblings on the number of alleles shared IBD between each sib pair at a given marker. Although the original HasemanElston method is simple, robust, and computational inexpensive, it may ignore information contained in the observed bivariate data. In fact, the squared mean corrected trait sum of sib pairs (S^{2} = (X_{1} + X_{2}  2μ)^{2}) may provide additional information on the genetic effect [7]. This observation has led to a number of methods to modify the original HasemanElston method by combining both D^{2} and S^{2} in linkage analysis to improve statistical power. One of the simplest methods was proposed by Elston et al. [1], which uses CP = (S^{2}  D^{2})/4 = (X_{1}  μ)(X_{2}  μ) as the response variable in the regression. It essentially averages both regression coefficient estimates (of D^{2} and S^{2}) with equal weights. An additional advantage of using CP as the response variable is that it may be more normally distributed than D^{2} and S^{2}.
Bayesian genome screening
Assume that we observe m markers along the genome. Among these m markers, some may display some evidence of tight linkage to genes with large effects and others only with weak effects. They may exhibit main effects and/or epistatic effects. In our regression setup, the CP s are treated as the responses and the IBD status at each marker X_{ j }is the candidate explanatory variable. Then the observed CP value of sibpair i can be described by the following linear model
where n is the number of sib pairs and the errors ε~N(0, σ^{2}) are assumed to be independent. A subset model is represented by a binary vector γ = (γ_{1}, ..., γ_{ m }), where γ_{ j }= 1 or γ_{ j }= 0 represents the presence or absence of variable j in the model. A large effect of α_{ j }indicates that marker j has the evidence for linkage and indicates the epistatic effect between markers j_{1} and j_{2}. The marker effects α_{ j }(j = 1,..., m) are given the prior of the mixture of normal and pointmass distributions conditional on the indicators γ_{ j }, α_{ j }γ_{ j }~γ_{ j }N(0,c^{2}) + (1  γ_{ j })δ_{ o }. Similarly, the epistatic effects (j_{1}, j_{2} = l,..., m) with j_{1} <j_{2} are given . Hence, if the effect is absent in the model choice (γ_{ j }= 0), we can simply omit that factor when building the model. If γ_{ j }= 1, the magnitude of the effect is large and then a nonzero estimate should be included in the model and its posterior distribution is largely determined by the data. The prior γ should represent beliefs about the relationship between factors. If only main effects but no epistatic effects are considered, the prior of γ can be simply set as prob(γ) = p^{k}, where k is the number of ones in γ. However, when epistatic effects are considered in model selection, the model space becomes enormous and the independent relationship may not be the proper assumption anymore. In this case, we incorporate Chipman's [8] hierarchical prior structure into our model space. Then the probability that the term is present may take on four different values, depending on the values of the pair ,
The choice of the values, (p_{00},p_{01},p_{10},p_{11}) represents the belief of the relationship between factors (main effects), and (epistatic effects). If we believe that an epistatic effect without main effects is quite unlikely, we can simply set the prior as (p_{00},p_{01},p_{10},p_{11}) = (0,0,0,p). Alternatively, we can use (p_{00},p_{01},p_{10},p_{11}) = (0,p_{1},p_{2},p_{3}) to relax these conditions.
With an appropriate prior distribution on σ^{2}, one can obtain the posterior distribution of γ using Gibbs sampling. Then, by examining the posterior distribution of γ, one can identify the optimal model with the highest posterior. However, since our interest is to find the markers having evidence for linkage, we proceeded to obtain the marginal posterior probabilities of different factors and rank their relative frequencies. That is, the marginal posteriors of all main effects and epistatic effects are obtained and used to rank these effects. We view the ranking as a measure of relative importance of all the factors, i.e., the degree of evidence for these factors linked to the disease genes. Our previous studies showed that the rankings of these factors are very robust [2], especially for those with the highest marginal probabilities, which are the focus of our study. Probabilities much less than 0.5 for both main and epistatic effects represent the belief that relatively few terms are active. This assumption is quite reasonable for our applications of gene screening, because there are only a few markers linked to the disease genes out of hundreds or thousands of markers. Even though the choice of p may be arbitrary, the results are robust to its setting. We should note that the priors for p are chosen partly to facilitate the computing.
In microsatellites and SNPs, there are 328 and 7,826 markers considered, respectively. Therefore, to consider epistatic effects, we needed to include 53,957 factors for microsatellites and about 30 millions factors for SNPs in the models. We set the prior as p = 0.05 for microsatellies and p = 0.005 for SNPs. There are about 1,499 sib pairs in the overall sample. We use 1,433 sib pairs for microsatellites analysis and 1,322 sib pairs for SNPs analysis in our study. The number of alleles shared IBD is obtained using GENEHUNTER (v.2.1) for microsatellites and MERLIN for SNPs for each sib pair at each marker. Because it is impractical to track the complete posterior of γ, only the marginal posterior of each marker is obtained. In our analysis, we use (p_{00},p_{01},p_{10},p_{11}) = (0,p,p,p) to give all the factors an equal opportunity to be included in the final model as long as one of main effects is in the model.
Results
In the full sample of cases, the quantitative trait M ranges from 0 to 160. Five individuals are in the highest threshold class (M > 128 drinks). The highest reported M is 160. The use of the log transformation minimizes their impact on the regression analysis, which can be inflated by selfreport [9]. As reported in Saccone [9], there is a close relationship between diagnosis and M.
In each analysis, the Markov chain Monte Carlo (MCMC) sampler was run for 100,000 cycles after discarding the first 2,000 cycles for the burnin period. Because MCMC samplers arise from recursive draws, they produce correlated samplers from the posteriors. Therefore, the chains are thinned (one iteration in every 10 cycles is saved) to reduce serial correlation in the stored samples. The total number of samples kept in the postBayesian analysis is 10,000. It takes ~2 hours for microsatellites and ~6 hours for SNPs to generate each sample with JAVA programs on a Linux cluster using 2.4GHz Intel processors. Table 1 displays the rankings of the marginal posteriors for the main effect and epistatic effect screening obtained from microsatellites and SNPs. Table 1 clearly shows that the high posteriors comparatively concentrate on chromosome 4 both for microsatellites and SNPs, which agrees with Saccone [9]. To localize the specific regions of the strong evidence for linkage on chromosome 4, we extract the relative frequencies of the marginal posteriors. Figure 1 shows the comparisons of SNPs vs. microsatellites in the regions having evidence of linkage, and they have similar patterns. For gene × gene interaction, microsatellites are able to locate epistatic effects between chromosomes 8 and 15, and between chromosomes 10 and 17, whereas SNPs locate epistatic effects only between sex and markers on chromosome 23 (sex chromosome). However, the marginal posterior probabilities are not strong enough to support the evidence of epistatic effects.
Discussion
In this study, we have compared the genomewide linkage analyses based on microsatellites and SNPs. Our methods located the main effects of markers both from microsatellites and SNPs and produced similar patterns between them. However, the results for epistatic effect screening are less consistent and revealing. This might be purely because these epistatic effects are weak in nature and further research in this area is warranted.
Conclusion
Bayesian genome screening methods provide a powerful and efficient tool in identifying potential markers and their epistatic effects. They are very effective because they are able to conduct searches over the entire model space; while the frequentist's best subset model selection procedure is constrained by computing power required to examine all candidate models. In addition, Bayesian genome screening methods can work on problems with many more candidate variables, which is essential to consider when epistatic effects are studied. When one tries to locate the epistatic effects, the number of covariates (factors) easily far outnumbers the sample size. Most traditional linkage methods do not work under this condition because they often assume a singlegene model and test effects one at a time. By using the prior structures that reflect the relationship among the candidate variables, our general approach can accommodate a large number of candidate markers as well as their epistatic effects by evaluating all factors simultaneously. We were able to locate markers on chromosome 4 that show the strong evidence of linkage with alcoholism related to quantitative phenotype, "maximum number of drinks consumed in a 24hour period", both from microsatellite and SNP scans and weak evidence for epistatic effects.
Abbreviations
 COGA:

Collaborative Study on the Genetics of Alcoholism
 GAW14:

Genetic Analysis Workshop 14
 IBD:

Identicalbydescent
 MCMC:

Markov chain Monte Carlo
 SNP:

Singlenucleotide polymorphism
 SSVS:

Stochastic search variable selection
References
Elston RC, Buxbaum S, Jacobs KB, Olson JM: Haseman and Elston revisited. Genet Epidemiol. 2000, 19: 117. 10.1002/10982272(200007)19:1<1::AIDGEPI1>3.0.CO;2E.
Oh C: Robust Bayesian variable selection. PhD dissertation. 2003, State University of New York at Stony Brook, Applied Math and Statistics Department
Miller A: Subset Selection in Regression. 2002, Boca Raton: Chapman & Hall/CRC
George EI, McCulloch RE: Variable selection via Gibbs sampling. Journal of American Statistical Association. 1993, 88: 881889. 10.2307/2290777.
Oh C, Ye KQ, He Q, Mendell NR: Locating disease genes using Bayesian variable selection with the HasemanElston method. BMC Genet. 2003, 4 (Suppl 1): S6910.1186/147121564S1S69.
Haseman JK: The investigation of linkage between a quantitative trait and a marker locus. Behav Genet. 1972, 2: 319. 10.1007/BF01066731.
Wright F: The phenotypic difference discards sibpair QTL linkage information. Am J Hum Genet. 1997, 60: 740774.
Chipman H: Bayesian variable selection with related predictors. Can J Stat. 1996, 24: 1736. 10.2307/3315687.
Saccone N: A genome screen of maximum number of drinks as an alcoholism phenotype. Am J Med Genet. 2000, 96: 632637. 10.1002/10968628(20001009)96:5<632::AIDAJMG8>3.0.CO;2#.
Acknowledgements
Supported in part by NIH grant R01 GM59507 and NSF grant DMS 0241160.
Author information
Authors and Affiliations
Corresponding author
Additional information
Authors' contributions
CO participated in the design of the study, performed the analysis, and drafted the manuscript. SW helped to obtain IBD values for linkage analysis. SW, NL, LC, and HZ participated in the design and the discussion of the study, and the preparation of the manuscript. All authors read and approved the final manuscript.
Rights and permissions
Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Oh, C., Wang, S., Liu, N. et al. A Bayesian genome screening of maximum number of drinks as an alcoholism phenotype with the new HasemanElston method. BMC Genet 6 (Suppl 1), S116 (2005). https://doi.org/10.1186/147121566S1S116
Published:
DOI: https://doi.org/10.1186/147121566S1S116