GCTA计算亲缘关系矩阵(两种不同方法)

亲缘关系矩阵

        亲缘关系矩阵被广泛用于全基因组关联分析(GWAS)和全基因组选择(GS)中,当然主要是混合线性模型。在GWAS中亲缘关系矩阵主要用于校正遗传背景,因为GWAS一般采用单点扫描的模型,这种模型极易形成假阳性。而遗传背景的本质是除当前扫描标记外,其他所有标记的效应(具体的我们后面再谈)。亲缘关系矩阵可分为A阵和G阵,A阵是通过系谱进行计算所得,G阵是通过基因标记计算所得。G阵又被称作realized relationship matrix(现实关系矩阵)。为什么G阵更realized呢?因为根据系谱推导除的A阵是无法区分全同胞的,全同胞的关系是完全一样的。但是我们知道即使是同卵双胞胎也会有不同,而通过基因标记就可以找到他们之间的不同的SNP。因此通过基因标记计算的亲缘关系矩阵更真实。

        关于G矩阵被应用最多的是Vanraden于2008在奶业科学上发表的文章。该文论述了几种G矩阵的计算方法,但目前主流用的是前两种。两种方法的区别在于校正标记期望方差的方法不同(具体我下次写一篇文章专门说公式推导)。`Vanraden`自己的模拟结果是第一种方法最好(使用所有标记期望方差的均值来校正),但是杨剑老师主推的是第二种方法(每个标记校正自己的期望方差)。目前在GS中主要用的是第一种方法来计算kinship。其实两种方法差别并不大,因为标记本来就是0,1,2三种分型,方差的区别也大不到哪去。

        以下正式讲解GCTA计算G矩阵的方法。


        今天我们来讲讲如何利用GCTA来计算亲缘关系矩阵。目前GCTA已经支持windows了,但作者仍建议我们在linux下使用。安装CCTA比较简单,直接下载解压,运行可执行文件即可。最好软链接到系统环境目录下(如~/bin),就可在任意目录下输入gcta64运行。

        进入GCTA官网,可以看到使用 ```--make-grm``` 就可以计算G矩阵:

    使用```--make-grm-alg``` 可以选择不同的算法,因为杨剑老师比较认同第二种算法,所以默认是第二种算法啦。

这里我们使用最常用的第一种G矩阵计算方法(设置make-grm-alg 1即可)。想用第二种算法的可以使用make-grm-alg 0。完整代码如下:

gcta --bfile test --make-grm --make-grm-alg 1 --out kinship

bfile参数用于指定待运行的plink二进制文件,这里用的数据是GCTA自带的测试数据。运行结果如下:


        运行后,会得到如下四个文件:

        1. test.grm.bin  含G阵下三角元素,是二进制文件

        2. test.grm.N.bin 记录计算G阵的SNP个数,是二进制文件

        3. test.grm.id 记录个体的family号和id号,即plink fam文件的前两列

        4. kinship.log 日志文件。

使用以上代码生成的是二进制文件的亲缘关系矩阵,不利于进行其他分析。如果仅仅为了计算亲缘关系矩阵,并需要用其他软件进一步分析,推荐使用以下代码:

gcta64 --bfile test --make-grm-gz --make-grm-alg 1 --out kinship

运行后,会得到如下四个文件:

1. kinship.grm.gz  个体间亲缘关系的gz压缩文件,解压后可获得四列数据结果,分别是 id1, id2, 用于计算id1和id2亲缘关系的SNP总数, id1和id2的亲缘关系系数。

2. kinship.grm.id 同上

3.kinship.log 同上

推荐阅读更多精彩内容