基本信息

原文：http://potatobreeding.cals.wisc.edu/wp-content/uploads/sites/161/2014/01/GS_tutorial.pdf

作者：Jeffrey Endelman

版本：4

更新：20130615

翻译：张敖

翻译更新：20221026 12:44:13

翻译内容

　　本文介绍如何使用第四版中新加入rrBLUP的特性(Endelman 2011)。

　　本包的基本核心仍然是mixed.solve，它求解了除残差之外的一个方差分量的混合模型。该函数估计两个方差分量，通过ML或REML模型估计，估计的执行使用Kang等(2008)描述的光谱分析算法。在这个过程中，很容易创建表型协方差矩阵的逆矩阵，逆矩阵随后用于固定效应和随机效应的BLUE和BLUP求解计算(Searle et al. 1992)。在(Endelman 2011)中，作者展示了mixed.solve如何用于全基因组预测，要么建模标记作为随机效应，要么家系随机效应。

　　第4版是围绕A.mat和kin.blup函数设计的。

A.mat

　　A.mat用标记估计真实的亲缘关系矩阵（A），对于没有缺失数据的高密度标记，A.mat的VanRaden(2008)建议的第一个公式：

　　这里，矩阵W中等位基因的表示数值通过群体均值进行中心化。Endelman and Jannink (2012)证明了该公式的平均对角元素为1+f，f是近交系数。

缺失标记数据

　　当基因型有缺失时，A.mat有两个补缺失的选项。一个是缺失值用标记的群体均值代替，这对SNP芯片类，只有少量缺失值是足够的。对GBS标记，缺失值的水平可能过高。这种情况下，A.mat可以使用基于多元正态分布的EM算法估计亲缘关系矩阵(Poland et al. 2012)。

　　为了证明EM算法，我下载了Poland et al. (2012) 的GBS数据https://www.crops.org/publications/tpg/supplements/5/tpg12-06-0006-dataset-s2.gz（链接已失效）。

　　下列代码会读取GBS数据，并转换为rrBLUP需要的{-1，0，1}格式。

GBS<-read.csv("tpg12-06-0006-dataset-s2",header=T,as.is=TRUE,row.names=1)
alleles <- setdiff(unique.x,union("H","N")){
    unique.x <- unique(x)
    y <- rep(0,length(x))
    y[which(x==alleles[1])] <- -1
    y[which(x==alleles[2])] <- 1
    y[which(x=="N")] <- NA
    return(y)
}
X <- apply(GBS[,-c(1:3)],1,parse.GBS)
dim(X) #lines by markers
frac.missing <- apply(X,2,function(z){length(which(is.na(z)))/length(z)})
length(which(frac.missing<0.5))
hist(frac.missing)

　　有1.6万的标记缺失率小于50%，足够估计出254个家系的关系矩阵。因为作者是在一个具有多个处理器的unix兼容的系统上运行这段代码的，所以我可以使用12个核来加速EM算法：

library(rrBLUP)
system.time(A1 <- A.mat(X,impute.method="EM",n.core=12,max.missing=0.5))

　　EM算法在其进行过程中显示收敛序列，它表示亲缘关系系数的根均方误差。默认停止标准为0.02，但是可以通过tol参数改变（输入?A.mat获得更多信息）。在本例中，当它达到0.0151时，计算停止，这仅仅需要1分钟的时间。如果只有一个核，可以简单的省略n.core选项，因为默认为1核。译者注：windows下指定核数量无效。

　　用均值进行补缺失，可以使用下列语法：

system.time(A2 <- A.mat(X, max.missing=0.5))

　　用均值补缺失肯定会更快，在很多情况下，在GEBV的预测精度上，均值的表现与EM算法和其他更先进的方法一样好。然而，与EM方法相比，均值的方法的育种值往往更有偏向性(Poland et al. 2012)。对两个矩阵的平均对角线元素比较表明，EM算法的结果更加接近给定1%杂合率的期望，表达式1+f≈2。提出

round(mean(diag(A2)),2)   # imputed with mean
round(mean(dia(A1)),2)   # imputed with EM

A矩阵的收缩估计

　　A.mat的另一个特性是收缩估计，它可以用于低密度的标记，例如来自384个SNP的芯片。当家系的数量与标记的数量相当或更多时，上面的方程可能是最优化的亲缘关系矩阵估计。Endelman and Jannink (2012)建议将估算值降低到(1+f)I，用收缩强度选择最小化均方误差。

　　为了说明这一点，我将使用BLR包中的一个数据集，该数据集由1279个DArT标记对599个小麦家系进行了基因分型。

library(BLR)
data(wheat)
M <- 2*X-1 #convert markers to {-1,1}
dim(M)   #=[1] 599 1279
A1 <- A.mat(M,shrink=TRUE)   #= [1] "Shrinkage intensity: 0.03"
A2 <- A.mat(M[,sample(1:1279,384)],shrink=TRUE)   #= "Shrinkage intensity: 0.1"
A3 <- A.mat(M[,sample(1:1279,192)],shrink=TRUE)   #= "Shrinkage intensity: 0.17"

　　如上例所示，收缩强度从0（无收缩）到1（完全收缩），标记密度减少，缩减强度增加。所有1279个标记，使用了非常小的缩减（3%），而384和192个标记的随机集合，缩减强度是10%和17%。

kin.blup

Endelman(2001)中，我介绍了一个mixed.solve的“包装器”，叫做kinship.BLUP，这是为家系的预测而设计的，使用加性遗传模型或高斯核。然而，“包装器”没有那么方便，所以我又设计了一个新的函数代替它，这个函数叫做kin.blup。该函数不需要用户创建设计矩阵，新函数会自动从数据框中完成创建部分。另一个不同是，用户用kin.blup传递的是亲缘关系矩阵，而不是标记，这样允许我们更好的计算A矩阵（参见A.mat）。

为了说明基础功能，我将继续以上述的小麦数据为例，其中，A1矩阵用所有的标记估计。这599个自交系已经在BLR包中分成了10个集合用于交叉验证。为了预测集合1的育种值，使用集合2至10的自交系表型，首先，表型和相应的基因型标识符必须组装成数据框。

test <- which(sets==1)
yNA <- Y[,1]   # grain yield in environment 1
yNA[test] <- NA   # mask yields for validation set
data1 <- data.frame(y=yNA, gid=1:599)

　　验证集合需要遮盖表型数据，需要将他们的值设成NA，如上面的例子。最小的数据集有两列，一列是表型值，一列是基因型的标签（注：基因型名称）。确保数据框中基因型标签（名称）与关系矩阵的行名对应，我们可以做全基因组预测：

rownames(A1) <- 1:599
ans1 <- kin.blup(data1, K=A1, geno="gid", pheno="y")
str(ans1)

　　正如上面的例子，亲缘关系矩阵传入参数“K”（给kinship）到kin.blup函数，以及数据框中的表型和基因型变量名。该函数返回方差组分（$Vg, $Ve）的REML估计，以及遗传值（$g）的BLUP，即本例中的育种值。同时，还将返回来自混合模型的残差。如上例所示，BLUP值在亲缘关系矩阵中返回599个条目（观测值），尽管这些系中10%的系没有表型（由于遮盖）。BLUP的顺序由K矩阵的顺序决定（array也如此命名）。

　　为了评估GEBV的预测精度，计算验证数据集的预测值与遮盖的表型值的相关性：

round(cor(ans1$g[test),Y[test,1],2)

　　对于具有高斯核的预测，而不是传递关系矩阵到K参数，使用欧式距离，并将GAUSS参数者只为TRUE。

D <- as.matrix(dist(M))   # Euclidean distance
system.time(ans2 <- kin.blup(data1, K=D, GAUSS=TRUE, geno="gid", pheno="y"))
system.time(ans2 <- kin.blup(data1, K=D, GAUSS=TRUE, geno="gid", pheno="y", n.core=10))
round(cor(ans$g[test],Y[test,1],2)

　　正如所见，多核可以提高GAUSS核的运算速度。高斯核的预测精度是0.63，高于使用亲缘关系矩阵。该结果与Endelman（2011）的结果不完全相同，用于确定最优尺度网络点（grid point）参数的并不相同；输入?kin.blup获得设置网络点的更多信息。

多环境试验

　　kin.blup函数能够处理不同环境中的重复测量。来自BLR包的小麦数据集平均在4个环境下测定599个自交系。环境2-4是相关的，本例中，将考虑为育种计划的一个目标环境。下面的代码将创建一个不平衡的数据集，使用跨环境的部分重复。

y <- c(Y[1:400,2],Y[101,400,3],Y[201:500,4])
env <- c(rep(2,400),rep(3,300),rep(4,300))
gid <- c(1:400,101:400,201:500)
data2 <- data.frame(y=y,env=env,gid=gid)
nrow(data2)

　　为了在预测模型中将环境效应作为固定效应，数据框的列被传递给函数：

system.time(ans <- kin.blup(data2,K=A1,geno="gid",pheno="y",fixed="env"))
round(cor(ans$g[501:599],rowMeans(Y[501:599,2:4])),3)

　　当有多个固定效应用于建模时，例如年份和地点，只需要传递一个列名数组，例如fixed=c(“year”,”location”)。

　　上面的例子是一步预测，这比首先计算线平均值的两步方法的计算要求更高。对于不平衡数据，kin.blup可以用两步法的速度做预测，同时保留关于不同重复水平的信息。这通过reduce=TRUE实现，转换混合模型的维度等于自交系的数量（参看使用手册获得更多细节）：

system.time(ans2<-kin.blup(data2,K=A1,geno="gid",pheno="y",fixed="env",reduce=TRUE))
round(cor(ans2$g[501:599],ans$g[501:599]),3)
round(cor(ans2$g[501:599],rowMeans(Y[501:599,2:4])),3)

　　在上面的例子中，采用约简方法的计算时间减少了近5倍。两种方法的预测非常相似（r = 0.96），但在本例中没有降低。

References

　　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

　　Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709–1723.

　　Pérez et al. 2010. Genomic-enabled prediction based on molecular markers and pedigree using the Bayesian Linear Regression package in R. Plant Genome 3:106–116.

　　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．

　　Searle et al. 1992. Variance Components. John Wiley & Sons, Hoboken.

　　VanRaden, P.M. 2008. Efficient methods to compute genomic predictions. J. Dairy Science
91:4414–4423.