Hierarchical likelihood opens a new way of estimating genetic values using genome-wide dense marker maps

Background Genome-wide 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 genome-wide 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. Marker-specific 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 marker-specific variances under the likelihoodist framework. Double hierarchical generalized linear models are powerful in localizing major QTL and computationally fast.

Results: Using double hierarchical generalized linear models, we analyzed the simulated dataset provided by the QTLMAS 2010 workshop. Marker-specific 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 marker-specific 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 genome-wide 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 [2][3][4][5]. Lee and Nelder developed the double hierarchical generalized linear model (DHGLM) in the likelihoodist framework [6]. DHGLM enables estimation of marker-specific 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 marker-specific variances.

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 where g~N(0, diag(l)) are the SNP effects, l = (λ 1 , λ 2 ,..., λ m )′ are the variances of the SNP effects, and the residuals e~N(0, s 2 I). The fixed effects b included an intercept and the sex effect in our application to reduce the residual errors. The SNP variances l are modeled as with an intercept a and normally distributed random effects b. The genomic estimated breeding value (GEBV) for individual i is computed as z g ij j j m  = ∑ 1 . QTL can be scanned using the marker-specific variances l. For a binary trait, the mean of y, is modeled by the same linear predictor Xb + Zg through a logit link function.
For the marker-specific variances, the correlated random effects, b, follow a multivariate normal distribution with a mean of zero and a variance-covariance matrix , where m is the number of SNPs and k, l are the SNP indices. When r = 1, all the SNPs have a constant variance (GLMM); when r = 0, the SNPs are assumed to be independent (DHGLM); and for 0 <r < 1, the correlation between two SNPs is a monomial function of r, which is referred to as the smoothed DHGLM [10]. We propose the use of smoothed DHGLMs since it reduces the noise in marker-specific variance estimates and highlights the signals of QTL. r, 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 s s s y z where s z. j 2 is the variance of z. j (the j-th 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 h-likelihood, fixed effects b through the marginal likelihood, and variance components l, s 2 and s b 2 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 b b  and ĝ, stand for individuals (1 to n) and SNPs (n + 1 to n + m), respectively.
• Solve the following WLS problem for â and b  ,

Estimation of SNP effects
The effect of each SNP was estimated by a smoothed DHGLM with spatial correlation parameter r = 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 r. r = 0.9 was specified in our analyses since it produced better shrinkage and smoothing results for this particular dataset.

QTL mapping
Moving from the mean part to the variance (dispersion) part of the models, marker-specific 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 ρ

Quantitative Trait
Marker index SNP effect = 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, joint-modeling both traits were not considered in this paper.

Marker index
Variance of SNP effect ρ = 0.9 ρ = 1 Additive QTL Chromosome borders Figure 2 QTL detection using estimated marker-specific variances The marker-specific variances were estimated using the smoothed DHGLM with spatial correlation parameter r = 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.  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.

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 10-20 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 R-Forge: https://r-forge.r-project.org/R/? group id=845.