1、不同类型个体间关系矩阵的构建
个体关系矩阵(Relationship Matrix)用于构建遗传随机效应各变量间的方差协方差,是线性混合模型的必要元素之一。HIBLUP可以使用系谱信息和基因型信息构建多种个体关系矩阵及其逆矩阵,如使用系谱信息构建的系谱加性关系矩阵(PA矩阵)和显性关系矩阵(PD矩阵)、使用基因型信息构建的基因组加性关系矩阵(GA矩阵)和显性关系矩阵(GD矩阵)以及同时使用系谱和基因型信息构建的加性关系矩阵(HA矩阵)和显性关系矩阵(HD矩阵)。以上使用不同信息构建的个体关系矩阵可以用于拟合单性状和多性状PBLUP、GBLUP、SSBLUP模型,计算个体的估计育种值。
在HIBLUP中,通过--make-xrm参数构建并输出个体关系矩阵。默认情况下,程序包含参数--add,仅输出加性关系矩阵(A矩阵),可以通过添加--dom构建显性关系矩阵(D矩阵),添加--add-inv(或--addinv)和--dom-inv(或--dominv)分别构建A矩阵和D矩阵的逆矩阵。个体关系矩阵储存在以.bin为后缀的二进制格式文件中,用户无法通过软件直接打开查看,因此HIBLUP提供了额外的参数--write-txt,命令行中添加该参数后,HIBLUP会额外输出方阵形式的文本格式关系矩阵,也可以在此基础上添加参数--triangle输出三联表形式的文本格式关系矩阵,可直接通过文本编辑器打开。
1.1 利用系谱构建个体关系矩阵(PRM)
通过参数--pedigree输入系谱信息即可计算系谱关系矩阵。如同时构建系谱加性关系矩阵及其逆矩阵的命令如下:
hiblup --make-xrm --pedigree demo.ped --add --add-inv --out demo
--pedigree:指定输入的系谱文件,HIBLUP接受的系谱文件的格式如demo.ped所示:
其中包含三列信息,第一列为个体ID,第二列和第三列分别是个体父亲和母亲的ID。第一列中个体ID可以按照出生日期排序,但这并不是强制要求的,HIBLUP能够对输入的系谱信息自动进行重排。注意:系谱文件中出现以下内容均被认为是缺失值:
NA Na . - NaN NAN nan na N/A n/a <NA>
程序运行效果如下:
程序运行后生成了demo.PA.bin、demo.PA.id、demo.PAinv.sparse.bin、demo.PAinv.sparse.id、demo.log共5个文件。
这里所构建的A矩阵是稠密矩阵,矩阵中大部分元素值不为零,矩阵中所有的下三角元素储存在后缀为.bin的文件中;而A矩阵的逆矩阵是稀疏矩阵,矩阵中只有少数元素不为零,如果仍然按照储存稠密矩阵的方式存储,会因储存大量零元素而浪费很多空间,同时在进行运算时也增加了许多无效计算。因此,HIBLUP仅将稀疏矩阵中所有非零元素以三联表的形式(即行索引、列索引、对应值)储存在后缀为.sparse.bin的文件中;个体ID储存在后缀为.id的文件中。
参数--add、--add-inv、--dom、--dom-inv可同时使用,也可以单独使用,生成对应的.id文件和.bin文件。
在命令行中添加参数--write-txt可以生成文本格式的关系矩阵,便于用户直接查看。命令行输入:
hiblup --make-xrm --pedigree demo.ped --add --write-txt --out demo
程序运行后,在用户的目录下会额外生成demo.txt文件,其格式如下:
部分软件接受的关系矩阵为三联表形式(如ASReml、DMU等),HIBLUP同样可通过参数--triangle实现,如:
hiblup --make-xrm --pedigree demo.ped --add-inv --write-txt --triangle --out demo
生成的三联表形式的A矩阵的逆矩阵文件demo.txt,其格式如下:
1.2构建基因组个体关系矩阵(GRM)
与构建系谱关系矩阵相似,构建基因组关系矩阵时只需要将参数--pedigree替换为--bfile并输入基因组文件即可。如同时构建基因组加性及显性关系矩阵的命令行输入:
hiblup --make-xrm --bfile demo --add --dom --out demo
--bfile:指定输入的二进制格式的基因型文件的文件名,不包含文件名后缀。HIBLUP接受的基因组文件格式可见第二期内容(
https://mp.weixin.qq.com/s/UT-w3DoPfg1a4_ybHz81JQ):
程序运行效果如下:
程序运行后生成了demo.GA.bin、demo.GA.id、demo.GD.bin、demo.GD.id、demo.log共5个文件。所构建的稠密矩阵的所有下三角元素储存在后缀为.bin的文件中,个体ID储存在后缀为.id的文件中。
HIBLUP构建基因组关系矩阵的公式如下:
其中,是数值化的基因型矩阵,维数为,为个体数,为标记数,是矩阵的迹。在HIBLUP中,共实现了4种方法获得矩阵,不同方法的基因型编码方式可参考下表:
其中和分别为A和a的等位基因频率。、和分别是AA、Aa和aa的基因型频率。
用户可赋值不同的编号给参数--code-method 自由切换对应的方法,HIBLUP默认使用--code-method 1。
不同方法的参考文献详见:
--code-method 1, (default) Su, Guosheng, et al. “Estimating additive and non-additive genetic variances and predicting genetic merits using genome-wide dense single nucleotide polymorphism markers.” Plos one (2012): e45293.
--code-method 2, Zeng, Zhao-Bang, Tao Wang, and Wei Zou. “Modeling quantitative trait loci and interpretation of models.” Genetics 169.3 (2005): 1711-1725.
--code-method 3, Yang, Jian, et al. “Common SNPs explain a large proportion of the heritability for human height.” Nature genetics 42.7 (2010): 565-569.
--code-method 4, Vitezica, Zulma G., et al. “Orthogonal estimates of variances for additive, dominance, and epistatic effects in populations.” Genetics 206.3 (2017): 1297-1307.
HIBLUP还可以通过参数--snp-weight输入包含不同SNP位点所占权重的文件来计算加权G矩阵。权重文件如demo.weight.txt所示:
其中,第一列为SNP位点的名称,第二列为该位点所占的权重,权重值均为正数。demo.weight.txt提供的是SNP位点的加性权重,可以用于构建加权基因组A矩阵。命令行输入如下:
hiblup --make-xrm --bfile demo --add --snp-weight snp.weight.txt --out demo_weight
需要注意的是,在构建加权基因组关系矩阵时,提供的权重值应与需要构建的矩阵类型相对应,即计算加权基因组A矩阵时,要提供SNP的加性权重;计算加权基因组D矩阵时,要提供SNP的显性权重;要同时计算A矩阵和D矩阵时,权重文件中应包含加性权重列和显性权重列,且加性权重在前显性权重在后。
在第二期中还提到,HIBLUP提供了进行SNP位点筛选的参数--extract和--exclude,以及个体筛选的参数--keep和--remove,在构建G矩阵时同样可以通过使用这些参数对标记和个体进行筛选。如:
hiblup --make-xrm --bfile demo --extract snp.filter.txt --add --out demo
在构建关系矩阵时添加--extract参数即可仅利用文件中包含的标记构建G矩阵,文件格式如snp.filter.txt所示,为一行一个的SNP位点;与之相反,--exclude则是将文件中包含的SNP位点删除,利用剩余的位点构建G矩阵,命令行输入如下:
hiblup --make-xrm --bfile demo --exclude snp.filter.txt --add --out demo
同理,使用个体筛选的参数--keep和--remove在构建G矩阵时可以过滤掉不需要的个体。命令行输入如下:
hiblup --make-xrm --bfile demo --keep id.filter.txt --add --out demo
在构建关系矩阵时添加--keep参数并输入包含个体ID的文件,即可利用文件中包含的个体构建G矩阵,文件格式如id.filter.txt所示,为一行一个的个体ID。
与--keep相反,--remove则是将文件中包含的个体删除,利用剩余的个体构建G矩阵。命令行输入如下:
hiblup --make-xrm --bfile demo --remove id.filter.txt --add --out demo
在构建GRM过程中,用户可通过参数--step控制计算过程中的内存消耗,默认为--step 10000,数值越大消耗内存越大,反之越小。
1.3构建SSGBLUP模型的H矩阵(HRM)
当同时赋值--bfile和--pedigree参数时,HIBLUP自动利用提供的基因型文件和系谱文件构建SSGBLUP模型的H矩阵,命令行输入:
hiblup --make-xrm --bfile demo --pedigree demo.ped --alpha 0.05 --add --add-inv --out demo
--alpha:指定构建H矩阵时系谱信息所占的权重,默认为0.05;
程序运行效果如下:
程序运行后生成了demo.HA.bin、demo.HA.id、demo.HAinv.sparse.bin、demo.HAinv.sparse.id、demo.log共5个文件。
这里构建的HA矩阵是稠密矩阵,矩阵中所有下三角元素储存在后缀为.bin的文件中;HA矩阵的逆矩阵是稀疏矩阵,行索引、列索引及对应值储存在后缀为.sparse.bin的文件中;个体ID储存在后缀为.id的文件中。
HIBLUP中,构建H矩阵的公式如下:
和分别是从由系谱构建的A矩阵中代表非基因型的个体和基因型个体的子矩阵,和是A矩阵中它们代表基因型和非基因型个体之间关系的子矩阵。由公式得到,其中是原始基因组关系矩阵G矩阵通过以下简单线性回归方程校正得到的:
调整G矩阵对角线的平均值与对角线平均值一致,G矩阵非对角线的平均值与非对角线平均值一致,计算出和的值,从而得到。
1.4构建环境随机效应矩阵
除了构建遗传随机效应矩阵以外,HIBLUP还能利用表型文件中的随机效应记录构建环境随机效应矩阵,命令行输入:
hiblup --make-xrm --pheno demo.phe --rand 6,7 --out demo
--pheno:指定输入的表型文件的路径和文件名,需要包含表头。HIBLUP的表型文件的第一列必须为个体ID,除此之外还可以包括表型记录,还应包括协变量、固定效应和随机效应,如demo.phe所示:
demo.phe中第一列为个体ID,第二列和第三列为固定效应性别和季节,第四列和第五列为协变量日龄和初生重,第六列和第七列为随机效应测定场和母本ID,第八列到第十列为三种性状的表型值。
--rand:随机效应在表型文件中列的位置,多个随机效应间用逗号隔开;
程序运行效果如下:
运行完成后每一个随机效应对应的环境随机效应矩阵都储存在.bin和.id文件中,同时生成日志文件。比如当使用测定场和母本作为随机效应构建环境随机效应矩阵时,生成demo.location.sparse.id、demo.location.sparse.bin、demo.dam.sparse.id、demo.dam.sparse.bin和demo.log文件。
2、关系矩阵的格式转换
2.1 二进制格式转换为文本格式
HIBLUP可将上述生成的二进制格式基因组关系矩阵转换为文本文件格式,也可以对GCTA、LDAK等软件输出的二进制格式基因组关系矩阵进行转换。新建一个包含基因型文件中6个个体ID的id.list.txt文件,构建其中个体的G矩阵,并将转换为文本格式,命令行输入:
1. #对HIBLUP生成的关系矩阵文件进行格式转换
2. hiblup --trans-xrm --xrm demo.GA --keep id.list.txt --out demo.GA
3. #对GCTA输出的关系矩阵文件进行格式转换
4. hiblup --trans-xrm --xrm GCTA.grm --keep id.list.txt --out demo.GA
--xrm:指定二进制文件的前缀,不包含文件名后缀.bin和.id;
程序运行效果如下:
从运行结果可以看到,HIBLUP将二进制文件demo.GA.bin中的矩阵数据转换成了文本格式的方阵并储存在demo.GA.txt中,6个个体的ID也储存到了demo.GA.id.txt文件中,同时生成日志文件demo.log。demo.GA.txt的格式如下图所示,矩阵的行和列与demo.GA.id.txt中的ID相对应。
另外,还可以通过添加--triangle参数将矩阵的下三角元素的位置和数值按行存储为三联表的形式,命令行输入:
hiblup --trans-xrm --xrm demo.GA --keep id.list.txt --triangle --out demo.GA
程序运行效果如下:
程序运行后生成demo.GA.txt、demo.GA.id.txt和demo.log文件。demo.GA.txt文件格式如下所示,其中第一列为行索引,第二列为列索引,第三列为对应值。对于稀疏矩阵,则省略零元素,只输出非零元素的行、列索引和元素值。
2.2文本格式转换为二进制格式
由于使用HIBLUP进行计算时需要输入二进制格式的文件,不接受文本格式文件,而有些时候我们只拥有文本格式的关系矩阵的数据,这种情况下可以使用--xrm-txt和--xrm-id参数分别输入文本格式的关系矩阵和个体ID,将其转换为二进制格式,命令行输入:
hiblup --trans-xrm --xrm-txt demo.GA.txt --xrm-id demo.GA.id.txt --out demo.GA
--xrm-txt:输入文本格式的关系矩阵文件;
--xrm-id:输入包含关系矩阵中所有个体ID的文件;
程序运行效果如下:
程序运行生成demo.GA.bin、demo.GA.id和demo.log,完成关系矩阵从文本格式到二进制格式的转换。
需要注意的是,如果输入的需要转换为二进制格式的矩阵文本文件为三联表的格式,则还需要添加参数--triangle将其转换为二进制格式,命令行输入如下:
hiblup --trans-xrm --xrm-txt demo.GA.txt --xrm-id demo.GA.id.txt --triangle –out demo.GA