这是之前写的DMU学习笔记，做一个汇总，更正了一些错误，添加了Linux版的安装包。 DMU遗传评估从入门到出家系列 1. 主要参数介绍 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)算法来最大化约束似然函数。被估计的（协）方差组分的渐近标准误是从平均信息矩阵中获得的。
想试试DMU处理一批数据, 发现这个软件, 竟然没有一个合适的操作说明文档, 我手头上有苏国生老师的PPT中文版DMU操作说明, 但看起来还是费劲.
刚好自己在学习这个软件, 用实际数据来演示如何使用这个软件进行数据分析.
「我想从四部分进行:」
其它内容, 包括测定日模型(随机回归模型), 母体效应模型, GBLUP模型, 显性上位性模型, 一步法GS模型等等以后再做总结.
说明文档是作者写的, 一般来说作者都想通过逻辑的构建, 让读者了解软件的方方面面, 但是读者一开始接触软件时, 迫切的是想解决问题, 不是来学理论, 不是来学知识, 只是想解决问题. 但是大多数文档无法满足这些迫切的需求. 所以, 最好的操作说明, 就是有数据, 有模型, 有结果说明, 可以很快上手. 我写此操作说明的目的就在于此.
「软件组成, 主要包括四类程序」
DMU1
这个主要是为了整理数据和模型, 相当于预处理程序, 其它三个程序都要经过它的处理才能分析. 类似BLUPF90的renumf90程序. 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$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
1 1 1 Vg11PAROUT-STD方差组分及标准误(计算遗传力)2. 单性状动物模型 本次主要是演示如何使用DMU分析单性状动物模型.
「数据使用learnasreml包中的数据」
learnasreml是我编写的辅助学习asreml的R包, 里面有相关的数据和代码, 这里我们用其中的animalmodel.dat和animalmodel.ped的数据.
如果没有软件包, 首先安装:
library(devtools)看一下数据:
> summary(dat)数据中,
有因子4个: 分别是ANIMAL, MOTHER, BYEAR, SEX
有变量2个: 分别是BWT和TARSUS
缺失值为0
系谱中,
有三列数据, 无出生时间一列, 缺失值为0
系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化 dir.create("d:/dmu-test")「数据文件:」 「系谱文件:」
想要分析的模型:
固定因子: BYEAR和SEX(第三列, 第四列) 所以这里编写DIR，第一部分, 是注释, 这里所写的东西会输出到结果文件, 基本上就是模型的解释, 这部分没有强制要求, 可以省略
$COMMENT第二部分是分析方法, 默认是AI
$ANALYSE 1 1 0 0第三部分是定义因子数和变量书, 以及文件位置:
$DATA ASCII (4,2,0) d:/dmu-test/animal-model.txt第四部分, 定义变量名, 也是为了方便结果输出, 相当于数据的行头名
$VARIABLE第五部分, 有6行, 定义模型
整体来说是:
第三行: 第1个y变量, 0无权重考虑,3个因子,第3列是第一个固定因子, 第4列第二个固定因子, 第1列是随机因子 $MODEL第六部分: 指定系谱
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt第七部分: 指定初始值(可以省略)
$PRIOR完整DIR文件, 放入model.txt中, 然后重命名为: Uni_animalmodel.DIR
$COMMENT这里运行的run_dmuai.bat, 将DMU安装路径下的文件run_dmuai.bat拷贝到d:/dmu-test文件夹, 在终端cmd界面键入:
run_dmuai.bat Uni_animalmodel执行结果:
D:\dmu-test>run_dmuai.bat Uni_animalmodel在文件*lst中有估算的方差组分, 结果如下:
SUMMARY OF MINIMIZATION PROCESS可以看到模型收敛
方差组分为:
Estimated (co)-variance components遗传力为:
Phenotypic co-variance matrix可以看出, 遗传力为0.536, 标准误为0.07
代码:
library(asreml)「方差组分:」
> summary(mod)$varcomp「遗传力:」
> pin(mod,h2 ~ V1/(V1+V2))两者方差组分和遗传力一致.
3. 单性状重复力模型 本次主要是演示如何使用DMU分析单性状重复力模型.
「重复力模型和动物模型的区别:」 不是所有的性状都可以分析重复力模型, 首先重复力模型是动物模型的拓展, 它适合一个个体多个观测值的情况.
只有这样的数据才可以将永久环境效应剖分出来.
「重复力是遗传力的上限:」 教科书上这样说, 这句话怎么理解呢?
首先, 我们认为
遗传力为
这里的Vg如果有重复测量的数据, 可以剖分为可以遗传的部分, 和不可以遗传的部分(永久环境效应), 那么遗传力的计算公式为:
重复力的计算公式为:
当Vpe为0时, 重复力=遗传力, 当Vpe>0时, 重复力>遗传力, 所以说重复力是遗传力的上限!
「数据使用learnasreml包中的数据」
learnasreml是我编写的辅助学习asreml的R包, 里面有相关的数据和代码, 这里我们用其中的repeatmodel.dat和repeatmodel.ped的数据.
如果没有软件包, 首先安装:
setwd("d:/dmu-test/")看一下数据:
> summary(dat)「数据介绍:」
有因子4个: 分别是 ANIMAL BYEAR AGE YEAR 系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化 dat = repeatmodel.dat想要分析的模型:
观测值: LAYDATE (第四列)
固定因子: 第二列, 第三列, 第四列
随机因子: ID, ide(ID)
所以这里编写DIR
第一部分, 是注释, 这里所写的东西会输出到结果文件, 基本上就是模型的解释, 这部分没有强制要求, 可以省略
$COMMENT第二部分是分析方法, 默认是AI
$ANALYSE 1 1 0 0第三部分是定义因子数和变量书, 以及文件位置:
$DATA ASCII (4,1,0) d:/dmu-test/repeat-model.txt第四部分, 定义变量名, 也是为了方便结果输出, 相当于数据的行头名
$VARIABLE第五部分, 有6行, 定义模型
整体来说是:
第一行: 单性状 # 1
第二行: 无吸收 # 0
第三行: 主要定义y变量, 固定因子, 随机因子
共五个因子(固定+随机, 固定写前面, 随机写后面) # 5 所以, 5个因子, 三个固定因子:2,3,4, 两个随机因子:1,1 #1 0 5 2 3 4 1 1 第四行: 有两个随机因子, 他们的关系是独立的, 所以是2 1 2
1第六部分: 指定系谱
$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt注意, 如果想要输出BLUP值, 定义:$DMUAI
$DMUAI完整DIR文件, 放入model.txt中, 然后重命名为: Uni_repeatmodel.DIR
$COMMENT这里运行的run_dmuai.bat, 将DMU安装路径下的文件run_dmuai.bat拷贝到d:/dmu-test文件夹, 在终端cmd界面键入:
run_dmuai.bat Uni_repeatmodel执行结果:
D:\dmu-test>run_dmuai.bat Uni_repeatmodel在文件*lst中有估算的方差组分, 结果如下:
SUMMARY OF MINIMIZATION PROCESS可以看到模型收敛
方差组分为:
Estimated (co)-variance components遗传力需要手动计算, 这里还没有找到解决方案.
代码:
library(asreml)「方差组分:」
> summary(mod)$varcomp「遗传力:」
> pin(mod,h2 ~ V1/(V1+V2+V3))两者方差组分一致.
4. 多性状动物模型 本次主要是演示如何使用DMU分析多性状动物模型.
「数据使用learnasreml包中的数据」
learnasreml是我编写的辅助学习asreml的R包, 里面有相关的数据和代码, 这里我们用其中的animalmodel.dat和animalmodel.ped的数据.
如果没有软件包, 首先安装:
setwd("d:/dmu-test/")看一下数据:
> summary(dat)「数据中：」
有因子4个: 分别是ANIMAL, MOTHER, BYEAR, SEX 「系谱中：」
系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化 想要分析的模型:
观测值: BWT(第五列), TARSUS (第六列） 固定因子: BYEAR和SEX(第三列, 第四列) 所以这里编写DIR
第一部分, 是注释, 这里所写的东西会输出到结果文件, 基本上就是模型的解释, 这部分没有强制要求, 可以省略
$COMMENT第二部分是分析方法, 默认是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第五部分, 有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第六部分: 指定系谱
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt注意, 如果想要输出BLUP值, 定义:$DMUAI
$DMUAI完整DIR文件, 放入model.txt中, 然后重命名为: mul-animalmodel.DIR
$COMMENT这里运行的run_dmuai.bat, 将DMU安装路径下的文件run_dmuai.bat拷贝到d:/dmu-test文件夹, 在终端cmd界面键入:
run_dmuai.bat mul_animalmodel执行结果:
D:\dmu-test>run_dmuai.bat mul_animalmodel在文件*lst中有估算的方差组分, 结果如下:
Eval Criterion !!Delta!! !!Gradient!! Parameters可以看到模型收敛
方差组分为:
Estimated (co)-variance components遗传力需要手动计算, 这里还没有找到解决方案.
代码:
library(asreml)「方差组分:」
> summary(mod2)$varcomp两者方差组分一致.
5. 单性状动物模型-母体效应 本次主要是演示如何使用DMU分析单性状动物模型-母体效应.
注意，这里的母体效应不太严谨，具体解释： 一文讲清楚动物模型中的母体效应
「数据使用learnasreml包中的数据」
learnasreml是我编写的辅助学习asreml的R包, 里面有相关的数据和代码, 这里我们用其中的animalmodel.dat和animalmodel.ped的数据.
如果没有软件包, 首先安装:
setwd("d:/dmu-test/")看一下数据:
> summary(dat)**数据中, **
有因子4个: 分别是ANIMAL, MOTHER, BYEAR, SEX **系谱中, **
有三列数据, 无出生时间一列, 缺失值为0
系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化 想要分析的模型:
固定因子: BYEAR和SEX(第三列, 第四列) 所以这里编写DIR
第一部分, 是注释, 这里所写的东西会输出到结果文件, 基本上就是模型的解释, 这部分没有强制要求, 可以省略
$COMMENT第二部分是分析方法, 默认是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第五部分, 有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第六部分: 指定系谱
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt注意, 如果想要输出BLUP值, 定义:$DMUAI
$DMUAI完整DIR文件, 放入model.txt中, 然后重命名为: Uni_animalmodel-maternal.DIR
$COMMENT这里运行的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在文件*lst中有估算的方差组分, 结果如下:
SUMMARY OF MINIMIZATION PROCESS可以看到模型收敛
方差组分为:
Estimated (co)-variance components「可以看到:」 加性方差组分为: 2.2778
母体效应方差组分为: 1.10404
残差方差组分为: 1.6569
代码:
library(asreml)「方差组分:」
> summary(mod)$varcomp「遗传力:」
> pin(mod,h2 ~ V2/(V1+V2+V3))两者方差组分一致.
6. DMU linux下语法高亮设置 用vim编程时, DMU的关键词没有语法高亮, 看着不舒服, 就进行一下设置, 并记录过程.
本次设置的比较简单, 将关键词分为:
不同组成part, 比如DATA, VARIATE, MODEL... 不同结构类型type, 比如PED, COR.... 新建DIR_suffix.vim文件, 设置后缀名读取 将上面内容, 保存为:DIR.vim文件, 放到:~/.vim/syntax文件夹中. 「如果没有syntax文件夹, 就新建一个.」
cp DIR.vim ~/.vim/syntax/au BufRead,BufNewFile *.DIR set filetype=DIR将上面内容保存到DIR_suffix.vim问价中, 放到:~/.vim/ftdetect文件夹中. 「如果没有ftdetect文件夹, 就新建一个.」
cp DIR_suffix.vim ~/.vim/ftdetect/使用下面代码, 新建文件test.DIR, 然后使用vim打开, 查看语法高亮是否成功:
$ANALYSE 1 1 0 0「效果如下:」
7. DMU在windows下安装测试 ❝ 官网好像挂掉了，可以公众号回复“DMU”获取软件和使用说明。
❞ 下载地址:http://dmu./
64为电脑安装「DMUv6-R5-2-EM64T.msi」 , 32为电脑安装「DMUv6-R5-2.msi」
** 7.3.1 保存测试数据**
#dir.create("d:/ddmu-test") "d:/ddmu-test/" )library (devtools)"dengfei2013/learnasreml" )library (learnasreml)"animalmodel.dat" )"animalmodel.ped" )2018 library (data.table)# write.table(dat,"animal-model.txt",row.names = F,col.names = F) "animal-model.txt" ,sep = " " ,col.names = F )"animal-ped.txt" ,sep = " " ,col.names = F )「7.3.2 查看数据」
「7.3.3 将下面内容保存为:Uni_animal.DIR」
$COMMENT「注意, 先建立txt, 然后将其重命名为:Uni_animal.DIR, 这里DIR是后缀, 而不是默认的txt」
我的默认的安装路径如下:C:\Program Files (x86)\QGG-AU\DMUv6\R5.2-EM64T\bin
将这些文件copy到你的数据文件夹下.
数据文件夹里面的内容有:
run_dmuai.bat Uni_animal
「系谱文件:」
