python通用读取vcf文件的类（复制粘贴即可用）
前言
处理vcf文件的时候，需要多种切割，正则匹配，如果要自己写其实会比较麻烦，并且每次还得根据vcf文件格式或者需要读取的值不同要修改相应的代码。因此很多人会选择一些python的vcf的库，但是首先你得安装这个库，并且有一些库它固定了能够读的内容，如果你的vcf的信息不在它固定的里面，就读不出来。比如最近我想读一个样本的AF，但是它放在最后样本的GT那列，不在INFO那一列，有一些库竟然无能为力。因此我写了这个通用的读vcf的类，直接复制粘贴这部分代码就可以方便的用这个类进行vcf文件的读取，过滤，写出等操作。
使用说明
首先复制类的代码，后面就可以直接用了
importsys
importosimportsubprocess
classRecord(object):
'''
Onelineinformationinvcffile
'''
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]iflen(pair_lst)>1else""}forpair_lstin[pair.split("=")forpairininfo[7].split(";")]]
self.FORMAT=info[8].split(":")
self.sample_num=len(info)-7
self.GT=[]
foriinrange(2):
GT_value=info[8+i+1].split(":")
GT_dict={}
forginrange(len(GT_value)):
GT_dict[self.FORMAT[g]]=GT_value[g]
self.GT.append(GT_dict)
classVCF(object):
'''
VCFclass,readVCF,writeVCF,getVCFinformation
'''
def__init__(self,uncompress_vcf):
self.header=[]
self.reader=open(uncompress_vcf,'r')
self.line=self.reader.readline().strip()
whileself.line.startswith('#'):
self.header.append(self.line)
self.line=self.reader.readline().strip()
self.record=Record(self.line)
def__iter__(self):
returnself
def__next__(self):
self.line=self.reader.readline().strip()
ifself.line!="":
self.record=Record(self.line)
returnself.record
else:
self.reader.close()
raiseStopIteration()
defreader_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:SB0/1:21,2:0.120:23:7,1:13,1:false:true:0.500:0.078:100.00:41.80:12,9,1,10/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]就行
把tumorAF大于0.5的line打印出来
forrecordingatk: #compareGATKtumorAFto0.05 iffloat(record.GT[0]['AF'])>0.05: print(record.line)
把FILTER为PASS的并且tumorAF>0.05写入列表并写出最后的VCF文件
snv="filter.vcf" result=gatk.header forrecordingatk: ifrecord.FILTER=="PASS"andfloat(record.GT[0]['AF'])>0.05: result.append(record.line) #writeoutresult withopen(snv,'w+')assnvf: forlineinresult: print(line,file=snvf)
查看gatk的下一个record，因为VCF类是可迭代的，因此除了for也支持next
record=next(gatk) print(record.line)
到此这篇关于python通用读取vcf文件的类（可以直接复制粘贴使用）的文章就介绍到这了,更多相关pythonvcf文件的类内容请搜索毛票票以前的文章或继续浏览下面的相关文章希望大家以后多多支持毛票票！
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。