Python通用VCF文件读取与处理类实现

前言

在处理VCF文件时，往往涉及多种切割和正则匹配操作，手动编写代码既繁琐又容易出错，且每次文件格式或读取需求变化时都需修改代码。虽然市面上有许多Python VCF库，但用户需先安装这些库，且部分库固定了可读取的字段，若所需信息（如位于样本GT列而非INFO列的AF值）不在固定范围内，则无法读取。为此，本文提供了一个通用的VCF读取类，用户只需复制粘贴代码即可方便地进行VCF文件的读取、过滤和写入操作。

使用说明

首先复制类的代码，后续即可直接使用。

import sys
import osimport subprocess
class Record(object):
 '''
 One line information in vcf file
 '''
 def __init__(self, line):
 info = line.split("\t")
 self.line = line
 self.CHROM = info[0] 
 self.POS = info[1]
 self.ID = info[2]
 self.REF = info[3]
 self.ALT = info[4]
 self.QUAL = info[5]
 self.FILTER = info[6]
 self.INFO = [{pair_lst[0]: pair_lst[1] if len(pair_lst)> 1 else ""} for pair_lst in [pair.split("=") for pair in info[7].split(";")]]
 self.FORMAT = info[8].split(":")
 self.sample_num = len(info) -7
 self.GT = []
 for i in range(2):
 GT_value = info[8 + i +1].split(":") 
 GT_dict = {}
 for g in range(len(GT_value)):
 GT_dict[self.FORMAT[g]] = GT_value[g] 
 self.GT.append(GT_dict) 
class VCF(object):
 '''
 VCF class, read VCF, write VCF, get VCF information
 '''
 def __init__(self, uncompress_vcf):
 self.header = []
 self.reader = open(uncompress_vcf, 'r')
 self.line = self.reader.readline().strip()
 while self.line.startswith('#'):
 self.header.append(self.line)
 self.line = self.reader.readline().strip()
 self.record = Record(self.line) 
 def __iter__(self): 
 return self 
 def __next__(self): 
 self.line = self.reader.readline().strip()
 if self.line != "":
 self.record = Record(self.line) 
 return self.record
 else:
 self.reader.close()
 raise StopIteration()
 def reader_close(self):
 self.reader.close() 

• 该类主要包含两个部分：一个是VCF类，负责存储VCF信息并对VCF文件进行操作；另一个是Record类，用于存储VCF文件中某一行的全部信息。
• 读入VCF文件

gatk_result = "realignment.vcf"
gatk = VCF(gatk_result)

查看VCF的Header

gatk.header

查看VCF当前行中存储的信息，初始状态为首行。信息以Record类对象保存。需要注意的是，VCF类是一个迭代器类，可以使用next函数和for循环来逐行读取信息。

record = gatk.record #这里record存储的是该Record类的地址

查看该Record对象的属性，包括line（行的内容，便于写出某行）、CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO（以字典形式存储）、FORMAT、sample_num（样本数量）以及GT（样本的基因型信息，通常在VCF中位于末尾用样本名表示的列）。

record.CHROM
record.line
record.ID #其他的属性同理

INFO字段的读取

以下是VCF中INFO字段的原始表示形式：

CONTQ=28;DP=38;ECNT=1;GERMQ=76;MBQ=20,37;MFRL=171,229;MMQ=60,60;MPOS=26;NALOD=1.16;NLOD=3.91; POPAF=6.00;RCNTS=0,0;ROQ=14;SEQQ=1;STRANDQ=11;TLOD=4.56

该信息在Record对象中的存储形式如下：

record.INFO [{'CONTQ': '28'}, {'DP': '38'}, {'ECNT': '1'}, {'GERMQ': '76'}, {'MBQ': '20,37'}, {'MFRL': '171,229'}, {'MMQ': '60,60'}, {'MPOS': '26'}, {'NALOD': '1.16'}, {'NLOD': '3.91'}, {'POPAF': '6.00'}, {'RCNTS': '0,0'}, {'ROQ': '14'}, {'SEQQ': '1'}, {'STRANDQ': '11'}, {'TLOD': '4.56'}]

GT字段的读取

以下是VCF中GT字段的存储形式，FORMAT字段对应着GT的值：

GT:AD:AF:DP:F1R2:F2R1:OBAM:OBAMRC:OBF:OBP:OBQ:OBQRC:SB 0/1:21,2:0.120:23:7,1:13,1:false:true:0.500:0.078:100.00:41.80:12,9,1,1 0/0:13,0:0.065:13:7,0:6,0:false:false:.:.:.:.:10,3,0,0 分别对应FORMAT、tumor样本GT、normal样本GT的值。

该信息在Record对象中的存储形式如下：

record.GT [{'GT': '0/1', 'OBQRC': '41.80', 'SB': '12,9,1,1', 'DP': '23', 'OBF': '0.500', 'OBAM': 'false', 'OBP': '0.078', 'AD': '21,2', 'F2R1': '13,1', 'F1R2': '7,1', 'AF': '0.120', 'OBQ': '100.00', 'OBAMRC': 'true'}, {'GT': '0/0', 'OBQRC': '.', 'SB': '10,3,0,0', 'DP': '13', 'OBF': '.', 'OBAM': 'false', 'OBP': '.', 'AD': '13,0', 'F2R1': '6,0', 'F1R2': '7,0', 'AF': '0.065', 'OBQ': '.', 'OBAMRC': 'false'}]

第一个字典代表tumor样本的GT，第二个字典代表normal样本的GT。根据你的样本数量，可能会有多个字典。可以通过索引取出特定样本的信息，例如要获取第一个样本的信息，只需使用record.GT[0]即可。

打印出tumor AF大于0.5的行

for record in gatk:
 # compare GATK tumor AF to 0.05
 if float(record.GT[0]['AF']) > 0.05:
 print(record.line)

将FILTER为PASS且tumor AF大于0.05的行写入列表，并最终写出VCF文件

snv = "filter.vcf"
result = gatk.header
for record in gatk:
 if record.FILTER == "PASS" and float(record.GT[0]['AF']) > 0.05:
 result.append(record.line)
# write out result
with open(snv, 'w+') as snvf:
 for line in result:
 print(line, file = snvf)

查看GATK的下一个Record，由于VCF类是可迭代的，因此除了for循环外，还支持next函数。

record = next(gatk)
print(record.line)

至此，关于Python通用读取VCF文件的类（可直接复制粘贴使用）的文章介绍完毕。更多相关Python VCF文件类的内容，请搜索本站以前的文章或继续浏览下面的相关文章。希望大家以后多多支持本站！