DMU-遺傳引數評估-學習筆記2
阿新 • • 發佈：2018-11-15
單性狀動物模型
本次主要是演示如何使用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比較
兩者方差組分和遺傳力一致.
