Abstract
Background
Genomewide dense markers have been used to detect genes and estimate relative genetic values. Among many methods, Bayesian techniques have been widely used and shown to be powerful in genomewide breeding value estimation and association studies. However, computation is known to be intensive under the Bayesian framework, and specifying a prior distribution for each parameter is always required for Bayesian computation. We propose the use of hierarchical likelihood to solve such problems.
Results
Using double hierarchical generalized linear models, we analyzed the simulated dataset provided by the QTLMAS 2010 workshop. Markerspecific variances estimated by double hierarchical generalized linear models identified the QTL with large effects for both the quantitative and binary traits. The QTL positions were detected with very high accuracy. For young individuals without phenotypic records, the true and estimated breeding values had Pearson correlation of 0.60 for the quantitative trait and 0.72 for the binary trait, where the quantitative trait had a more complicated genetic architecture involving imprinting and epistatic QTL.
Conclusions
Hierarchical likelihood enables estimation of markerspecific variances under the likelihoodist framework. Double hierarchical generalized linear models are powerful in localizing major QTL and computationally fast.
Background
Genetic analyses in livestock studies are generally based on information from pedigrees and molecular markers. Traditionally, a kinship matrix can be calculated using the pedigree data, which can be used in a generalized linear mixed model (GLMM) to estimate breeding values. By including genetic marker information, genomic estimated breeding values (GEBV) can be obtained taking into account the information from these markers, and also quantitative trait loci (QTL) can be mapped by associating genotypes at a certain locus to the phenotype observations.
Dense marker genotypes along genome can now be affordably obtained due to new and efficient methods for typing single nucleotide polymorphism (SNP) markers. The dense SNP maps have made genomewide association (GWA) studies popular for gene detection. Classic GWA methods [1], commonly applied to study genetic diseases in humans, are based on simple repeated single marker tests across the genome. To achieve more powerful mapping and better prediction, a unified model including all the SNPs in the genome is preferred. Such models have been estimated using Bayesian methods, implemented by Markov chain Monte Carlo (MCMC) techniques that are computationally demanding [25]. Lee and Nelder developed the double hierarchical generalized linear model (DHGLM) in the likelihoodist framework [6]. DHGLM enables estimation of markerspecific variances using a fast iterative algorithm without specifying any prior distributions [7]. The likelihoodist way of estimation is conducted through a likelihood function named hierarchical likelihood (hlikelihood) [8].
The aim of this paper is to map QTL and report GEBV for the simulated dataset provided by QTLMAS 2010 workshop. We employ a unified analysis via the hlikelihood and model the data using DHGLM. GEBV are calculated from the estimated marker effects, and QTL are mapped by the estimated markerspecific variances.
Methods
Data
The dataset used in this paper was simulated for the QTLMAS 2010 workshop (Poznań, Poland). A pedigree consisting of 3226 individuals in 5 generations (F_{0}  F_{4}) was simulated, where F_{0} contains 5 males and 15 females. Each female was mated once and gave birth to about 30 progeny. Two traits were simulated, where one is quantitative (QT), and the other is binary (BT). Young individuals in F_{4} (individuals 2327 to 3226) had no phenotypic records. The genome was assumed to be about 5 × 10^{8} bp long, consisting of 5 chromosomes, each of which contained about 1 × 10^{8} bp. Each individual was genotyped for 10031 biallelic SNPs in the genome.
Models
DHGLM provides a unified analysis for both QTL mapping and genomic breeding value estimation. Similar to BayesA, the data are modeled on two levels, i.e. both the phenotypic mean and the variance are modeled with random effects. For a quantitative trait, the phenotype y (n × 1 vector) is postulated as a random effect model
y = Xβ + Zg + e (1)
where g ~ N(0, diag(λ)) are the SNP effects, λ = (λ_{1}, λ_{2},..., λ_{m})′ are the variances of the SNP effects, and the residuals e ~ N(0, σ^{2}I). The fixed effects β included an intercept and the sex effect in our application to reduce the residual errors. The SNP variances λ are modeled as
log λ = 1a + b (2)
with an intercept a and normally distributed random effects b. The genomic estimated breeding value (GEBV) for individual i is computed as . QTL can be scanned using the markerspecific variances λ. For a binary trait, the mean of y, is modeled by the same linear predictor Xβ + Zg through a logit link function.
For the markerspecific variances, the correlated random effects, b, follow a multivariate normal distribution with a mean of zero and a variancecovariance matrix , where m is the number of SNPs and k, l are the SNP indices. When ρ = 1, all the SNPs have a constant variance (GLMM); when ρ = 0, the SNPs are assumed to be independent (DHGLM); and for 0 <ρ < 1, the correlation between two SNPs is a monomial function of ρ, which is referred to as the smoothed DHGLM [10]. We propose the use of smoothed DHGLMs since it reduces the noise in markerspecific variance estimates and highlights the signals of QTL. ρ, regarded as a spatial correlation parameter, was chosen to be 0.9 in this paper, which nicely shrank the SNPs with zero effect.
The overall phenotypic variance can be expressed as
where is the variance of z._{j} (the jth column of Z) across individuals. These variance values can be directly calculated from the data. The contribution (heritability) of a particular SNP is expressed by [4].
Fitting algorithm
According to the extended likelihood principle, inference of the random SNP effects g should be drawn through the hlikelihood, fixed effects β through the marginal likelihood, and variance components λ, σ^{2} and through the adjusted profile likelihood [11]. However, for efficient estimation, we propose to initialize variance components and iterate the following steps until convergence [7],
• Solve the following WLS problem for and ĝ,
Where and . The subscript M stands for ‘mean’.
• Update σ^{2} by fitting the deviance residuals using an interceptonly gamma GLM and prior weight w_{M} = (1 – q_{M})/2, where are the residuals of (4), and are the diagonal elements of The subscript 1 and 2 stand for individuals (1 to n) and SNPs (n + 1 to n + m), respectively.
• Solve the following WLS problem for â and ,
where , , z = log λ + (d_{M2} – λ)/λ is linearized λ in a gamma GLM with a log link, and L satisfies LL′ = A. The subscript D stands for ‘dispersion’.
• Update by fitting the deviance residuals using an interceptonly gamma GLM and prior weight w_{D} = (1 – q_{D})/2, where ê_{D} are the last m residuals of (5), and q_{D} are the last m diagonal elements of .
Results and Discussion
Estimation of SNP effects
The effect of each SNP was estimated by a smoothed DHGLM with spatial correlation parameter ρ = 0.9 for both traits (Figure 1). For both traits, DHGLM shrank the estimated SNP effects for the loci not linked to main QTL towards zero; meanwhile, the SNPs linked to QTL were highlighted. Note that the extent of shrinkage depends on the spatial correlation parameter ρ. ρ = 0.9 was specified in our analyses since it produced better shrinkage and smoothing results for this particular dataset.
Figure 1. Estimated SNP effects The SNP effects were estimated using the smoothed DHGLM with spatial correlation parameter ρ = 0.9. The dashed vertical lines indicate the chromosome borders.
QTL mapping
Moving from the mean part to the variance (dispersion) part of the models, markerspecific variances were estimated and used to detect QTL (Figure 2). The overall variance component estimate from GLMM can be regarded as a reference value (smoothed DHGLM with ρ = 1), which was estimated using the hglm package [12] in R [13]. The 6 peaks for QT, corresponding to SNP number 163, 952, 2719, 3957, 4493 and 5492, were QTL which had values greater than the overall variance component estimate. The two strong QTL for BT had similar positions as two for QT. Other small peaks lower than the reference line were suggestive QTL. Simulated main QTL were precisely mapped. The two main epistatic QTL pairs for QT were detected as two single QTL due to the very short distance between interacting SNPs. Heritability for QT and BT was calculated for detected QTL and suggestive QTL (Table 1). 30.35% and 33.42% of the phenotypic variance were explained for QT and BT, respectively. Phenotypes of QT and BT are significantly correlated with a Spearman’s rank correlation coefficient of 0.2431. However, jointmodeling both traits were not considered in this paper.
Figure 2. QTL detection using estimated markerspecific variances The markerspecific variances were estimated using the smoothed DHGLM with spatial correlation parameter ρ = 0.9. The dashed horizontal line is the overall variance of SNP effects estimated by GLMM. The peaks higher than this line were detected as QTL, and other small peaks below were suggestive QTL. Simulated QTL are also shown as vertical bars with their heights proportional to the variances they explained. For nice visualization, simulated variances are 1/50 magnified for QT and 1/1500 magnified for BT.
Table 1. Estimated heritability of the detected QTL and suggestive QTL for QT and BT.
GEBV
GEBV were estimated for all the 3226 individuals in the pedigree. Examining outsample prediction, we compare the GEBV with the true breeding values (TBV) for the young individuals (23273226) without phenotypic records (Figure 3). The correlation coefficients between GEBV and TBV were 0.60 for QT and 0.72 for BT. The linear regression slopes were 0.41 for QT and 0.62 for BT. Accuracy of GEBV was worse for QT than for BT mainly because three imprinted QTL were simulated only for QT, and QT had a more complicated genetic architecture.
Figure 3. Scatterplots of GEBV against TBV for the young individuals without phenotypic records The GEBV were estimated using the smoothed DHGLM with spatial correlation parameter ρ = 0.9. The values are not scaled on the same mean.
Conclusions
DHGLM were shown to be an efficient and reliable approach for both QTL mapping and genomic selection. Since DHGLM can be estimated by iterating interlinked GLMs, the execution time is greatly shortened comparing to the Bayesian computation. On a Macintosh laptop with a 2 GHz processor and 4 GB memory (1067 MHz), it took about 1020 minutes, depending on starting values, to obtain our results using our implementation in R. No priors are required for parameters in DHGLM. Main QTL mapped via DHGLM showed very good accuracy though some QTL with small effects were shrunk or smoothed down. An R package iQTL has been implemented and is available on RForge: https://rforge.rproject.org/R/?group id=845 webcite.
List of abbreviations used
bp: base pair; DHGLM: double hierarchical generalized linear model; DNA: deoxyribonucleic acid; GEBV: genomic estimated breeding values; GLM: generalized linear model; GLMM: generalized linear mixed model; GWA: Genomewide association; hlikelihood: hierarchical likelihood; HGLM: hierarchical generalized linear model; MCMC: Markov chain Monte Carlo; QTL: quantitative trait locus/loci; QTLMAS: quantitative trait loci and marker assisted selection; REML: restricted maximum likelihood; SNP: single nucleotide polymorphism; TBV: true breeding values; WLS: weighted least squares.
Competing interests
No competing interest to declare by any of the authors.
Authors contributions
XS, LR and ÖC initiated the study. XS analyzed the simulated common dataset of the QTLMAS 2010 workshop and drafted the paper. LR initiated the smoothed version of double hierarchical generalized linear models. XS, LR and ÖC worked on the revision together and approved the final manuscript.
Acknowledgements
Xia Shen is funded by a Future Research Leaders grant from the Swedish Foundation for Strategic Research (SSF) to Örjan Carlborg. Lars Rönnegård is funded by the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (FORMAS). François Besnier is acknowledged for sharing his IBD calculation program to validate our results by variance component methods.
This article has been published as part of BMC Proceedings Volume 5 Supplement 3, 2011: Proceedings of the 14th QTLMAS Workshop. The full contents of the supplement are available online at http://www.biomedcentral.com/17536561/5?issue=S3.
References

Cantor RM, Lange K, Sinsheimer JS: Prioritizing GWAS results: A review of statistical methods and recommendations for their application.
Am. J. Hum. Genet. 2010, 86:622. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Meuwissen THE, Hayes BJ, Goddard ME: Prediction of total genetic value using genomewide dense marker maps.
Genetics 2001, 157(4):18191829. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Xu S: Estimating polygenic effects using markers of the entire genome.
Genetics 2003, 163:789801. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Xu S: An empirical Bayes method for estimating epistatic effects of quantitative trait loci.
Biometrics 2007, 63:513521. PubMed Abstract  Publisher Full Text

Yi N, Xu S: Bayesian LASSO for quantitative trait loci mapping.
Genetics 2008, 179:10451055. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lee Y, Nelder JA: Double hierarchical generalized linear models (with discussion).

Lee Y, Nelder JA, Pawitan Y: Generalized linear models with random effects: unified analysis via hlikelihood. Chapman & Hall/CRC; 2006.

Lee Y, Nelder JA: Hierarchical generalized linear models (with discussion).

Yi N, Banerjee S: Hierarchical generalized linear models for multiple quantitative trait locus mapping.
Genetics 2009, 181(3):11011113. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rönnegård L, Lee Y: Hierarchical generalized linear models have a great potential in genetics and animal breeding. In Proc. WCGALP. Leipzig, Germany; 2010.

Lee Y, Nelder JA, Noh M: Hlikelihood: problems and solutions.
Statistical Computing 2007, 17:4955. Publisher Full Text

Rönnegård L, Shen X, Alam M: hglm: a package for fitting hierarchical generalized linear models.

R Development Core Team: R: A Language and Environment for Statistical Computing. [http://www.Rproject.org] webcite
R Foundation for Statistical Computing, Vienna, Austria; 2009.