rrBLUP使用手册

1.rrBLUP包简介

R 包rrBLUP能方便快速地实现RRBLUP(Ridge Regression Best Linear Unbiased Prediction)模型，该包主要有A.mat和mixed.solve两个核心函数。

1.1 A.mat函数

该函数主要用来过滤和填充基因型，返回加性效应关系矩阵（即Kinship）。

用法

A.mat(X,min.MAF=NULL,max.missing=NULL,impute.method="mean",tol=0.02,
n.core=1,shrink=FALSE,return.imputed=FALSE)

---

参数

X	n个品系与m个双等位基因标记的非相基因型矩阵(n × m)，编码为{-1，0，1}。允许使用小数点(表明已填充过)和缺失值(NA)。
min.MAF	最小等位基因频率。A矩阵对稀有等位基因不敏感，因此默认情况下只删除单态标记。
max.missin	最大缺失值比例；默认完全删除丢失的标记
impute.method	有两种选择。默认值为“Mean”，表示每个标记的平均值。“EM”选项使用EM算法(请参阅详细信息)。
tol	指定EM算法的收敛标准(请参阅详细信息).
n.core	指定用于EM算法的并行核心数量(仅在UNIX命令行中使用)。
shrink	设置shrink=FALSE以禁用收缩估计。有关如何启用收缩估算的详细信息，请参阅。
return.imputed	如果为True，则返回计算的标记矩阵。

---

细节

• 在高密度标记矩阵中，关系矩阵被估计为 $A=\frac{W{W}^{\prime }}{c}$,其中$W_{ik}=X_{ik}+1-2p_k$，$P_k$是标记k的1个等位基因的频率。通过使用归一化常数$c=2{\sum }_k{p}_k(1-p_k)$，对角线元素的平均值为$1+f$¹。

• EM补偿算法基于多变量正态分布，设计用于GBS(genotyping-by-sequencing)标记，GBS倾向于高密度，但有大量缺失的数据²。当$RMSerror=n^{-1}||A_t-A_{t-1}|{|}_2<tol$时，EM算法在迭代 $t$ 处停止。

• 收缩估计可以提高全基因组标记辅助选择的准确性，特别是在低标记密度的情况下¹。收缩强度的范围从$0$(无收缩)到$1(A=(1+f)I)$。有两种算法可用于估计收缩强度。

  • 第一个是Endelman和Jannink(2012)中描述的方法，由$Shrink=list(method=“EJ”)$指定¹。
  • 第二种方法是将随机的标记样本指定为模拟QTL，然后基于QTL将A矩阵与基于剩余标记的A矩阵进行回归。³⁴ 回归方法由$shrink=list(method=“reg”，n.qtl=100，n.iter=5)$指定，其中参数$n.qtl$和$n.iter$可以分别调整模拟QTL的数量和迭代次数。

• 收缩和EM补偿选项是为相反的场景(低密度与高密度)设计的，不能同时使用。当使用EM算法时，推定的等位基因可以位于区间[-1，1]之外。不满足min.MAF和max.Missing标准的多态标记不被补偿。

---

值

• 如果 $return.imputed=FALSE$,返回 $n\times n$ 的加性关系矩阵
• 如果 $return.imputed=TRUE$ ,该函数返回包含以下内容的列表
  • $A – A矩阵
  • $imputed – 补偿标记矩阵

---

范例

#random population of 200 lines with 1000 markers
library(rrBLUP)
X <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  X[i,] <- ifelse(runif(1000)<0.5,-1,1)
}

A <- A.mat(X)

---

---

1.2 GWAS

基于混合模型进行全基因组关联分析⁵：

$$y=X\beta +Zg+S\tau +ϵ$$
其中$\beta$是固定效应向量，既可以对环境因素进行建模，也可以对群体结构进行建模。变量$g$将为每个株系的遗传背景建立以下具有随机效应的模型：$V_{ar}[g]=K{\sigma }^2$。变量$\tau$将加性SNP效应建模为固定效应。残差是：$V_{ar}[ϵ]=K{\sigma }_e^2$

用法

GWAS(pheno, geno, fixed=NULL, K=NULL, n.PC=0,min.MAF=0.05, n.core=1, P3D=TRUE, plot=TRUE)

---

参数

选项	意义
pheno	数据框，其中第一列是品种名(gid)。其余的列可以是表型，也可以是固定效果的级别。任何没有被指定为固定效应的列都被假定为表型。
geno	第一列中具有标记名称的数据框。第二列和第三列分别包含染色体和定位位置(bp或cm)，只有当plot=true时才使用它们来绘制曼哈顿图。如果标记未映射，只需为这两列使用占位符。列4和更高的列包含每行的分数，编码为{-1，0，1}={aa，Aa，AA}。允许使用小数(填充)和缺失(NA)值。列名必须与pheno数据框中的行名相匹配。
fixed	一个字符串数组，其中包含应作为(分类)固定效果包含在混合模型中的列的名称。
K	由于多基因效应导致的行间协方差的血缘关系矩阵。如果未通过，则使用A.mat从标记计算。
n.PC	要包含为固定效果的主分量数。默认值为0(等于K模型)。
min.MAF	指定最小次要等位基因频率(MAF)。如果有标记的MAF小于min.MAF，则为其分配零分。
n.core	设置n.core>1将在具有多核的计算机上启用并行执行(仅在UNIX命令行使用)。
P3D	当P3D=TRUE时，REML只估计方差分量一次，模型中没有任何标记。当P3D=FALSE时，通过REML分别估计每个标记的方差分量。
plot	当plot=TRUE时，生成QQ和曼哈顿图。

---

细节

• 对于表型来自不同环境的不平衡设计，可以使用fixed选项对环境平均值进行建模(例如，如果表型数据框中的列称为“env”，则fixed=“env”)。当包括主分量时(P+K模型)，载荷由K矩阵的特征值分解确定。
• 术语P3D(先前确定的总体参数)是由Zhang等人引入的⁶。当P3D=FALSE时，此函数等同于具有REML的EMMA⁷.当P3D=TRUE时，它等同于EMMAX⁸。与P3D=FALSE相比，P3D=TRUE选项速度更快，但可能低估了重要性。
• 曼哈顿图像中的虚线对应的是FDR=0.05，并且是使用qvalue包计算的⁹。通过插值确定对应于q-valu =0.05的p值。如果没有小于0.05的q值，则省略虚线。

值

返回一个数据框，其中前三列是标记名称、染色体和位置，后续列是性状的标记分数($-lo{g}_{10}p$)。

#random population of 200 lines with 1000 markers
library(rrBLUP)
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) #选择5个QTL
u <- rep(0,1000)
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 = TRUE)

---

---

1.3 kin.blup

描述

用G-BLUP预测基因型值，其中基因型协方差G可以是加性的，也可以是基于高斯核的。

用法

kin.blup(data,geno,pheno,GAUSS=FALSE,K=NULL,fixed=NULL,covariate=NULL,
PEV=FALSE,n.core=1,theta.seq=NULL)

参数

选项	意义
data	包含表型、基因型标识符和任何环境变量列的数据框。
geno	数据框中包含基因型标识符的列名的字符串。
pheno	数据框中包含表型的列名的字符串。
GAUSS	要使用高斯核对遗传协方差进行建模，请设置GAUSS=TRUE并传递K的欧几里德距离(见下文)。
K	有三个选项可用于指定亲属关系:1)如果 $K=NULL$,基因型被认为是独立的($G=I{V}_g$)；2)对于育种值预测，设置$GAUSS=FALSE$并使用$K$的相加关系矩阵来创建模型($G=K{V}_g$)；3）对于高斯内核，设置$GAUSS=TRUE$并传递$K$的欧几里得距离矩阵来创建模型$G_{ij}=e^{-(\frac{K_{ij}}{\theta }{)}^2}{V}_g$
fixed	一个字符串数组，其中包含应作为(分类)固定效果包含在混合模型中的列的名称。
covariate	包含应作为协变量包含在混合模型中的列名的字符串数组。
PEV	当PEV=TRUE时，该函数返回遗传型值的预测误差方差（$PE{V}_i=V_{ar}[g_i^{\ast }-g_i]$）
n.core	指定用于并行执行高斯内核方法的核心数量(仅用于UNIX命令行)。
theta.seq	高斯核的比例参数是通过最大化值网格上的受限对数似然率来设置的。默认情况下，栅格是通过将间隔$(0，max(K)]$分割为 10个点来构建的。将数值数组传递给此变量(theta.seq=“theta Sequence”)将指定一组不同的网格点(例如，对于较大的问题，您可能希望少于10个)。

细节

• 此函数是Mixed.solve的包装器，从而求解以下形式的混合模型：
  $$y=X\beta +[Z0{]}_g+ϵ$$
  其中，$\beta$是固定效应向量，$g$是具有协变量$G=V_{ar}[g]$的随机表型值向量，残差为$V_{ar}[{ϵ}_i]=R_i{\sigma }_e^2$,默认$R_i=1$。遗传值的设计矩阵已经被划分，以说明并不是所有的品系都需要表型(例如，对于基因组选择来说)。与mixed.solve不同，这一函数并不会返回固定效应的估计值，只返回基因型值的BLUP解。它的设计目的是取代kinship.BLUP，并使用户不必显式构造设计矩阵。方差分量由REML估计，并且为K中的每个条目返回BLUP值，无论它是否已有表型。K的行名必须与表现型系数据帧中的基因型标签相匹配;
  丢失的表型(NA)被简单地省略了。
• 与它的前身不同，这个函数不直接处理标记数据。对于育种值的预测，用户必须提供一个关系矩阵，该矩阵可以从带有A.mat的标记中计算出来。对于高斯核预测，通过K的欧氏距离矩阵，它可以用dist来计算。
• 在混合模型的术语中，“固定”变量和“协变量”变量都是固定效应( 上式中的$\beta$):前者被视为具有不同层次的因子，而后者是连续的，每个变量有一个系数。总体均值作为一个固定效应自动包括在内。预测误差方差(PEV)是BLUPs的SE的平方(见mix .solve)，可用于根据$r_i^2=1-\frac{PE{V}_i}{V_g{K}_{ii}}$估计BLUP预测的预期精度。

值

函数总是返回以下值：

$$V_g$	遗传方差的REML估计
$$V_e$	误差方差的REML估计
$$g$	遗传价值的BLUP值
$$resid$	残差
$$pred$	预测的遗传值，平均于固定效应
IF $PEV=TRUE$,返回值也会包括$$PEV$
$$PEV$	遗传值的预测误差方差
if $GAUSS=TRUE$,返回值也会包括 $$profile$
$$profile$	高斯核中尺度参数的对数似然分布

范例

#加载rrBLUP包
library(rrBLUP)
#随机生成具有1000个标记，由200个株系组成的群体
M = matrix(rep(0,200),200,1000)
for (i in 1:200) {
  M[i,] = ifelse(runif(1000)<0.5,-1,1)
}
rownames(M) = 1:200
A <- A.mat(M)

#随机生成表型值
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5 #遗传力
y <- g + rnorm(200, mean = 0,sd=sqrt((1-h2)/h2*var(g)))
data = data.frame(y=y,gid=1:200)

#预测育种值
ans <- kin.blup(data = data,geno = "gid",pheno = "y",K=A)
accuracy = cor(g,ans$g)

Kinship.BLUP

描述

这一功能已被$kin.blup$取代。请参考它的帮助页面。

1.6 mixed.solve函数

该函数主要用于求解模型参数。

描述

计算方程的混合模型的最大似然(ML/REML)解
$$y=X\beta +Zu+ϵ$$
其中$\beta$是固定效应向量，$u$是是具有协方差$V_{ar}[u]=K{\sigma }_u^2$的随机效应向量；残差是$V_{ar}[ϵ]=I{\sigma }_e^2$。这类混合模型除了残差外，还有一个单一的方差成分，与岭回归有密切的关系(岭回归参数：$\lambda =\frac{{\sigma }_e^2}{{\sigma }_u^2}$)

使用方法

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

参数

y	观测值向量$(n\times 1)$,遗漏的值(NA)被省略，同时也是X和Z对应的行。
Z	随机效应的涉及矩阵$(n\times m)$，如果没有通过，假定为单位矩阵
K	随机效应的协方差矩阵$(m\times m)$,一定是半正定义的。如果没有通过，则假定为单位矩阵。
X	固定效应的设计矩阵$(m\times m)$,如果没有通过，则使用一个向量1来建模截距。X必须是全列秩(意味着是可估计的)。
method	指定是否使用完整(“ML”)或限制(“REML”)最大似然方法。
bounds	包含两个元素的数组，它们指定岭回归参数的上下边界。
SE	如果为真，则计算标准误差
return.Hinv	如果为真，函数返回$H=ZK{Z}^{\prime }+\lambda I$的逆函数。这对GWAS很有用。

细节

这个函数可以用来预测标记效应或育种值(见例子)。数值方法基于$ZK{Z}^{\prime }$和$SZK{Z}^{\prime }S$的谱分解，其中$S=I-X(X^{\prime }X{)}^{-1}{X}^{\prime }$是X零空间的投影算子(Kang et al.， 2008)。该算法生成逆表型协方差矩阵$V^{-1}$，然后可以使用标准公式分别计算固定效应和随机效应的BLUE和BLUP解(Searle et al. 1992):
$$BLUE(\beta )={\beta }^{\ast }=(X^{\prime }{V}^{-1}X{)}^{-1}{X}^{\prime }{V}^{-1}y$$
$$BLUP(u)=u^{\ast }={\sigma }_u^{2}K{Z}^{\prime }{V}^{-1}(y-X{\beta }^{\ast })$$

标准误差计算为下列矩阵对角线元素的平方根(Searle等，1992):
$$V_{ar}[{\beta }^{\ast }]=(X^{\prime }{V}^{-1}X{)}^{-1}$$

$$V_{ar}[u^{\ast }-u]=K{\sigma }_u^2-{\sigma }_u^{4}K{Z}^{\prime }{V}^{-1}ZK+{\sigma }_u^{4}K{Z}^{\prime }{V}^{-1}X{V}_{ar}[{\beta }^{\ast }]X^{\prime }{V}_{-1}ZK$$
对于K = I的标记效果，如果不传递K，函数会比用户传递单位矩阵运行得更快。

值

如果 $SE=FALSE$,函数返回值包括:

$$V_u$	${\sigma }_u^2$的估计值
$$V_e$	${\sigma }_e^2$的估计值
$$beta$	BLUE($\beta$)
$$u$	BLUEP(u)
$$LL$	最大化对数似然(完全或有限，取决于方法)

如果 $SE=TRUE$,函数返回值包括:

$$beta.SE$	$\beta$的标准误差
$$u.SE$	$u^{\ast }-u$的标准误差

如果 $return.Hinv=TRUE$,函数返回值包括:

$$Hinv$	H的逆函数

范例

#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)

---

1. 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 ↩︎ ↩︎ ↩︎

2. Mueller et al. 2015. Shrinkage estimation of the genomic relationship matrix can improve genomic estimated breeding values in the training set. Theor Appl Genet doi: 10.1007/s00122-015-2464-6 ↩︎

3. Poland, J., J. Endelman et al. 2012. Genomic selection in wheat breeding using genotyping-bysequencing.Plant Genome 5:103-113. doi: 10.3835/plantgenome2012.06.0006 ↩︎

4. Yang et al. 2010. Common SNPs explain a large proportion of the heritability for human height.Nat. Genetics 42:565-569. ↩︎

5. Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208. ↩︎

6. Zhang et al. 2010. Mixed linear model approach adapted for genome-wide association studies. Nat.Genet. 42:355-360. ↩︎

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

8. Kang et al. 2010. Variance component model to account for sample structure in genome-wide association studies. Nat. Genet. 42:348-354. ↩︎

9. Storey and Tibshirani. 2003. Statistical significance for genome-wide studies. PNAS 100:9440-9445. ↩︎