【原】DMU遺傳評估從入門到出家系列

育種數(shù)據(jù)分析 2021-11-18

展開全文

這是之前寫的DMU學(xué)習(xí)筆記，做一個(gè)匯總，更正了一些錯(cuò)誤，添加了Linux版的安裝包。

DMU遺傳評估從入門到出家系列

1. 主要參數(shù)介紹

1.1 DMU軟件介紹

DMU是一個(gè)數(shù)量遺傳學(xué)工具包，主要功能包括估計(jì)方差組分和固定效應(yīng)，預(yù)測育種值。DMU的開發(fā)歷史可以追溯到25年前，大部分功能基于數(shù)量遺傳學(xué)研究的需求而開發(fā)。在丹麥動(dòng)物育種研究中，DMU是一個(gè)主要的統(tǒng)計(jì)研究工具（估計(jì)和預(yù)測）。此外，DMU也應(yīng)用于丹麥牛，羊，貂和馬等常規(guī)遺傳評估研究。因此，DMU不但在一些特定的項(xiàng)目中具備高性能優(yōu)勢，也適用于常規(guī)數(shù)量遺傳學(xué)研究?！癉MU”名稱最初來自于程序包中用來進(jìn)行初始化的過程名字縮寫。這些過程利用約束最大似然法（REML），通過Derivative-free方式執(zhí)行MUltivariate analysis，因此得名DMU。但是，在當(dāng)前的DMU版本中，并不包括DF-REML模塊，現(xiàn)在D僅代表DJF（丹麥農(nóng)業(yè)科學(xué)學(xué)院的縮寫）。

DMU安裝包有很多模塊，如DMU1、DMU4、DMU5、DMUAI和RJMC。DMUAI模塊可利用平均信息限制最大似然(AI-REML) (Jensen et al. 1997)算法進(jìn)行（協(xié)）方差組分的估計(jì)。AI是通過平均觀察和預(yù)期信息的信息矩陣得到的。該模塊還可以使用期望最大化(EM)算法來最大化約束似然函數(shù)。被估計(jì)的（協(xié)）方差組分的漸近標(biāo)準(zhǔn)誤是從平均信息矩陣中獲得的。

1.2 學(xué)習(xí)DMU初衷

想試試DMU處理一批數(shù)據(jù), 發(fā)現(xiàn)這個(gè)軟件, 竟然沒有一個(gè)合適的操作說明文檔, 我手頭上有蘇國生老師的PPT中文版DMU操作說明, 但看起來還是費(fèi)勁.

剛好自己在學(xué)習(xí)這個(gè)軟件, 用實(shí)際數(shù)據(jù)來演示如何使用這個(gè)軟件進(jìn)行數(shù)據(jù)分析.

「我想從四部分進(jìn)行:」

1, DMU語法介紹
2, 單性狀動(dòng)物模型
3, 單性狀重復(fù)力模型
4, 多性狀動(dòng)物模型

其它內(nèi)容, 包括測定日模型(隨機(jī)回歸模型), 母體效應(yīng)模型, GBLUP模型, 顯性上位性模型, 一步法GS模型等等以后再做總結(jié).

說明文檔是作者寫的, 一般來說作者都想通過邏輯的構(gòu)建, 讓讀者了解軟件的方方面面, 但是讀者一開始接觸軟件時(shí), 迫切的是想解決問題, 不是來學(xué)理論, 不是來學(xué)知識(shí), 只是想解決問題. 但是大多數(shù)文檔無法滿足這些迫切的需求. 所以, 最好的操作說明, 就是有數(shù)據(jù), 有模型, 有結(jié)果說明, 可以很快上手. 我寫此操作說明的目的就在于此.

1.3 DMU語法介紹

「軟件組成, 主要包括四類程序」

DMU1 這個(gè)主要是為了整理數(shù)據(jù)和模型, 相當(dāng)于預(yù)處理程序, 其它三個(gè)程序都要經(jīng)過它的處理才能分析. 類似BLUPF90的renumf90程序.
DMUAI 這個(gè)主要估算方差組分的程序
DMU4和DMU5 DMU4主要是求解混合線性方程組, 它不估算方差組分, 只求解. 類似BLUPF90包中的blupf90程序.DMU5功能和DMU4類似, 也是求解方程組, 適用于大數(shù)據(jù)
RGMC 主要是貝葉斯抽樣, 估算方差組分, 計(jì)算育種值.

「數(shù)據(jù)和系譜及逆矩陣格式」

全部數(shù)據(jù), 不要有行頭
數(shù)據(jù)中不能含有字符, 字母, 都必須是數(shù)字
逆矩陣可以是下三角或者上三角矩陣的三列形式
系譜數(shù)據(jù)包括四列: ID, Sire, Dam, Birth
數(shù)據(jù)中, 因子(ID, Sex...)放在前面, 觀測值(y1, y2, y3)放在后面, 因子用整數(shù)表示, 不能含有字母

因此, 在進(jìn)行分析之前, 首先需要對數(shù)據(jù)進(jìn)行轉(zhuǎn)化, 比如系譜要變?yōu)檎麛?shù), 要有第四列信息出生信息, 如果沒有, 就寫成2018年就行. 數(shù)據(jù)中也要重新編號(hào), 特別是某些因子含有字母, 需要轉(zhuǎn)化為數(shù)字. 可以使用R語言進(jìn)行轉(zhuǎn)化, 將系譜的所有水平編號(hào)為1...n, 然后替換. 將數(shù)據(jù)的所有水平, 重新編碼.

「參數(shù)文件」文件名為name.DIR, 其中name為程序名稱, DIR必須要有, 并保持大寫.

$COMMENT文件注釋, 一般是解釋你所使用的模型
$ANALYSIS你分析所使用的模型, 如果你需要估算方差組分, 那么簡單寫為:$ANALYSIS 1 1 0 0
$DATA指定數(shù)據(jù)格式,因子數(shù)目, 觀測值數(shù)目, 缺失值, 和數(shù)據(jù)位置如果是txt文件, 有5個(gè)因子, 4個(gè)觀測值, 缺失值-999, 在D盤根目錄$DATA ASCII(5,4,-999) d:/dat.txt
$VARIABLE寫出因子和變量的名稱, 第一行為因子, 第二行為變量

ID Loc Year Herd Sex Hy
y1 y2 y3 y4

$MODEL指定分析模型中, 觀測值個(gè)數(shù), 固定因子, 隨機(jī)因子比如單性狀, 正態(tài)數(shù)據(jù)

1 1 0 0 0

比如二性狀, 正態(tài)數(shù)據(jù)

2 2 0 0 0

固定因子: 每個(gè)性狀一行, 包含若干整數(shù) 單性狀中, y1 = Loc + Year + Herd + Sex, random = ID

1 0 5 1 2 3 4 5

隨機(jī)因子: 每個(gè)性狀一行, 包含若干整數(shù) 1

$VAR_STR定義方差協(xié)方差結(jié)構(gòu) 可以支持系譜, 和自定義關(guān)系矩陣inv 定義系譜文件:$VAR_STR 2 PED 2 ASCII ped.txt定義逆矩陣:$VAR_STR 1 COR ASCII ginv
$PRIOR定義初始值, 不過不定義, 默認(rèn)是方差組分為1, 協(xié)方差組分為0, 定義格式, 下三角行列形式. 比如兩性狀, Vg和Ve

$VAR_REST(可選項(xiàng), 主要是固定初始值)

1.4 文件輸出

lst描述統(tǒng)計(jì), 模型迭代, 方差組分估計(jì)
PAROUT方差組分估計(jì)(行列形式顯示)
PAROUT-STD方差組分及標(biāo)準(zhǔn)誤(計(jì)算遺傳力)
LLIK最后一次迭代情況

1.5 命令行文件執(zhí)行

run_dmuai運(yùn)行dmuai程序
run_dmu4運(yùn)行dmu4程序
run_dmu5
run_rjmc

2. 單性狀動(dòng)物模型

本次主要是演示如何使用DMU分析單性狀動(dòng)物模型.

2.1 示例數(shù)據(jù)

「數(shù)據(jù)使用learnasreml包中的數(shù)據(jù)」

learnasreml是我編寫的輔助學(xué)習(xí)asreml的R包, 里面有相關(guān)的數(shù)據(jù)和代碼, 這里我們用其中的animalmodel.dat和animalmodel.ped的數(shù)據(jù).

如果沒有軟件包, 首先安裝:

library(devtools)
install_github("dengfei2013/learnasreml")
library(learnasreml)
data("animalmodel.dat")
data("animalmodel.ped")

dat = animalmodel.dat
ped = animalmodel.ped

summary(dat)
summary(ped)

看一下數(shù)據(jù):

> 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

數(shù)據(jù)中, 有因子4個(gè): 分別是ANIMAL, MOTHER, BYEAR, SEX 有變量2個(gè): 分別是BWT和TARSUS 缺失值為0

系譜中, 有三列數(shù)據(jù), 無出生時(shí)間一列, 缺失值為0

2.2 數(shù)據(jù)預(yù)處理

系譜增加第四列出生時(shí)間, 因?yàn)閿?shù)據(jù)都是數(shù)字, 沒有字符串, 不需要轉(zhuǎn)化
在保存數(shù)據(jù)時(shí), 去掉行頭
編輯DIR文件

2.3 使用R語言清洗數(shù)據(jù)并保存數(shù)據(jù)到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)

2.4 文件類型

「數(shù)據(jù)文件:」「系譜文件:」

2.5 編寫DIR文件

想要分析的模型:

觀測值: BWT(第五列)
固定因子: BYEAR和SEX(第三列, 第四列)
隨機(jī)因子: ID

所以這里編寫DIR，第一部分, 是注釋, 這里所寫的東西會(huì)輸出到結(jié)果文件, 基本上就是模型的解釋, 這部分沒有強(qiáng)制要求, 可以省略

$COMMENT
Estimate variance components for BWT
Model
y: BWT
fixed: BYEAR + SEX
random: ANIMAL

第二部分是分析方法, 默認(rèn)是AI

$ANALYSE 1 1 0 0

第三部分是定義因子數(shù)和變量書, 以及文件位置:

$DATA  ASCII (4,2,0) d:/dmu-test/animal-model.txt

第四部分, 定義變量名, 也是為了方便結(jié)果輸出, 相當(dāng)于數(shù)據(jù)的行頭名

$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS

第五部分, 有6行, 定義模型整體來說是:

第一行: 單性狀
第二行: 無吸收
第三行: 第1個(gè)y變量, 0無權(quán)重考慮,3個(gè)因子,第3列是第一個(gè)固定因子, 第4列第二個(gè)固定因子, 第1列是隨機(jī)因子
第四行:1個(gè)隨機(jī)因子
第五行: 無回歸項(xiàng)
第六行: 無約束

$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

2.6 執(zhí)行DIR文件

這里運(yùn)行的run_dmuai.bat, 將DMU安裝路徑下的文件run_dmuai.bat拷貝到d:/dmu-test文件夾, 在終端cmd界面鍵入:

run_dmuai.bat Uni_animalmodel

執(zhí)行結(jié)果:

D:\dmu-test>run_dmuai.bat Uni_animalmodel

D:\dmu-test>Echo OFF
Starting DMU using Uni_animalmodel.DIR as directive file

2.7 查看結(jié)果

在文件*lst中有估算的方差組分, 結(jié)果如下:

                             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, 標(biāo)準(zhǔn)誤為0.07

2.8 對比asreml的結(jié)果:

代碼:

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

2.9 DMU和asreml比較

兩者方差組分和遺傳力一致.

3. 單性狀重復(fù)力模型

本次主要是演示如何使用DMU分析單性狀重復(fù)力模型.

3.1 概念解析

「重復(fù)力模型和動(dòng)物模型的區(qū)別:」不是所有的性狀都可以分析重復(fù)力模型, 首先重復(fù)力模型是動(dòng)物模型的拓展, 它適合一個(gè)個(gè)體多個(gè)觀測值的情況.

比如豬的產(chǎn)仔數(shù), 一個(gè)母豬有多個(gè)胎次
比如雞的產(chǎn)蛋, 不同時(shí)間段, 雞都有產(chǎn)蛋量
牛的產(chǎn)奶量, 不同的測定日, 產(chǎn)奶量不同
豬的飼料消耗, 也是重復(fù)測量的數(shù)據(jù)

只有這樣的數(shù)據(jù)才可以將永久環(huán)境效應(yīng)剖分出來.

「重復(fù)力是遺傳力的上限:」教科書上這樣說, 這句話怎么理解呢? 首先, 我們認(rèn)為

遺傳力為

這里的Vg如果有重復(fù)測量的數(shù)據(jù), 可以剖分為可以遺傳的部分, 和不可以遺傳的部分(永久環(huán)境效應(yīng)), 那么遺傳力的計(jì)算公式為:

重復(fù)力的計(jì)算公式為:

當(dāng)Vpe為0時(shí), 重復(fù)力=遺傳力, 當(dāng)Vpe>0時(shí), 重復(fù)力>遺傳力, 所以說重復(fù)力是遺傳力的上限!

3.2 使用數(shù)據(jù)

「數(shù)據(jù)使用learnasreml包中的數(shù)據(jù)」

learnasreml是我編寫的輔助學(xué)習(xí)asreml的R包, 里面有相關(guān)的數(shù)據(jù)和代碼, 這里我們用其中的repeatmodel.dat和repeatmodel.ped的數(shù)據(jù).

如果沒有軟件包, 首先安裝:

setwd("d:/dmu-test/")
library(devtools)
# install_github("dengfei2013/learnasreml")
library(learnasreml)
data("repeatmodel.dat")
data("repeatmodel.ped")

dat = repeatmodel.dat
ped = repeatmodel.ped

summary(dat)
summary(ped)

看一下數(shù)據(jù):

> summary(dat)
     ANIMAL         BYEAR      AGE          YEAR         LAYDATE
 1      :   5   1000   : 109   2:308   1004   :  79   Min.   : 0.00
 3      :   5   1001   :  98   3:322   1005   :  78   1st Qu.:20.00
 9      :   5   999    :  86   4:339   1003   :  69   Median :24.00
 17     :   5   1002   :  85   5:315   1006   :  64   Mean   :23.54
 42     :   5   987    :  70   6:323   1002   :  60   3rd Qu.:27.00
 50     :   5   989    :  66           988    :  54   Max.   :41.00
 (Other):1577   (Other):1093           (Other):1203
> 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

「數(shù)據(jù)介紹:」

有因子4個(gè): 分別是 ANIMAL BYEAR AGE YEAR
有變量1個(gè): LAYDATE
缺失值用0表示
系譜中有三列數(shù)據(jù), 無出生時(shí)間一列, 缺失值為0

3.3 需要做的處理

系譜增加第四列出生時(shí)間, 因?yàn)閿?shù)據(jù)都是數(shù)字, 沒有字符串, 不需要轉(zhuǎn)化
在保存數(shù)據(jù)時(shí), 去掉行頭
編輯DIR文件

3.4 使用R語言清洗數(shù)據(jù), 并保存數(shù)據(jù)到D盤dmu-test

dat = repeatmodel.dat
ped = repeatmodel.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,"repeat-model.txt",sep = " ",col.names = F)
fwrite(dmuped,"repeat-ped.txt",sep = " ",col.names = F)

3.5 編寫DIR文件

想要分析的模型:

觀測值: LAYDATE (第四列)
固定因子: 第二列, 第三列, 第四列
隨機(jī)因子: ID, ide(ID)

所以這里編寫DIR 第一部分, 是注釋, 這里所寫的東西會(huì)輸出到結(jié)果文件, 基本上就是模型的解釋, 這部分沒有強(qiáng)制要求, 可以省略

$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)

第二部分是分析方法, 默認(rèn)是AI

$ANALYSE 1 1 0 0

第三部分是定義因子數(shù)和變量書, 以及文件位置:

$DATA  ASCII (4,1,0) d:/dmu-test/repeat-model.txt

第四部分, 定義變量名, 也是為了方便結(jié)果輸出, 相當(dāng)于數(shù)據(jù)的行頭名

$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT

第五部分, 有6行, 定義模型整體來說是: 第一行: 單性狀 # 1 第二行: 無吸收 # 0 第三行: 主要定義y變量, 固定因子, 隨機(jī)因子

分析的是第一個(gè)變量 # 1
無權(quán)重考慮 # 0
共五個(gè)因子(固定+隨機(jī), 固定寫前面, 隨機(jī)寫后面) # 5
第一個(gè)固定因子是第二列因子 #2 #BYEAR
第二個(gè)固定一致是第三列因子 #3 #AGE
第三個(gè)固定因子是第四列 #4 #YEAR
第四個(gè)隨機(jī)因子是第一列 #1 #ANIMAL
第五個(gè)隨機(jī)因子是第一列 #1 #ANIMAL
所以, 5個(gè)因子, 三個(gè)固定因子:2,3,4, 兩個(gè)隨機(jī)因子:1,1 #1 0 5 2 3 4 1 1

第四行: 有兩個(gè)隨機(jī)因子, 他們的關(guān)系是獨(dú)立的, 所以是2 1 2

1
0
1 0 5 2 3 4 1 1
2 1 2
0
0

第六部分: 指定系譜

$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt

注意, 如果想要輸出BLUP值, 定義:$DMUAI

$DMUAI
10
1D-7
1D-6
1

完整DIR文件, 放入model.txt中, 然后重命名為: Uni_repeatmodel.DIR

$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)

$ANALYSE 1 1 0 0
$DATA  ASCII (4,1,0) d:/dmu-test/repeat-model.txt
$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT
$MODEL
1
0
1 0 5 2 3 4 1 1
2 1 2
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt


$DMUAI
10
1D-7
1D-6
1

3.6 執(zhí)行DIR文件

這里運(yùn)行的run_dmuai.bat, 將DMU安裝路徑下的文件run_dmuai.bat拷貝到d:/dmu-test文件夾, 在終端cmd界面鍵入:

run_dmuai.bat Uni_repeatmodel

執(zhí)行結(jié)果:

D:\dmu-test>run_dmuai.bat Uni_repeatmodel

D:\dmu-test>Echo OFF
Starting DMU using Uni_repeatmodel.DIR as directive file

3.7 查看結(jié)果

在文件*lst中有估算的方差組分, 結(jié)果如下:

                             SUMMARY OF MINIMIZATION PROCESS
                                                            
 Eval Criterion    !!Delta!!   !!Gradient!!                   Parameters
 ---- ---------     ---------  ------------ |------------------------------------------
   1   12629.2     0.8574        4.330      |    1.8100        1.8898        1.8705
  
   2   8234.59      1.370        6.822      |    2.9917        3.3812        3.2879
  
   3   6444.28      1.776        8.529      |    4.2397        5.4642        5.1761
  
   4   5857.47      1.566        6.869      |    4.9013        7.4662        6.8832
  
   5   5736.36     0.6798        2.497      |    4.9407        8.3737        7.6324
  
   6   5727.01     0.7387E-01   0.2311      |    4.9325        8.4634        7.7233
  
   7   5726.91     0.1399E-02   0.1596E-02  |    4.9341        8.4621        7.7245
  
   8   5726.91     0.7706E-04   0.5119E-04  |    4.9340        8.4622        7.7245
  
   9   5726.91     0.4564E-05   0.2610E-05  |    4.9340        8.4622        7.7245
  
  10   5726.91     0.2695E-06   0.1558E-06  |    4.9340        8.4622        7.72

可以看到模型收斂

方差組分為:

                   Estimated (co)-variance components
                          ----------------------------------
                                                            
             Parameter vector for L at convergence
             Asymptotic SE based on AI-information matrix
                                                            
               No          Parameter             Asymp. S.E.
                                                            
                1           4.93404               1.76364
                2           8.46217               1.63818
                3           7.72445              0.329943

遺傳力需要手動(dòng)計(jì)算, 這里還沒有找到解決方案.

3.8 對比asreml的結(jié)果:

代碼:

library(asreml)

head(dat)
dat[dat$LAYDATE==0,]$LAYDATE=NA
ainv = asreml.Ainverse(ped)$ginv
mod = asreml(LAYDATE ~ BYEAR + AGE + YEAR, random = ~ ped(ANIMAL)+ide(ANIMAL), ginverse = list(ANIMAL=ainv),data=dat)
summary(mod)$varcomp
pin(mod,h2 ~ V1/(V1+V2+V3))

「方差組分:」

> summary(mod)$varcomp
                    gamma component std.error   z.ratio constraint
ped(ANIMAL)!ped 0.6387559  4.934041 1.7636385  2.797649   Positive
ide(ANIMAL)!id  1.0955038  8.462169 1.6381812  5.165588   Positive
R!variance      1.0000000  7.724454 0.3299432 23.411466   Positive

「遺傳力:」

> pin(mod,h2 ~ V1/(V1+V2+V3))
   Estimate         SE
h2 0.233612 0.07907261

3.9 DMU和asreml比較

兩者方差組分一致.

4. 多性狀動(dòng)物模型

4.1 使用數(shù)據(jù)

本次主要是演示如何使用DMU分析多性狀動(dòng)物模型.

「數(shù)據(jù)使用learnasreml包中的數(shù)據(jù)」

learnasreml是我編寫的輔助學(xué)習(xí)asreml的R包, 里面有相關(guān)的數(shù)據(jù)和代碼, 這里我們用其中的animalmodel.dat和animalmodel.ped的數(shù)據(jù).

如果沒有軟件包, 首先安裝:

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)

看一下數(shù)據(jù):

> 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

「數(shù)據(jù)中：」

有因子4個(gè): 分別是ANIMAL, MOTHER, BYEAR, SEX
有變量2個(gè): 分別是BWT和TARSUS
缺失值為0

「系譜中：」

有三列數(shù)據(jù), 無出生時(shí)間一列, 缺失值為0

4.2 需要做的處理

系譜增加第四列出生時(shí)間, 因?yàn)閿?shù)據(jù)都是數(shù)字, 沒有字符串, 不需要轉(zhuǎn)化
在保存數(shù)據(jù)時(shí), 去掉行頭
編輯DIR文件

4.3 編寫DIR文件

想要分析的模型:

觀測值: BWT(第五列), TARSUS (第六列）
固定因子: BYEAR和SEX(第三列, 第四列)
隨機(jī)因子: ID

所以這里編寫DIR 第一部分, 是注釋, 這里所寫的東西會(huì)輸出到結(jié)果文件, 基本上就是模型的解釋, 這部分沒有強(qiáng)制要求, 可以省略

$COMMENT

Model
y: BWT  TARSUS
fixed: BYEAR + SEX
random: ANIMAL

第二部分是分析方法, 默認(rèn)是AI

$ANALYSE 1 1 0 0

第三部分是定義因子數(shù)和變量數(shù), 以及文件位置:

$DATA  ASCII (4,2,0) d:/dmu-test/animal-model.txt

?
上面的意思是, 數(shù)據(jù)是ASCII格式, 有4個(gè)固子, 2個(gè)變量, 缺失值用0表示, 數(shù)據(jù)的絕對路徑是: d:/dmu-test/animal-model.txt
?

第四部分, 定義變量名, 也是為了方便結(jié)果輸出, 相當(dāng)于數(shù)據(jù)的行頭名

$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS

第五部分, 有6行, 定義模型整體來說是:

第一行: 兩性狀 # 2
第二行: 1性狀無吸收 # 0
第三行: 2性狀無吸收 # 0
第四行: 1性狀, 是由3個(gè)因子, 兩個(gè)固定因子:3,4, 一個(gè)隨機(jī)因子:1 # 1 0 3 3 4 1
第五行: 2性狀, 是由3個(gè)因子, 兩個(gè)固定因子:3,4, 一個(gè)隨機(jī)因子:1 # 2 0 3 3 4 1
第六行: 性狀1, 1個(gè)隨機(jī)因子 # 1 1
第七行: 性狀2, 1個(gè)隨機(jī)因子 # 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

4.4 執(zhí)行DIR文件

這里運(yùn)行的run_dmuai.bat, 將DMU安裝路徑下的文件run_dmuai.bat拷貝到d:/dmu-test文件夾, 在終端cmd界面鍵入:

run_dmuai.bat mul_animalmodel

執(zhí)行結(jié)果:

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>

4.5 查看結(jié)果

在文件*lst中有估算的方差組分, 結(jié)果如下:

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

遺傳力需要手動(dòng)計(jì)算, 這里還沒有找到解決方案.

4.6 對比asreml的結(jié)果:

代碼:

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   Positive

4.7 DMU和asreml比較

兩者方差組分一致.

5. 單性狀動(dòng)物模型-母體效應(yīng)

本次主要是演示如何使用DMU分析單性狀動(dòng)物模型-母體效應(yīng).

注意，這里的母體效應(yīng)不太嚴(yán)謹(jǐn)，具體解釋：一文講清楚動(dòng)物模型中的母體效應(yīng)

5.1 示例數(shù)據(jù)

「數(shù)據(jù)使用learnasreml包中的數(shù)據(jù)」

learnasreml是我編寫的輔助學(xué)習(xí)asreml的R包, 里面有相關(guān)的數(shù)據(jù)和代碼, 這里我們用其中的animalmodel.dat和animalmodel.ped的數(shù)據(jù).

如果沒有軟件包, 首先安裝:

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)

看一下數(shù)據(jù):

> 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

**數(shù)據(jù)中, **

有因子4個(gè): 分別是ANIMAL, MOTHER, BYEAR, SEX
有變量2個(gè): 分別是BWT和TARSUS
缺失值為0

**系譜中, ** 有三列數(shù)據(jù), 無出生時(shí)間一列, 缺失值為0

5.3 需要做的處理

系譜增加第四列出生時(shí)間, 因?yàn)閿?shù)據(jù)都是數(shù)字, 沒有字符串, 不需要轉(zhuǎn)化
在保存數(shù)據(jù)時(shí), 去掉行頭
編輯DIR文件

5.4 編寫DIR文件

想要分析的模型:

觀測值: BWT(第五列)
固定因子: BYEAR和SEX(第三列, 第四列)
隨機(jī)因子: ID + MOTHER

所以這里編寫DIR 第一部分, 是注釋, 這里所寫的東西會(huì)輸出到結(jié)果文件, 基本上就是模型的解釋, 這部分沒有強(qiáng)制要求, 可以省略

$COMMENT

Model
y: BWT  TARSUS
fixed: BYEAR + SEX
random: ANIMAL+MOTHER

第二部分是分析方法, 默認(rèn)是AI

$ANALYSE 1 1 0 0

第三部分是定義因子數(shù)和變量數(shù), 以及文件位置:

$DATA  ASCII (4,2,0) d:/dmu-test/animal-model.txt

?
上面的意思是, 數(shù)據(jù)是ASCII格式, 有4個(gè)固子, 2個(gè)變量, 缺失值用0表示, 數(shù)據(jù)的絕對路徑是: d:/dmu-test/animal-model.txt
?

第四部分, 定義變量名, 也是為了方便結(jié)果輸出, 相當(dāng)于數(shù)據(jù)的行頭名

$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT

第五部分, 有6行, 定義模型整體來說是: 第一行: 單性狀 # 1 第二行: 1性狀無吸收 # 0 第三行: 1個(gè)性狀, 是由3個(gè)因子, 兩個(gè)固定因子:3,4, 一個(gè)隨機(jī)因子:1 # 1 0 3 3 4 1 第四行: 2個(gè)隨機(jī)因子, 他們的關(guān)系是獨(dú)立, 沒有協(xié)方差, 1 2 # 2 1 2

?
注意, 如果兩個(gè)隨機(jī)因子之間, 是有協(xié)方差的, 那么寫作 2 1 1, 即表示隨機(jī)因子有: 加性+ 母體 + 加性:母體協(xié)方差
?

第五行: 無回歸項(xiàng) # 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

5.5 執(zhí)行DIR文件

這里運(yùn)行的run_dmuai.bat, 將DMU安裝路徑下的文件run_dmuai.bat拷貝到d:/dmu-test文件夾, 在終端cmd界面鍵入:

run_dmuai.bat Uni_animalmodel-maternal

執(zhí)行結(jié)果:

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>

5.6 查看結(jié)果

在文件*lst中有估算的方差組分, 結(jié)果如下:

                            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 母體效應(yīng)方差組分為: 1.10404 殘差方差組分為: 1.6569

4.7 對比asreml的結(jié)果:

代碼:

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

5.8 DMU和asreml比較

兩者方差組分一致.

6. DMU linux下語法高亮設(shè)置

用vim編程時(shí), DMU的關(guān)鍵詞沒有語法高亮, 看著不舒服, 就進(jìn)行一下設(shè)置, 并記錄過程.

6.1 設(shè)置的效果如下

6.2 設(shè)置流程

本次設(shè)置的比較簡單, 將關(guān)鍵詞分為:

模型model, 比如DMU1, DMU2...
不同組成part, 比如DATA, VARIATE, MODEL...
不同結(jié)構(gòu)類型type, 比如PED, COR....
新建DIR.vim文件, 里面設(shè)置相關(guān)參數(shù)
新建DIR_suffix.vim文件, 設(shè)置后綴名讀取

6.3 DIR.vim文件:


"------------------------------------------------------------------------------"
" Description
"------------------------------------------------------------------------------"

" vim syntax highlighting file for DMU programs

" Author:   Deng Fei <dengfei_2013@163.com>
" Created:  Unknown
" Modified: 2018-11-18
" License:  GPLv2

syn keyword model DMU1 DMU4 DMU5 DMUAI RJMC

syn keyword part COMMENT ANALYSE DATA VARIABLE MODEL GLMM GLMM_PRED REDUCE MIXTURE VAR_STR VAR_REST PRECOND SOLUTION PRIOR RESIDUALS TRAITS ABSORB RANDOM REGRES NOCOV

syn keyword type PED DOM COR GRE PGMIX ABS_QTL GROUP VAR COV COR V_RATIO ASCII

hi model             ctermfg=Yellow
hi part              ctermfg=red
hi type              ctermfg=Green

將上面內(nèi)容, 保存為:DIR.vim文件, 放到:~/.vim/syntax文件夾中. 「如果沒有syntax文件夾, 就新建一個(gè).」

cp DIR.vim ~/.vim/syntax/

6.4 DIR_suffix.vim文件:

au BufRead,BufNewFile *.DIR            set filetype=DIR

將上面內(nèi)容保存到DIR_suffix.vim問價(jià)中, 放到:~/.vim/ftdetect文件夾中. 「如果沒有ftdetect文件夾, 就新建一個(gè).」

cp DIR_suffix.vim ~/.vim/ftdetect/

6.5 測試

使用下面代碼, 新建文件test.DIR, 然后使用vim打開, 查看語法高亮是否成功:

$ANALYSE 1 1 0 0
$DATA  ASCII (8,15,0) dat_dmu.txt
$VARIABLE
IDF1F2F3F4F5
y1
$MODEL
1
0
1 0 5 2 3 4 5 1
1
0
0
$VAR_STR 1 PED 2 ASCII ped_dmu.txt
$DMUAI
10
1d-7
1d-6
1
0

「效果如下:」

7. DMU在windows下安裝測試

7.1 下載地址

?
官網(wǎng)好像掛掉了，可以公眾號(hào)回復(fù)“DMU”獲取軟件和使用說明。
?

下載地址:http://dmu./

64為電腦安裝「DMUv6-R5-2-EM64T.msi」, 32為電腦安裝「DMUv6-R5-2.msi」

7.2 安裝DMU

7.3 測試DMU是否安裝成功

** 7.3.1 保存測試數(shù)據(jù)**

#dir.create("d:/ddmu-test")
setwd("d:/ddmu-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
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)

「7.3.2 查看數(shù)據(jù)」

「7.3.3 將下面內(nèi)容保存為:Uni_animal.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