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比较
两者方差组分和遗传力一致.
标签：性状,animalmodel,dat,DMU,dmu,ped,动物模型,Qu,BWT 来源： https://blog.51cto.com/yijiaobani/2867179