软件学习-DMU

DMU软件使用方法。

官网

https://dmu.ghpc.au.dk/dmu/

模块介绍

• dmu1: 预处理模块
• dmuai : 估计方差组分模块
• dmu4 : BLUP 模块
• dmu5 ：通过数据迭代求解的 BLUP模块
• rjmc : 贝叶斯分析模块

输入数据

所有输入文件只能含有数值，不能含有字符（中文或英文字符）

表型数据

• 整数变量在前（如个体ID，固定因子），实数变量在后（性状，协变量等）
• 整数变量缺失值表示为 0 ，实数变量用一个极小值表示（如 -99，实际是等于小于参数卡中 miss 值的任意值）

系谱数据

• 4列，第4列是排序列（出生日期；如果忽略近交，最好使用行号，避免报错）
• 旧版本的DMU软件要求系谱要按照出生日期排序，现在版本的DMU软件不需要。

基因型数据

使用gmatrix软件构建G阵，这里我遇到报错，好像是 root 安装的 gcc ，g++ 版本太低，通过Conda环境解决GMatrix依赖问题

export LD_LIBRARY_PATH=$CONDA_PREFIX/lib:\$LD_LIBRARY_PATH

运行 gmatrix 如下，这里 renew_chipid.fam 中的个体号是已经重编码之后的，与上面的表型系谱文件保持一致。

gmatrix --bfile renew_chipid --grm agrm --out renew_chipid

生成了 renew_chipid.id 和 renew_chipid.agrm.id_fmt 。这里 renew_chipid.id 是一列基因型个体号， renew_chipid.agrm.id_fmt 为G阵内容, 为三列:个体号、个体号、值，示例如下

1 1 1.02875
2 1 0.333399
2 2 0.8651

一步法参数卡说明

生长参数卡-DMU4

DMU4 参数卡如下，这里

• $ANALYSE 11 9 0 0 : 表示是运行DMU4，方法是速度优先的 FSPAK算法（计算se），最后两个用不上。

• $DATA ASCII (3,2,-99) dmu.dat : 这里 3 是整型变量数目，2是实型变量数目，-99是缺失值，dmu.dat 是表型数据

• $VARIABLE ：每一列变量的名称缩写

• $MODEL :

  • 2 2 0 0 0

    几个性状就写 n n 0 0 0 即可（用一行5个整数描述分析性状的数目和性状类型，即1)总性
    状数, 2)正态分布性状数,3)正态截尾分布性状数, 4)多类阈性状数，5) 二类阈性状数），但是我看官方说明书，这里只需要写一个数字 n 即可

  • 0 0

    只有几个性状就放几个0的行，只用于 DMU5 模块，描述在分析过程中是否采用吸收法吸收一个固定效应。

  • 1 0 3 2 3 1 2 0 3 2 3 1

  用于指定每一个性状的模型，以 1 0 3 2 3 1 为例，逐个数字解释，1表示目标性状在实型数据中的序号，0表示权重变量在实型数据中的序号（0表示没有权重），3表示因子总数（固定+随机），之后便是具体的因子所在的列号，2表示第二列（固定因子），3表示第三列（窝效应），1表示第一列（加性效应）(官方文档上说，0表示截距，应该就是值为1的列)

  • 2 1 2 2 1 2

    用于指定每一个性状的随机因子组的设置。第一个数字 2 是这个性状随机因子总数，之后是每一个随机因子所在的随机因子组（从1开始编号，相同数字表示属于同一个因子组，例如母体效应模型）

  • 0 0

    用于指定每一个性状的协变量，第一个数字是协变量数目（0表示没有），如果存在协变量，则之后的数字就是指定的协变量在实型数据中的序号。

  • 0

    用于指定哪些残差协方差不存在（例如当两个性状是在不同的个体中获得的情况），第一行的数字表示不存在的残差协方差的数目，之后的行就是描述不存在的残差协方差的两个表型所在实型数据中的序号。举个例子假如第一个性状和第二个性状的残差协方差不存在，则第一行为1，第二行为1 2。

• $VAR_STR 2 PGMIX 2 ASCII dmu.ped renew_chipid.id renew_chipid.agrm.id_fmt 0.05 G-ADJUST

  这里指定随机效应的方差协方差结构。

  这里第一个数字 2 指定随机因子组的编号，这里我们加性效应是第2个随机因子，因此设置为2

  PGMIX 表示使用 H 阵，即一步法。

  2 表示构建 A逆的方式，2表示公畜-母畜系谱，不考虑近交

  ASCII 表示文本格式

  dmu.ped 系谱文件

  renew_chipid.id 基因型个体的ID

  renew_chipid.agrm.id_fmt G阵元素

  0.05 ：A阵权重

  G-ADJUST ：将G阵的对角线和非对角线元素均值校正到与A阵相同

  常规评估则使用 $VAR_STR 2 PED 2 ASCII dmu.ped 即可，将 PGMIX 改为 PED ，后面的内容删掉即可。

• $SOLUTION ：输出solution结果

• $RESIDUALS ASCII ：输出残差结果

• $PRIOR ：方差组分的值，对于 DMUAI 则为初始值

  以 1 1 1 11.281 为例，第一个数字1表示随机因子(组)编号，第二和第三个数字1表示这个随机因子组的方差组分的行号和列号，最后是方差组分的值（可以按照上三角或者下三角的形式描述）。

• $DMU4 : 可以用于对于某些解计算标准误，最多100个（每一行表示一个需要的解的序号）。

$COMMENT
This is a test example for estimating cow DRP variance component
Model drp=u+a+e

$ANALYSE 11 9 0 0

$DATA  ASCII (3,2,-99) dmu.dat
$VARIABLE
#1   2  3   
id  Gu wo
#1  2  
age bf  

$MODEL
2 2 0 0 0                                                                            
0 
0                                                                                    
#Trait Weight #effects u id
1  0   3  2 3 1
2  0   3  2 3 1 	
#Random  id  
2  1 2 
2  1 2                                                                                                  
#Regressions
0 
0
#NoCov
0                                                                             

$VAR_STR 2 PGMIX 2 ASCII dmu.ped  renew_chipid.id renew_chipid.agrm.id_fmt 0.05 G-ADJUST  

$SOLUTION

$RESIDUALS ASCII

$PRIOR
1  1  1    11.281    
1  2  1   0.87081
1  2  2   4.3223
2  1  1    17.450    
2  2  1    1.0489    
2  2  2    5.3171    
3  1  1    30.568    
3  2  1   0.89017    
3  2  2    14.207    

$DMU4

生长参数卡-DMUAI

这里只需要将上面的参数卡改2个地方

• $ANALYSE 1 1 0 0 ： 表示采用DMUAI，采用AI+EM混合的方法（应该是说如果AI迭代失败就用 EM算法的意思）

• $DMU4 改为 $DMUAI， 默认值如下

  $DMUAI
  10 ! Emstep = Number of steps before full weight on EM
  1.0d-7 ! Conv_ndelta = Conv. criteria for norm of the update vector
  1.0d-6 ! Conv_gnorm = Conv. criteria for norm of the gradient vector (AI)
  1 ! Printout = 1 -> Solution vector is printed/written to file SOL
  0 ! Fspopt = 0 -> time = 1 memory optimized FSPAK
  0 ! For development should always be zero (0)

一步法繁殖参数卡

这里只展示 DMU4 的参数卡，DMUAI的参数卡同上进行修改。

这里模型是2个固定因子（场年季和胎次），2个随机因子（永久环境效应和加性效应）

$COMMENT
This is a test example for estimating cow DRP variance component
Model drp=u+a+e

$ANALYSE 11 9 0 0

$DATA  ASCII (3,1,-99) dmu.dat
$VARIABLE
#1   2  3   
id  Gu taici
#1  
tnb  

$MODEL
1 1 0 0 0
0
#Trait Weight #effects id
1  0   4  2 3 1 1
#Random  id  
2  1 2 
#Regressions
0 
#NoCov
0                                                                             

$VAR_STR 2 PGMIX 2 ASCII dmu.ped  renew_chipid.id renew_chipid.agrm.id_fmt 0.05 G-ADJUST  

$SOLUTION

$RESIDUALS ASCII

$PRIOR
  1  1  1   0.95886    
  2  1  1   0.70218    
  3  1  1    10.566    

$DMU4

输出文件

结果文件 .SOL

每一列说明如下，正常来说我们只要提取第1，2，5，8，9 列，只需要第1列是4的行（加性效应）

列	说明	备注
1	效应类型：1表示用于回归的协变量；2表示固定因子；3表示除了加性效应外的随机因子；4表示加性效应；5表示DMU5中的吸收效应	正常分析只需要提取第一列是4的结果
2	性状	单性状模型则均为1
3	随机效应在协方差矩阵中的编号（0表示固定因子）	这应该是用于随机因子组的情况，如母体效应模型
4	因子编号	第几个固定因子/第几个随机因子(组)
5	水平原始编号	对应育种值就是个体号
6	这个水平的观测值数目
7	重编码的水平编号
8	估计值/预测值
9	SE

结果文件 .PAROUT

方差组分结果，示例如下

1  1  1    10.958817    
1  2  1   0.77859959    
1  2  2    4.1533056    
2  1  1    17.123556    
2  2  1   0.80002186    
2  2  2    5.7508502    
3  1  1    30.272560    
3  2  1    1.1302146    
3  2  2    13.807863    

GBLUP

基因型数据

此时这里我们需要使用 G逆 矩阵，命令如下

gmatrix --bfile renew_chipid --grm agrm --inv --out renew_chipid

输出 renew_chipid.id 为一列基因型个体号，renew_chipid.agiv.id_fmt 为G逆矩阵的元素，为三列:个体号、个体号、值。

G逆矩阵首行需要提供行列式的值（不清楚是什么行列式），实际操作只要首行添加0 0 1即可。不然估计方差组分时会报错。

# G逆文件首行添加 0 0 1 
sed -i '1i\0 0 1' renew_chipid.agiv.id_fmt

修改参数卡

方差组分结果修改如下， 这里 chip_id 是所有基因型个体的个体号文件，fmt 是 G逆 矩阵的元素（两个个体号，G逆元素，上/下三角元素）

$VAR_STR 2 GREL ASCII renew_chipid.agiv.id_fmt

这里第一个数字 2 指定随机因子组的编号，这里我们加性效应是第2个随机因子，因此设置为2

GREL 表示使用GBLUP。

ASCII 表示文本格式

renew_chipid.agiv.id_fmt 系谱文件

以生长性状举例，完整 DMU4 参数卡如下

$COMMENT
This is a test example for estimating cow DRP variance component
Model drp=u+a+e

$ANALYSE 11 9 0 0

$DATA  ASCII (3,2,-99) dmu.dat
$VARIABLE
#1   2  3   
id  Gu wo
#1  2  
age bf  

$MODEL
2 2 0 0 0                                                                            
0 
0                                                                                    
#Trait Weight #effects u id
1  0   3  2 3 1
2  0   3  2 3 1 	
#Random  id  
2  1 2 
2  1 2                                                                                                  
#Regressions
0 
0
#NoCov
0                                                                             

$VAR_STR 2 GREL ASCII renew_chipid.agiv.id_fmt

$SOLUTION

$RESIDUALS ASCII

$PRIOR
1  1  1    11.281    
1  2  1   0.87081
1  2  2   4.3223
2  1  1    17.450    
2  2  1    1.0489    
2  2  2    5.3171    
3  1  1    30.568    
3  2  1   0.89017    
3  2  2    14.207    

$DMU4

结果文件

注意这里结果中没有第一列是4的结果，第一列是3，并且第4列是加性效应编号（如2）的结果就是育种值结果。

结果举例如下

...
3  1  1   1         881           2         881       -0.821932             2.76516    
3  2  2   1         881           2         881        0.181913             1.74172    
3  1  1   1         882           2         882        0.528774             2.76009    
3  2  2   1         882           2         882         1.68531             1.73936    
3  1  1   1         883           1         883        -1.20906             3.00678    
3  2  2   1         883           1         883        0.931078E-01         1.88361    
3  1  1   2           1           1           1         4.23958             2.55934    
3  2  2   2           1           1           1         1.01104             1.53087    
3  1  1   2           2           1           2        -2.48123             2.50535    
3  2  2   2           2           1           2        -1.16639             1.49150    
3  1  1   2           3           1           3         1.01817             2.57584    
3  2  2   2           3           1           3         1.69694             1.53178    
3  1  1   2           4           1           4        -2.85359             2.62452    
3  2  2   2           4           1           4         1.28764             1.56361   
...