DMU-单性状动物模型-学习笔记2

单性状动物模型

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

数据使用learnasreml包中的数据

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

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

library(devtools)
install_github("dengfei2013/learnasreml")
library(learnasreml)
data("animalmodel.dat")
data("animalmodel.ped")

dat = animalmodel.dat
ped = animalmodel.ped

summary(dat)
summary(ped)

看一下数据:

> 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文件

使用R语言清洗数据, 并保存数据到D盘dmu-test

dir.create("d:/dmu-test")
setwd("d:/dmu-test/")
library(devtools)
install_github("dengfei2013/learnasreml")
library(learnasreml)
data("animalmodel.dat")
data("animalmodel.ped")
dat = animalmodel.dat
ped = animalmodel.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)

文件类型

数据文件:
系谱文件:

编写DIR文件

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

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

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

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

$ANALYSE 1 1 0 0

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

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

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

$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS

第五部分, 有6行, 定义模型
整体来说是:
第一行: 单性状
第二行: 无吸收
第三行: 第1个y变量, 0无权重考虑,3个因子,第3列是第一个固定因子, 第4列第二个固定因子, 第1列是随机因子
第四行:1个随机因子
第五行: 无回归项
第六行: 无约束

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

第六部分: 指定系谱

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

第七部分: 指定初始值(可以省略)

$PRIOR
1 1 1 0.3
2 1 1 0.7

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

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

$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 3 3 4 1
1
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt
$PRIOR
1 1 1 0.3
2 1 1 0.7

执行DIR文件

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

run_dmuai.bat Uni_animalmodel

执行结果:

D:\dmu-test>run_dmuai.bat Uni_animalmodel

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

查看结果

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

                             SUMMARY OF MINIMIZATION PROCESS
                                                            
 Eval Criterion    !!Delta!!   !!Gradient!!          Parameters
 ---- ---------     ---------  ------------ |----------------------------
   1   2656.56     0.6019        3.038      |    1.6342        1.5678    
  
   2   2279.44     0.5828        2.916      |    2.2850        2.0736    
  
   3   2194.38     0.2913        1.424      |    2.6342        2.2923    
  
   4   2186.51     0.4243E-01   0.2060      |    2.6859        2.3227    
  
   5   2186.39     0.9753E-03   0.3300E-02  |    2.6857        2.3241    
  
   6   2186.39     0.1209E-03   0.6814E-04  |    2.6858        2.3240    
  
   7   2186.39     0.1714E-04   0.9622E-05  |    2.6858        2.3240    
  
   8   2186.39     0.2431E-05   0.1365E-05  |    2.6858        2.3240    
  
   9   2186.39     0.3449E-06   0.1936E-06  |    2.6858        2.3240    

可以看到模型收敛

方差组分为:

                          Estimated (co)-variance components
                          ----------------------------------
                                                            
             Parameter vector for L at convergence          
             Asymptotic SE based on AI-information matrix   
                                                            
               No          Parameter             Asymp. S.E.
                                                            
                1           2.68577              0.443729
                2           2.32401              0.347584    
                                                            
                                                            
               Asymp. correlation matrix of parameter vector

遗传力为:

   Phenotypic co-variance matrix           
  
               1
     1        5.0097857
  
  
  
          Intra Class
   Trait  correlation     V(t)         SE(t)         SD-A          SD-P
  
     1     0.53611       0.00552       0.07431       1.63

可以看出, 遗传力为0.536, 标准误为0.07

对比asreml的结果:

代码:

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

方差组分:

> summary(mod)$varcomp
                   gamma component std.error  z.ratio constraint
ped(ANIMAL)!ped 1.155638  2.685531 0.4437949 6.051288   Positive
R!variance      1.000000  2.323852 0.3475778 6.685846   Positive

遗传力:

> pin(mod,h2 ~ V1/(V1+V2))
    Estimate         SE
h2 0.5361002 0.07432624

DMU和asreml比较

两者方差组分和遗传力一致.