python通讯录文件读取_python通用读取vcf文件的类（可以直接复制粘贴使用）

前言

处理vcf文件的时候，需要多种切割，正则匹配，如果要自己写其实会比较麻烦，并且每次还得根据vcf文件格式或者需要读取的值不同要修改相应的代码。因此很多人会选择一些python的vcf的库，但是首先你得安装这个库， 并且有一些库它固定了能够读的内容，如果你的vcf的信息不在它固定的里面，就读不出来。比如最近我想读一个样本的AF，但是它放在最后样本的GT那列，不在INFO那一列，有一些库竟然无能为力。

因此我写了这个通用的读vcf的类，直接复制粘贴这部分代码就可以方便的用这个类进行vcf文件的读取，过滤，写出等操作。

使用说明

首先复制类的代码，后面就可以直接用了

import sys

import os

import 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的读取

这是GT在vcf的存储形式，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的line打印出来

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)