rrBLUP

可用于genomic selection基因组选择和association mapping关联映射

rrBLUP包中的核心函数：

1.mixed.solve

把标记效应建模为随机效应

或

把行数据的基因型（genotypic）值用于A.mat函数（计算加性关系矩阵，预测育种值）

2.kinship.BLUP：包含上位效应在基因型值预测中

3.GWA：关联映射

 

处理缺失的数据：

函数mixed.solve，kinship.BLUP和GWA都允许缺失表型数据，那些observations只是简单地不被使用就行。

当表型数据缺失时,kinship.BLUP（选项‘RR’）和GWA都依赖于A.mat中的EM算法，也可以与mixed.solve联合使用。

kinship.BLUP中的非加性核基于dist，它使用成对的完整observations

 

并行计算：

对MAC，Linux和UNIX用户，R包multicore可与rrBLUP联合使用，以便使用单台机器多处理器的优势。对大的数据集，尤其是具有缺失数据时，推荐试用这一特性。对kinship.BLUP，A.mat和GWA都有效。

需要R>=24.1才能正常工作，而且必须在命令行使用R而不是GUI中。

 

函数使用：

1.A.mat：计算加性关系矩阵

A.mat(G,min.MAF=NULL,max.missing=NULL,impute=TRUE,tol=0.02,n.core=1,return.G=FALSE)

其中：

G：n*m unphased基因型矩阵，n个observations，每个observation有m个biallelic标记，编码为{-1,0,1}={aa,Aa,AA}，允许有部分的（impute）和缺失的（NA）值。

min.MAF：最小等位基因频率下限。默认时移除monomorphic markers。

max.missing：缺失值比例的上限，默认时移除完全缺失的标记。

impute：若为真，则缺失的基因型数据会被估算；否则，从完整的observations对中计算A，不保证半正定性。

tol：指定用EM算法估计缺失值时的收敛标准。若tol<0，用每个标记的种群均值估算缺失值。

n.core：n.core>1允许在多核机器上并行执行。（需要安装R包multicore，只能用命令行执行，不能在GUI中运行）

return.G：若为真，返回估算的标记矩阵。当使用了EM算法时，估算的等位基因可能位于区间[-1,1]之外。

多态基因不满足min.MAF，也不估算max.missing标准。

 

A.mat函数的值：

return.G为假时，返回n*n的加性关系矩阵

return.G为真时，返回包含以下信息的list：

（1）A：加性关系矩阵A

（2）G.imp：估算的G矩阵

 

2.GWA关联映射

GWA(y,G,Z=NULL,X=NULL,min.MAF=0.05,n.core=1,check.rank=FALSE)

其中：

y：n*1向量，表示observations，不包含NA缺失值。

G：t*m的基因型矩阵，t行observations（？），每行包含m个biallelic标记。

Z：n*t的0-1矩阵，n行observations，如不给出，则使用单位矩阵。

X：n*p的固定效应设计矩阵。若不给出，由1组成的向量对截距建模。

check.rank：若为真，将检查每一个标记的增广设计矩阵的rank。设计矩阵奇异的标记，将赋予零值。

 

GWA函数的值：m*1向量，表示标记的分数，值等于-log10(p-value)

 

该函数实现了2010年Kang等人的迭代广义最小二乘法iterative, generalized least-squares method，使用mixed.solve估计方差组分。

使用min.MAF足够保证问题是well-posed。然而，如果一个错误消息显示问题是奇异的，则将check.rank设为TRUE，这将使算法变慢，但能fix the error。

 

3.kinship.BLUP：使用mixed.solve，基于行之间的kinship进行基因组预测。

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)

其中：

G.train：n.train*m个训练群体的unphased基因型值，n.train行*m个biallelic标记。

Z.train：n.obs*n.train的0-1矩阵，在训练集中将observations与行关联。

K.method：默认为“RR”，岭回归。其中的K是用A.mat计算出的加性关系矩阵。其它选项：

（1）GAUSS：高斯核

（2）EXP：指数核

其中的D是用dist计算出的欧几里得距离。

n.profile：用于k.method=“GAUSS”或“EXP”，对尺度参数，在对数似然profile中使用的点的数量。

mixed.method：为“REML”（默认）或“ML”

 

kinship.BLUP函数的值：

（1）$g.train：训练集的BLUP解

（2）$g.pred：预测集的BLUP解（当G.pred不为空时）

（3）$beta：固定效应的ML估计

对于GAUSS或EXP核，还返回$profile：尺度参数的对数似然profile

 

4.mixed.solve：计算混合模型的最大似然ML或REML

mixed.solve(y,Z=NULL,K=NULL,X=NULL,method="REML",bounds=c(1e-9,1e+9),SE=FALSE,return.Hinv=FALSE)

其中：

y：n*1的observations向量，忽略缺失的值，和X与Z中相应的列对应。

Z：n*m的随机效应设计矩阵。若不传递，默认为单位矩阵。

K：m*m的随机效应协方差矩阵，必须为半正定。若不传递，则设为单位矩阵。

X：固定效应的设计矩阵。若不传递，设单位向量作为截距建模。

method：为完全（ML）或受限（REML）最大似然方法。

bounds：有2个元素的数组，指定为岭参数的下限和上限。

SE：若为真，则计算标准差。

return.Hinv：若为真，则函数返回H的inverse（对GWA有用）。

 

mixed.solve函数的值：

（1）$Vu

（2）$Ve

（3）$beta

（4）u：BLUP(u)

（5）$LL：最大对数似然

若SE=TRUE，返回的结果list还包括：

（1）$beta.SE

（2）$u.SE

若return.Hinv为真，返回的结果list还包括：$Hinv（H的逆矩阵）