AFEchidna示例11--批量进行计算基因组关系矩阵和运行Gblup分析

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)
 }

运行过程：

>  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)
+  }
For matrix:  1 -----Generating G matrix is under going.
For matrix:  2 -----Generating G matrix is under going.
For matrix:  3 -----Generating G matrix is under going.
For matrix:  4 -----Generating G matrix is under going.
For matrix:  5 -----Generating G matrix is under going.

通过上述代码生成5种基因组关系矩阵，并批量输出为相应文件。

生成.es0文件，如下：

 get.es0.file(dat.file='G.data.csv')
 get.es0.file(es.file='G.data.es')

通过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)

运行过程：

// process
>  Gblup.mG<-update(Ablup, random=c(G1~grm1(ID),
+                                    G2~grm2(ID),
+                                    G3~grm3(ID),
+                                    G4~grm4(ID),
+                                    G5~grm5(ID)),
+                    batch.G=T)

Program runs for 2 more G-structure at one time. ------ 

run 1 -- random effects: grm1(ID)
Running Echidna for analysis:  t1 

 Sat Jul  3 20:32:24 2021 
  Iteration      LogL eSigma NEDF
1         1 -19010.39  193.5 6028
2         2 -18972.54  189.9 6028
3         3 -18934.72  185.5 6028
4         4 -18921.93  183.2 6028
5         5 -18919.44  182.1 6028
6         6 -18919.32  181.9 6028
7         7 -18919.32  181.9 6028
Sat Jul  3 20:32:25 2021 LogL Converged 


run 2 -- random effects: grm2(ID)
Running Echidna for analysis:  t1 

 Sat Jul  3 20:32:25 2021 
  Iteration      LogL eSigma NEDF
1         1 -19008.43  193.4 6028
2         2 -18972.15  189.9 6028
3         3 -18935.40  185.6 6028
4         4 -18922.76  183.3 6028
5         5 -18920.23  182.2 6028
6         6 -18920.11  181.9 6028
7         7 -18920.11  181.9 6028
Sat Jul  3 20:32:26 2021 LogL Converged 


run 3 -- random effects: grm3(ID)
Running Echidna for analysis:  t1 

 Sat Jul  3 20:32:27 2021 
  Iteration      LogL eSigma NEDF
1         1 -19107.82  201.9 6028
2         2 -19045.64  196.6 6028
3         3 -18999.52  192.5 6028
4         4 -18946.33  187.0 6028
5         5 -18925.25  184.0 6028
6         6 -18919.76  182.4 6028
7         7 -18919.33  181.9 6028
8         8 -18919.32  181.9 6028
Sat Jul  3 20:32:28 2021 LogL Converged 


run 4 -- random effects: grm4(ID)
Running Echidna for analysis:  t1 

 Sat Jul  3 20:32:28 2021 
  Iteration      LogL eSigma NEDF
1         1 -19029.56  195.2 6028
2         2 -18985.68  191.2 6028
3         3 -18940.18  186.2 6028
4         4 -18923.40  183.6 6028
5         5 -18919.56  182.2 6028
6         6 -18919.33  181.9 6028
7         7 -18919.32  181.9 6028
Sat Jul  3 20:32:29 2021 LogL Converged 


run 5 -- random effects: grm5(ID)
Running Echidna for analysis:  t1 

 Sat Jul  3 20:32:30 2021 
  Iteration      LogL eSigma NEDF
1         1 -18980.91  190.7 6028
2         2 -18953.28  187.8 6028
3         3 -18927.54  184.4 6028
4         4 -18920.07  182.5 6028
5         5 -18919.34  182.0 6028
6         6 -18919.32  181.9 6028
7         7 -18919.32  181.9 6028
Sat Jul  3 20:32:31 2021 LogL Converged

结果查看：

// 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

上述结果显示，不同基因组关系矩阵，方差分量结果存在差异。

AFEchidna包免费对外开放，仅用于学术研究，不可用于商业研究，否则造成的后果自负。

感兴趣者，可发邮件到yzhlinscau@163.com免费索取程序包。

参考文献：
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.