通过比对bim文件查找有问题基因型文件代码实现

当你用 plink 软件合并同一款芯片的多对plink文件时,出现了下述报错,提示有些位点不是二等位基因型。

1
Error: Variant '1_242598' is not biallelic.

这时候你就需要找到有问题的plink文件

前言

这种情况比较罕见,正常情况下不会发生,最可能的原因有两个

  1. 虽然是同一款芯片,但是可能存在更换了注释文件,导致同一个位点分型出来的两个碱基不同。
  2. 弄混了 Top 链 和 Forward 链的 plink 文件

使用软件

  • python 3.8 及以上版本
  • plink 软件 v1.90 版本

输入文件格式

格式说明见根据芯片号拆分plink文件的代码实现

通过 plink 软件的 --make-bed 命令生成的二进制基因型文件,具体格式说明见plink学习笔记。这里我们需要的就是后缀为.bim 的文件,这个文件在map文件的基础上新增了两列,为该SNP的两个碱基。

举个例子,下表为 bim 文件中的内容,第一行中的SNP的第一个碱基为 0 ,说明这个SNP均为一种分型(GG)。

1 CNCB10000545 0 119120907 0 G
1 CNC10012269 0 119133853 A G
1 CNC10012270 0 119144804 G C

比对多对plink文件代码

运行代码

代码文件见: cmp_bims.py

将所有需要检验的plink文本文件和本程序放在一个文件夹内,运行命令示范如下:

1
python cmp_bims.py in_path

参数说明

in_path :需要检验的多对plink文本文件所在的文件夹路径,注意一对ped文件和map文件的前缀必须相同。

代码说明

每对plink文件生成同名二进制文件

首先导入需要的模块

1
2
import sys,os
from collections import defaultdict

唯一的一个参数命名为 work_path,即工作路径

1
work_path = sys.argv[1]

创建一个报错提示逻辑值total_error_status ,初始化为 False 。如果值为 True ,表示存在报错。

1
total_error_status = False

切换目录

1
2
# 切换目录
os.chdir(work_path)

每对plink文件生成 bim 文件。这里是遍历输入路径下的所有文件,找到后缀为.ped 的文件,查找是否存在相应前缀的.map 文件,如果不存在则报错;之后使用 plink 软件生成二进制文件,将退出码赋给 exit_code ,如果退出码不为0,说明plink 软件生成二进制文件存在报错。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 每对plink文件生成 bim 文件
for i in os.listdir("."):
if os.path.isfile(i): #判断是文件
if i[-4:] == ".ped": #判断后缀是.ped
prefix = i[:-4] #去除.ped
if not os.path.exists(f"{prefix}.map"):
print(f"Error: have {prefix}.ped but not have {prefix}.map\n")
total_error_status = True
else:
exit_code = os.system(f"plink --allow-extra-chr --chr-set 95 --file {prefix} --out {prefix} &> /dev/null")
if exit_code != 0:
print(f"Error: plink file {prefix} has problems and can't go through plink software, please see {prefix}.log\n")
total_error_status = True

在上面一步中,只要存在报错,total_error_status 的值就会改为True 。如果有问题,程序中止。

1
2
if total_error_status:
sys.exit(1)

遍历所有bim文件,查找有问题的plink文件

这里首先是遍历输入文件的 bim 文件,将第一个bim文件所有SNP的碱基信息存入到dict1 字典中,键为SNP名称,值为该SNP的碱基的集合(包含 “0”)。

然后以第一个bim文件为模块,比对其他bim 文件。如果其他bim 文件的位点不再第一个bim 文件中,或者某个SNP的两个碱基与第一个bim文件不匹配,则打印报错信息,否则就打印***.bim is well ,说明与第一个bim文件相比较,这个bim文件没问题。

注意,由于这里第一个bim文件是随机选择的,因此如果恰好第一个bim文件就是有问题的,那么可能出现其他bim 文件均报错的情况。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#遍历所有的bim文件
j = 0
for i in os.listdir("."):
if (os.path.isfile(i)) and (i[-4:] == ".bim"): #判断是文件,后缀为bim
# prefix = i[:-4] #去除.bim
j += 1
if j == 1:
print(f"{i} is first bim") # 打印第一个文件
bim_file = open(i,'r')
for i2 in bim_file:
f2 = i2.split()
temp_alleles_set = dict1[f2[1]] # 该SNP的碱基集合
temp_alleles_set.add("0") # 先添加0
temp_alleles_set.add(f2[-2])
temp_alleles_set.add(f2[-1])
bim_file.close()
else:
# 以第一个bim为模板,比对其他bim文件。
bim_file = open(i,'r')
for i2 in bim_file:
f2 = i2.split()
if f2[1] in dict1:
temp_alleles_set = dict1[f2[1]].copy() # 第一个文件该SNP碱基集合的拷贝
temp_alleles_set.add(f2[-2])
temp_alleles_set.add(f2[-1])
if len(temp_alleles_set) > 3: # 如果超过了两个碱基
print(f"Error: SNP {f2[1]} in {i} is not consistent with first bim\n") #打印错误文件的错误碱基
break
else:
print(f"Error: SNP {f2[1]} in {i} not in first bim\n")
break
else:
print(f"{i} is well") #没有问题则打印没问题
bim_file.close()

比对两个bim文件代码

前言

假如你在上一步找到有问题的plink文件,你先看看这个plink文件与所谓的正常的plink文件到底是哪些位点的等位基因不一样,你就可以接着运行这个脚本。

或者你总共就只有两对plink文件,合并报错,你也可以直接运行这个脚本。

运行代码

代码文件见: cmp_two_bims.py

将需要检验的bim文件和本程序放在一个文件夹内,运行命令示范如下:

1
python cmp_two_bims.py 1.bim 2.bim cmp_two_bims.txt

参数说明

1.bim: 第一个bim文件名称

2.bim: 第二个bim文件名称

cmp_two_bims.txt : 结果文件名称

注意:这里直接用的是bim文件,如果没有生成就先用plink软件生成一下。

输出文件格式

输出文件为不一致位点的碱基信息,总共5列,第一列为SNP名称,第二列和第三列为该SNP在第一个bim文件的两个碱基,第四列和第五列为该SNP在第二个bim文件的两个碱基。示例如下

CNC10012269 A G A C
CNC10012270 G C A T

代码说明

读取文件

1
2
3
4
import sys
bim1 = open(sys.argv[1],'r')
bim2 = open(sys.argv[2],'r')
out_file = open(sys.argv[3],'w')

输出文件写入标题

1
out_file.write("\t".join(["snp","bim1_A1","bim1_A2","bim2_A1","bim2_A2"])+'\n')

遍历第一个bim文件

将第一个bim文件SNP的两个碱基信息写入一个字典中

1
2
3
4
dick = {}
for i in bim1:
f = i.split()
dick[f[1]] = f[-2:] #snp名称:倒数第二列

遍历第二个bim文件,比对SNP

首先判断第二个bim文件的SNP是否均在第一个bim文件中,如果不是则报错中止程序。

之后检查每个SNP在两个bim文件中的碱基数目是否超过 2 (加上 “0” 超过3) ,超过了就说明有问题,将该SNP数据的信息写入输出文件,不一致位点数 error_base_num 加1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
error_base_num = 0
for i in bim2:
f = i.split()
if f[1] not in dick:
print("Error: {} in bim2 not in bim1".format(f[1]))
sys.exit()
else:
bim1_alleles_set = set(dick[f[1]])
bim1_alleles_set.add("0")
bim1_alleles_set.add(f[-2])
bim1_alleles_set.add(f[-1])
if len(bim1_alleles_set) > 3: #如果包括0,两个bim合并的等位基因数目大于3,就有问题
error_base_num += 1
bim1_alleles = '\t'.join(dick[f[1]])
out_file.write(f"{f[1]}\t{bim1_alleles}\t{f[-2]}\t{f[-1]}\n")
# print("Error: {} in bim2 file has more than 2 alleles".format(f[1]))
# sys.exit()

打印提示信息

如果不一致位点数为0,那么打印祝贺信息;反之提示有几个不一致位点。

1
2
3
4
if error_base_num == 0:
print("Congradulation: two bim fils is consistent")
else:
print(f"Error: {error_base_num} bases has more than 2 alleles")

关闭文件

用完记得关闭文件,这是个好习惯

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

请我喝杯茶吧~

支付宝
微信