Don't BLUP Twice
Year: 2024
Type: article
Abstract: Large and complex data sets can be difficult to model in a single comprehensive analysis. In many crop experiments, phenotypic data are collected on lines replicated within and across environments. Optimal analyses of these experiments account for replication and experimental design factors along with genomic relationships among lines. In the context of genome-wide association study (GWAS) scans, fitting complex models accounting for multiple levels of random variation (both genetic and nongenetic) in the data can result in prohibitively large computation time requirements. Ignoring these complexities, however, results in increased residual variance, reduced statistical power, and potential confounding of marker–trait associations with population structure. A useful analytical approach in this situation is a 2-stage analysis, in which a first-stage analysis is conducted to account for the complexities of experimental design and generate summary values for each line. For data from a single trial, this will be a summary value across replicates for each line. For multienvironment data, the summary will be across environments. Importantly, a first-stage model can include many covariance parameters and incur a relatively high computational burden and time cost because it needs to be fit to the data only once per trait, no matter how many markers have been scored. The summary values of each line obtained from the first stage can then be used as phenotypic inputs to the second-stage analysis to test each marker effect. Because the replicated line phenotypes have been reduced to a single summary value already adjusted for any covariates or random nongenetic factors, a computationally efficient GWAS model can be fit to test millions of markers in reasonable time. The 2-stage analysis, therefore, is a good compromise that accounts for complex experimental designs and nongenetic factors while keeping computational time constrained. However, confusion remains about how to model line effects in 2-stage analyses for GWAS, quantitative trait locus (QTL) mapping, and genomic selection (GS) analyses. Line effects can be modeled as fixed effects (with no distributional assumption), as independent random effects (assumed to be random uncorrelated samples from a normal distribution whose variance can be estimated in the analysis), or as random effects with covariances proportional to genomic relationships. Typical GWAS models fit line effects as random with genomic relationships to account for population structure and pedigree relationships when testing individual marker–trait associations (Yu et al. 2006). This is appropriate for a single-stage analysis or in the second-stage GWAS when using a 2-stage approach. But lines should not be treated as random effects in the first-stage analysis! Fitting lines as random effects in the first-stage analysis when data are unbalanced results in best linear unbiased predictors (BLUPs) of lines that are affected by the data structure in addition to genetic effects. Treating lines as random with genomic relationships in the first stage has the additional problem of reducing or eliminating the effects of individual genetic variants that the GWAS is searching for. The problems with fitting line effects as random in the first of a 2-stage analysis arise due to the properties of BLUPs of random effects in linear mixed models. Because of their distributional assumptions, line BLUPs are shrunk back toward the population mean compared to analogous estimators of line estimated marginal means, also called best linear unbiased estimators (BLUEs; Robinson 1991). The shrinkage effect on BLUPs depends on their estimated reliability, which in turn depends on the number of observations of each line, the trait heritability, and the amount of information available on relatives weighted by the genetic similarity (if lines are modeled with genetic relationships). If line effects are not modeled with kinship relationships, such that the genetic values are considered to be identically and independently distributed, and if data are balanced (no missing data, complete blocks, all lines measured in every replication and every environment), BLUEs and BLUPs will have identical rankings. In this scenario, the deviation between a line's BLUP value and the overall population mean will be equal to the corresponding deviation of its BLUE and the overall mean multiplied by the broad-sense heritability of the line means. Thus, the only difference is the scaling of the line value deviations from the overall mean due to the shrinkage effect. In real-world experiments, however, data are very often unbalanced due to missing data or incomplete block experimental designs. Consider an unbalanced data set, wherein some lines are missing in some replication blocks; in this case, the lines vary in terms of how much information on their effects exists in the data set (Isik et al. 2017). If lines are considered independent random effects, their BLUP values are shrunk back toward the population mean by different amounts corresponding to the number of replications and environments in which they are measured. Lines measured in more environments will have less shrinkage in their genetic value predictions, reflecting their greater precision due to more information. At the extreme, we can predict the BLUP of a line with no information: it will simply be the population mean. This is due to the assumption of the distribution of random effects: they are drawn from some larger reference population, and a priori, we expect their mean value to be 0. In contrast, if lines are considered fixed effects, a line with no data will have a nonestimable BLUE value: we have no distributional assumption about those effects, they do not have to be normally distributed, and there is no larger reference population whose central value is meaningful. For the purposes of selection, this differential shrinkage effect that reflects the precision of line BLUPs is a good thing (Robinson 1991): it downweights unreliable predictions; but it is a bad thing for any second-stage analysis because the line BLUPs reflect not only genetics but also the data structure and possibly the experimental design! Further, consider the situation of using BLUPs predicted from a first-stage analysis that incorporates the genomic relationships among lines. The resulting genomic BLUPs (GBLUPs) from this first-stage analysis will have differential shrinkage depending on the amount of direct measurements and also on the frequency of closely related lines in the data set. Further, the GBLUP of a particular line will reflect not only the actual observations of that line but also the observations on other lines, weighted by their genomic relationships. If the GBLUP is doing something useful, we expect GBLUPs to have different rankings than the BLUEs. For selection, this is desirable: GBLUPs provide the best predictions of the polygenic effects, which typically correspond to breeding values. But the purpose of GWAS is to identify individual genetic variants that have measurable effects after accounting for polygenic background effects and population structure. So, if the input values we feed into the GWAS analysis are predictions of the polygenic background effects, then we have already taken out most of the signal that the GWAS is searching for! Additional problems with the use of BLUPs in the context of evolutionary biology studies were identified previously by Hadfield et al. (2010). It may seem contradictory to treat line effects as statistically fixed in stage 1 followed by treating them as random in stage 2, but in fact, this is necessary to recover results close to what would be obtained from a single-stage analysis (Smith et al. 2001; Piepho et al. 2012; Damesa et al. 2017; Verbyla 2023). This apparent contradiction probably has contributed to lingering confusion over how to properly handle genetic effects in a 2-stage GWAS. Furthermore, crop geneticists may hesitate to treat lines as fixed effects in the first stage of a 2-stage GWAS or GS analysis because they are aware that for the purposes of selection and breeding, BLUPs are better than BLUEs (Piepho et al. 2008), and GBLUPs are even better (Isik et al. 2017). Therefore, it is not surprising that researchers may want to obtain BLUPs or GBLUPs in the first stage of the analysis to get the “best” prediction of the line's value to use as input for second-stage GWAS. But this is not correct if the goal of the analysis is to obtain summary values to use as input “phenotypes” to a subsequent analysis; instead, lines should be fit as fixed rather than random effects in the first-stage analysis. Then, in the second-stage analysis, lines should be fit as random with their polygenic effects modeled with covariances proportional to the realized kinship coefficients. This recommendation holds for 2-stage GWAS, GS, or QTL analyses (Damesa et al. 2019; Feldmann et al. 2023; George and Cavanagh 2015). In each case, the first stage models function to boil down complex data sets to generate a single summary phenotype value for each line; the differences arise in the second stage of these analyses. It is inappropriate to treat genetic effects as random in both stages, as this results in a BLUP analysis of BLUPs, hence the recommendation “don't BLUP twice.” Animal breeders have used deregression (“unshrinking”) of first-stage BLUPs for 2-stage GS analyses (Garrick et al. 2009), and Verbyla (2023) demonstrated that estimating BLUEs in first-stage analysis as recommended here is closely related to deregression (Garrick et al. 2009). In summary, a best practice for 2-stage GWAS is to first conduct an analysis of phenotype data alone without using genotypic information. One should fit a model with lines treated as fixed effects to obtain adjusted line means (BLUEs), considering other environmental effects as random or fixed factors or covariates as appropriate. The adjusted mean values from that mixed model are BLUEs, and these BLUEs can be used in the second-stage analysis. In the second stage, the line effects can be considered random effects with covariances proportional to their genomic relationships to account for polygenic background effects in GWAS (Fig. 1; Supplementary Files 1–3). R code and data to demonstrate two-stage GWAS using data from Choquette et al. (2023) is provided in Supplementary Files 1–3. This general procedure holds for both single-environment (Fig. 1a) and multienvironment trials (Fig. 1b) if the goal is to estimate marker effects averaged across environments. If data are from multienvironment trials and the goal of GWAS is to estimate marker-by-environment interactions, however, the second-stage model will necessarily be more complex, and those details will not be pursued here. Two-stage analysis workflow for (A) single-environment trials and (B) multienvironment trials. In the first stage of a single-environment trial A), a model Yjk = µ + Gj + Bk + εjk is fit to the data with Gj line effects treated as fixed and Bk block effects treated as random. BLUEs of line values are estimated and used as dependent variables in a second-stage analysis testing each marker's effect: Y¯j = µ + xjpβp + Aj + εj*, where xjp is a numeric indicator variable for allele dosage at marker p in line j, βp is the fixed effect of marker p, Aj is the polygenic background effect line j with covariances equal to the genomic relationship matrix G times the additive genetic variance (σ^A2), and εj* is the residual effect, distinguished from plot-level experimental error effects (εjk in stage 1), as it includes both estimation error of the line phenotypic mean value and genetic residual reflecting genotypic values not capture by the genomic relationship matrix. For multienvironment trials B), random terms for environment main effects (Ei) and genotype-by-environment interactions (GEij) are added to the first-stage model, but the second-stage model for marker main effects remains unchanged. A refinement of this general approach may be beneficial when the accuracy of variance parameter estimates for nongenetic effects is critical because of sparseness of the data and can be improved by modeling genotypes as random. The suggested refined approach in the first stage involves 2 steps (George and Cavanagh 2015; Smith et al. 2001; Verbyla 2023): (i) fit the mixed model taking genotypes as random and store all nongenetic variance component estimates and (ii) fit the same model taking genotypes as fixed, plugging in the nongenetic variance component estimates from (i) and compute the line BLUEs (Fig. 2; Supplementary Files 1–3). The second-stage GWAS would proceed as outlined in the other cases (Fig. 2). A refined approach for optimizing BLUEs in first stage of a 2-stage analysis. The first stage is split into 2 iterations. First, lines are treated as random effects with genetic covariances to obtain optimal estimates of random covariance components (random block effects are illustrated here but often more complex models involving random spatial effects are considered in the demonstration R code provided in Supplementary File 1). Next, the first-stage model is refit treating lines as fixed effects but forcing the random nongenetic covariance components to equal their estimated values from the initial model and line BLUEs are estimated. The second-stage model proceeds as before. Going forward, some details remain to be considered: if the data are unbalanced, the 2-stage analysis loses some information if it treats line BLUEs as independent “observations” having equal precision. For genomic prediction, differences in precision and covariances among BLUEs introduced by data unbalance can be accounted for using various weighted analysis schemes (Piepho et al. 2012; Damesa et al. 2017; Endelman 2023). The most accurate method carries the full variance–covariance matrix of the line BLUEs forward and integrates it as the residual variance–covariance structure for the line BLUEs as the response in the second-stage analysis. This method, which is also standard procedure in network meta-analysis (Schwarzer et al. 2015; Madden et al. 2016), may also be computationally demanding, especially with a large number of lines. Computational savings are afforded by diagonal approximations, but these come at the expense of some loss of efficiency (Damesa et al. 2017; Verbyla 2023). If data are balanced with respect to complete blocks, fitting complete blocks as fixed in stage 1 will not impact the line BLUEs (Piepho et al. 2003) and may reduce computational burden of weighted stage 2 analysis by removing the contribution of complete block variation to covariances among line BLUEs (Möhring and Piepho 2009). This trick is not always useful, however, because incomplete block designs and spatial analyses are favored for field trials with many entries, and properly modeling random incomplete blocks and spatial trends as random effects unavoidably results in covariances among entry BLUEs (Möhring and Piepho 2009). There does not seem to be a method for data weighting to account for variable precision (and possibly covariances) among input values for popular GWAS software, however. Developing such methods will require adding an extra random term to the stage 2 model to capture genotypic residuals (genetic effects not captured by the genomic relationship matrix but not influenced by residual weights/covariances; Feldmann et al. 2023). Perhaps, this genetic residual variance component could be estimated once and then fixed for subsequent marker tests to speed up the GWAS scan. Further software development could help improve GWAS to account for deviations from identical and independent residual errors that can arise in 2-stage analyses. The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, and supplementary files. Supplemental material available at G3 online. This research was supported by the United States Department of Agriculture – Agriculture Research Service and German Research Foundation (DFG), project PI 377/24-1.
Source: G3 Genes Genomes Genetics
Authors: James B Holland, Hans-Peter Piepho
Language: English
Subfield: Genetics
Domain: Life Sciences
Open Access status: gold
Estimated APC paid: $2,790
Award: G139225745
OpenAlex ID: W4404506419