单性状动物模型-母体效应
本次主要是演示如何使用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(第五列)
固定因子: 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
1
执行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>
查看结果
在文件*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
对比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.08842455
DMU和asreml比较
两者方差组分一致.
