博文
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比较
两者方差组分和遗传力一致.
本次分析所使用数据:
https://blog.sciencenet.cn/blog-2577109-1146234.html
上一篇：DMU-遗传参数评估-学习笔记1
下一篇：Python学生物统计-数据可视化-学习笔记5
