统计SNP基因型频率及MAF代码实现

最近准备一些材料时需要计算每个SNP三种分型的基因型频率和MAF,plink 软件只给出 MAF,因此写了一个简单的代码实现它。

使用软件

  • python 3.8 及以上版本,事先安装好 Numpy 模块

背景知识

假设一个SNP,它有两种碱基 A 和 T,那么就有三种可能的基因型 AA, AT, TT。

假设某个群体样本数目为 100 ,在这个SNP 中有 5 个样本分型失败,有 20 个样本分型为 AA, 有 30 个样本分型为 AT,有45 个样本分型为 TT,计算三种分型的基因型频率如下:

注意,分母计算总数时需要剔除分型失败的样本,三种基因型频率相加为1。

现在再看MAF计算方式,MAF 是最小等位基因频率的缩写。因此,我们可以将 P(A) 和 P(T) 都计算出来,看看哪个数值更小,那就是MAF

从上面三种基因型频率结果可以看出,A 是较小的等位基因,而 T 是较大的等位基因。判断理由很简单,因为 P(AT) 不影响 P(A) 和 P(T) 的相对大小,只要看 P(AA) 和 P(TT) 哪个更小,哪个就是最小等位基因了。

输入文件

使用 plink 软件的 recodeA 命令生成的 raw 文件。这里每个SNP占据一列,值为最小等位基因的数目(0 1 2),缺失值记为 NA。

输出文件

输出文件共5列,第一列为SNP名称,第二列为较小等位基因纯和子的基因型频率,第三列为杂合子基因型频率,第四列为另一类纯和子基因型频率,第五列为MAF。

seq-rs80960919 0.119152 0.497401 0.383447 0.367853
seq-rs81350349 0.195922 0.530588 0.273491 0.461216
seq-rs81349911 0.2088 0.5384 0.2528 0.478

运行代码

代码文件见:genotype_frequency.py

将输入文件和本程序放在同一文件夹下,运行命令示范如下。

1
python genotype_frequency.py A.raw A_frequency.txt

参数说明:

A.raw: 输入文件名称

A_frequency.txt: 结果文件名称

运行结束后,生成以第二个参数为名称的结果文件

代码说明

替换分型缺失值

首先我用 linux 的 sed 命令来精准查找 NA,替换为 3,生成中间文件 temp.raw。

1
2
3
4
5
6
import os,sys
import numpy as np

raw_file_name = sys.argv[1]

os.system(f"sed 's/\<NA\>/3/g' {raw_file_name} > temp.raw") # 改为精准匹配

读取raw文件

首先从标题中提取 SNP 名称

1
2
3
4
5
raw_file = open("temp.raw",'r')
raw_file_title = raw_file.readline().split() #剔除标题
title_list = []
for i in raw_file_title[6:]:
title_list.append(i[:-2]) # snp 名称列表

之后将基因型信息存在到一个二维数组 array1 中,每一行表示一个样本,每一列表示一个SNP

1
2
3
4
5
6
7
8
raw_list=[]
for i in raw_file:
f = i.split()
raw_list.append(f[6:])

raw_file.close()

array1 = np.array(raw_list, dtype=np.int8)

统计基因型频率和MAF

首先统计每列SNP中 2,1, 0 三种分型的数目,各自除以三种分型总数,便得到了三种基因型频率。

1
2
3
4
5
6
7
8
9
10
# 统计三种基因型数据

p_array = (array1==2).sum(axis=0) # 按列求和, 较小等位基因的纯合子
h_array = (array1==1).sum(axis=0) # 按列求和
r_array = (array1==0).sum(axis=0) # 按列求和

p_rate = p_array/(p_array+h_array+r_array)
h_rate = h_array/(p_array+h_array+r_array)
r_rate = r_array/(p_array+h_array+r_array)

根据基因型频率,计算MAF

1
maf_rate = (2*p_rate + h_rate)/2

写入结果文件

1
2
3
4
5
6
7
8
# 写入结果文件
out_file = open(sys.argv[2],'w')

snp_num = len(title_list)
for i in range(snp_num):
out_file.write(title_list[i]+"\t"+str(p_rate[i])+"\t"+str(h_rate[i])+"\t"+str(r_rate[i])+"\t"+str(maf_rate[i])+"\n")

out_file.close()

最后,剔除中间文件

1
os.system("rm temp.raw")

小结

这个代码中计算的MAF 和 Plink 软件计算的 MAF 比对过,是一样的。

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2022 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信