本文主要是給大家介紹了pythonPyVCF文檔處理VCF文件類型范例詳細說明，感興趣的小伙伴值得借鑒參考一下，希望可以有一定的幫助，祝愿大家多多的不斷進步，盡快工作上得到晉升

　　前言

　　vcf文件的全名是是variantcallfile，即突變性鑒別文檔，這是基因工作內容過程中產生的一類文檔，存放的是基因里的突變性信息內容。根據對vcf文件展開分析，可以獲得自我的基因變異信息內容。嗯，總而言之，這是非常重要的文檔，因此如何處理它也變得極為重要。它文件信息如下：

　　文檔開頭是一堆以“##”開始注解行，包括了文檔的相關信息。然后就是以“#”開頭一列，共9+若干個一部分，前九一部分注明的是后邊行每一部分意味著的數據，等同于頁眉。后邊部分是樣本名稱，可以有多個。注解行完成后是實際的突變性信息內容，每一行分成9+若干個一部分，每一部分間用分隔符（‘t’）隔開。

　　一般解決vcf文件時，在載入，解決環節總會寫許多重復代碼，關鍵任務編碼非常少。自然，如果只是找結構域的CHROM，POS，ID，REF，ALT，QUAL這些主要參數時，這么做還可以。由于vcf格式標準，這些參數構造較為簡單。但如果解決庫函數信息內容，或是解決INFO，FORMAT主要參數時，寫較為復雜的正則匹配，這么做不但繁雜，還很容易出差錯。

　　Python的PyVCF庫克服了這種情況，主要是通過正則匹配把vcf文件信息內容轉化成結構型的數據，優化了vcf文件的處理方式，便捷后面獲取主要參數及解決。

　　PyVCF庫安裝

　　cmd界面

　　pipinstallPyVCF

　　或者從https://github.com/jamescasbon/PyVCF網站上下載安裝包，自行安裝。

　　PyVCF庫的導入

　　importvcf

　　PyVCF庫詳細介紹

　　使用范例：

　　>>>importvcf
　　>>>vcf_reader=vcf.Reader(filename=r'D:testexample.hc.vcf.gz')
　　>>>forrecordinvcf_reader:
　　printrecordRecord(CHROM=chr1,POS=10146,REF=AC,ALT=[A])Record(CHROM=chr1,POS=10347,REF=AACCCT,ALT=[A])Record(CHROM=chr1,POS=10439,REF=AC,ALT=[A])Record(CHROM=chr1,POS=10492,REF=C,ALT=[T])Record(CHROM=chr1,POS=10583,REF=G,ALT=[A])

　　調用vcf.Reader類解決vcf文件，vcf文件信息內容就被保存到vcf_reader中了。它是一個可迭代對象，它迭代元素都是一個_Record對象的范例，保存著非注解行的一列信息內容，即基因變異結構域的具體信息。通過它，我們可以很輕易地得到結構域的詳細信息。

　　_Record對象------結構域信息的儲存形式

　　classvcf.model._Record(CHROM,POS,ID,REF,ALT,QUAL,FILTER,INFO,FORMAT,sample_indexes,samples=None)

　　_Record是vcf.model中的一個對象，除了它還有_Call，_AltRecord等對象。它基本屬性為CHROM，POS，ID，REF，ALT，QUAL，FILTER，INFO，FORMAT，也就是vcf中的一列結構域信息內容。接下來對這些屬性一一說明：

　　CHROM：染色體名稱，類型為str。

　　POS：結構域在染色體里的位置，種類為int。

　　ID：一般是突變性的rs號，類型為str。如果是.’，則為None。

　　REF：參考基因組在這個結構域里的堿基，類型為str。

　　ALT：在這個結構域的測序結果。是_AltRecord類的派生類范例的目錄。種類為list。_AltRecord類有4個子類，代表著突變幾類種類：如snp，indel，structualvariants等。每一個范例都可以直接相對比較（僅限相等相對比較，并沒有大于小于之說），一部分派生類沒有實現str方法，也就是說不能轉成字符串。

　　QUAL：該結構域的測序質量，種類為int或float。

　　FILTER：過濾信息。將FILTER列按分號隔開形成的字符串目錄，種類為list。如果未給出主要參數則為None。

　　INFO：該結構域的一些測試指標。將‘=’前的主要參數作為鍵，后面的主要參數作為值，構建成的字典。種類為dict。

　　FORMAT：基因型信息內容。保存vcf的FORMAT列的原始形式

　　>>>for record in vcf_reader:
　　print type(record.CHROM),record.CHROM
　　print type(record.POS),record.POS
　　print type(record.ID),record.ID
　　print type(record.REF),record.REF
　　print type(record.ALT),record.ALT
　　print type(record.QUAL),record.QUAL
　　print type(record.FILTER),record.FILTER
　　print type(record.INFO),record.INFO
　　print type(record.INFO['BaseQRankSum']),record.INFO['BaseQRankSum']
　　print type(record.FORMAT),record.FORMAT
　　<type'str'>chr1<type'int'>
　　234481<type'NoneType'>None<type'str'>
　　T<type'list'>[A]<type'float'>2025.77<type'NoneType'>
　　None<type'dict'>{'ExcessHet':3.0103,'AC':[1],
　　'BaseQRankSum':-2.743,'MLEAF':[0.5],'AF':[0.5],
　　'MLEAC':[1],'AN':2,'FS':2.371,'MQ':42.83,
　　'ClippingRankSum':0.0,'SOR':0.972,'MQRankSum':-2.408,
　　'ReadPosRankSum':1.39,'DP':156,'QD':13.07}<type'float'>
　　-2.743<type'str'>GT:AD:DP:GQ:PL

　　除了這些基本屬性之外，_Record對象還有一些其他屬性：

　　samples：把FORMAT信息作為鍵，后面對應的信息做為值，構建成的字典（CallData對象），以及sample名稱，這兩個值組成一個Call對象，共同構成samples的一個元素。這樣就把sample和基因型信息給關聯起來，按下標訪問每一個Call對象。samples類型為list。

　　start：突變開始的位置

　　end：突變結束的位置

　　alleles：該位點所有的可能情況，由REF和ALT參數組成的列表（REF類型是str，ALT參數是_AltRecord對象的子類實例），類型是list。

　　>>>for record in vcf_reader:
　　print record.samples,
　　'n',record.samples[0].sample,
　　'n',record.samples[0]['GT']
　　#按下標訪問Call，按.sample訪問sample，按鍵訪問FORMAT對應信息
　　print record.start,record.POS,record.end
　　print record.REF,record.ALT,record.alleles
　　#注意G沒有引號，它是_AltRecord對象
　　[Call(sample=192.168.1.1,CallData(GT=0/1,AD=[39,14],
　　DP=53,GQ=99,PGT=0|1,PID=13116_T_G,PL=[449,0,2224]))]
　　192.168.1.10/113115 13116 13116T[G]['T',G]

　　_Record對象方法：

　　對象之間比較大小方法：根據染色體名稱和位置信息比較。Python3只實現了‘=’和‘<’的比較。

　　迭代方法：對samples里的元素進行迭代。

　　字符串方法：只返回CHROM，POS，REF，ALT四列信息。

　　genotype(name)方法，和samples按下標訪問不同，這個方法提供按sample名稱進行訪問的功能。

　　add_format(fmt)，add_filter(flt)，add_info(info,value=True)：給相應的屬性添加元素。

　　get_hom_refs()：拿到samples中該位點未突變的所有sample，返回列表。

　　get_hom_alts()：拿到samples中該位點100%突變的所有sample，返回列表。

　　get_hets()：拿到samples中該位點基因型為雜合的所有sample，返回列表。

　　get_unknown()：拿到samples中該位點基因型未知的所有sample，返回列表。

　　>>>record=next(vcf_reader)
　　>>>record2=next(vcf_reader)
　　>>>print record>record2
　　#按染色體名稱和位置進行比較False
　　>>>for i in record:
　　#按samples列表進行迭代
　　print i
　　Call(sample=192.168.1.1,
　　CallData(GT=0/1,AD=[18,11],
　　DP=29,GQ=99,PL=[280,0,528]))
　　>>>print str(record)
　　#字符串方法Record(CHROM=chr1,POS=10492,REF=C,ALT=[T])
　　>>>print record.genotype('192.168.1.1')
　　#按sample名字進行訪問
　　Call(sample=192.168.1.1,
　　CallData(GT=0/1,AD=[39,14],DP=53,GQ=99,PGT=0|1,PID=13116_T_G,PL=[449,0,2224]))
　　_Record對象還有很多有用的方法屬性：
　　num_called：該位點已識別的sample數目。
　　call_rate：已識別的sample數目占sample總數的比例。
　　num_hom_ref，num_hom_alt，num_het，num_unknown：四種基因型的數量
　　aaf：所有sample等位基因的頻率（即除開REF），返回列表。
　　heterozygosity：該位點的雜合度，0.5為雜合突變，0為純合突變。
　　var_type：突變類型，包括‘snp’，‘indel’，‘sv’（structural variant），‘unknown’。
　　var_subtype：更加細化的突變類型，如‘indel’包括‘del’，‘ins’，‘unknown’。
　　is_snp，is_indel，is_sv，is_transition，is_deletion：判斷突變是不是snp，indel，sv，transition，deletion等等。
　　>>>record=next(vcf_reader)
　　>>>print recordRecord(CHROM=chr1,POS=13118,REF=A,ALT=[G])
　　>>>print record.samples
　　#只有一個sample
　　[Call(sample=192.168.1.1,CallData(GT=0/1,AD=[41,13],DP=54,GQ=99,PGT=0|1,PID=13116_T_G,PL=[449,0,2224]))]
　　>>>record.num_called1
　　>>>record.call_rate1.0
　　>>>record.num_hom_ref0
　　>>>record.aaf[0.5]
　　>>>record.num_het1
　　>>>record.heterozygosity0.5
　　>>>record.var_type'snp'
　　>>>record.var_subtype'ts'
　　>>>record.is_snpTrue
　　>>>record.is_indelFalse

　　Reader對象------處理vcf文件，構建結構化信息

　　class Reader(fsock=None,filename=None,compressed=None,prepend_chr=False,strict_whitespace=False,encoding='ascii')

　　在讀vcf文件時，總共有六個參數可供選擇，如上圖所示。

　　fsock：目標文件的文件對象，可以用open(文件名)得到這個文件對象。

　　filename：文件名，當fsock和filename同時存在時，優先考慮fsock。

　　compressed：是否要解壓，不提供參數時由程序自行判斷（以文件名是否以.gz結尾判斷是否需要解壓）。

　　prepend_chr：在保存染色體名稱時，是否加前綴‘chr’，默認不加，如果vcf文件的染色體名稱本來沒有前綴‘chr’，可設置為True，自動加上。

　　strict_whitespace：是否嚴格以制表符‘t’分隔數據。True則表示嚴格按制表符分，False表示可以夾雜空格。

　　encoding：文件編碼。

　　>>>vcf_reader=vcf.Reader(open('vcf/test/example-4.0.vcf','r'))
　　#fsock
　　>>>vcf_reader=vcf.Reader(filename=r'D:testexample.hc.vcf.gz')
　　#filename

　　頭文件信息主要保存在Reader對象的屬性中，包括alts，contigs，filters，formats，infos，metadata。

　　alts使用實例：

　　>>>vcf_reader=vcf.Reader(filename=r'D:testexample.hc.vcf.gz')
　　>>>vcf_reader.altsOrderedDict([('NON_REF',Alt(id='NON_REF',desc='Represents any possible alternative allele at this location'))])

　　#字典類型

　　>>>vcf_reader.alts['NON_REF'].id'NON_REF'
　　>>>vcf_reader.alts['NON_REF'].desc'Represents any possible alternative allele at this location'

　　其他的屬性用法類似。

　　Reader對象實現了兩個方法：

　　next()：獲得下一行的數據，也就是返回下一個_Record對象。可以顯式調用next()得到下一行數據，也可以直接迭代Reader對象，它會自動調用next()函數以獲得下一行數據。

　　fetch(chrom，start=None，end=None)：返回chrom染色體從start+1到end坐標的所有突變位點。不給end，就返回chrom染色體從start+1到末尾的所有突變位點；

　　start和end都不給，就返回chrom染色體所有的突變位點。這個方法需要用另一個第三方Python模塊pysam來建立文件索引，如果沒有安裝這個模塊，會導致錯誤。

　　另外，使用這個方法之后，它會將對象的可迭代范圍改成fetch()得到的突變位點，所以用這個方法，原來的迭代進度就失效了。

　　>>>vcf_reader=vcf.Reader(filename=r'D:testexample.hc.vcf.gz')
　　>>>vcf_reader.next()<vcf.model._Record object at 0x0000000003ED8780
　　>>>>record=vcf_reader.next()
　　>>>print recordRecord(CHROM=chr1,POS=10347,REF=AACCCT,ALT=[A])
　　>>>for record in vcf_reader:
　　print recordRecord(CHROM=chr1,POS=10439,REF=AC,ALT=[A])Record(CHROM=chr1,POS=10492,REF=C,ALT=[T])

　　這個庫還有一個Writer對象，在此就不詳細介紹了，因為大部分對vcf文件的處理都可以用上面兩個對象的知識搞定。

　　綜合使用：

　　#!/usr/bin/env python
　　#-*-coding:utf-8-*-import vcf
　　#導入PyVCF庫
　　filename=r'D:testexample.hc.vcf.gz'
　　vcf_reader=vcf.Reader(filename=filename)
　　#調用Reader對象處理vcf文件
　　for record in vcf_reader:
　　#迭代Reader對象,返回的是_Record對象
　　#record是_Record對象
　　print record.CHROM,record.POS,record.ID,record.ALT
　　if record.is_snp:
　　#判斷是否是snp
　　print"I'm a snp"
　　elif record.var_type!='sv':
　　#和elif record.is_sv:等價
　　print"I'm not a sv"
　　if record.heterozygosity==0.5:
　　#判斷是否為雜合突變
　　print"I'm a heterozygous mutation"
　　......

　　這個庫實現的所有功能，都可以自己寫代碼實現，而且實現方法比較簡單。之所以要用這個庫來處理vcf文件，是因為這個庫考慮的東西可能比我們自己了解的更多，其實現也可能比我們自己的代碼更加完備合理。

　　還有，重復造車總歸是不好的。

　　綜上所述，這篇內容就給大家介紹到這里了，希望可以給各位讀者帶來幫助。