DMU遗传评估从入门到出家系列

邓飞

发布于 2020-05-17 09:00:29

3.8K1

发布于 2020-05-17 09:00:29

文章被收录于专栏：育种数据分析之放飞自我

这是之前写的DMU学习笔记，做一个汇总，更正了一些错误，添加了Linux版的安装包。

DMU遗传评估从入门到出家系列

1. 主要参数介绍

1.1 DMU软件介绍

DMU是一个数量遗传学工具包，主要功能包括估计方差组分和固定效应，预测育种值。DMU的开发历史可以追溯到25年前，大部分功能基于数量遗传学研究的需求而开发。在丹麦动物育种研究中，DMU是一个主要的统计研究工具（估计和预测）。此外，DMU也应用于丹麦牛，羊，貂和马等常规遗传评估研究。因此，DMU不但在一些特定的项目中具备高性能优势，也适用于常规数量遗传学研究。“DMU”名称最初来自于程序包中用来进行初始化的过程名字缩写。这些过程利用约束最大似然法（REML），通过Derivative-free方式执行MUltivariate analysis，因此得名DMU。但是，在当前的DMU版本中，并不包括DF-REML模块，现在D仅代表DJF（丹麦农业科学学院的缩写）。

DMU安装包有很多模块，如DMU1、DMU4、DMU5、DMUAI和RJMC。DMUAI模块可利用平均信息限制最大似然(AI-REML) (Jensen et al. 1997)算法进行（协）方差组分的估计。AI是通过平均观察和预期信息的信息矩阵得到的。该模块还可以使用期望最大化(EM)算法来最大化约束似然函数。被估计的（协）方差组分的渐近标准误是从平均信息矩阵中获得的。

1.2 学习DMU初衷

想试试DMU处理一批数据, 发现这个软件, 竟然没有一个合适的操作说明文档, 我手头上有苏国生老师的PPT中文版DMU操作说明, 但看起来还是费劲.

刚好自己在学习这个软件, 用实际数据来演示如何使用这个软件进行数据分析.

「我想从四部分进行:」

1, DMU语法介绍
2, 单性状动物模型
3, 单性状重复力模型
4, 多性状动物模型

其它内容, 包括测定日模型(随机回归模型), 母体效应模型, GBLUP模型, 显性上位性模型, 一步法GS模型等等以后再做总结.

说明文档是作者写的, 一般来说作者都想通过逻辑的构建, 让读者了解软件的方方面面, 但是读者一开始接触软件时, 迫切的是想解决问题, 不是来学理论, 不是来学知识, 只是想解决问题. 但是大多数文档无法满足这些迫切的需求. 所以, 最好的操作说明, 就是有数据, 有模型, 有结果说明, 可以很快上手. 我写此操作说明的目的就在于此.

1.3 DMU语法介绍

「软件组成, 主要包括四类程序」

DMU1 这个主要是为了整理数据和模型, 相当于预处理程序, 其它三个程序都要经过它的处理才能分析. 类似BLUPF90的renumf90程序.
DMUAI 这个主要估算方差组分的程序
DMU4和DMU5 DMU4主要是求解混合线性方程组, 它不估算方差组分, 只求解. 类似BLUPF90包中的blupf90程序.DMU5功能和DMU4类似, 也是求解方程组, 适用于大数据
RGMC 主要是贝叶斯抽样, 估算方差组分, 计算育种值.

「数据和系谱及逆矩阵格式」

全部数据, 不要有行头
数据中不能含有字符, 字母, 都必须是数字
逆矩阵可以是下三角或者上三角矩阵的三列形式
系谱数据包括四列: ID, Sire, Dam, Birth
数据中, 因子(ID, Sex...)放在前面, 观测值(y1, y2, y3)放在后面, 因子用整数表示, 不能含有字母

因此, 在进行分析之前, 首先需要对数据进行转化, 比如系谱要变为整数, 要有第四列信息出生信息, 如果没有, 就写成2018年就行. 数据中也要重新编号, 特别是某些因子含有字母, 需要转化为数字. 可以使用R语言进行转化, 将系谱的所有水平编号为1...n, 然后替换. 将数据的所有水平, 重新编码.

「参数文件」文件名为name.DIR, 其中name为程序名称, DIR必须要有, 并保持大写.

$COMMENT文件注释, 一般是解释你所使用的模型
$ANALYSIS你分析所使用的模型, 如果你需要估算方差组分, 那么简单写为:$ANALYSIS 1 1 0 0
$DATA指定数据格式,因子数目, 观测值数目, 缺失值, 和数据位置如果是txt文件, 有5个因子, 4个观测值, 缺失值-999, 在D盘根目录$DATA ASCII(5,4,-999) d:/dat.txt
$VARIABLE写出因子和变量的名称, 第一行为因子, 第二行为变量

ID Loc Year Herd Sex Hy
y1 y2 y3 y4

$MODEL指定分析模型中, 观测值个数, 固定因子, 随机因子比如单性状, 正态数据

1 1 0 0 0

比如二性状, 正态数据

2 2 0 0 0

固定因子: 每个性状一行, 包含若干整数单性状中, y1 = Loc + Year + Herd + Sex, random = ID

1 0 5 1 2 3 4 5

随机因子: 每个性状一行, 包含若干整数 1

$VAR_STR定义方差协方差结构可以支持系谱, 和自定义关系矩阵inv 定义系谱文件:$VAR_STR 2 PED 2 ASCII ped.txt定义逆矩阵:$VAR_STR 1 COR ASCII ginv
$PRIOR定义初始值, 不过不定义, 默认是方差组分为1, 协方差组分为0, 定义格式, 下三角行列形式. 比如两性状, Vg和Ve

$VAR_REST(可选项, 主要是固定初始值)

1.4 文件输出

lst描述统计, 模型迭代, 方差组分估计
PAROUT方差组分估计(行列形式显示)
PAROUT-STD方差组分及标准误(计算遗传力)
LLIK最后一次迭代情况

1.5 命令行文件执行

run_dmuai运行dmuai程序
run_dmu4运行dmu4程序
run_dmu5
run_rjmc

2. 单性状动物模型

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

2.1 示例数据

「数据使用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

2.2 数据预处理

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

2.3 使用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)

2.4 文件类型

「数据文件:」

「系谱文件:」

2.5 编写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

2.6 执行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

2.7 查看结果

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

2.8 对比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

2.9 DMU和asreml比较

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

3. 单性状重复力模型

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

3.1 概念解析

「重复力模型和动物模型的区别:」不是所有的性状都可以分析重复力模型, 首先重复力模型是动物模型的拓展, 它适合一个个体多个观测值的情况.

比如猪的产仔数, 一个母猪有多个胎次
比如鸡的产蛋, 不同时间段, 鸡都有产蛋量
牛的产奶量, 不同的测定日, 产奶量不同
猪的饲料消耗, 也是重复测量的数据

只有这样的数据才可以将永久环境效应剖分出来.

「重复力是遗传力的上限:」教科书上这样说, 这句话怎么理解呢? 首先, 我们认为

遗传力为

这里的Vg如果有重复测量的数据, 可以剖分为可以遗传的部分, 和不可以遗传的部分(永久环境效应), 那么遗传力的计算公式为:

重复力的计算公式为:

当Vpe为0时, 重复力=遗传力, 当Vpe>0时, 重复力>遗传力, 所以说重复力是遗传力的上限!

3.2 使用数据

「数据使用learnasreml包中的数据」

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

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

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

dat = repeatmodel.dat
ped = repeatmodel.ped

summary(dat)
summary(ped)

看一下数据:

> summary(dat)
     ANIMAL         BYEAR      AGE          YEAR         LAYDATE
 1      :   5   1000   : 109   2:308   1004   :  79   Min.   : 0.00
 3      :   5   1001   :  98   3:322   1005   :  78   1st Qu.:20.00
 9      :   5   999    :  86   4:339   1003   :  69   Median :24.00
 17     :   5   1002   :  85   5:315   1006   :  64   Mean   :23.54
 42     :   5   987    :  70   6:323   1002   :  60   3rd Qu.:27.00
 50     :   5   989    :  66           988    :  54   Max.   :41.00
 (Other):1577   (Other):1093           (Other):1203
> 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 BYEAR AGE YEAR
有变量1个: LAYDATE
缺失值用0表示
系谱中有三列数据, 无出生时间一列, 缺失值为0

3.3 需要做的处理

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

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

dat = repeatmodel.dat
ped = repeatmodel.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,"repeat-model.txt",sep = " ",col.names = F)
fwrite(dmuped,"repeat-ped.txt",sep = " ",col.names = F)

3.5 编写DIR文件

想要分析的模型: 观测值: LAYDATE (第四列) 固定因子: 第二列, 第三列, 第四列随机因子: ID, ide(ID)

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

$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)

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

$ANALYSE 1 1 0 0

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

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

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

$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT

第五部分, 有6行, 定义模型整体来说是: 第一行: 单性状 # 1 第二行: 无吸收 # 0 第三行: 主要定义y变量, 固定因子, 随机因子

分析的是第一个变量 # 1
无权重考虑 # 0
共五个因子(固定+随机, 固定写前面, 随机写后面) # 5
第一个固定因子是第二列因子 #2 #BYEAR
第二个固定一致是第三列因子 #3 #AGE
第三个固定因子是第四列 #4 #YEAR
第四个随机因子是第一列 #1 #ANIMAL
第五个随机因子是第一列 #1 #ANIMAL
所以, 5个因子, 三个固定因子:2,3,4, 两个随机因子:1,1 #1 0 5 2 3 4 1 1

第四行: 有两个随机因子, 他们的关系是独立的, 所以是2 1 2

1
0
1 0 5 2 3 4 1 1
2 1 2
0
0

第六部分: 指定系谱

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

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

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

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

$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)

$ANALYSE 1 1 0 0
$DATA  ASCII (4,1,0) d:/dmu-test/repeat-model.txt
$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT
$MODEL
1
0
1 0 5 2 3 4 1 1
2 1 2
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt


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

3.6 执行DIR文件

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

run_dmuai.bat Uni_repeatmodel

执行结果:

D:\dmu-test>run_dmuai.bat Uni_repeatmodel

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

3.7 查看结果

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

                             SUMMARY OF MINIMIZATION PROCESS
                                                            
 Eval Criterion    !!Delta!!   !!Gradient!!                   Parameters
 ---- ---------     ---------  ------------ |------------------------------------------
   1   12629.2     0.8574        4.330      |    1.8100        1.8898        1.8705
  
   2   8234.59      1.370        6.822      |    2.9917        3.3812        3.2879
  
   3   6444.28      1.776        8.529      |    4.2397        5.4642        5.1761
  
   4   5857.47      1.566        6.869      |    4.9013        7.4662        6.8832
  
   5   5736.36     0.6798        2.497      |    4.9407        8.3737        7.6324
  
   6   5727.01     0.7387E-01   0.2311      |    4.9325        8.4634        7.7233
  
   7   5726.91     0.1399E-02   0.1596E-02  |    4.9341        8.4621        7.7245
  
   8   5726.91     0.7706E-04   0.5119E-04  |    4.9340        8.4622        7.7245
  
   9   5726.91     0.4564E-05   0.2610E-05  |    4.9340        8.4622        7.7245
  
  10   5726.91     0.2695E-06   0.1558E-06  |    4.9340        8.4622        7.72

可以看到模型收敛

方差组分为:

                   Estimated (co)-variance components
                          ----------------------------------
                                                            
             Parameter vector for L at convergence
             Asymptotic SE based on AI-information matrix
                                                            
               No          Parameter             Asymp. S.E.
                                                            
                1           4.93404               1.76364
                2           8.46217               1.63818
                3           7.72445              0.329943

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

3.8 对比asreml的结果:

代码:

library(asreml)

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

「方差组分:」

> summary(mod)$varcomp
                    gamma component std.error   z.ratio constraint
ped(ANIMAL)!ped 0.6387559  4.934041 1.7636385  2.797649   Positive
ide(ANIMAL)!id  1.0955038  8.462169 1.6381812  5.165588   Positive
R!variance      1.0000000  7.724454 0.3299432 23.411466   Positive

「遗传力:」

> pin(mod,h2 ~ V1/(V1+V2+V3))
   Estimate         SE
h2 0.233612 0.07907261

3.9 DMU和asreml比较

两者方差组分一致.

4. 多性状动物模型

4.1 使用数据

本次主要是演示如何使用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

4.2 需要做的处理

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

4.3 编写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

4.4 执行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>

4.5 查看结果

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

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

4.6 对比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

4.7 DMU和asreml比较

两者方差组分一致.

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

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

5.1 示例数据

「数据使用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

5.3 需要做的处理

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

5.4 编写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

5.5 执行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>

5.6 查看结果

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

4.7 对比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

5.8 DMU和asreml比较

两者方差组分一致.

6. DMU linux下语法高亮设置

用vim编程时, DMU的关键词没有语法高亮, 看着不舒服, 就进行一下设置, 并记录过程.

6.1 设置的效果如下

在这里插入图片描述

6.2 设置流程

本次设置的比较简单, 将关键词分为:

模型model, 比如DMU1, DMU2...
不同组成part, 比如DATA, VARIATE, MODEL...
不同结构类型type, 比如PED, COR....
新建DIR.vim文件, 里面设置相关参数
新建DIR_suffix.vim文件, 设置后缀名读取

6.3 DIR.vim文件:

"------------------------------------------------------------------------------"
" Description
"------------------------------------------------------------------------------"

" vim syntax highlighting file for DMU programs

" Author:   Deng Fei <dengfei_2013@163.com>
" Created:  Unknown
" Modified: 2018-11-18
" License:  GPLv2

syn keyword model DMU1 DMU4 DMU5 DMUAI RJMC

syn keyword part COMMENT ANALYSE DATA VARIABLE MODEL GLMM GLMM_PRED REDUCE MIXTURE VAR_STR VAR_REST PRECOND SOLUTION PRIOR RESIDUALS TRAITS ABSORB RANDOM REGRES NOCOV

syn keyword type PED DOM COR GRE PGMIX ABS_QTL GROUP VAR COV COR V_RATIO ASCII

hi model             ctermfg=Yellow
hi part              ctermfg=red
hi type              ctermfg=Green

将上面内容, 保存为:DIR.vim文件, 放到:~/.vim/syntax文件夹中. 「如果没有syntax文件夹, 就新建一个.」

cp DIR.vim ~/.vim/syntax/

6.4 DIR_suffix.vim文件:

au BufRead,BufNewFile *.DIR            set filetype=DIR

将上面内容保存到DIR_suffix.vim问价中, 放到:~/.vim/ftdetect文件夹中. 「如果没有ftdetect文件夹, 就新建一个.」

cp DIR_suffix.vim ~/.vim/ftdetect/

6.5 测试

使用下面代码, 新建文件test.DIR, 然后使用vim打开, 查看语法高亮是否成功:

$ANALYSE 1 1 0 0
$DATA  ASCII (8,15,0) dat_dmu.txt
$VARIABLE
ID	F1	F2	F3	F4	F5
y1
$MODEL
1
0
1 0 5 2 3 4 5 1
1
0
0
$VAR_STR 1 PED 2 ASCII ped_dmu.txt
$DMUAI
10
1d-7
1d-6
1
0

「效果如下:」

7. DMU在windows下安装测试

7.1 下载地址

❝官网好像挂掉了，可以公众号回复“DMU”获取软件和使用说明。 ❞

下载地址:http://dmu.agrsci.dk/

在这里插入图片描述

64为电脑安装「DMUv6-R5-2-EM64T.msi」, 32为电脑安装「DMUv6-R5-2.msi」

7.2 安装DMU

在这里插入图片描述

7.3 测试DMU是否安装成功

** 7.3.1 保存测试数据**

#dir.create("d:/ddmu-test")
setwd("d:/ddmu-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
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)

「7.3.2 查看数据」

在这里插入图片描述

「7.3.3 将下面内容保存为:Uni_animal.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

在这里插入图片描述

「注意, 先建立txt, 然后将其重命名为:Uni_animal.DIR, 这里DIR是后缀, 而不是默认的txt」

在这里插入图片描述

7.4 找到DMU的安装路径, copy文件

我的默认的安装路径如下:C:\Program Files (x86)\QGG-AU\DMUv6\R5.2-EM64T\bin

将这些文件copy到你的数据文件夹下.

数据文件夹里面的内容有:

7.5 打开cmd

进入D盘
进入文件夹:ddmu-test
查看当前文件内容

在这里插入图片描述

7.6 执行命令

run_dmuai.bat Uni_animal

在这里插入图片描述

7.7 查看结果

在这里插入图片描述

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-05-04，如有侵权请联系 cloudcommunity@tencent.com 删除

网站

编程算法

本文分享自育种数据分析之放飞自我微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

网站

编程算法

DMU遗传评估从入门到出家系列

DMU遗传评估从入门到出家系列

这是之前写的DMU学习笔记，做一个汇总，更正了一些错误，添加了Linux版的安装包。

DMU遗传评估从入门到出家系列

1. 主要参数介绍

1.1 DMU软件介绍

1.2 学习DMU初衷

1.3 DMU语法介绍

1.4 文件输出

1.5 命令行文件执行

2. 单性状动物模型

2.1 示例数据

2.2 数据预处理

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

2.4 文件类型

2.5 编写DIR文件

2.6 执行DIR文件

2.7 查看结果

2.8 对比asreml的结果:

2.9 DMU和asreml比较

3. 单性状重复力模型

3.1 概念解析

3.2 使用数据

3.3 需要做的处理

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

3.5 编写DIR文件

3.6 执行DIR文件

3.7 查看结果

3.8 对比asreml的结果:

3.9 DMU和asreml比较

4. 多性状动物模型

4.1 使用数据

4.2 需要做的处理

4.3 编写DIR文件

4.4 执行DIR文件

4.5 查看结果

4.6 对比asreml的结果:

4.7 DMU和asreml比较

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

5.1 示例数据

5.3 需要做的处理

5.4 编写DIR文件

5.5 执行DIR文件

5.6 查看结果

4.7 对比asreml的结果:

5.8 DMU和asreml比较

6. DMU linux下语法高亮设置

6.1 设置的效果如下

6.2 设置流程

6.3 DIR.vim文件:

6.4 DIR_suffix.vim文件:

6.5 测试

7. DMU在windows下安装测试

7.1 下载地址

7.2 安装DMU

7.3 测试DMU是否安装成功

7.4 找到DMU的安装路径, copy文件

7.5 打开cmd

7.6 执行命令

7.7 查看结果

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐