DMU-多性状动物模型-学习笔记4

最新推荐文章于 2023-11-24 15:18:21 发布

育种数据分析之放飞自我

最新推荐文章于 2023-11-24 15:18:21 发布

阅读量2.7k

点赞数

分类专栏：统计分析数量遗传学大数据时代的动植物育种数据分析

互秀朋友圈呀, 微信: yijiaobani, 备注: 姓名+ 单位

本文链接：https://blog.csdn.net/yijiaobani/article/details/83959539

版权

数量遗传学同时被 3 个专栏收录

90 篇文章 55 订阅

订阅专栏

大数据时代的动植物育种数据分析

33 篇文章 37 订阅

订阅专栏

统计分析

31 篇文章 2 订阅

订阅专栏

多性状动物模型

本次主要是演示如何使用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(第五列), TARSUS (第六列）
固定因子: BYEAR和SEX(第三列, 第四列)
随机因子: ID

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

$COMMENT

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

第二部分是分析方法, 默认是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 TARSUS

第五部分, 有6行, 定义模型
整体来说是:
第一行: 两性状 # 2
第二行: 1性状无吸收 # 0
第三行: 2性状无吸收 # 0
第四行: 1性状, 是由3个因子, 两个固定因子:3,4, 一个随机因子:1 # 1 0 3 3 4 1
第五行: 2性状, 是由3个因子, 两个固定因子:3,4, 一个随机因子:1 # 2 0 3 3 4 1
第六行: 性状1, 1个随机因子 # 1 1
第七行: 性状2, 1个随机因子 # 1 1
第八行: 性状1,无回归 # 0
第九行: 性状2,无回归 # 0
第十行: 残差0

$MODEL
2  
0
0
1 0 3 3 4 1
2 0 3 3 4 1
1 1
1 1
0
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中, 然后重命名为: mul-animalmodel.DIR

$COMMENT
Model
y: BWT  TARSUS 
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
2  
0
0
1 0 3 3 4 1
2 0 3 3 4 1
1 1
1 1
0
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt
$DMUAI
10
1D-7
1D-6
1

执行DIR文件

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

run_dmuai.bat mul_animalmodel

执行结果:

D:\dmu-test>run_dmuai.bat mul_animalmodel

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

D:\dmu-test>

查看结果

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

Eval Criterion    !!Delta!!   !!Gradient!!                        Parameters
 ---- ---------     ---------  ------------ |--------------------------------------------------------
   1   12028.8     0.6039        4.096      |    1.5877       0.73966E-01    1.8936        1.4327    
                                            |   0.12929        1.9136    
  
   2   7774.73     0.9673        6.170      |    2.1162       0.31777        3.4204        1.7356    
                                            |   0.49187        3.5631    
  
   3   5909.74      1.510        8.930      |    2.3621       0.82080        5.6988        1.9228    
                                            |    1.1827        6.3310    
  
   4   5161.67      1.984        10.91      |    2.4806        1.4486        8.3095        2.1217    
                                            |    2.1515        10.257    
  
   5   4917.50      1.785        9.047      |    2.5638        1.8830        10.081        2.3066    
                                            |    3.0591        14.120    
  
   6   4867.84     0.7835        3.603      |    2.5821        1.9975        10.541        2.3927    
                                            |    3.4651        15.932    
  
   7   4864.20     0.8472E-01   0.3898      |    2.5817        2.0041        10.586        2.4033    
                                            |    3.5105        16.129    
  
   8   4864.17     0.1682E-02   0.4107E-02  |    2.5819        2.0049        10.590        2.4033    
                                            |    3.5102        16.128    
  
   9   4864.17     0.4621E-04   0.6606E-04  |    2.5819        2.0049        10.590        2.4032    
                                            |    3.5102        16.128    
  
  10   4864.17     0.7679E-05   0.1041E-04  |    2.5819        2.0049        10.590        2.4032    
                                            |    3.5102        16.128    
  
  11   4864.17     0.1192E-05   0.1748E-05  |    2.5819        2.0049        10.590        2.4032    
                                            |    3.5102        16.128    
  
  12   4864.17     0.1937E-06   0.3123E-06  |    2.5819        2.0049        10.590        2.4032    
                                            |    3.5102        16.128

可以看到模型收敛

方差组分为:

                          Estimated (co)-variance components
                          ----------------------------------
                                                            
             Parameter vector for L at convergence          
             Asymptotic SE based on AI-information matrix   
                                                            
               No          Parameter             Asymp. S.E.
                                                            
                1           2.58189              0.437110    
                2           2.00491              0.857216    
                3           10.5895               2.68116    
                4           2.40324              0.348455    
                5           3.51022              0.727723    
                6           16.1280               2.36436

遗传力需要手动计算, 这里还没有找到解决方案.

对比asreml的结果:

代码:

library(asreml)
dat[dat$BWT==0,]$BWT=NA
dat[dat$TARSUS==0,]$TARSUS=NA
ainv = asreml.Ainverse(ped)$ginv
mod2 = asreml(cbind(BWT,TARSUS) ~ trait + trait:(BYEAR + SEX),
             random = ~ us(trait):ped(ANIMAL), rcov = ~ units:us(trait),ginverse = list(ANIMAL=ainv),data=dat)
summary(mod2)$varcomp

方差组分:

> summary(mod2)$varcomp
                                          gamma component std.error  z.ratio constraint
trait:ped(ANIMAL)!trait.BWT:BWT        2.581883  2.581883 0.4371085 5.906732   Positive
trait:ped(ANIMAL)!trait.TARSUS:BWT     2.004949  2.004949 0.8572152 2.338910   Positive
trait:ped(ANIMAL)!trait.TARSUS:TARSUS 10.589430 10.589430 2.6811944 3.949520   Positive
R!variance                             1.000000  1.000000        NA       NA      Fixed
R!trait.BWT:BWT                        2.403246  2.403246 0.3484542 6.896879   Positive
R!trait.TARSUS:BWT                     3.510189  3.510189 0.7277219 4.823531   Positive
R!trait.TARSUS:TARSUS                 16.128117 16.128117 2.3643446 6.821390   Positive