来自Michael Goddard 和Ruidong Xiang 培训。
介绍
大多数经济性状是复杂或数量性状，被多基因与环境控制。
传统计算的EBVs来自表型和系谱，现在增加了基因型数据。
2007 发布了SNPs芯片，使GWAS变为常用手段。理论是基于SNP可与引起变异的位点处于连锁不平衡（LD）
同时预测基因育种值：
image.png
使用贝叶斯方法：
image.png
image.png
需要假设SNP的先验分布：
其最推荐的分布为 mixture of N (Bayes R)
image.png
例子：
image.png
使用MCMC（Gibbs sampling）:
image.png
k为多变量分布，
image.png
image.png
人类的例子：
image.png
预测基因育种值：
image.png
牛的产奶量， BayesR的结果更好一些：
image.png
我们想SNP突因果突变变更近
image.png
怎么找到好的SNP：
image.png
统计方法：应该同时使用所有的SNP比较好：
image.png
Bayes R中 效应分布:长尾分布
image.png
更准确的预测：
image.png
例子：
image.png
羊：
image.png
练习
两部分
MCMC 抽样器在线性回归使用
贝叶斯在GS使用
Gibbs sampler虚假代码
image.png
Gibbs抽样器在线性回归：
image.png
a,b, V的后验分布，也是需要估计的参数：
image.png
###############################################################################
#====================code Gibbs lm====================================
###############################################################################
.libPaths("C:/Users/Public/Rlib")
setwd("C:/Users/rx01/teach/stats_lecture1_Nov2021")
library(ggplot2)
set.seed(2)
#simulate a trait
phe <- rnorm(100,1000,500)
head(phe,10)
SNP <- sample(rep(c(0,1,2),500),100)
head(SNP,50)
#----set up function
mcmc.lr <- function(Y,X,Intcept,Slope,s2,mu0,alpha.0,beta.0,s20,nit){
n <- length(Y)
sampl <- matrix(0,nit,3)
colnames(sampl) <- c("Intcept","Slope","sigma2")
sampl[1,] <- c(Intcept,Slope,s2)
for(iter in 2:nit){
#sample Intcept
V <- n/s2+mu0/s20
M <- sum(Y-X*Slope)/s2+1/s20
Intcept <- rnorm(1,M/V,1/sqrt(V))
#sample Slope
V <- sum(X^2)/s2+mu0/s20
M <- sum(X*(Y-Intcept))/s2+1/s20
Slope <- rnorm(1,M/V,1/sqrt(V))
#sample s2|mu,Y,Z
A <- n/2 + alpha.0
B <- sum((Y-Intcept-X*Slope)^2)/2 + beta.0
s2 <- 1/rgamma(1,A,B)
#keep track of the results
sampl[iter,] <- c(Intcept,Slope,s2)
}
return(sampl)
}
### Priors
mu0 <- 0
s20 <- 1000
alpha.0 <- 0.01
beta.0 <- 0.01
# Initial values
Intcept <- 1
Slope <- 1
s2 <- 1
#----number of iterations
nit <- 5000
nb <- 2000
res <- mcmc.lr(phe,SNP,Intcept,Slope,s2,mu0,alpha.0,beta.0,s20,nit)
head(res)
data.frame(Mean=colMeans(res[nb:nit,]),SD=apply(res[nb:nit,],2,sd))
#---frequentist approach
lm(phe ~ SNP)
#---histograms
hist(res[,1],breaks = 100,main="Intercept")
abline(v=mean(res[,1]),col="blue")
hist(res[,2],breaks = 100,main="Slope")
abline(v=mean(res[,2]),col="blue")
hist(res[,3],breaks = 100,main="sigma2")
abline(v=mean(res[,3]),col="blue")
#---iterations
plot(res[,1]/(1:nit), type="l", main="Intercept", ylab="", xlab='')
plot(res[,2]/(1:nit), type="l", main="Slope", ylab="", xlab='')
plot(res[,3]/(1:nit), type="l", main="sigma2", ylab="", xlab='')
gBLUP vs BayesR 在GS使用
数据：
image.png
软件方法：
image.png
image.png
GCTB命令：
image.png
GCTA：
gBLUP
image.png
SNP-BLUP
image.png
image.png
使用PLINK 对验证群体（SNP—BLUP结果）：
image.png
结果：
image.png
使用PLINK 对验证群体（BayeR结果）：
image.png
image.png
比较在R
image.png
image.png
BayeR需要更长的计算时间，国际中国实际使用目前是ssGBLUP
网友评论