||
这篇是基因组选择的理论加实践,因为我看到一句话,Talk is cheap. Show me the code,很有感触,有感而写。使用的包是R的sommer和asreml,其实强健的还是成熟的软件,比如DMU,BLUPF90,PIBLUP,ASreml等,但sommer作为基本功能的演示,非常合适。
基因组选择(Genomic Selection, GS), 利用覆盖全基因组的高密度分子遗传标记进行的标记辅助选择.
选择进展的定义
$$ 选择进展 = \frac{选择强度 \times 选择准确性\times 遗传标准差}{世代间隔}$$
基于系谱的动物模型
局限:
分子标记辅助育种(MAS)
局限:
全基因组选择
优点:
三种方法的区别
别说话, 上代码
数据描述
软件使用
# 载入软件包
if (!requireNamespace("devtools")) install.packages("devtools")
library(devtools)
if (!requireNamespace("learnasreml")) install_github("dengfei2013/learnasreml")
library(learnasreml)
if (!requireNamespace("sommer")) install.packages("sommer")
library(sommer)
# 载入数据
data("gs_geno")
data("gs_phe")
# 动物模型animal model
gphe = gs_phe[gs_phe$G %in%c(3:5), ] # 这里提取测序个体
gphe$Progeny = as.factor(as.character(gphe$Progeny)) # 将Progeny转化为因子
ainv = asreml.Ainverse(gphe[,1:3])$ginv # 计算A逆矩阵
moda_as = asreml(Phen ~ Sex + G, random = ~ ped(Progeny),ginverse= list(Progeny=ainv),data = gphe) # 构建模型
summary(moda_as)$varcomp #结果
# GBLUP sommer的操作代码
gdat = gs_geno # 重命名数据
gdat[1:10,1:10] # 预览数据
G = A.mat(gdat-1) # sommer计算G矩阵时, 支持的是-1,0,1, 所以需要转化
time_sommer = system.time({
modg_mm = mmer2(Phen ~ Sex + G, random = ~ g(Progeny),G = list(Progeny=G),data = gphe) # 构建模型
})
# GBLUP asreml的操作代码
library(asreml)
diag(G) = diag(G)+0.01 # 矩阵奇异, 对角线+0.01防止奇异
ginv = write_relation_matrix(G,type="ginv") # 转化为三元组形式
attr(ginv,"rowNames") = rownames(G) # 对ginv进行rowNames的设置
time_asreml = system.time({
modg_as = asreml(Phen ~ Sex + G, random = ~ giv(Progeny),ginverse= list(Progeny=ginv),data = gphe,workspace=1e9) # 构建模型
})
# 比较sommer和asreml结果
summary(modg_mm)
summary(modg_as)$varcomp
# 比较运行时间
time_sommer
time_asreml
# 计算准确性:asreml
blup_as = as.data.frame(coef(modg_as)$random)
re1 = cbind(gphe,blup_as)
cor(re1$Polygene,re1$effect)
# 计算准确性:sommer
blup_mm = as.data.frame(randef(modg_mm))
head(blup_mm)
names(blup_mm) = "blup_mm"
re2 = cbind(gphe,blup_mm)
cor(re1$Polygene,re2$blup_mm)
方差组分结果比较
> # 比较sommer和asreml结果
> summary(modg_mm)$var.comp.table
VarComp VarCompSE Zratio
g(Progeny).Phen-Phen 0.1947930 0.04421259 4.405827
units.Phen-Phen 0.8076997 0.04320582 18.694236
> summary(modg_as)$varcomp
gamma component std.error z.ratio constraint
giv(Progeny).giv 0.2417449 0.1947875 0.04518119 4.311252 Positive
R!variance 1.0000000 0.8057563 0.04384776 18.376227 Positive
可以看出, sommer和asreml方差组分结果基本一致, Vg=0.194, Ve=0.805
运行时间比较
> # 比较运行时间
> time_sommer
用户 系统 流逝
33.53 3.88 9.47
> time_asreml
用户 系统 流逝
19.00 3.90 23.03
可以看到,sommer运行时间很快,sommer支持的是G矩阵类型,如果使用pedigree,有点麻烦。asreml-r版慢一点,个人体会,asreml-w更快一点。
准确性比较
> # 计算准确性:asreml
> blup_as = as.data.frame(coef(modg_as)$random)
> re1 = cbind(gphe,blup_as)
> cor(re1$Polygene,re1$effect)
[1] 0.6338289
>
> # 计算准确性:sommer
> blup_mm = as.data.frame(randef(modg_mm))
Returning object of class 'list' where each element correspond to one random effect.
> names(blup_mm) = "blup_mm"
> re2 = cbind(gphe,blup_mm)
> cor(re1$Polygene,re2$blup_mm)
[1] 0.6327671
因为这是模拟数据,有真值(True breeding value),如果是真实数据,需要用交叉验证的手段评价准确性。
Meuwissen T H , Hayes B J , Goddard M E . Prediction of total genetic value using genome-wide dense marker maps[J]. Genetics, 2001, 157(4):1819-29.
Vanraden P M. Efficient methods to compute genomic predictions[J]. Journal of Dairy Science, 2008, 91(11):4414-4423.
Legarra A, Aguilar I, Misztal I. A relationship matrix including full pedigree and genomic information[J]. Journal of Dairy Science, 2009, 92(9):4656-4663.
Aguilar I, Misztal I, Johnson D L, et al. Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score.[J]. Journal of Dairy Science, 2010, 93(2):743-752.
Ishizaka A, Labib A. Review of the main developments in the analytic hierarchy process[J]. Expert Systems with Applications, 2011, 38(11):14336-14345.
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-19 22:36
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社