统计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

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

python genotype_frequency.py A.raw A_frequency.txt

参数说明：

A.raw: 输入文件名称

A_frequency.txt: 结果文件名称

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

代码说明

替换分型缺失值

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

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 名称

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

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 三种分型的数目，各自除以三种分型总数，便得到了三种基因型频率。

# 统计三种基因型数据

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

maf_rate = (2*p_rate + h_rate)/2

写入结果文件

# 写入结果文件
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()

最后，剔除中间文件

os.system("rm temp.raw")

小结

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