引子： 本人于2019年主译出版的《动植物育种遗传数据分析》，原书主要使用ASReml进行分析。本博客将陆续演示使用AFEchidna包达到类似ASReml的结果。
第4章的代码如下：
library(AFEchidna)
### chapter 4 Breeding Values
m1<-echidna(height~1+Prov,
random=~Female*Block,
es0.file="pine_provenance.es0")
Var(m1)
pin(m1,mulp=c(Va~4*V3,
VP~V1+V3+V4,
H2i~4*V3/(V1+V3+V4)))
## Variation Among Family Means
pin(m1,mulp=c(VP.fm~V1/(5*5.2)+V3+V4/5,
H2.fs~V3/(V1/(5*5.2)+V3+V4/5)))
## Within-Family Variation
pin(m1,mulp=c(VP.wf~V1*(5*5.2-1)/(5*5.2)+V4*(5-1)/5,
H2.wf~3*V3/(V1*(5*5.2-1)/(5*5.2)+V4*(5-1)/5)))
## 4.2 calculate BV and its Accuracy
GCA <- coef(m1)$random %>% filter(Term=='Female')
names(GCA)[3]="GCA"
head(GCA)
GCA$BV=2*GCA$GCA+coef(m1)$fixed[5,3] # BV
##GCA variance
GCA.var=Var(m1)[3,'Sigma']
GCA.se=GCA$SE
## Accurancy of BV
Corr=sqrt(1-GCA.se^2/GCA.var)
GCA$Corr=Corr
head(GCA)
## 4.3 tree model
tm1<-echidna(height~1+Prov,
random=~nrm(Treeid)+Block,
residual=~idv(units),
es0.file="pine_provenance.es0")
Var(tm1)
bv<-coef(tm1)$random
head(bv)
tbv<- bv %>% filter(Term=='nrm(Treeid)')
mu=coef(tm1)$fixed[5,3]
BV.se=tbv$SE
BV.var=Var(tm1)[3,'Sigma']
Corr=sqrt(1-BV.se^2/BV.var)
tbv$Corr=Corr
head(tbv)
## 4.4 animal model for pig dataset
pm1<-echidna(weanwt~year+sex+weanage,
random=~nrm(pig),
es0.file='pig_data.es0')
Var(pm1)
pin(pm1)
pin(pm1,mulp=c(Va~V2,
VP~V1+V2,
h2i~V2/(V2+V1)) )
pm2<-echidna(weanwt~year+sex+weanage,
random=~str(~nrm(pig)+nrm(dam),~us(2):nrm(pig)),
#delf=F,foldN='pm2',
es0.file='pig_data.es0')
Var(pm2)
pin(pm2,mulp=c(VP~V1+V2+V3+V4,
h2i.d~V2/(V1+V2+V3+V4), # Direct heritability
h2i.m~V4/(V1+V2+V3+V4), # Maternal heritability
dm.cor~V3/sqrt(V2*V4)
))
pm3<-update(pm2,maxit=35,
random=~str(~nrm(pig)+nrm(dam),~corh(2):nrm(pig))+ide(dam))
Var(pm3)
summary(pm3)
waldT(pm3)
### Genetic Groups Effect
# method 1: same to code 4.1
# method 2
#file.edit('pine_provenance_pedigree.txt')
#file.copy('pine_provenance.es0','pine_provenance2.es0')
#file.edit('pine_provenance2.es0')
ggm<-echidna(height~1,random=~nrm(Treeid)+Block,
es0.file='pine_provenance2.es0')
Var(ggm)
tbv<-coef(ggm)$random %>% filter(Term=='nrm(Treeid)')
head(tbv)
##### 4.5 Effect of Self-Fertilization
# file.edit('pine_provenance.es0')
sfm<-echidna(height~1+Prov,
random=~ nrm(Treeid)+Block,
es0.file='pine_provenance.es0')
Var(sfm)
sfm.s1<-update(sfm,selfing=0.1)
Var(sfm.s1)
参考文献
1 . Zhang WH, Wei RY, Liu Y, Lin YZ. AFEchidna is a R package for genetic evaluation of plant and animal breeding datasets. BioRxiv. DOI: 10.1101/2021.06.24.449740.
2. 林元震 丁昌俊 主译，动植物育种遗传数据分析, 科学出版社，2019.9
