【原】基因型数据012及-1,0,1计算基因频率

G矩阵计算时, 有不同的编码形式, 我们看一下计算的区别.

第一种方法按照012, 即AA是0, 表示major基因, 1 表示杂合, 2表示aa(minor).

第二种方法按照-1, 0, 1, 即-1是AA, 表示major基因型, 0表示杂合, 1表示aa(minor).

下面我们模拟一下数据, 看一下两者计算的区别.

数据

dat = data.frame(snpid = paste0("snp",1:7),type = c("GG","GG","AG","GG","AG","GG","AA"))dat

snpid	type
snp1	GG
snp2	GG
snp3	AG
snp4	GG
snp5	AG
snp6	GG
snp7	AA

上面是7个SNP的基因型, 可以看出GG是major, AA为minor.

1, 手动计算

• GG个数: 4

• AG个数: 2

• AA个数: 1

计算方法, 所有G的个数, 除以总配子的个数, 即为G的基因频率

p_g = (4*2+2*1)/(7*2);p_g

0.714285714285714

p_a = (2*1+1*2)/(7*2);p_a

0.285714285714286

2, 转化为012

如何计算基因频率(major为0, minor为2)

dat$type012 = 0 dat[dat$type=="GG",]$type012 =0 dat[dat$type %in% c("AG","GA"),]$type012 =1 dat[dat$type=="AA",]$type012 =2dat

snpid	type	type012
snp1	GG	0
snp2	GG	0
snp3	AG	1
snp4	GG	0
snp5	AG	1
snp6	GG	0
snp7	AA	2

计算方法:

将type012所有列之和相加, 为minor的基因个数, 除以2*n, 即为minor的基因频率

sum(dat$type012)/(length(dat$type012)*2)

0.285714285714286

3, 转化为-1,0,1

计算方法(-1为major, 0为杂合,1为minor)

dat$type_101 = dat$type012 -1dat

snpid	type	type012	type_101
snp1	GG	0	-1
snp2	GG	0	-1
snp3	AG	1	0
snp4	GG	0	-1
snp5	AG	1	0
snp6	GG	0	-1
snp7	AA	2	1

计算方法:

分子为: 列之和 + 行数

分母为: 行数 * 2

(sum(dat$type_101) + length(dat$type_101))/(2*length(dat$type_101))

0.285714285714286

结论:

1, 如果转化为0, 1, 2, 那么他们之和即为minor的个数, 除以2n即为minor的频率. 因为major的值为0, 1为杂合, 同时表示有一个minor, 2为纯合, 表示2个minor, 计算比较方便.

2, 如果基因型转化为-1, 0, 1的话, 那么将其进行求和, -1和1抵消. 推断过程如下:

假设A的频率为(1-q), a的频率为q, n为个数.

AA的个数为x

Aa的个数为y

aa的个数为z

x + y + z =n

q = (2*z + y)/(2*n)

-1*x +0*y + 1*z = z-x

2*z +y = 2*z + (n - x -z) = n + z - x #分子

因此:

q = (n + z -x)/(2*n)

这种方法, 计算量比较小, 因为-1 和1 进行了抵消.

这种方法不太直观, 需要用笔画一下即可理解.

所以, 我们计算G矩阵时, 将SNP分型转化为-1, 0, 1的形式.

如何利用系谱计算近交系数和亲缘关系系数

如何构建G矩阵-基因组亲缘关系矩阵（Genomic relationships matrix）

育种中一般配合力和特殊配合力的计算方法