Title: | Ridge Regression and Other Kernels for Genomic Selection |
---|---|
Description: | Software for genomic prediction with the RR-BLUP mixed model (Endelman 2011, <doi:10.3835/plantgenome2011.08.0024>). One application is to estimate marker effects by ridge regression; alternatively, BLUPs can be calculated based on an additive relationship matrix or a Gaussian kernel. |
Authors: | Jeffrey Endelman |
Maintainer: | Jeffrey Endelman <[email protected]> |
License: | GPL-3 |
Version: | 4.6.3 |
Built: | 2024-10-31 20:53:13 UTC |
Source: | CRAN |
This package has been developed primarily for genomic prediction with mixed models (but it can also do genome-wide association mapping with GWAS
). The heart of the package is the function mixed.solve
, which is a general-purpose solver for mixed models with a single variance component other than the error. Genomic predictions can be made by estimating marker effects (RR-BLUP) or by estimating line effects (G-BLUP). In Endelman (2011) I made the poor choice of using the letter G to denotype the genotype or marker data. To be consistent with Endelman (2011) I have retained this notation in kinship.BLUP
. However, that function has now been superseded by kin.blup
and A.mat
, the latter being a utility for estimating the additive relationship matrix (A) from markers. In these newer functions I adopt the usual convention that G is the genetic covariance (not the marker data), which is also consistent with the notation in Endelman and Jannink (2012).
Vignettes illustrating some of the features of this package can be found at https://potatobreeding.cals.wisc.edu/software/.
Endelman, J.B. 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4:250-255. <doi:10.3835/plantgenome2011.08.0024>
Endelman, J.B., and J.-L. Jannink. 2012. Shrinkage estimation of the realized relationship matrix. G3:Genes, Genomes, Genetics 2:1405-1413. <doi:10.1534/g3.112.004259>
Calculates the realized additive relationship matrix
A.mat( X, min.MAF = NULL, max.missing = NULL, impute.method = "mean", tol = 0.02, n.core = 1, shrink = FALSE, return.imputed = FALSE )
A.mat( X, min.MAF = NULL, max.missing = NULL, impute.method = "mean", tol = 0.02, n.core = 1, shrink = FALSE, return.imputed = FALSE )
X |
matrix ( |
min.MAF |
Minimum minor allele frequency. The A matrix is not sensitive to rare alleles, so by default only monomorphic markers are removed. |
max.missing |
Maximum proportion of missing data; default removes completely missing markers. |
impute.method |
There are two options. The default is "mean", which imputes with the mean for each marker. The "EM" option imputes with an EM algorithm (see details). |
tol |
Specifies the convergence criterion for the EM algorithm (see details). |
n.core |
Specifies the number of cores to use for parallel execution of the EM algorithm |
shrink |
set shrink=FALSE to disable shrinkage estimation. See Details for how to enable shrinkage estimation. |
return.imputed |
When TRUE, the imputed marker matrix is returned. |
At high marker density, the relationship matrix is estimated as , where
and
is the frequency of the 1 allele at marker k. By using a normalization constant of
, the mean of the diagonal elements is
(Endelman and Jannink 2012).
The EM imputation algorithm is based on the multivariate normal distribution and was designed for use with GBS (genotyping-by-sequencing) markers, which tend to be high density but with lots of missing data. Details are given in Poland et al. (2012). The EM algorithm stops at iteration
when the RMS error =
< tol.
Shrinkage estimation can improve the accuracy of genome-wide marker-assisted selection, particularly at low marker density (Endelman and Jannink 2012). The shrinkage intensity ranges from 0 (no shrinkage) to 1 (
). Two algorithms for estimating the shrinkage intensity are available. The first is the method described in Endelman and Jannink (2012) and is specified by
shrink=list(method="EJ")
. The second involves designating a random sample of the markers as simulated QTL and then regressing the A matrix based on the QTL against the A matrix based on the remaining markers (Yang et al. 2010; Mueller et al. 2015). The regression method is specified by shrink=list(method="REG",n.qtl=100,n.iter=5)
, where the parameters n.qtl
and n.iter
can be varied to adjust the number of simulated QTL and number of iterations, respectively.
The shrinkage and EM-imputation options are designed for opposite scenarios (low vs. high density) and cannot be used simultaneously.
When the EM algorithm is used, the imputed alleles can lie outside the interval [-1,1]. Polymorphic markers that do not meet the min.MAF and max.missing criteria are not imputed.
If return.imputed = FALSE, the additive relationship matrix is returned. If return.imputed = TRUE, the function returns a list containing
the A matrix
the imputed marker matrix
Endelman, J.B., and J.-L. Jannink. 2012. Shrinkage estimation of the realized relationship matrix. G3:Genes, Genomes, Genetics. 2:1405-1413. <doi:10.1534/g3.112.004259>
Mueller et al. 2015. Shrinkage estimation of the genomic relationship matrix can improve genomic estimated breeding values in the training set. Theor Appl Genet 128:693-703. <doi:10.1007/s00122-015-2464-6>
Poland, J., J. Endelman et al. 2012. Genomic selection in wheat breeding using genotyping-by-sequencing. Plant Genome 5:103-113. <doi:10.3835/plantgenome2012.06.0006>
Yang et al. 2010. Common SNPs explain a large proportion of the heritability for human height. Nat. Genetics 42:565-569. <doi:10.1038/ng.608>
Performs genome-wide association analysis based on the mixed model (Yu et al. 2006):
where is a vector of fixed effects that can model both environmental factors and population structure.
The variable
models the genetic background of each line as a random effect with
.
The variable
models the additive SNP effect as a fixed effect. The residual variance is
.
GWAS(pheno, geno, fixed=NULL, K=NULL, n.PC=0, min.MAF=0.05, n.core=1, P3D=TRUE, plot=TRUE)
GWAS(pheno, geno, fixed=NULL, K=NULL, n.PC=0, min.MAF=0.05, n.core=1, P3D=TRUE, plot=TRUE)
pheno |
Data frame where the first column is the line name (gid). The remaining columns can be either a phenotype or the levels of a fixed effect. Any column not designated as a fixed effect is assumed to be a phenotype. |
geno |
Data frame with the marker names in the first column. The second and third columns contain the chromosome and map position (either bp or cM), respectively, which are used only when plot=TRUE to make Manhattan plots. If the markers are unmapped, just use a placeholder for those two columns. Columns 4 and higher contain the marker scores for each line, coded as {-1,0,1} = {aa,Aa,AA}. Fractional (imputed) and missing (NA) values are allowed. The column names must match the line names in the "pheno" data frame. |
fixed |
An array of strings containing the names of the columns that should be included as (categorical) fixed effects in the mixed model. |
K |
Kinship matrix for the covariance between lines due to a polygenic effect. If not passed, it is calculated from the markers using |
n.PC |
Number of principal components to include as fixed effects. Default is 0 (equals K model). |
min.MAF |
Specifies the minimum minor allele frequency (MAF). If a marker has a MAF less than min.MAF, it is assigned a zero score. |
n.core |
Setting n.core > 1 will enable parallel execution on a machine with multiple cores (use only at UNIX command line). |
P3D |
When P3D=TRUE, variance components are estimated by REML only once, without any markers in the model. When P3D=FALSE, variance components are estimated by REML for each marker separately. |
plot |
When plot=TRUE, qq and Manhattan plots are generated. |
For unbalanced designs where phenotypes come from different environments, the environment mean can be modeled using the fixed option (e.g., fixed="env" if the column in the pheno data.frame is called "env"). When principal components are included (P+K model), the loadings are determined from an eigenvalue decomposition of the K matrix.
The terminology "P3D" (population parameters previously determined) was introduced by Zhang et al. (2010). When P3D=FALSE, this function is equivalent to EMMA with REML (Kang et al. 2008). When P3D=TRUE, it is equivalent to EMMAX (Kang et al. 2010). The P3D=TRUE option is faster but can underestimate significance compared to P3D=FALSE.
The dashed line in the Manhattan plots corresponds to an FDR rate of 0.05 and is calculated using the qvalue package (Storey and Tibshirani 2003). The p-value corresponding to a q-value of 0.05 is determined by interpolation. When there are no q-values less than 0.05, the dashed line is omitted.
Returns a data frame where the first three columns are the marker name, chromosome, and position, and subsequent columns are the marker scores for the traits.
Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.
Kang et al. 2010. Variance component model to account for sample structure in genome-wide association studies. Nat. Genet. 42:348-354.
Storey and Tibshirani. 2003. Statistical significance for genome-wide studies. PNAS 100:9440-9445.
Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208.
Zhang et al. 2010. Mixed linear model approach adapted for genome-wide association studies. Nat. Genet. 42:355-360.
#random population of 200 lines with 1000 markers M <- matrix(rep(0,200*1000),1000,200) for (i in 1:200) { M[,i] <- ifelse(runif(1000)<0.5,-1,1) } colnames(M) <- 1:200 geno <- data.frame(marker=1:1000,chrom=rep(1,1000),pos=1:1000,M,check.names=FALSE) QTL <- 100*(1:5) #pick 5 QTL u <- rep(0,1000) #marker effects u[QTL] <- 1 g <- as.vector(crossprod(M,u)) h2 <- 0.5 y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) pheno <- data.frame(line=1:200,y=y) scores <- GWAS(pheno,geno,plot=FALSE)
#random population of 200 lines with 1000 markers M <- matrix(rep(0,200*1000),1000,200) for (i in 1:200) { M[,i] <- ifelse(runif(1000)<0.5,-1,1) } colnames(M) <- 1:200 geno <- data.frame(marker=1:1000,chrom=rep(1,1000),pos=1:1000,M,check.names=FALSE) QTL <- 100*(1:5) #pick 5 QTL u <- rep(0,1000) #marker effects u[QTL] <- 1 g <- as.vector(crossprod(M,u)) h2 <- 0.5 y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) pheno <- data.frame(line=1:200,y=y) scores <- GWAS(pheno,geno,plot=FALSE)
Genotypic value prediction by G-BLUP, where the genotypic covariance G can be additive or based on a Gaussian kernel.
kin.blup(data,geno,pheno,GAUSS=FALSE,K=NULL,fixed=NULL,covariate=NULL, PEV=FALSE,n.core=1,theta.seq=NULL)
kin.blup(data,geno,pheno,GAUSS=FALSE,K=NULL,fixed=NULL,covariate=NULL, PEV=FALSE,n.core=1,theta.seq=NULL)
data |
Data frame with columns for the phenotype, the genotype identifier, and any environmental variables. |
geno |
Character string for the name of the column in the data frame that contains the genotype identifier. |
pheno |
Character string for the name of the column in the data frame that contains the phenotype. |
GAUSS |
To model genetic covariance with a Gaussian kernel, set GAUSS=TRUE and pass the Euclidean distance for K (see below). |
K |
There are three options for specifying kinship:
(1) If K=NULL, genotypes are assumed to be independent |
fixed |
An array of strings containing the names of columns that should be included as (categorical) fixed effects in the mixed model. |
covariate |
An array of strings containing the names of columns that should be included as covariates in the mixed model. |
PEV |
When PEV=TRUE, the function returns the prediction error variance for the genotypic values ( |
n.core |
Specifies the number of cores to use for parallel execution of the Gaussian kernel method (use only at UNIX command line). |
theta.seq |
The scale parameter for the Gaussian kernel is set by maximizing the restricted log-likelihood over a grid of values. By default, the grid is constructed by dividing the interval (0,max(K)] into 10 points. Passing a numeric array to this variable (theta.seq = "theta sequence") will specify a different set of grid points (e.g., for large problems you might want fewer than 10). |
This function is a wrapper for mixed.solve
and thus solves mixed models of the form:
where is a vector of fixed effects,
is a vector of random genotypic values with covariance
, and the residuals follow
, with
by default. The design matrix for the genetic values has been partitioned to illustrate that not all lines need phenotypes (i.e., for genomic selection). Unlike
mixed.solve
, this function does not return estimates of the fixed effects, only the BLUP solution for the genotypic values. It was designed to replace kinship.BLUP
and to relieve the user of having to explicitly construct design matrices. Variance components are estimated by REML and BLUP values are returned for every entry in K, regardless of whether it has been phenotyped. The rownames of K must match the genotype labels in the data frame for phenotyped lines; missing phenotypes (NA) are simply omitted.
Unlike its predecessor, this function does not handle marker data directly. For breeding value prediction, the user must supply a relationship matrix, which can be calculated from markers with A.mat
. For Gaussian kernel predictions, pass the Euclidean distance matrix for K, which can be calculated with dist
.
In the terminology of mixed models, both the "fixed" and "covariate" variables are fixed effects ( in the above equation): the former are treated as factors with distinct levels while the latter are continuous with one coefficient per variable. The population mean is automatically included as a fixed effect.
The prediction error variance (PEV) is the square of the SE of the BLUPs (see mixed.solve
) and can be used to estimate the expected accuracy of BLUP predictions according to .
The function always returns
REML estimate of the genetic variance
REML estimate of the error variance
BLUP solution for the genetic values
residuals
predicted genetic values, averaged over the fixed effects
If PEV = TRUE, the list also includes
Prediction error variance for the genetic values
If GAUSS = TRUE, the list also includes
the log-likelihood profile for the scale parameter in the Gaussian kernel
Endelman, J.B. 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4:250-255. <doi:10.3835/plantgenome2011.08.0024>
#random population of 200 lines with 1000 markers M <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { M[i,] <- ifelse(runif(1000)<0.5,-1,1) } rownames(M) <- 1:200 A <- A.mat(M) #random phenotypes u <- rnorm(1000) g <- as.vector(crossprod(t(M),u)) h2 <- 0.5 #heritability y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) data <- data.frame(y=y,gid=1:200) #predict breeding values ans <- kin.blup(data=data,geno="gid",pheno="y",K=A) accuracy <- cor(g,ans$g)
#random population of 200 lines with 1000 markers M <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { M[i,] <- ifelse(runif(1000)<0.5,-1,1) } rownames(M) <- 1:200 A <- A.mat(M) #random phenotypes u <- rnorm(1000) g <- as.vector(crossprod(t(M),u)) h2 <- 0.5 #heritability y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) data <- data.frame(y=y,gid=1:200) #predict breeding values ans <- kin.blup(data=data,geno="gid",pheno="y",K=A) accuracy <- cor(g,ans$g)
***This function has been superseded by kin.blup
; please refer to its help page.
kinship.BLUP(y, G.train, G.pred=NULL, X=NULL, Z.train=NULL, K.method="RR", n.profile=10, mixed.method="REML", n.core=1)
kinship.BLUP(y, G.train, G.pred=NULL, X=NULL, Z.train=NULL, K.method="RR", n.profile=10, mixed.method="REML", n.core=1)
y |
Vector ( |
G.train |
Matrix ( |
G.pred |
Matrix ( |
X |
Design matrix ( |
Z.train |
0-1 matrix ( |
K.method |
"RR" (default) is ridge regression, for which K is the realized additive relationship matrix computed with |
n.profile |
For K.method = "GAUSS" or "EXP", the number of points to use in the log-likelihood profile for the scale parameter |
mixed.method |
Either "REML" (default) or "ML". |
n.core |
Setting n.core > 1 will enable parallel execution of the Gaussian kernel computation (use only at UNIX command line). |
BLUP solution for the training set
BLUP solution for the prediction set (when G.pred != NULL)
ML estimate of fixed effects
For GAUSS or EXP, function also returns
log-likelihood profile for the scale parameter
Endelman, J.B. 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4:250-255.
#random population of 200 lines with 1000 markers G <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { G[i,] <- ifelse(runif(1000)<0.5,-1,1) } #random phenotypes g <- as.vector(crossprod(t(G),rnorm(1000))) h2 <- 0.5 y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) #split in half for training and prediction train <- 1:100 pred <- 101:200 ans <- kinship.BLUP(y=y[train],G.train=G[train,],G.pred=G[pred,],K.method="GAUSS") #correlation accuracy r.gy <- cor(ans$g.pred,y[pred])
#random population of 200 lines with 1000 markers G <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { G[i,] <- ifelse(runif(1000)<0.5,-1,1) } #random phenotypes g <- as.vector(crossprod(t(G),rnorm(1000))) h2 <- 0.5 y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) #split in half for training and prediction train <- 1:100 pred <- 101:200 ans <- kinship.BLUP(y=y[train],G.train=G[train,],G.pred=G[pred,],K.method="GAUSS") #correlation accuracy r.gy <- cor(ans$g.pred,y[pred])
Calculates maximum-likelihood (ML/REML) solutions for mixed models of the form
where is a vector of fixed effects and
is a vector of random effects with
. The residual variance is
. This class
of mixed models, in which there is a single variance component other than the residual error,
has a close relationship with ridge regression (ridge parameter
).
mixed.solve(y, Z=NULL, K=NULL, X=NULL, method="REML", bounds=c(1e-09, 1e+09), SE=FALSE, return.Hinv=FALSE)
mixed.solve(y, Z=NULL, K=NULL, X=NULL, method="REML", bounds=c(1e-09, 1e+09), SE=FALSE, return.Hinv=FALSE)
y |
Vector ( |
Z |
Design matrix ( |
K |
Covariance matrix ( |
X |
Design matrix ( |
method |
Specifies whether the full ("ML") or restricted ("REML") maximum-likelihood method is used. |
bounds |
Array with two elements specifying the lower and upper bound for the ridge parameter. |
SE |
If TRUE, standard errors are calculated. |
return.Hinv |
If TRUE, the function returns the inverse of |
This function can be used to predict marker effects or breeding values (see examples). The numerical method
is based on the spectral decomposition of and
, where
is
the projection operator for the nullspace of
(Kang et al., 2008). This algorithm generates the inverse phenotypic covariance matrix
, which can then be used to calculate the BLUE and BLUP solutions for the fixed and random effects, respectively, using standard formulas (Searle et al. 1992):
The standard errors are calculated as the square root of the diagonal elements of the following matrices (Searle et al. 1992):
For marker effects where K = I, the function will run faster if K is not passed than if the user passes the identity matrix.
If SE=FALSE, the function returns a list containing
estimator for
estimator for
BLUE()
BLUP()
maximized log-likelihood (full or restricted, depending on method)
If SE=TRUE, the list also contains
standard error for
standard error for
If return.Hinv=TRUE, the list also contains
the inverse of
Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.
Endelman, J.B. 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4:250-255.
Searle, S.R., G. Casella and C.E. McCulloch. 1992. Variance Components. John Wiley, Hoboken.
#random population of 200 lines with 1000 markers M <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { M[i,] <- ifelse(runif(1000)<0.5,-1,1) } #random phenotypes u <- rnorm(1000) g <- as.vector(crossprod(t(M),u)) h2 <- 0.5 #heritability y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) #predict marker effects ans <- mixed.solve(y,Z=M) #By default K = I accuracy <- cor(u,ans$u) #predict breeding values ans <- mixed.solve(y,K=A.mat(M)) accuracy <- cor(g,ans$u)
#random population of 200 lines with 1000 markers M <- matrix(rep(0,200*1000),200,1000) for (i in 1:200) { M[i,] <- ifelse(runif(1000)<0.5,-1,1) } #random phenotypes u <- rnorm(1000) g <- as.vector(crossprod(t(M),u)) h2 <- 0.5 #heritability y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g))) #predict marker effects ans <- mixed.solve(y,Z=M) #By default K = I accuracy <- cor(u,ans$u) #predict breeding values ans <- mixed.solve(y,K=A.mat(M)) accuracy <- cor(g,ans$u)