Comparison of linear mixed model analysis and genealogy-based haplotype clustering with a Bayesian approach for association mapping in a pedigreed population
- Golam R Dashab†1, 2,
- Naveen K Kadri†1Email author,
- Mohammad M Shariati1, 2 and
- Goutam Sahana1Email author
© Dashab et al.; licensee BioMed Central Ltd. 2012
Published: 21 May 2012
Despite many success stories of genome wide association studies (GWAS), challenges exist in QTL detection especially in datasets with many levels of relatedness. In this study we compared four methods of GWA on a dataset simulated for the 15th QTL-MAS workshop. The four methods were 1) Mixed model analysis (MMA), 2) Random haplotype model (RHM), 3) Genealogy-based mixed model (GENMIX), and 4) Bayesian variable selection (BVS). The data consisted of phenotypes of 2000 animals from 20 sire families and were genotyped with 9990 SNPs on five chromosomes.
Out of the eight simulated QTL, these four methods MMA, RHM, GENMIX and BVS identified 6, 6, 8 and 7 QTL respectively and 4 QTL were common across the methods. GENMIX had the highest power to detect QTL however it also produced 4 false positives. BVS was the second best method in terms of power, detecting all QTL except the one on chromosome 5 with epistatic interaction. Two spurious associations were obtained across methods. Though all the methods considered the full pedigree in the analyses, it was not sufficient to avoid all the spurious associations arising due to family structure.
Using several methods with divergent approaches for GWAS can be useful in gaining confidence on the QTL identified. In our comparison, GENMIX was found to be the best method in terms of power but it needs appropriate correction for multiple testing to avoid the false positives. This study shows that the issues of multiple testing and the relatedness among study samples need special attention in GWAS.
Despite many successes, genome-wide association studies (GWAS) still present major challenges. This is particularly true for samples drawn from a population with multiple levels of relatedness, such as population structure and/or familial relatedness. The efficiency of a GWAS method to detect a quantitative trait locus (QTL) depends on several factors, for example, the genetic architecture, allele frequency and heritability of the QTL, and the linkage disequilibrium with the marker. The population structure and relatedness of the samples may result in spurious associations. We applied a range of GWAS methods to map quantitative trait loci (QTL) in the simulated dataset provided by the 15th QTL-MAS workshop  and compared their efficiency in QTL detection with respected to this particular dataset.
We compared four different methods of GWAS, 1) Mixed model analysis (MMA); 2) Random haplotype model (RHM); 3) Genealogy-based mixed model (GENMIX) and 4) Bayesian variable selection method (BVS). The mixed model approach  utilizes the full relationship matrix and is the method of choice when the samples are drawn from a complex pedigreed population. The haplotype-based association methods using mixed models are generally regarded as more powerful than methods based on single markers [3, 4] since they fully exploit LD information from multiple markers. On the other hand, genealogy based clustering of haplotypes in GENMIX not only consider the local LD but also takes the history of the origin of these haplotypes . Contrary to the above three methods which analyze single markers or a few markers at a time, Bayesian variable selection  simultaneously fits multiple marker effects and avoids the problem of multiple testing. Therefore, it is useful to compare such Bayesian methods with more standard frequentist approaches where a single or a few SNPs are fitted at a time. The above-mentioned methods were compared for power, precision of location estimate, and type I error rate.
The simulated population consisted of 20 sire families, each sire was mated to 10 dams and each full-sib family had 15 progeny. The phenotype was available for 10 progeny per full-sib family i.e. a total of 2000 individuals. There were five chromosomes each with 1998 SNPs at equal distance of 0.05 cM. The four GWAS method used for association mapping are described below.
Mixed model analysis (MMA)
Where y is the vector of 2,000 phenotypes, 1 is a vector of 1s of length 2,000, μ is the general mean, g is the additive effect of the SNP and X is a vector with genotypic indicators (0, 1, or 2) associating records to the marker effect, u is the random polygenic effect with the normal distribution , where A is the additive relationship matrix and is the polygenic variance. Z is an incidence matrix relating phenotypes to the corresponding random polygenic effect, and e is a vector of random environmental deviates with the normal distribution , where is the error variance and I is the identity matrix. Testing was done using a Wald test against a null hypothesis of H0:g i =0. The significance threshold was determined using a Bonferroni correction for the number of markers tested to obtain an experiment-wise P-value of 0.05.
Random haplotype model (RHM)
Where y i is the phenotype of animal i, μ is the population mean, ui is the random polygenic effect, and are the random effects of the maternal and paternal haplotypes carried by individual i, and e i is the random residual effect as defined for MMA. The other random effect q was assumed to be normally distributed with mean zero and variances (assuming equal variance for paternal and maternal haplotypes). The significance of the haplotype substitution effect was assessed with a likelihood ratio test comparing the RHM model with a null-model containing mean, polygenic effect and random error terms but no haplotype effects. Analysis was performed using the DMU software package . Significant threshold was fixed at genome wide 5% level after Bonferroni correction and the mid-point of significant haplotypes were considered as the putative QTL positions.
Genealogy based mixed-model (GENMIX)
where y i is the phenotype of individual i, μ is the population mean, ui is as described above in the MMA; and are two haplotype effects of individual i, where h1 i and h2 i can take values 11, 12, 13, 22, 23, and 33 and , is the haplotype variance, and e i is a random residual as defined for MMA. The local genealogies were constructed using the software Blossoc (http://www.daimi.au.dk/~mailund/Blossoc/) and variance component analysis was carried out using the software DMU . The significance of the SNP association was tested using likelihood ratio test and the significant threshold was fixed at genome-wide 5% level after Bonferroni correction for multiple testing for the total number of markers.
Bayesian variable selection (BVS)
The method is based on specifying a mixture distribution for SNP effects while all SNP are fitted simultaneously in the model . It was assumed that most markers had very small effects on the trait (98% of SNP in this analysis) and only few markers (2%) had large effects. The allocation of each SNP to either of these two distributions is done using an indicator variable in Gibbs sampling. The averaged mixture indicator estimates a posterior probability for that SNP to come from the distribution with large effects, which is interpreted as the probability for presence of an associated marker or QTL. The analysis was performed using BAYZ software  and the variances of the two mixture components were estimated. The SNP with posterior probability of the mixture indicator higher than 0.10; that corresponds to a Bayes factor of 5.5 were reported as QTL. In cases where adjacent markers showed a decreasing or increasing posterior probability of association due to linkage disequilibrium, only the SNP with highest probability was reported as QTL.
Positions (cM) of identified QTL with the four methods
QTL effects estimated by single marker analysis based on linear mixed model; fitting all the detected QTL simultaneously
Allele substitution effect¤
Precision of the methods was assessed by the average of absolute differences between the positions of the simulated and the detected QTL, whenever it was identified. The QTL with the biggest effect on chromosome 1 was detected with high precision (on average ± 0.3 cM from the simulated QTL) and the epistatic QTL on chromosome 5 was detected with least precision (on an average 8 cM from the simulated QTL). In general the MMA identified the QTL with higher precision.
In our study we used additive models without considering the genetic architecture of the simulated QTL; however the methods performed well in localising the true simulated QTL. Out of the four methods employed, GENMIX performed comparatively better in QTL detection. It detected all the 8 simulated QTL and 6 were mapped accurately within a 2 cM region of the QTL region. However, it also identified 4 false positives (FP). The number of tests carried out in GENMIX was approximately 7 times the number of markers and we used Bonferroni multiple testing correction for the number of total marker but not for the total number of tests (i.e. ~7 times the number of markers) which could have resulted in increased number of false positives. Besides the number of haplotypes in a lineage goes down as we moved down the tree  which can give numerical instability. Out of these four FPs in GENMIX, two (on chromosome 4 and 5) were identified by other methods at the same location (Figure 1). Divergent approaches of GWAS picking up the same FP could be due to insufficient correction for family structure. A likely explanation is that some SNPs in these two regions were positively correlated (in linkage disequilibrium) with the QTL because of linkage (within family). It is thus not straightforward to distinguish true associations from spurious, regardless of the correction for the pedigree structure. This underlines the importance of replication study before a follow-up study can be taken up for identifying causal mutation underlying a QTL.
BVS was the second best method in terms of power to identify QTL and it had less FP compared to GENMIX. It detected all the simulated QTL except the one on chromosome 5 with epistatic interactions. BVS fits all the SNPs simultaneously and given that the first epistatic QTL was fitted in the model, there was a little chance for the second one to be significant in the model. In other words, the first QTL explains most of the variation induced by both QTL because of their dependency. Especially, this can happen if the epistasis is of additive by additive nature, where most of the epistatic variance is converted to additive . In order to confirm this, we ran the MMA for all SNP on chromosome 5 where the first epistatic QTL was already in the model. As a result, the second epistatic QTL was not detected (results not shown).
The MMA identified six QTL. The two linked QTL on chromosome 2 were both identified by MMA but only the first one (the most significant) was reported in the workshop as the second QTL was not significant when fitted along with first one in the model. On the other hand RHM detected both of them but the first QTL was mapped 10 cM downstream the true QTL.
The highest significant SNP for the multi-allelic QTL on chromosome 1 (largest QTL) in MMA was 0.7 cM away from the true position, while the other methods mapped it closer to its position. No individual SNP (bi-allelic) can be in perfect LD with this QTL (multi-allelic) which might have resulted in poor precision for this QTL in MMA.
The imprinted QTL on chromosome 4 was only detected by GENMIX and BVS. The power of detection of the QTL will decrease if the model does not reflect the true genetic architecture of the QTL. However, GENMIX and BVS methods were sensitive enough to identify the imprinted QTL, though both of them model its effect as additive.
Sahana et al.  observed very high false positives when haplotypes were considered as fixed effects in the model. Because the frequency of some haplotypes can be very low, this could result in low accuracy of estimates and result in false positive when haplotypes are fitted as fixed effect. We expected this problem can be taken care by fitting haplotypes as random where the effects of the low frequent haplotypes will be regressed towards zero. However, RHM still had very high false positive rate.
Using several methods in analysing GWA can be useful in gaining confidence on the QTL identified. Though, genealogy-based mixed model can be a powerful approach for GWAS, appropriate multiple testing correction is necessary to avoid false positives. Our study also shows that correction for pedigree relationship is not always enough to avoid spurious association arising due to family structure.
NKK and GS were supported by a grant No. 3405-10-0137, funded jointly by the Danish Ministry of Food, Agriculture and Fisheries, The Milk Levy Fund, Viking Genetics, and Nordic Cattle Genetic Evaluation. GRD acknowledges financial support by The Ministry of Science, Research and Technology of Iran.
This article has been published as part of BMC Proceedings Volume 6 Supplement 2, 2012: Proceedings of the 15th European workshop on QTL mapping and marker assisted selection (QTL-MAS). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcproc/supplements/6/S2.
- Yu J, Pressoir G, Briggs WH, Bi IV, Yamasaki M, Doebley JF, McMullen MD, Gaut BS, Nielsen DM, Holland JB: A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nature genetics. 2005, 38: 203-208.View ArticlePubMedGoogle Scholar
- Akey J, Jin L, Xiong M: Haplotypes vs single marker linkage disequilibrium tests: what do we gain?. European Journal of Human Genetics. 2001, 9: 291-300. 10.1038/sj.ejhg.5200619.View ArticlePubMedGoogle Scholar
- Morris RW, Kaplan NL: On the advantage of haplotype analysis in the presence of multiple disease susceptibility alleles. Genetic epidemiology. 2002, 23: 221-233. 10.1002/gepi.10200.View ArticlePubMedGoogle Scholar
- Sahana G, Mailund T, Lund M, Guldbrandtsen B: A New Powerful Method for Genome-wide Association Mapping Using Local Genealogies in a Mixed Model. Local Genealogies in a Linear Mixed Model for Genome-Wide Association Mapping in Complex Pedigreed Populations. PLoS ONE. 2011, 11 (6): e27061-Google Scholar
- George EI, McCulloch RE: Variable selection via Gibbs sampling. Journal of the American Statistical Association. 1993, 88: 881-889. 10.2307/2290777.View ArticleGoogle Scholar
- Madsen P, Jensen J: DMU-A user's guide, A package for analyzing multivariate mixed models. Version 6, DJF, Foulum, Denmark. 2011, [http://www.dmu.agrsci.dk/dmuv6_guide.5.0.pdf] . Release 5Google Scholar
- Scheet P, Stephens M: A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. The American Journal of Human Genetics. 2006, 78: 629-644. 10.1086/502802.View ArticlePubMedGoogle Scholar
- Mailund T, Besenbacher S, Schierup MH: Whole genome association mapping by incompatibilities and local perfect phylogenies. BMC bioinformatics. 2006, 7: 454-10.1186/1471-2105-7-454.PubMed CentralView ArticlePubMedGoogle Scholar
- Hudson RR, Kaplan NL: Statistical properties of the number of recombination events in the history of a sample of DNA sequences. Genetics. 1985, 111: 147-PubMed CentralPubMedGoogle Scholar
- Janss L: BAYZ Manual version 2.02. Janss Biostatistics, Leiden, Netherlands. 2011Google Scholar
- Hill WG, Goddard ME, Visscher PM: Data and theory point to mainly additive genetic variance for complex traits. PLoS Genetics. 2008, 4: e1000008-10.1371/journal.pgen.1000008.PubMed CentralView ArticlePubMedGoogle Scholar
- Sahana G, Guldbrandtsen B, Janss L, Lund MS: Comparison of association mapping methods in a complex pedigreed population. Genetic epidemiology. 2010, 34: 455-462. 10.1002/gepi.20499.View ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.