前言

處理vcf文件的時(shí)候，需要多種切割，正則匹配，如果要自己寫(xiě)其實(shí)會(huì)比較麻煩，并且每次還得根據(jù)vcf文件格式或者需要讀取的值不同要修改相應(yīng)的代碼。因此很多人會(huì)選擇一些python的vcf的庫(kù)，但是首先你得安裝這個(gè)庫(kù)， 并且有一些庫(kù)它固定了能夠讀的內(nèi)容，如果你的vcf的信息不在它固定的里面，就讀不出來(lái)。比如最近我想讀一個(gè)樣本的AF，但是它放在最后樣本的GT那列，不在INFO那一列，有一些庫(kù)竟然無(wú)能為力。因此我寫(xiě)了這個(gè)通用的讀vcf的類(lèi)，直接復(fù)制粘貼這部分代碼就可以方便的用這個(gè)類(lèi)進(jìn)行vcf文件的讀取，過(guò)濾，寫(xiě)出等操作。

使用說(shuō)明

首先復(fù)制類(lèi)的代碼，后面就可以直接用了

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() 

• 主要有兩個(gè)類(lèi)，一個(gè)是VCF類(lèi),存儲(chǔ)的是vcf的信息，及對(duì)vcf文件的操作，一個(gè)是Record類(lèi)，它包括vcf某一行存儲(chǔ)的全部信息
• 讀入vcf文件

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

查看vcf的header

gatk.header

查看vcf當(dāng)前行中儲(chǔ)存的信息，一開(kāi)始是首行。它以Record這個(gè)類(lèi)保存的。注意VCF類(lèi)是個(gè)迭代器類(lèi)，可以用next和for循環(huán)來(lái)讀入每一行的信息

record = gatk.record #這里record存儲(chǔ)的是該Record類(lèi)的地址

查看該record的屬性，包括line(行的內(nèi)容，方便寫(xiě)出某行), CHROM, POS, ID，REF,ALT, QUAL, FILTER, INFO(字典的形式存儲(chǔ)), FORMAT, sample_num(多少個(gè)樣本),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中的存儲(chǔ)形式

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的存儲(chǔ)形式，F(xiàn)ORMAT對(duì)應(yīng)著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對(duì)應(yīng)的值

這是在record中的存儲(chǔ)形式

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'}]

第一個(gè)字典就是tumor的GT，第二個(gè)字典就是normal的GT，當(dāng)然，根據(jù)你的樣本數(shù)量會(huì)有多個(gè)字典，這里可以按索引取出比如要取出第一個(gè)樣本的，只需要record.GT[0]就行

把tumor AF大于0.5的line打印出來(lái)

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寫(xiě)入列表并寫(xiě)出最后的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的下一個(gè)record， 因?yàn)閂CF類(lèi)是可迭代的，因此除了for也支持next

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

到此這篇關(guān)于python通用讀取vcf文件的類(lèi)（可以直接復(fù)制粘貼使用）的文章就介紹到這了,更多相關(guān)python vcf文件的類(lèi)內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

最新評(píng)論