R包-rrBLUP

使用 rrBLUP 包计算 SNP 效用的方法。

表型和基因型处理

1. 表型文件中必须剔除存在缺失的行。

2. 基因型质控填充，将位点编码必须为 {-1,0,1} 。

3. 表型文件和基因型文件中的个体顺序必须保持一致（只能使用既有表型又有基因型的个体数据）。

求解

这里我们主要用到的就是 rrBLUP 包中的 mixed.solve 函数，这里还可以手动填充固定效应（必须列满秩），用法如下

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

参数：

y: 的表型向量，不能有缺失 （NA），必须与 X 和 Z 的行顺序保持一致。

Z: 的随机效应矩阵，如果不设置，默认为单位矩阵。

K: 的随机效应协方差矩阵，必须为半正定矩阵。如果不设置，默认为单位矩阵。

X: 的固定效应设计矩阵，如果不设置，则为 1 向量。X 必须为列满秩序矩阵。

method: “ML”, “REML”

bounds: 两个元素的数组，表明 ridge 参数的下限和上限。

SE ：如果 TRUE , 则计算标准误

return.Hinv: 如果 TRUE, 返回 reverse of ，这对 GWAS 有用。

正常来说，如果固定效应为均值，那么就只需要设置 y 和 Z (填充后的基因型矩阵，由**{-1,0,1}**构成)即可。

out <- mixed.solve(y, Z=marker)

返回值如下， 其中我们需要的 snp 效应就是$u。

举例

下面是一个实际应用的脚本例子。

这里的两个输入文件 rrblup.phe 与 rrblup.raw 是事先处理好的，二者中的个体顺序保持一致。 而 rrblup.raw 中的内容是通过 plink 软件的 --recodeA 选项生成的。

输出文件夹 rrblup_results 中的每一个文件表示一个结果，其中 u.csv 就是我们需要的 SNP 效应值。

经过手动比较，rrBLUP 和 GBLUP 计算得到的个体 GEBV 相关系数为 1， 没有问题。

# 运行 rrblup
library("rrBLUP")

raw_file = read.table("rrblup.raw", header=T)
phe_file = read.table("rrblup.phe", header=F)

# 表型
y = phe_file[,2]

# 基因型
X = as.matrix(raw_file[,7:ncol(raw_file)])

X = X -1

# rrblup
out <- mixed.solve(y, Z=X)

dir.create("rrblup_results", showWarnings = FALSE) 
setwd("rrblup_results")
for (i in 1:length(out)){
    list_name = names(out)[i]
    out_file_name = paste0(list_name,".csv")   
    write.csv(out[i],file=out_file_name,quote = FALSE, row.names=FALSE)
}