单性状动物模型-母体效应

本次主要是演示如何使用DMU分析单性状动物模型-母体效应.

数据使用learnasreml包中的数据

learnasreml是我编写的辅助学习asreml的R包, 里面有相关的数据和代码, 这里我们用其中的animalmodel.dat和animalmodel.ped的数据.

如果没有软件包, 首先安装:

setwd("d:/dmu-test/")
library(devtools)
# install_github("dengfei2013/learnasreml")
library(learnasreml)
data("animalmodel.dat")
data("animalmodel.ped")

dat = animalmodel.dat
ped = animalmodel.ped

summary(dat)
summary(ped)

dmuped = ped
dmuped$Birth = 2018

head(dat)
library(data.table)
# write.table(dat,"animal-model.txt",row.names = F,col.names = F)
fwrite(dat,"animal-model.txt",sep = " ",col.names = F)
fwrite(dmuped,"animal-ped.txt",sep = " ",col.names = F)

看一下数据:

> summary(dat)
     ANIMAL         MOTHER         BYEAR     SEX          BWT             TARSUS     
 1      :   1   96     :   8   998    : 53   1:470   Min.   : 0.000   Min.   : 0.00  
 2      :   1   541    :   8   994    : 47   2:614   1st Qu.: 2.730   1st Qu.: 0.00  
 3      :   1   581    :   8   983    : 45           Median : 6.385   Median :16.27  
 5      :   1   584    :   8   987    : 45           Mean   : 5.802   Mean   :12.93  
 6      :   1   1302   :   8   991    : 45           3rd Qu.: 8.660   3rd Qu.:21.94  
 7      :   1   12     :   7   997    : 44           Max.   :15.350   Max.   :39.66  
 (Other):1078   (Other):1037   (Other):805                                           
> summary(ped)
       ID           FATHER           MOTHER      
 Min.   :   1   Min.   :   0.0   Min.   :   0.0  
 1st Qu.: 328   1st Qu.:   0.0   1st Qu.: 135.0  
 Median : 655   Median :   0.0   Median : 538.0  
 Mean   : 655   Mean   : 261.5   Mean   : 547.4  
 3rd Qu.: 982   3rd Qu.: 458.0   3rd Qu.: 932.0  
 Max.   :1309   Max.   :1304.0   Max.   :1306.0

数据中,
有因子4个: 分别是ANIMAL, MOTHER, BYEAR, SEX
有变量2个: 分别是BWT和TARSUS
缺失值为0

系谱中,
有三列数据, 无出生时间一列, 缺失值为0

需要做的处理

• 系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化
• 在保存数据时, 去掉行头
• 编辑DIR文件

编写DIR文件

想要分析的模型:
观测值: BWT(第五列)
固定因子: BYEAR和SEX(第三列, 第四列)
随机因子: ID + MOTHER

所以这里编写DIR
第一部分, 是注释, 这里所写的东西会输出到结果文件, 基本上就是模型的解释, 这部分没有强制要求, 可以省略

$COMMENT

Model
y: BWT  TARSUS 
fixed: BYEAR + SEX
random: ANIMAL+MOTHER

第二部分是分析方法, 默认是AI

$ANALYSE 1 1 0 0

第三部分是定义因子数和变量数, 以及文件位置:

$DATA  ASCII (4,2,0) d:/dmu-test/animal-model.txt

上面的意思是, 数据是ASCII格式, 有4个固子, 2个变量, 缺失值用0表示, 数据的绝对路径是: d:/dmu-test/animal-model.txt

第四部分, 定义变量名, 也是为了方便结果输出, 相当于数据的行头名

$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT

第五部分, 有6行, 定义模型
整体来说是:
第一行: 单性状 # 1
第二行: 1性状无吸收 # 0
第三行: 1个性状, 是由3个因子, 两个固定因子:3,4, 一个随机因子:1 # 1 0 3 3 4 1
第四行: 2个随机因子, 他们的关系是独立, 没有协方差, 1 2 # 2 1 2

注意, 如果两个随机因子之间, 是有协方差的, 那么写作 2 1 1, 即表示随机因子有: 加性+ 母体 + 加性:母体协方差

第五行: 无回归项 # 0
第六行: 无约束 # 0

$MODEL
1
0
1 0 4 3 4 1 2
2 1 2
0 
0

第六部分: 指定系谱

$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt

注意, 如果想要输出BLUP值, 定义:$DMUAI

$DMUAI
10
1D-7
1D-6
1

完整DIR文件, 放入model.txt中, 然后重命名为: Uni_animalmodel-maternal.DIR

$COMMENT
Estimate variance components for BWT
Model
y: BWT
fixed: BYEAR + SEX
random: ANIMAL + MOTHER

$ANALYSE 1 1 0 0
$DATA  ASCII (4,2,0) d:/dmu-test/animal-model.txt
$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS
$MODEL
1
0
1 0 4 3 4 1 2
2 1 2
0 
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt

$SOLUTION
$DMUAI
10
1D-7
1D-6
1

执行DIR文件

这里运行的run_dmuai.bat, 将DMU安装路径下的文件run_dmuai.bat拷贝到d:/dmu-test文件夹, 在终端cmd界面键入:

run_dmuai.bat Uni_animalmodel-maternal

执行结果:

D:\dmu-test>run_dmuai.bat Uni_animalmodel-maternal

D:\dmu-test>Echo OFF
Starting DMU using Uni_animalmodel-maternal.DIR as directive file

D:\dmu-test>

查看结果

在文件*lst中有估算的方差组分, 结果如下:

                            SUMMARY OF MINIMIZATION PROCESS

 Eval Criterion    !!Delta!!   !!Gradient!!                   Parameters
 ---- ---------     ---------  ------------ |------------------------------------------
   1   2311.63     0.4281        3.796      |    1.6003        1.1819        1.3955    

   2   2170.63     0.3174        2.343      |    2.1014        1.1469        1.6188    

   3   2150.86     0.1004       0.5529      |    2.2655        1.1053        1.6588    

   4   2150.13     0.7153E-02   0.2589E-01  |    2.2777        1.1040        1.6570    

   5   2150.13     0.8765E-04   0.8250E-04  |    2.2778        1.1040        1.6569    

   6   2150.13     0.8343E-06   0.1020E-05  |    2.2778        1.1040        1.6569    

   7   2150.13     0.5385E-07   0.4235E-07  |    2.2778        1.1040        1.6569

可以看到模型收敛

方差组分为:

                          Estimated (co)-variance components
                          ----------------------------------

             Parameter vector for L at convergence          
             Asymptotic SE based on AI-information matrix   

               No          Parameter             Asymp. S.E.

                1           2.27778              0.497101    
                2           1.10404              0.239802    
                3           1.65690              0.373448    


               Asymp. correlation matrix of parameter vector

可以看到:
加性方差组分为: 2.2778
母体效应方差组分为: 1.10404
残差方差组分为: 1.6569

对比asreml的结果:

代码:

library(asreml)
head(dat)
dat[dat$BWT==0,]$BWT=NA
ainv = asreml.Ainverse(ped)$ginv
mod = asreml(BWT ~ BYEAR + SEX, random = ~ ped(ANIMAL) + MOTHER, ginverse = list(ANIMAL=ainv),data=dat)
summary(mod)$varcomp
pin(mod,h2 ~ V2/(V1+V2+V3))

方差组分:

> summary(mod)$varcomp
                     gamma component std.error  z.ratio constraint
MOTHER!MOTHER.var 0.666325  1.104035 0.2398025 4.603936   Positive
ped(ANIMAL)!ped   1.374724  2.277783 0.4971012 4.582131   Positive
R!variance        1.000000  1.656902 0.3734477 4.436772   Positive

遗传力:

> pin(mod,h2 ~ V2/(V1+V2+V3))
    Estimate         SE
h2 0.4520558 0.08842455

DMU和asreml比较

两者方差组分一致.