全基因组选择介绍及实践-1

2019-06-13 20:23:44 浏览数 (2)

1, 序言

这篇是基因组选择的理论加实践，因为我看到一句话，Talk is cheap. Show me the code，很有感触，有感而写。使用的包是R的sommer和asreml，其实强健的还是成熟的软件，比如DMU，BLUPF90，PIBLUP，ASreml等，但sommer作为基本功能的演示，非常合适。

2, 定义

基因组选择(Genomic Selection, GS), 利用覆盖全基因组的高密度分子遗传标记进行的标记辅助选择.

  • 表型选择, 这对于遗传力高的性状选择有效, 即表型可以遗传比例较大的性状, 比如植物里面千粒重, 比如动物里面的初生重等.
  • EBV selection, 应用系谱的动物模型
  • Marker assisted selection(MAS), 分子标记辅助选择, 对于有主效QTL或者主效基因的性状, 有优势.
  • Genomic selection, 全基因组选择

选择进展的定义

  • 选择强度: Intensity of selection
  • 选择准确性: Accuracy of selection
  • 遗传标准差: Genetic standard deviation
  • 世代间隔: Generation interval

3, 不同选择方法的比较

基于系谱的动物模型 局限:

  • 对低遗传力性状效果较差
  • 不容易度量的性状(胴体性状, 肉质性状)效果较差
  • 不能早期度量的性状, 效果较差

分子标记辅助育种(MAS) 局限:

  • 需要先对主效基因或者QTL进行检测
  • 不同群体变化较大
  • 标记可解释的遗传变异百分比较低
  • 在动物育种中的应用非常有限

全基因组选择 优点:

  • 无需进行主效基因或者QTL的检测
  • 不依赖于表型信息(候选群)
  • 能够捕获基因组中的全部变异
  • 对于低遗传力, 难以度量的性状提升效果明显

4, 基因组选择流程

  • 1, 建立参考群体
    • 获得每个个体的性状表型值
    • 测定每个个体的SNP基因型(使用SNP芯片)
    • 估计SNP效应值
  • 2, 在候选群体中进行基因组选择
    • 测定候选个体的SNP基因型
    • 计算个体的GEBV
    • 依据GEBV进行选择

动物模型, GBLUP方法和Single-setp方法

三种方法的区别

  • 动物模型是利用的系谱构建的A矩阵
  • GBLUP是利用基因组信息构建的G矩阵
  • 一步法(single-setp)是利用系谱和基因组信息构建的H矩阵 5，其它方法 除了GBLUP和Single-step, 还有其它方法用于基因组选择, 不过应用最广最强建最有优势的还是GBLUP和Single-step.
  • RR-BLUP(Rirdge Regression BLUP)
  • Bayesian method
    • Bayes A
    • Bayes B
    • Bayes C
    • Bayesian Lasso
  • Principal component regression
  • Semiparametric procedures
  • Machine learning methods

6，Talk is cheap, show me the code

别说话, 上代码

数据描述

  • 数据使用QMSim软件模拟而成
  • 共有2000个体
  • 1200个有基因型
  • 共有5个世代, 3,4为参考群, 5为候选群, 每个世代400个体
  • 系谱数据为50k
  • 遗传力为0.3

软件使用

  • QMSim模拟数据
  • R包: learnasreml
  • R包: sommer
  • R包：asreml
代码语言：javascript复制
# 载入软件包
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)

方差组分结果比较

代码语言：javascript复制
> # 比较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

运行时间比较

代码语言：javascript复制
> # 比较运行时间
> time_sommer
 用户  系统  流逝
33.53  3.88  9.47
> time_asreml
 用户  系统  流逝
19.00  3.90 23.03

可以看到，sommer运行时间很快，sommer支持的是G矩阵类型，如果使用pedigree，有点麻烦。asreml-r版慢一点，个人体会，asreml-w更快一点。

准确性比较

代码语言：javascript复制
> # 计算准确性: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），如果是真实数据，需要用交叉验证的手段评价准确性。

网站 量子密钥管理服务

0 人点赞