python通用读取vcf文件的类(复制粘贴即可用)

前言

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

这是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文件的类内容请搜索我们以前的文章或继续浏览下面的相关文章希望大家以后多多支持我们!

(0)

相关推荐

  • python通用读取vcf文件的类(复制粘贴即可用)

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

  • Python实现读取TXT文件数据并存进内置数据库SQLite3的方法

    本文实例讲述了Python实现读取TXT文件数据并存进内置数据库SQLite3的方法.分享给大家供大家参考,具体如下: 当TXT文件太大,计算机内存不够时,我们可以选择按行读取TXT文件,并将其存储进Python内置轻量级splite数据库,这样可以加快数据的读取速度,当我们需要重复读取数据时,这样的速度加快所带来的时间节省是非常可观的,比如,当我们在训练数据时,要迭代10万次,即要从文件中读取10万次,即使每次只加快0.1秒,那么也能节省几个小时的时间了. #创建数据库并把txt文件的数据存进

  • Python实现的Excel文件读写类

    本文实例讲述了Python实现的Excel文件读写类.分享给大家供大家参考.具体如下: #coding=utf-8 ####################################################### #filename:ExcelRW.py #author:defias #date:2015-4-27 #function:read or write excel file #################################################

  • Java实现的读取资源文件工具类ResourcesUtil实例【可动态更改值的内容】

    本文实例讲述了Java实现的读取资源文件工具类ResourcesUtil.分享给大家供大家参考,具体如下: package com.gcloud.common; import java.io.Serializable; import java.text.MessageFormat; import java.util.ArrayList; import java.util.Iterator; import java.util.List; import java.util.Locale; impor

  • python简单读取大文件的方法

    本文实例讲述了python简单读取大文件的方法.分享给大家供大家参考,具体如下: Python读取大文件(GB级别)采用的办法很简单: with open(...) as f: for line in f: <do something with line> 例如: with open(filepath,'r') as infile: for line in infile: print line 一切都交给python解释器处理,读取效率很高,且占用资源少. stackoverflow参考链接:

  • Python实现读取json文件到excel表

    本文实例为大家分享了Python实现读取json文件到excel表,供大家参考,具体内容如下 一.需求 1.'score.json' 文件内容: { "1":["小花",99,100,98.5], "2":["小王",90,30.5,95], "3":["小明",67.5,49.6,88] } 2.读取json文件保存到数据库,并计算出每个人的总分和平均分 二.实现代码 import j

  • Python简单读取json文件功能示例

    本文实例讲述了Python简单读取json文件功能.分享给大家供大家参考,具体如下: read_json.json: { "rule":{ "namespace":"strategy", "name":"test_exp_1496234234223400", "version":0, "last_modify_time":1434234236819000, "

  • Python实现读取txt文件并画三维图简单代码示例

    记忆力差的孩子得勤做笔记! 刚接触python,最近又需要画一个三维图,然后就找了一大堆资料,看的人头昏脑胀的,今天终于解决了!好了,废话不多说,直接上代码! #由三个一维坐标画三维散点 #coding:utf-8 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d.axes3d import Axes3D x = [] y = [] z = [] f = open("data\\record.

  • python批量读取txt文件为DataFrame的方法

    我们有时候会批量处理同一个文件夹下的文件,并且希望读取到一个文件里面便于我们计算操作.比方我有下图一系列的txt文件,我该如何把它们写入一个txt文件中并且读取为DataFrame格式呢? 首先我们要用到glob模块,这个python内置的模块可以说是非常的好用. glob.glob('*.txt') 得到如下结果: all.txt是我最后得到的结果文件.可以见返回的是一个包含txt文件名称的列表,当然如果你的文件夹下面只有txt文件,那么你用os.listdir()可以得到一个一样的列表 然后

  • Python之读取TXT文件的方法小结

    方法一: <span style="font-size:14px;">#read txt method one f = open("./image/abc.txt") line = f.readline() while line: print line line = f.readline() f.close() </span> 方法二: #read txt method two f = open("./image/abc.txt&q

随机推荐