||
更好的示范请见链接:https://blog.csdn.net/yzhlinscau/article/details/118444615
AFEchidna包可以批量进行计算基因组关系矩阵,简单示例如下:
library(AFEchidna) setwd("D:/Rdata") G.marker<-read.csv(file="G.marker.csv",header=T) G.ped<-read.csv(file="G.ped.csv",header=T) GmN<-paste0(c('GOF','GD','G05','GMF','Greg'),'.grm',sep='') Gmat<-vector("list", 5) for(i in 1:5){ cat('For matrix: ',i,'-----') Gmat[[i]]<-GenomicRel( G.marker,i, ped=G.ped, Gres=T) write.csv(Gmat[[i]],file=GmN[[i]],row.names=F,quote=F) }
通过file.edit('G.data.es0')修改.es0文件,添加上述5种基因组关系矩阵文件,示例如下:
// .es0 temple #!WORK 2 !REN !ARG TITLE: G.data #!DOPART $1 # "ID","Female","Male","Year","Site","t1","t2" ... # "26","1","12","2001","6",NA,NA ... ID !P # 26 Female !I # 1 Male !I # 12 Year !I # 2001 Site !I # 6 t1 # * t2 # * # Verify data fields are correctly classified ... G.ped.csv !SKIP 1 GOF.grm #G1 GD.grm #G2 G05.grm #G3 GMF.grm #G4 Greg.grm !ND #G5 G.data.csv !SKIP 1
运行普通的个体(动物)模型:
// animal model Ablup<-echidna(fixed=t1~1+Site, random=~ nrm(ID), es0.file='G.data.es0')
基于上述结果,进一步进行Gblup的批量分析:
// batch--Gblup Gblup.mG<-update(Ablup, random=c(G1~grm1(ID), G2~grm2(ID), G3~grm3(ID), G4~grm4(ID), G5~grm5(ID)), batch.G=T)
结果查看:
// results > Gblup.mG2<-b2s(Gblup.mG) > lapply(Gblup.mG2, Var) $G1 Term Sigma SE Z.ratio 1 Residual 181.870 3.3543 54.219956 2 grm1(ID) 98.327 13.0270 7.547939 $G2 Term Sigma SE Z.ratio 1 Residual 181.920 3.3557 54.212236 2 grm2(ID) 99.225 13.2310 7.499433 $G3 Term Sigma SE Z.ratio 1 Residual 181.87 3.3543 54.219956 2 grm3(ID) 214.87 28.4560 7.550956 $G4 Term Sigma SE Z.ratio 1 Residual 181.87 3.3543 54.219956 2 grm4(ID) 116.58 15.4430 7.549051 $G5 Term Sigma SE Z.ratio 1 Residual 181.870 3.3543 54.219956 2 grm5(ID) 73.167 9.6942 7.547503
上述结果显示,不同基因组关系矩阵,方差分量结果存在差异。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-21 01:11
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社