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

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

Error: Variant '1_242598' is not biallelic.

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

前言

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

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

使用软件

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

输入文件格式

plink 文本文件

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

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文本文件和本程序放在一个文件夹内，运行命令示范如下：

python cmp_bims.py in_path

参数说明：

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

代码说明

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

首先导入需要的模块

import sys,os
from collections import defaultdict

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

work_path = sys.argv[1]

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

total_error_status = False

切换目录

# 切换目录
os.chdir(work_path) 

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

# 每对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 。如果有问题，程序中止。

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 文件均报错的情况。

#遍历所有的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文件和本程序放在一个文件夹内，运行命令示范如下：

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

代码说明

读取文件

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

输出文件写入标题

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

遍历第一个bim文件

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

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。

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，那么打印祝贺信息；反之提示有几个不一致位点。

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")

关闭文件

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

bim1.close()
bim2.close()
out_file.close()