Abstract
In single-step analyses, missing genotypes are explicitly or implicitly imputed, and this requires centering the observed genotypes using the means of the unselected founders. If genotypes are only available for selected individuals, centering on the unselected founder mean is not straightforward. Here, computer simulation is used to study an alternative analysis that does not require centering genotypes but fits the mean 
of unselected individuals as a fixed effect. Starting with observed diplotypes from 721 cattle, a five-generation population was simulated with sire selection to produce 40,000 individuals with phenotypes, of which the 1000 sires had genotypes. The next generation of 8000 genotyped individuals was used for validation. Evaluations were undertaken with (J) or without (N) 
when marker covariates were not centered; and with (JC) or without (C) 
when all observed and imputed marker covariates were centered. Centering did not influence accuracy of genomic prediction, but fitting 
did. Accuracies were improved when the panel comprised only quantitative trait loci (QTL); models JC and J had accuracies of 99.4%, whereas models C and N had accuracies of 90.2%. When only markers were in the panel, the 4 models had accuracies of 80.4%. In panels that included QTL, fitting 
in the model improved accuracy, but had little impact when the panel contained only markers. In populations undergoing selection, fitting 
in the model is recommended to avoid bias and reduction in prediction accuracy due to selection.
- centering genotype covariates
- estimated breeding value
- genomic prediction
- selection
- single-step
- GenPred
- Shared Data Resources
- Genomic Selection
In pedigree-based analyses, the expected value of breeding values is zero. In order to achieve similar properties in whole-genome analyses, marker genotype covariates are often transformed. When all individuals are genotyped, it has been shown that inference on genotype effects does not depend on how the covariates are transformed (Strandén and Christensen 2011). However, when data includes genotyped and nongenotyped individuals, inference on marker effects from single-step analyses may depend on how the covariates are transformed (Fernando et al. 2014). In single-step analyses using marker effects models, the breeding values of nongenotyped individuals are partitioned into components representing the prediction of nongenotyped individuals conditional on their genotyped relatives and an independent deviation (Fernando et al. 2014). The prediction of nongenotyped individuals conditional on their genotyped relatives is done based on best linear prediction, which requires the first moments to be known without error. This is straightforward if the mean of the genomic breeding value is zero in the absence of selection. Centering the observed genotype covariates using what their means would be in the absence of selection would result in genomic breeding values with null means. However, such genotype covariate means are typically unavailable. Fernando et al. (2014) proposed a solution for the marker effects model that involves fitting an additional fixed covariate that estimates the mean 
of the linear component of the genotypic value, which is denoted by 
in Equation 1 below, in a population where selection is absent. Using this approach, even when there is selection, the selection process can be ignored, because the analysis is conditional on the data used for selection (Goffinet 1983; Gianola and Fernando 1986; Im et al. 1989; Sorensen et al. 2001). In Markov chain Monte Carlo (MCMC) analyses, centering results in better mixing (Strandén and Christensen 2011), reducing the number of iterations required to obtain converged genomic predictions. In practice, centering the entire matrix of genotype covariates, including the observed and imputed genotypes, using their sample means can be done in addition to the Fernando et al. (2014) approach of fitting 
This type of centering of the entire genotype matrix does not affect inference on marker effects.
The same issue with centering of the observed genotype covariates that we described above for the marker effects model is also implicit for the single-step breeding value model (single-step GBLUP), and a similar solution was proposed by Vitezica et al. (2011). In their proposed solution, the observed genotype covariates are centered using their means; in addition, the genomic covariance matrix is corrected for the change in the mean breeding value of the genotyped individuals (Vitezica et al. 2011). It was shown in that paper that this is equivalent to fitting the change in breeding value due to selection as a random effect.
Here, we use simulated data to compare the accuracy and bias in genomic prediction applied to populations under selection with and without centering the entire matrix of genotype covariates, and with and without fitting 
as a fixed effect. Further, we show that when the observed genotype covariates are centered using means calculated from selected individuals rather than means from all individuals, the meaning of 
changes from the mean of unselected individuals to become the mean breeding value in selected individuals, as claimed by Vitezica et al. (2011).
Methods
Theory
To simplify the presentation of the genetic model without loss of generality, we will assume that the unconditional expectation of the phenotypic value for all individuals is the same. Let 
denote the row vector of genotypes for individual i, which is often coded as −1,0,1. Then, under additive gene action, the genotypic value, 
which is the expected phenotypic value of an individual with genotypes 
, can be written as
(1)where β is the value of 
when 
and 
is the vector of substitution effects. The scalar β and the vector 
are constants, but 
will be a random variable because of randomness in 
owing to the randomness in the genotypes for a randomly sampled individual. Note that the expected value of the linear component 
of the genotypic value in Equation 1 is 
where 
which may not be equal to zero. Thus, it is customary to write the model for the genotypic value, as can be derived from Equation 1, as follows:
(2)where 
is a constant, representing 
and 
is a random variable that has null expectation, which is the breeding value predicted in a pedigree-based BLUP (PBLUP) evaluation. When genotypes are observed and used in a genomic analysis, they may be transformed or coded by subtracting their expectations, 
from the observed values, 
In both Equations 1 and 2, 
is the same substitution effect for locus j. The intercepts in these models, however, are different. In Equation 1, the intercept is β, and it is the value of 
when 
In Equation 2, on the other hand, the intercept is 
and it is the value of 
when 
More generally, 
is not known, so genotypes are coded by subtracting a different vector 
from the observed genotypes as 
Still, 
is the substitution effect for locus j, but the intercept will change to become 
which is the value of 
when 
Thus, as more rigorously shown in (Strandén and Christensen 2011), inference about 
does not depend on how the genotypes are coded. A simpler but rigorous proof is given in the Appendix of this paper.
In single-step analyses, where some individuals are not genotyped, the missing genotypes are imputed either implicitly (Legarra et al. 2009) or explicitly (Fernando et al. 2014) using best linear prediction. Let 
denote the matrix of genotypes for individuals that were genotyped. Then, the genotypes of the individuals with missing genotypes are imputed as
where 
is the matrix of pedigree-based additive relationships between the nongenotyped and genotyped individual, and 
is the matrix of additive relationships among genotyped individuals. Now, the model for the genotypic values, when genotypes are coded as in Equation 2, becomes
where 
is that part of 
that cannot be imputed from knowledge of the breeding values of genotyped relatives. In practice, the true value of 
is not known, and data for its estimation may not be available. Rearranging these equations in terms of the uncentered 
rather than the centered matrix of genotype covariates 
results in
and substituting 
as previously defined, and letting 
results in
(3)which suggests that 
could be treated as an unknown constant and estimated as a fixed effect from the data (Fernando et al. 2014). The covariate vector for 
is denoted by 
for nongenotyped individuals, and by 
for genotyped individuals. So, Equation 3 becomes
which can be combined as
(4)where 
and the 0 matrix in the last term of Equation 4 is required because 
does not appear in the model for genotyped individuals.
When the vector 
represents the substitution effects of a large number of loci containing positive and negative effects, 
will tend to have a value close to zero. Accordingly, we have simulated some scenarios with positive 
so that the entire 
vector is positive to exacerbate the impact of 
Nevertheless, when marker rather than causal alleles are fitted in the model, the sign of the substitution effects depends on the phase relationship between marker and causal allele, which is equally likely to be positive or negative.
Even if 
in a population undergoing selection, it is expected that 
in nonfounders. Suppose 
is the mean of the observed genotype covariates in such a population undergoing selection, and these means are used to center the matrix 
of observed genotypes. Then, the model for the genotypic values can be written in terms of the matrix 
of centered covariates as
(5)and using 
for the covariate corresponding to 
Equation 5 can be written as
which can be combined as
(6)Note that the regression coefficients for 
in Equation 4 and 
in Equation 6, must be equal. This implies that 
Similarly, the intercepts of these two models must be equal too, and this implies 
Simulations
Phenotypic and genotypic data were simulated using XSim (Cheng et al. 2015), based on haplotypes from 10 genomic regions of 721 US Hereford beef cattle that were genotyped with the Illumina 770K BovineHD BeadChip, and reported in terms of the number of copies of the A allele at each locus. The selected regions came from choosing every 10th single nucleotide polymorphism (SNP), starting at SNP 5001, to get a low-density panel of 200 SNPs from each of chromosomes 1–10 (BTA1–BTA10), after eliminating SNPs with MAF < 0.01. These 2000 SNPs represent 10 0.1 M chromosomes. Average linkage disequilibrium (LD) between adjacent SNPs was 0.300. 50 SNPs from each chromosome were randomly chosen to represent QTL. The QTL effects were sampled from a normal distribution with mean 
= 0.2 and multiplied by the number of copies of the A allele to produce the true breeding value (TBV). The TBVs were added to a normally distributed residual term scaled by the sample variance of the TBV to simulate a trait with a heritability of 0.5. The first 20 SNPs from each of the 10 chromosomes were also used to simulate a smaller panel, with 5 QTL and 15 markers per chromosome. The average LD between adjacent SNPs was 0.289. TBV were simulated in the same way as for the 2000 SNP scenario, then scaled to simulate traits with heritabilities (
) of 0.1, 0.3, or 0.5. An additional scenario with 
= 0 was used to simulate TBV for a trait with heritability 0.5.
Half the observed diplotypes from US Hereford cattle were assigned to represent males and the remainder to represent females. Those 360 males and 361 females were sampled in pairs, with replacement, to produce 4000 male and 4000 female offspring representing generation G4. There were no mutations. Four more nonoverlapping generations of random mating were carried out, with one male and one female offspring per dam mated to randomly chosen sires to produce the founder population (G0).
The G1 generation was produced by mass phenotypic selection of the top 200 G0 males, and this was repeated for five more generations. Each female was randomly mated twice to selected males to produce one offspring of each sex each generation. Across nonoverlapping generations G0–G5, a total of 48,000 individuals with phenotypes, genotypes, and TBVs were simulated for each scenario.
The training data included phenotypes from all individuals in G0–G4 (n = 40,000), and genotypes from all 1000 sires and all 8000 G5 animals. Fixed loci, if any, were filtered from the panel before genomic prediction analyses. The genetic and residual variances used in genomic prediction were the sample variance of the TBV in G0 and the corresponding residual variance used to define the desired heritability in the founder population.
Models
Five statistical models were compared for differences in accuracy and bias of prediction. These included models with or without 
, and with or without centering of observed and imputed marker covariates, and a model that used pedigree relationships but not marker covariates.
Mixed linear model:
Accuracy of PBLUP was quantified using the correlation of TBV and estimated breeding values (EBV), where TBV was as simulated and EBV were obtained by fitting the mixed linear model (Henderson 1973, 1984):
where 
is a vector of phenotypic observations, 1 is a vector of 1s, 
is a general mean, 
is a vector of random additive genetic effects, 
is a vector of random residual effects, and 
is a known incidence matrix relating observations to 
In this model, 
so that 
Further, 
with 
being the numerator relationship matrix, and 
so that 
Single-step Bayesian regression model:
All genomic EBVs (GEBV) were obtained by single-step Bayesian regression (Fernando et al. 2014), using BayesC priors for marker effects with π = 0, which gives predictions that are identical to those from single-step GBLUP. The model was implemented in Julia (http://julialang.org) based on the SSBR package (http://QTL.rocks) to construct an MCMC chain of 50,000 samples. Individuals were separated into two groups designated with subscripts g or n according to whether or not simulated genotypes were assumed to be observed or missing. The single-step Bayesian regression model including a covariate J for 
(model J) was:
where 
and 
are vectors of phenotypes for nongenotyped and genotyped individuals, 1 is a vector of 1s, μ is a general mean, 
is the expected value of the linear component 
of the genotypic value if selection was absent, 
is a vector of random substitution effects of markers, 
is a vector of imputation residuals, 
and 
are incidence matrices relating the breeding values of nongenotyped and genotyped individuals to their phenotypes, 
which is defined for genotyped individuals, is a vector of −1s, 
which is defined for nongenotyped individuals, is a vector computed as 
is the matrix of imputed marker covariates, 
is the matrix of observed marker covariates, and 
and 
are vectors of random residual effects for nongenotyped and genotyped individuals. This model can be represented as
(7)where 
U = 
and 
is a vector of random residual effects.
This model was used for four variants of the single-step Bayesian regression analysis, depending on whether or not the covariate J corresponding to the mean 
was in the model, and whether or not the columns in the marker covariate matrix 
were centered using means of the observed and imputed genotype covariates. The analyses with J or without J are denoted as J or N when covariates were not centered, and as JC or C when the entire matrix of imputed and observed genotype covariates were centered, respectively (Table 1).
Accuracy of genomic prediction was quantified using the correlation of TBV and GEBV (
), where GEBVs were obtained from each of the four analyses described above. Bias of genomic prediction was quantified using the deviation from unity of the coefficient of regression of TBV on GEBV (
). In models JC and J, GEBVs are obtained using Equation 24 in (Fernando et al. 2014):
where 
is the GEBV, 
is the best linear unbiased estimate of the mean of breeding values, 
is the BLUP of the vector of random substitution effects of all markers, and 
is the BLUP of the imputation residual.
In model J, the matrix 
contains the uncentered number of copies of the A allele at each locus, and the uncentered version is used to impute 
In model JC, the entire matrix of imputed, 
and observed, 
genotype covariates is centered. In models C and N, the GEBVs are computed in a corresponding manner, except that the covariate J and its coefficient 
are not included in the model.
The four analyses J, N, JC, and C were all applied to three different genotype panels, comprising the causal QTL plus markers, just the causal QTL, or just the markers. All 12 combinations of four analyses and three genotype panels were applied to data simulated with 200 loci comprising 50 QTL whose effects were sampled from a normal distribution with 
to construct phenotypes with 
The four analyses JC, J, C, and N were repeated using only genotype panels comprising QTL plus markers for three other scenarios: 200 loci, 
200 loci, 
and 200 loci, 
Fernando et al. (2014) had observed that including 
in the model was necessary only when 
and the number of observations exceeds the number of markers. So, to examine the impact of 
the scenario with 200 loci and 
was repeated with 
and to examine the impact of the number of loci, the scenario with 
and 
was repeated with 2000 loci.
Every scenario was replicated 10 times with each replicate having been constructed starting from the sampling of G5 which represented simulated offspring from the haplotypes of real animals. Every phenotypic dataset was also fitted using the PBLUP model. All reported correlations and regression coefficients are the means of 10 replicates. These are presented along with the SEs of those means.
In single-step GBLUP (Legarra et al. 2009; Aguilar et al. 2010; Christensen and Lund 2010), the missing genotypes are not explicitly imputed, and only the observed genotype covariates are centered using their means. We argued earlier that when the number of loci is large, 
will be close to zero, especially when the panel consists of only marker loci. However, in a population undergoing selection, 
which is the mean breeding value of genotyped individuals, is not expected to be zero even when the panel only includes marker loci and 
So, in addition to the above, single-step Bayesian regression analyses, using the model given in Equation 7 with and without J (models JC* and C*), were applied to a marker panel with 200 loci, 
and 
when the matrix of observed genotype covariates were centered using their means as: 
where 
and 
is the number of individuals with genotypes. Recall that even when the means of the observed genotypes are used for centering, the model for the genotypic values can be written in terms of the matrix 
of uncentered covariates as shown by Equation 6. We will compare the estimates of 
in this model with those of 
from Equation 4 where the means of observed genotype covariates are not used for centering.
Data availability
The genotypes representing G0 from each replicate and scenario are available at: https://figshare.com/s/d7798b811a9a6a4172fc. These genotypes and the methodology described previously are sufficient to reproduce the simulations used in this study.
Results and Discussion
Effect of fitting a genotypic mean and centering of observed and imputed marker covariates
Accuracy:
The accuracies of genomic prediction as assessed by validation in G5 after training using G0–G4 for a trait with 50 QTL whose effects were sampled from a normal distribution with 
= 0.2 and 
= 0.5 are in Table 2. The accuracy of PBLUP in predicting breeding values for individuals without phenotypes was 41.5%, accounting for <20% of genetic variance. All analyses using genotypes resulted in more accurate predictions than using pedigree alone. Centering the entire matrix of imputed and observed genotype covariates had no effect on the accuracy of the genomic analyses regardless of the nature of the marker panel. This is because, as shown in the Appendix, this type of centering is equivalent to adding and subtracting 
from the model equation, and this has no effect on the mixed model solutions for 
) between TBV and (G)EBV for alternative analysesIn this study, selection resulted in the successive advances in mean TBV from G0 to G5 being 10.14, 10.63, 11.18, 11.68, 12.15, and 12.57. The mean genotypic value was not zero in G0 because QTL genotypes were not centered and the mean QTL effect was 
= 0.2. Since the QTL effects do not change with selection, the advance in TBV reflects changes in the frequencies of the favorable alleles of the 50 QTL. So centering using the allele frequency means of the genotyped sires in G1–G4 and all individuals in G5 does not closely approximate the centering that would have occurred if the allele frequency means had been obtained from the unselected population. By contrast, fitting 
in the model estimates the relevant mean from the data.
In panels that included causal variants (QTL), fitting 
in the model substantially improved the accuracy to being near perfect. This is not surprising, given that there were only 50 QTL, the heritability was 0.5, and there were 40,000 phenotyped ancestors, including 200 genotyped sires per generation in the training. However, in the panel that contained only markers with no causal variants, fitting 
in the model had little impact. This is because in the population that was simulated here the founders were not selected, and thus 
was close to zero for the panel with only markers. This would not be the case in a population where even the founders are descendants of selected parents. Thus, for traits that have been subject to selection, fitting 
in the model is expected to improve accuracy.
Using one replicate as an example, for the panel including both QTL and markers, the estimate of μ was ∼9.90 for the analyses using J and N. The estimate of μg was 4.98 for the analysis using J. For the genotyped individuals, the covariate values in 
are all −1, so 
+ 
is a vector of values equal to 9.90−4.98 = 4.92. For nongenotyped individuals, the covariate values in 
can vary widely, but many are close to −1 while others are close to 0. This means that 
will include values that range from 9.90 to 4.92, accounting for variation in accuracy of imputation. When 
is not included in the model, these effects are ignored, which can reduce the accuracy of predicting nongenotyped individuals. Failing to account for these effects will propagate errors in 
and 
the latter impacting the accuracy of predicting genotyped individuals. Collectively, these errors reduced accuracy from 97 to 96% for the panel including QTL and markers, and from 99 to 90% for the panel including only QTL. However, when the panel comprised only markers, the estimates 
will include both positive and negative values because the phase of markers and QTL are equally likely to take either sign, in which case 
will be close to zero as confirmed in the above mentioned replicate where the estimate was −1.18.
Here, results from the analysis with only QTL on the panel are used to show that 
is the mean of 
in the founder population and not the mean of the breeding value 
of the genotyped individuals. Recall that the QTL model was used with an intercept of 
to simulate the data. Thus, when only QTL are on the panel, the true value of β is zero. In analysis J because 
both 
and 
are estimates of 
and could be pooled, which for the replicate above would be 
In that replicate, the actual mean of 
in G0 was 10.16, which was estimated in the analysis to be 9.05. On the other hand, the mean of the breeding value 
in the 9000 genotyped individuals was 2.4, which is clearly not near the pooled estimate of 9.05 for 
These genotyped individuals included 1000 selected sires, of which 200 were genotyped in each generation from G0 to G4, and 8000 offspring from G5. The mean values of 
for the selected sires were 1.00, 1.59, 2.08, 2.50, and 2.88, respectively, for G0 through G4, and 2.45 for the offspring in G5. It is apparent that the 
parameter corresponding to the covariate 
is the mean of 
in the founder population and not the mean breeding value 
in the selected individuals. In analysis JC with the covariates centered, the intercept β is the value of 
when 
which is the case when 
The estimate 
was 16.11 in this analysis, but 
remained about the same value, namely 7.80. This shows that 
has the same interpretation whether the entire matrix of observed and imputed genotypes is centered or not. In neither case does it represent the mean breeding value of selected individuals.
Accuracy of PBLUP increased with heritability, as expected (Table 3). Further, genomic predictions using panels that include causal mutations were near perfect when 
was included in the model. These high accuracies are a reflection of these phenotypes being influenced by only 50 QTL and there being a large training dataset. Accuracy was reduced when 
was not fitted in the model. There was no advantage in terms of accuracy to centering the covariates, but MCMC mixing may have been improved, although this was not investigated.
) between TBV and (G)EBV for alternative analyses for different heritabilities
) of TBV on (G)EBVBias:
Table 4 shows the regression coefficients of TBV on (G)EBV for 
and 
= 0.2, the same scenarios represented in Table 2. The regression coefficients of TBV on GEBV for each scenario were close to 1 with very low SE, which indicates that the genomic predictions exhibited almost no bias. The differences in regression coefficients between analyses were very small, but the marker panel comprising only markers was biased upward, whereas the marker panels that included causal mutations were biased slightly downward.
Effect of mean QTL effect (
vs.)
We had hypothesized that the impact of omitting 
from the model would be greatest when 
departs significantly from 0, which is more likely to occur when 
departs from 0. For that reason, our base simulation used 
= 0.2. Results are shown in Table 5 for the panel including QTL and markers with 
= 0.5 for 
= 0.2 compared with 
= 0. These results confirmed that the benefit of fitting 
was greatest when 
= 0.2, but there was still an advantage to fitting 
when 
= 0. That advantage is likely to erode as the number of QTL increases. Changing the mean QTL effect had no impact on bias, except for a slight influence on PBLUP.
) for alternative QTL distributions and analysesEffect of more QTL and markers (200 SNP vs. 2000 SNP)
We had hypothesized that the improvement of accuracy from adding an extra covariate for 
would reduce as the number of QTL increases, because 
is likely to be closer to zero for a trait that is more polygenic. Table 6 shows that PBLUP was largely unaffected by changes to genetic architecture, but the accuracy of genomic prediction declined slightly as the number of QTL increased. This reflects the fact that the precision of estimating QTL effects is greater when the effects are large, and polygenic traits with more QTL must have on average smaller effects when compared at the same genetic variance. The benefit of fitting 
in the model was virtually eliminated when the number of substitution effects to estimate increased from 150 to 1500. In contrast to the results for accuracy, there was no impact of QTL number on bias. Centering had no impact on accuracy or bias.
) for different numbers of QTL and alternative analysesCentering using the entire matrix of genotype covariates or only the observed genotype covariates
Table 7 shows the accuracies and regression coefficients of TBV on (G)EBV for the genotype panel with 150 markers, 
, and 
= 0. The analyses were performed after centering: the entire matrix of imputed and observed genotype covariates (
); only observed genotype covariates (
), which is the type of centering done in single-step GBLUP; or not centering the covariates (
). The accuracy of the genomic analysis with covariates centered as 
but without J (model C*) was ∼17% lower than the other genomic analyses. However, when J was included in the model with covariates centered as 
the accuracy of prediction was markedly improved.
) when centering for all genotypes or observed genotypesAs explained previously, 
where 
is the expected value of the covariates in the founders, will tend to zero for the marker panel that does not include QTL, even with 
However, even if 
in a population undergoing selection when selected individuals are genotyped, 
where 
is the expected value of the observed genotype covariates. In this study, selection was used to increase the mean of the trait. Thus, 
is expected to be negative because most of the genotyped individuals were from G5, whereas 
is expected to be zero. The negative estimate of 
from 10 replicates of the JC* analysis, −2.69, confirms that 
On the other hand, the mean of 
from 10 replicates of the JC analysis was −0.75. This explains why fitting J in the model improved the accuracy of genomic prediction when covariates were centered as in single-step GBLUP.
Fernando et al. (2014) found that centering using 
improved the accuracy without J in the model when the population was not under selection and the genotyped individuals were unselected. In that study, mating was random with no selection, so the allele frequency means of the genotyped individuals were a reasonable approximation of the allele frequency means in the founder population. By contrast, our simulation here shows that centering using 
can reduce the accuracy when the population is under selection, unless J is fitted in the model.
In single-step GBLUP, the observed genotypes are commonly centered by subtracting their mean and used to construct a genomic relationship matrix, for example, by using the first method proposed by VanRaden (2008). Using that genomic relationship matrix in the single-step GBLUP formula in Aguilar et al. (2010) does not account for J. This was recognized by Vitezica et al. (2011), who proposed a modification for populations under selection that involved adding a constant to all elements of the genomic relationship matrix that they derived by equating the sum of the elements of the genomic relationship matrix to the sum of the elements of the numerator relationship matrix. In the appendix of that paper they showed that this modification is equivalent to fitting a covariate Q = −J and treating −
as a random effect. In addition to this modification, Christensen et al. (2012) proposed a multiplicative scaling to the genomic relationship matrix such that its diagonals have the same mean as the diagonals of the numerator relationship matrix. Vitezica et al. (2011) claimed that −
represents the mean breeding value of selected individuals, and we have confirmed here that this is true provided the observed genotype covariates are centered by their mean.
Most populations are under natural or artificial selection. In many cases, genotypes are only available on selected individuals. In single-step genomic analysis that combine genotyped and nongenotyped individuals in a joint analysis, the mean of observed genotypes are available for centering. If the observed genotypes include QTL, the accuracy of genomic prediction can be severely compromised, unless the J covariate is fitted in the model. If the observed genotypes are only markers, the accuracy of genomic prediction may not necessarily be improved by fitting J in the model, but it doesn’t do any harm. However, if centering is applied only to the observed genotypes, which is the type of centering used in single-step GBLUP, accuracy could be severely compromised, unless the J covariate is fitted in the model or an equivalent approach is adopted.
Acknowledgments
The authors are grateful to Bruce L. Golden and Hao Cheng for assistance in the implementation of SSBR models, to Jack C. M. Dekkers for his constructive comments on design of the simulation, and to the very helpful anonymous reviewers. R.L.F. acknowledges useful discussions with Andres Legarra and Zulma Vitezica. This work was supported by the US Department of Agriculture’s Agriculture and Food Research Initiative National Institute of Food and Agriculture competitive grant 2015-67015-22947.
Appendix
Here we show that inference about 
does not depend on how the genotypes are coded. The marker effects model can be described by the following general model:
(8)where y is a vector of observed phenotypes, 1 is a vector of 1s, μ is a general mean, M is a matrix of marker covariates, coded 0, 1, 2, which might represent the number of copies of the A allele, 
is a vector of random substitution effects of markers, and e is a vector of residuals. Henderson’s mixed model equations (MME) that correspond to Equation 8 are:
where 
is the best linear unbiased estimate of the mean, and 
is the best linear unbiased predictor of the vector of random substitution effects of all markers. Now we can eliminate 
from the equations for 
by subtracting from those equations the equation for 
premultiplied by 
Then, the MME are transformed to
and substituting 
, and 
and its transpose, the transformed MME become
(9)where 
is the row vector of column means of 
as in 
Now, consider the coding obtained by centering the marker genotypes as 
Then the model can be written as
(10)where 
The MME that correspond to Equation 10 are
but 
and, similarly, its transpose is 
Then the MME become
Expanding 
gives 
but because 
and 
as previously shown, the second term, 
which is equal to the last term in the expansion. Thus, the MME become
(11)The equations for 
in (9) and (11) are identical, and this proves that centering with 
doesn’t affect inference about 
Now suppose an arbitrary vector 
is used to transform the genotypes as 
Then the model becomes
Adding and subtracting 
the above equation can be written a:
which is identical to Equation 10, proving that inference about 
does not depend on how the genotypes are coded.
Footnotes
Communicating editor: D. J. de Koning
- Received December 14, 2016.
- Accepted June 9, 2017.
- Copyright © 2017 Hsu et al.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
