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

这是之前写的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和Ve1 1 1 Vg11 1 2 1 Vg12 1 2 2 Vg22 2 1 1 Ve11 2 2 1 Ve21 2 2 2 Ve22$VAR_REST(可选项, 主要是固定初始值)1.4 文件输出lst描述统计, 模型迭代, 方差组分估计PAROUT方差组分估计(行列形式显示)PAROUT-STD方差组分及标准误(计算遗传力)LLIK最后一次迭代情况1.5 命令行文件执行run_dmuai运行dmuai程序run_dmu4运行dmu4程序run_dmu5run_rjmc2. 单性状动物模型

本次主要是演示如何使用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-testdir.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.72.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 file2.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.074326242.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表示系谱中有三列数据, 无出生时间一列, 缺失值为03.3 需要做的处理系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化在保存数据时, 去掉行头编辑DIR文件3.4 使用R语言清洗数据, 并保存数据到D盘dmu-testdat = 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 13.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 file3.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.079072613.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

「系谱中：」

有三列数据, 无出生时间一列, 缺失值为04.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 14.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 Positive4.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 15.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.088424555.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-AUDMUv6R5.2-EM64Tbin

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

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

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

在这里插入图片描述

7.6 执行命令run_dmuai.bat Uni_animal

在这里插入图片描述

7.7 查看结果

在这里插入图片描述

---来自腾讯云社区的---邓飞