av一区二区在线观看_亚洲男人的天堂网站_日韩亚洲视频_在线成人免费_欧美日韩精品免费观看视频_久草视

您的位置:首頁技術文章
文章詳情頁

python通用讀取vcf文件的類(復制粘貼即可用)

瀏覽:4日期:2022-08-05 08:42:10

前言

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

使用說明

首先復制類的代碼,后面就可以直接用了

import sysimport osimport subprocessclass 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.CHROMrecord.linerecord.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.headerfor record in gatk: if record.FILTER == 'PASS' and float(record.GT[0][’AF’]) > 0.05: result.append(record.line)# write out resultwith 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文件的類內容請搜索好吧啦網以前的文章或繼續瀏覽下面的相關文章希望大家以后多多支持好吧啦網!

標簽: Python 編程
相關文章:
主站蜘蛛池模板: 国产视频日韩 | 久久综合一区 | 欲色av| 狠狠综合久久av一区二区小说 | 成人免费视频观看视频 | 欧美亚洲视频在线观看 | 欧美日韩国产在线观看 | 国产欧美在线观看 | 欧美 日韩 国产 在线 | 欧美激情久久久 | 成人免费高清 | 国产无人区一区二区三区 | 午夜免费电影院 | 欧美小视频在线观看 | 你懂的国产 | 中文精品视频 | 亚洲一区二区不卡在线观看 | 91精品久久久久久久久 | 欧美视频一区 | 亚洲精品国产第一综合99久久 | 免费小视频在线观看 | 日本五月婷婷 | 久久久久久久av | 亚洲码欧美码一区二区三区 | 欧美日韩综合 | 久久精品中文字幕 | 成人超碰| 欧美精品综合在线 | 国产精品久久久久久久久久久久冷 | 色婷婷激情综合 | 欧美性一区二区三区 | 久色网| www.日韩系列 | 成人影院网站ww555久久精品 | 亚洲一区二区久久 | h片在线看 | 天堂资源最新在线 | 精品亚洲一区二区 | 日韩精品一区二区三区中文在线 | 国产一二区免费视频 | 欧美一区二区三区在线看 |