python如何读取vcf文件的类-创新互联
这篇文章主要介绍python如何读取vcf文件的类,文中介绍的非常详细,具有一定的参考价值,感兴趣的小伙伴们一定要看完!
创新互联长期为千余家客户提供的网站建设服务,团队从业经验10年,关注不同地域、不同群体,并针对不同对象提供差异化的产品和服务;打造开放共赢平台,与合作伙伴共同营造健康的互联网生态环境。为老城企业提供专业的成都网站设计、网站建设,老城网站改版等技术服务。拥有十多年丰富建站经验和众多成功案例,为您定制开发。使用说明
首先复制类的代码,后面就可以直接用了
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的读取
这是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)
以上是“python如何读取vcf文件的类”这篇文章的所有内容,感谢各位的阅读!希望分享的内容对大家有帮助,更多相关知识,欢迎关注创新互联行业资讯频道!
标题名称:python如何读取vcf文件的类-创新互联
网站链接:http://scyanting.com/article/dcdcgj.html