Genome-wide association analyses of the 15th QTL-MAS workshop data using mixed model based single locus regression analysis

Background The mixed model based single locus regression analysis (MMRA) method was used to analyse the common simulated dataset of the 15th QTL-MAS workshop to detect potential significant association between single nucleotide polymorphisms (SNPs) and the simulated trait. A Wald chi-squared statistic with df =1 was employed as test statistic and the permutation test was performed. For adjusting multiple testing, phenotypic observations were permutated 10,000 times against the genotype and pedigree data to obtain the threshold for declaring genome-wide significant SNPs. Linkage disequilibrium (LD) in term of D' between significant SNPs was quantified and LD blocks were defined to indicate quantitative trait loci (QTL) regions. Results The estimated heritability of the simulated trait is approximately 0.30. 82 genome-wide significant SNPs (P < 0.05) on chromosomes 1, 2 and 3 were detected. Through the LD blocks of the significant SNPs, we confirmed 5 and 1 QTL regions on chromosomes 1 and 3, respectively. No block was detected on chromosome 2, and no significant SNP was detected on chromosomes 4 and 5. Conclusion MMRA is a suitable method for detecting additive QTL and a fast method with feasibility of performing permutation test. Using LD blocks can effectively detect QTL regions.


Background
Recently, the high-density single nucleotide polymorphism (SNP) arrays have been developed for almost all domestic animals, which offer the prerequisite of genome-wide association study (GWAS), a more powerful approach for high-resolution mapping of loci controlling phenotypic traits in domestic animals [1]. In GWAS, two basic designs of resource population have been widely used for association analysis, one is the case-control design with unrelated individuals, and the other is the family-based design with pedigree structure. Corresponding to these two designs, different approaches for association analysis have been proposed. However, there is no clear evidence showing general superiority of one approach over others. In farm animals, family based design is more relevant because of complex pedigree structure in almost all animal populations. In our previous GWAS study [2], we employed a mixed model based single locus regression analysis (MMRA) to test the association between SNPs and milk production traits in dairy cattle. We found this method was more powerful than the TDT-based single locus regression analysis. To further verify its performance in terms of power and type I error, we applied it to the common dataset provided in the 15th QTL-MAS workshop.

Methods
The simulated population consisted of 3,220 individuals in two generations. The first generation consisted of 20 sires and 200 dams, which were assumed to be unrelated. Each sire mated with 10 dams and each dam produced 15 progenies, leading to a total of 3,000 individuals in the second generation. Of the 15 progenies of each dam, 10 were phenotyped for a continuous trait. All of the 3,220 individuals were genotyped for 9,990 SNP markers distributed on 5 chromosomes without missing. Each chromosome had a size of 1 Morgan (M) and carried 1,998 evenly distributed SNPs.

Variance component estimation
We applied the software DMU (Version 6, release 5.0) [3] to estimate the variance components of the simulated trait, which would be used in the subsequent association analysis, based on the following model Where y is the vector of phenotypes of the 2,000 phenotyped individuals, μ is the overall mean, a is the vector of the residual polygenic effect with a ∼ N(0, Aσ 2 a ) (where A is the additive genetic relationship matrix and σ 2 a is the additive genetic variance), Z is the incidence matrix of a, and e is the vector of residual errors with e ∼ N(0, Iσ 2 e ) (where I is a unit matrix and σ 2 e is the residual error variance).

Genotype quality control
We removed the 1,000 progenies without phenotypes off the genotype data, and we calculated the minor allele frequency (MAF) for each SNP for the remained 2,220 individuals (2,000 progenies and 220 parents). We found that 2,879 SNPs were homozygous (MAF = 0) for all the tested individuals and additionally 715 SNPs had a MAF less than 0.03. These SNPs were removed and 6,396 SNPs remained for the subsequent analyses.

Association analysis
The mixed model based single locus analysis [2,4] was performed based on the following linear mixed model: where y is the vector of phenotypes of the 2000 phenotyped individuals, μ is the overall mean, × is the vector of the SNP genotype indicators which takes values 0, 1 or 2 corresponding to the three genotypes 11, 12 and 22 (assuming 2 is the allele with a minor frequency), b is the regression coefficient of phenotypes on SNP genotypes (i.e., the substitution effect of the SNP), a is the vector of the residual polygenic effect with a ∼ N(0, Aσ 2 a ) , Z is the incidence matrix of a, and e is the vector of residual errors with e ∼ N(0, Iσ 2 e ). For each SNP, the estimate of b and the corresponding sampling variances Var( ∧ b) can be obtained via mixed model equations (MME), and a Wald chi-squared statistic b 2 /Var( b) with df =1 was constructed to examine whether the SNP is associated with the trait.

Statistical inference
For the analyses above, the permutation method was adopted to adjust for multiple testing from the number of SNP loci detected. In our method, the phenotypes were permuted 10,000 times against the genotype and pedigree data and the empirical distribution of the Wald chisquared statistic under the null hypothesis (no association existed between any SNP and the trait in genome-wide level) was obtained using the largest Wald chi-squared statistic value across all SNPs from each permuted dataset. The threshold value for declaring a significant association was determined by choosing the 95th percentile of the empirical distribution, i.e., we declared a significant SNP at a 0.05 genome-wide significance level if its raw value of the Wald chi-squared statistic was larger than the empirical threshold value.
For the significant SNPs, linkage disequilibrium (LD) in term of D' between them was quantified using Haploview [5] and the LD blocks were defined by the criteria of Gabriel et al. [6] with default parameters.

Association analysis results
The estimates of σ 2 a and σ 2 e are 24.82 and 58.65, respectively, so that the heritability estimate is 0.30 approximately. The profile of the raw p values (from the chidistribution and in terms of -log10 p) of all tested SNPs is shown in Figure 1. By using simply Bonferroni correction, we detected 119 significant SNPs of 0.05 genomewide significance level (raw p values < 7.82E-6). However, by using permutation test, we detected 82 significant SNPs of 0.05 genome-wide significance level (raw To further pinpoint the relationship among the detected SNPs, we analysed the LD levels in terms of D' between the significant SNPs (Figures 2, 3 and 4) for chromosomes 1-3, respectively. Through the criteria of Gabriel et al. [6] with default parameters in Haploview [5], we defined 5 LD blocks on chromosome 1, which harbour 4 to 10 significant SNPs, and 1 LD block on chromosome 3, which harbour 10 significant SNPs with 6 outside. No block was detected on chromosome 2. The LD patterns show that these significant SNPs links to each other in different LD levels.

Comparison of the significant SNPs with the simulated QTN
On chromosome 1, there is one simulated QTN located at 2.85cM (No.57), which had the largest effect among all  simulated QTNs. We detected 63 significant SNPs on this chromosome. However, the true QTN at 2.85cM has a MAF of 0 and was discarded after quality control, and the adjacent SNP at 2.90cM (No.58), which has the largest estimated effect among all significant SNPs, is accordingly considered as the putative QTN. Although a large number of pseudo significant SNPs were identified on this chromosome, the LD levels between the most effective SNP and other 62 significant SNPs (Figure 2) showed that 47 of them are in strong LD (D'>0.5) with it. This suggests that the simulated QTN may be surrogated by a suite of "ghost" QTNs nearby due to high LD level.
On chromosome 2, there are two simulated QTNs in coupling linkage phase located at 81.90cM (No.3638) and 93.75cM (No.3875), respectively. We detected 3 significant SNPs, the first is exactly at 81.90cM, and the second (No.3662, at 83.10cM) is in strong LD (D' = 0.97) with the first one ( Figure 3). But the third one No.3916 is at 95.80cM and is 2.05cM away from the second simulated QTN, while the LD level between them is strong (D'=0.69, Figure 3).
On chromosome 3, there are two simulated QTN in repulsion linkage phase located at 5cM (No.4100) and 15cM (No.4300), respectively. However, the first simulated SNP on this chromosome also has a MAF of 0 and was discarded after quality control. Of the 16 significant SNPs detected, 10 are harboured in the LD block covering the interval between 4.75cM and 5.65cM with an average LD level of 0.97 (D'), in which SNP No.4101 is just adjacent to the first true QTN. The second simulated QTN is 1.10cM away from the significant SNP (No.4322) and the LD level between them is strong (D'=0.93).
The one simulated imprinting QTN on chromosome 4, and 2 simulated epistatic QTNs on chromosome 5 were not detected by our analysis. This is because our method dose not account for both imprinting effect and epistatic effect. Our method needs to be further improved to account for interaction effects between SNPs and imprinting effects from parents.

Comparison of the significant SNPs with the those with high effects estimated via Bayesian approaches
To further validate significant SNPs identified herein, we compared the most promising SNPs detected with those with highest effects estimated via Bayesian approaches (BayesA, BayesB and BayesCπ) reported in our another analysis on prediction of genomic breeding values for the same data set [7]. Since the results from the three Bayesian approaches are similar and BayesCπ performed best, we only compare with BayesCπ here. Specifically, on chromosome 1, the most effective SNP (No.58) identified by MMRA is exactly the same as that by BayesCπ. In all, most of findings herein are largely consistent with those with highest effect estimates via BayesCπ. This further demonstrates that the Bayesian approaches (particularly BayesCπ) could also sever as tools for QTL mapping, as suggested by Fan et al. [8].

Computing time
All analyses were implemented through Fortran programs and performed on an octal-core Linux Server (Intel Xeon E5504 2.00GHz; 48.00GB RAM). The time needed was about 1.5 minutes for one permutation analysis. The 10,000 permutations were performed through 8 threads, each was assigned 1,250 permutations. So, the total computing time was about 31 hours. This shows that MMRA is a fast method with feasibility of performing a large number of permutations.

Conclusion
Our results herein show that the MMRA method is suitable for detecting additive QTL, and it is a fast method with feasibility of performing permutation test. And the LD region on chromosome 3 can effectively integrate significant SNPs for QTL region detection. However, we detects only one true additive QTN (No.3638), two SNPs (No.58 and No.4101) close to two true additive QTNs (No.57 and No.4100) with many false positives, which remains to be further investigated and the MMRA method needs to be further improved to account for other nonadditive effects. 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.