遗传学常用分析笔记

【GS模型】全基因组选择之rrBLUP

2020-12-09  本文已影响0人  生物信息与育种

[toc]

1. 理论

rrBLUP是基因组选择最常用的模型之一，也是间接法模型的代表。回顾一下，所谓间接法是指：在参考群中估计标记效应，再结合预测群的基因型信息将标记效应累加，最终获得预测群的个体估计育种值。而直接法则是指：将个体作为随机效应，参考群体和预测群体遗传信息构建的亲缘关系矩阵作为方差协方差矩阵，通过迭代法估计方差组分，进而求解混合模型获取待预测个体的估计育种值。简言之，直接法是通过构建A/G/H等矩阵求解育种值，间接法是通过计算标记效应来获得育种值。

RRBLUP全称“ridge regression best linear unbiased prediction”，即岭回归最佳线性无偏预测。光从长串名字看，里面涉及到大量的统计学基本概念，可分为两部分看。一是Ridge Regression，岭回归是一种改良的最小二乘估计法，专用于解决共线性数据分析的有偏估计回归方法，通过在损失函数中加上一个正则化项，放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更符合实际；二是BLUP，对一个随机效应（如个体育种值）的预测具有线性（预测量是样本观察值的线性函数）、无偏（预测量的数学期望等于随机效应本身的数学期望）和预测误差方差最小等统计学性质。

最佳线性无偏预测（ best linear unbiased prediction，BLUP）是线性模型中用来评估随机效应的，等同于固定效应中的最佳线性无偏估计（best linear unbiased estimates, BLUE）。

越解释越混乱，这些基础知识需要好好巩固巩固。

先看间接法的模型公式：

image.png

关键在于求解标记效应：

image.png

间接法重难点在于如何对参数的先验分布（即对gi及其方差服从的分布）进行合理的假设。RRBLUP模型假设所有标记都具有效应，且来源于同一个分布，即σgi2相等（理论上RRBLUP与GBLUP方法是等价的）。但实际上所有标记不会都具有效应，标记方差也会不同，因而出现了各种假设的Bayes方法，这样就带来更多待估参数，提高准确性的同时也增加了计算量。

2. 实操

采用rrBLUP Package 来进行实操学习。

2.1 rrBLUP包简介

rrBLUP包主要两个函数：

A.mat(
    X, #编码为{-1,0,1}的标记矩阵，可包含小数点（表明已填充过）和NA
    min.MAF=NULL, #最小等位基因频率
    max.missing=NULL, #最大缺失值比例
    impute.method="mean", #填充方法，mean/EM
    tol=0.02, #EM算法的收敛标准
    n.core=1, #EM算法的并行核心数（仅unix）
    shrink=FALSE, #收缩估计
    return.imputed=FALSE #返回填充标记矩阵
)

返回加性效应关系矩阵（即kinship），填充后的分子标记矩阵。

mixed.solve(
       y,  #观测值
       Z=NULL, #随机效应矩阵
       K=NULL, #随机效应协变量矩阵
       X=NULL,  #固定效应矩阵
       method="REML",  #最大似然估计方法，ML/REML
       bounds=c(1e-09, 1e+09),   #岭回归参数两元素边界
       SE=FALSE, #是否计算标准误
       return.Hinv=FALSE #是否H取逆，一般在GWAS中用，忽略之
)

mixed.solve函数返回值（列表）：

如果设了标准误参数，则会返回：

2.2 实操

#load in the sample files
library(rrBLUP)

#load the Markers and Phenotypes
Markers <- as.matrix(read.table(file="snp.txt"), header=F)
dim(Markers)
Markers[1:10,1:10]

Pheno <-as.matrix(read.table(file ="traits.txt", header=TRUE))
dim(Pheno)
head(Pheno)
image.png
image.png 
#####
#what if markers are NA?
#impute with A.mat
#remove markers with more than 50% missing data
impute=A.mat(Markers,max.missing=0.5,impute.method="mean",return.imputed=T) 
#均值填充导致标记不止0，-1，1三类数
Markers_impute2=impute$imputed

dim(Markers)
dim(Markers_impute2)
image.png 
#######
#define the training and test populations
#training-60% validation-40%
train= as.matrix(sample(1:96, 58))
test<-setdiff(1:96,train)
Pheno_train=Pheno[train,]
m_train=Markers_impute2[train,]
Pheno_valid=Pheno[test,]
m_valid=Markers_impute2[test,]

########
yield=(Pheno_train[,1])
yield_answer<-mixed.solve(yield, Z=m_train, K=NULL, SE = FALSE, return.Hinv=FALSE)
# yield_answer<-mixed.solve(yield, Z=m_train, K=NULL, SE = TRUE, return.Hinv=FALSE)
YLD = yield_answer$u
e = as.matrix(YLD)
pred_yield_valid =  m_valid %*% e
pred_yield=(pred_yield_valid[,1])+yield_answer$beta
pred_yield
yield_valid = Pheno_valid[,1]
YLD_accuracy <-cor(pred_yield_valid, yield_valid, use="complete" )
YLD_accuracy
image.png
image.png 
#### cross validation for many cycles for yield only
traits=1
cycles=500
accuracy = matrix(nrow=cycles, ncol=traits)
for(r in 1:cycles){
  train= as.matrix(sample(1:96, 58))
  test<-setdiff(1:96,train)
  Pheno_train=Pheno[train,]
  m_train=Markers_impute2[train,]
  Pheno_valid=Pheno[test,]
  m_valid=Markers_impute2[test,]
  
  yield=(Pheno_train[,1])
  yield_answer<-mixed.solve(yield, Z=m_train, K=NULL, SE = FALSE, return.Hinv=FALSE)
  YLD = yield_answer$u
  e = as.matrix(YLD)
  pred_yield_valid =  m_valid %*% e
  pred_yield=(pred_yield_valid[,1])+yield_answer$beta
  pred_yield
  yield_valid = Pheno_valid[,1]
  accuracy[r,1] <-cor(pred_yield_valid, yield_valid, use="complete" )
}
mean(accuracy)
image.png

3. 补充说明

关于模型

需要提前清洗数据，并编码基因型数据。


image.png

关于交叉验证

rrBLUP预测准确性与训练群和验证群大小、标记数目以及遗传力等有关。可以看到，上面循环500次的重复性并不是很好。


image.png

使用K折交叉验证的方法可能会好些。关于交叉验证，一般有三种方法：

参考资料

上一篇下一篇

猜你喜欢

热点阅读