博文

自编程序包AAfun的一些功能（之一）

已有 3895 次阅读 2014-5-14 10:40 |个人分类:ASReml|系统分类:科研笔记

本人在ASReml-R程序包的使用中，编写了一个程序包AAfun，增加一些额外的函数，作为ASReml-R的辅助。主要功能如下：

功能一、asreml.pin，用于计算遗传力、遗传相关等参数，并得到遗传相关的显著性水平。

library(AAfun)
 
### function 1 asreml.pin ####
## pin(): 1. count standard error(se) for heritability(h2) and corr value;
##        2. output significent level for corr value
## work for asreml package
 
data(PrSpa)
df<-PrSpa
 
library(asreml)
# exmaple 1 for sigle trait model
fm<-asreml(h5~1+Rep, random=~Fam, subset=Spacing=='3',data=df)
summary(fm)$varcomp[,1:3]
pin(fm, h2 ~4*V1/(V1+V2)) # asreml result, h2 formlua
 
# exmaple 2 for us model for bi-trait
fm2<-asreml(cbind(dj,h5)~ trait+trait:Rep,
            random=~ us(trait):Fam, rcov=~units:us(trait),
            subset=Spacing=='3',data=df,maxit=40)
 
summary(fm2)$varcomp[,1:3]
 
pin(fm2, h2_A ~ 4 * V1/(V1+V5)) # heritability for trait A
pin(fm2, h2_B ~ 4 * V3/(V3+V7)) # heritability for trait B
 
# genetic corr: asreml result, corr formlua, 1 for signif
pin(fm2, gCORR ~ V2/sqrt(V1*V3),1)
pin(fm2, pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)),1) # phenotype corr
 
# exmaple 3 for corr model for multi-trait
fm3<-asreml(cbind(dj,h3,h5)~ trait+trait:Rep,
            random=~ corgh(trait):Fam, rcov=~units:us(trait),
            subset=Spacing=='3',data=df,maxit=40)
 
summary(fm3)$varcomp[,1:3]
pin(fm3,N=3) #  asreml result, corr coef number
pin(fm3) # corr coef default N is 1, for the first corr

运行结果如下：

> pin(fm, h2 ~4*V1/(V1+V2)) # asreml result, h2 formlua
 
   Estimate    SE
h2    0.306 0.124
 
> pin(fm2, h2_A ~ 4 * V1/(V1+V5)) # heritability for trait A
 
     Estimate    SE
h2_A    0.455 0.145
 
> pin(fm2, gCORR ~ V2/sqrt(V1*V3),1)
 
      Estimate    SE sig.level
gCORR    0.442 0.257         *
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
 
> pin(fm3,N=3)
                                 Estimate    SE sig.level
trait:Fam!trait.h3:!trait.dj.cor    0.751 0.233       ***
trait:Fam!trait.h5:!trait.dj.cor    0.448 0.257         *
trait:Fam!trait.h5:!trait.h3.cor    0.798 0.123       ***
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

功能二、asreml.met，用于辅助ASReml在多地点试验(MET)的因子分析。

### function 2 asreml.met ####
##  met.plot(): plots MET data
##  met.corr(): calculate var/cov/corr from asreml MET factor analytic results
##  met.biplot(): biplots MET factor analytic results from asreml
##  work for asreml package
 
data(MET)
 
##  plot MET data -- example 1
# variable order: genotype,yield,site,row,col
MET2<-MET[,c(1,9,2,4:5)]
met.plot(MET2) 
 
## plot MET data -- example 2
MET3<-MET[,c(1,9,2,4:7)] # add variable order on MET2: Rep, Block
met.plot(MET3,"My met trials")
 
MET$yield <- 0.01* MET$yield
met.asr<-asreml(yield~Loc, random=~ Genotype:fa(Loc,2),
                rcov=~ at(Loc):ar1(Col):ar1(Row),
                data=MET, maxiter=50)
 
## count var/cov/corr matrix, etc.
met.corr(met.asr,MET$Loc,2)
met.corr(met.asr,MET$Loc,kn=2) # asreml met-result, Site, Site cluster number
 
## biplot asreml-met results
met.biplot(met.asr,6,36,2)
met.biplot(met2.asr,6,36,3)
 
met.biplot(met.asr,siteN=6,VarietyN=36,faN=2,dSco.u=1.8,dLam.u=1.0)
# dLam.u -- least distance from center
# dSco.u -- least Variety breeding value
# if can not draw fig 3, try multiplying or being devided by 10 for aim trait data.

运行结果如下：

1. 绘制MET数据的图形

2. 输出MET分析的相关矩阵

> met.corr(met.asr,MET$Loc,2)
Site cluster results:
S1 S2 S3 S4 S5 S6
 1  1  1  2  1  2
 
CovVarCorr matrix
 
     S1    S2     S3    S4    S5     S6
S1 5.38 0.816  0.875 0.520 0.979  0.130
S2 3.90 4.246  0.688 0.472 0.811  0.161
S3 4.97 3.469  5.988 0.042 0.759 -0.366
S4 2.84 2.288  0.243 5.535 0.682  0.914
S5 4.31 3.169  3.525 3.044 3.598  0.328
S6 0.32 0.353 -0.954 2.291 0.662  1.135

3. 绘制MET分析结果的biplot图

> met.biplot(met.asr,6,36,2)
 
Scores.u is: 4
       Fa1   Fa2 Scores
V29 -1.164 -5.99   6.11
V30  0.470 -4.49   4.51
V31 -0.739 -5.02   5.08
V32  0.541 -4.03   4.06
V33 -0.218 -4.96   4.97
 
distFC.u is: 1.2
     FA1   FA2 distFC
S3 -2.03  1.37   1.68
S4 -1.39 -1.90   1.65
S6 -0.23 -1.04   1.57