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）。
- best：估計誤差最小，估計與真實育種值相關最大；
- linear：估計基於線性模型，假設個體性能是遺傳和非遺傳效應的總和；
- unbiased：估計過程中估計值是無偏的，估計與真實育種值間平均差為零，即所有可能估計值的平均值等於真實育種值；
- prediction：一個體將來作為親本的預測育種值。
越解釋越混亂，這些基礎知識需要好好鞏固鞏固。
先看間接法的模型公式：
- y為表型向量；
- X為固定效應系數矩陣，其長度為訓練群個體數目，元素值均為1的向量；
- b為固定效應，即為訓練群表型均值；
- Zi為第i個位點數字化基因型向量（如：{0,1,2}或{-1,0,1}之類的編碼，不同編碼方式有無差異？但rrBLUP必須為{-1,0,1}編碼方式），之和即為標記編碼關聯矩陣；
- gi為第i個位點效應值，需要根據模型估計出的分子標記效應向量；
- e為殘余誤差，服從分布N~(0, Iσe2)。
關鍵在於求解標記效應：
- σgi2為第i個標記方差，直接與性狀遺傳構建相關。
間接法重難點在於如何對參數的先驗分布（即對gi及其方差服從的分布）進行合理的假設。RRBLUP模型假設所有標記都具有效應，且來源於同一個分布，即σgi2相等（理論上RRBLUP與GBLUP方法是等價的）。但實際上所有標記不會都具有效應，標記方差也會不同，因而出現了各種假設的Bayes方法，這樣就帶來更多待估參數，提高准確性的同時也增加了計算量。
2. 實操
采用rrBLUP Package 來進行實操學習。
2.1 rrBLUP包簡介
rrBLUP包主要兩個函數：
A.mat函數用來過濾和填充基因型。
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函數直接求解模型參數。
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函數返回值（列表）：
- Vu：σ^2_u的估計值
- Ve：σ^2_e的估計值
- beta：BLUE(β)，即訓練群表型均值
- u：BLUP(u)，即預測值
- LL：maximized log-likelihood（ML/REML值）
如果設了標准誤參數，則會返回：
- beta.SE：β標准誤
- u.SE：u^*-u
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)
- 過濾和填充
#####
#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)
- 數據划分訓練集和驗證集
#######
#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
- 循環n次
#### 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)
3. 補充說明
關於模型
需要提前清洗數據，並編碼基因型數據。
關於交叉驗證
rrBLUP預測准確性與訓練群和驗證群大小、標記數目以及遺傳力等有關。可以看到，上面循環500次的重復性並不是很好。
使用K折交叉驗證的方法可能會好些。關於交叉驗證，一般有三種方法：
- 簡單交叉驗證：就是上面的方法，隨機將樣本分為一定比例的訓練集和測試集，然后訓練和驗證模型，再把樣本打亂，重新選擇訓練集和測試集，繼續訓練數據和驗證模型。一般用於初步模型的建立。
- K折交叉驗證：先將數據集隨機划分為K個大小相近的互斥子集，每次隨機選擇K-1份作為訓練集，剩下一份做測試集。當這一輪完成后，下一輪重新隨機選擇K-1份來訓練數據，最后多輪結果取均值。交叉驗證法評估結果的穩定性和保真性在很大程度上取決於K的取值。
- 留一交叉驗證：是K折交叉驗證的特例，即K等於樣本數N。每次N-1樣本訓練，留一個樣本驗證。一般用於樣本量很少的情況（如小於50）。
- booststrapping自助法也算是一種比較特殊的交叉驗證法：在m個樣本中隨機采一個樣放入訓練集，再將樣本放回，重復采集m次得到m個樣本組成的訓練集（可能有重復樣本），同時用原始的m個樣本做測試集。一般用於樣本量非常少的情況（如小於20）。