博文
DMU-多性状动物模型-学习笔记4
||
多性状动物模型
本次主要是演示如何使用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
需要做的处理
- 系谱增加第四列出生时间, 因为数据都是数字, 没有字符串, 不需要转化
- 在保存数据时, 去掉行头
- 编辑DIR文件
编写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执行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>查看结果
在文件*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遗传力需要手动计算, 这里还没有找到解决方案.
对比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 PositiveDMU和asreml比较
两者方差组分一致.
https://blog.sciencenet.cn/blog-2577109-1147133.html
上一篇：DMU-单性状重复力模型-学习笔记3
下一篇：DMU-单性状动物模型-母体效应--学习笔记5
