HIBLUP教程|第3期:不同类型个体间关系矩阵的构建及格式转换

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 215,794评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,050评论 3 391
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 161,587评论 0 351
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,861评论 1 290
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,901评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,898评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,832评论 3 416
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,617评论 0 271
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,077评论 1 308
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,349评论 2 331
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,483评论 1 345
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,199评论 5 341
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,824评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,442评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,632评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,474评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,393评论 2 352

推荐阅读更多精彩内容