python?PyVCF文件處理VCF文件格式實例詳解
引言
vcf文件的全稱是variant call file，即突變識別文件，它是基因組工作流程中產(chǎn)生的一種文件，保存的是基因組上的突變信息。通過對vcf文件進行分析，可以得到個體的變異信息。嗯，總之，這是很重要的文件，所以怎么處理它也顯得十分重要。它的文件信息如下：
文件的開頭是一堆以“##”開始的注釋行，包含了文件的基本信息。然后是以“#”開頭的一行，共9+n個部分，前九部分標注的是后面行每部分代表的信息，相當于表頭。后面部分是樣本名稱，可以有多個。注釋行結(jié)束后是具體的突變信息，每一行分為9+n個部分，每部分之間用制表符（‘\t’）分隔。
通常處理vcf文件時，在讀取，處理階段總是會寫很多重復代碼，核心的任務代碼很少。當然，如果僅僅是找位點的CHROM，POS，ID，REF，ALT，QUAL這幾個參數(shù)時，這樣做也可以。因為vcf格式規(guī)范，這幾個參數(shù)的結(jié)構(gòu)相對簡單。但是如果處理頭文件信息，或者處理INFO，F(xiàn)ORMAT參數(shù)時，要寫比較復雜的正則表達式，這樣做不僅繁瑣，而且容易出錯。
Python的PyVCF庫解決了這個問題，它通過正則表達式把vcf文件信息轉(zhuǎn)換成結(jié)構(gòu)化的信息，簡化了vcf文件的處理過程，方便后續(xù)提取相關參數(shù)及處理。
PyVCF庫的安裝
cmd界面
pip install PyVCF
或者從https://github.com/jamescasbon/PyVCF網(wǎng)站上下載安裝包，自行安裝。
PyVCF庫的導入
import vcf
PyVCF庫的名字為vcf，導入之后可以使用其方法對vcf文件做處理。
PyVCF庫詳細介紹
使用實例：
>>> import vcf >>> vcf_reader = vcf.Reader(filename=r'D:\test\example.hc.vcf.gz') >>> for record in vcf_reader: print recordRecord(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])
調(diào)用vcf.Reader類處理vcf文件，vcf文件信息就被保存到vcf_reader中了。它是一個可迭代對象，它的迭代元素都是一個_Record對象的實例，保存著非注釋行的一行信息，即變異位點的具體信息。通過它，我們可以很輕易地得到位點的詳細信息。
_Record對象------位點信息的儲存形式
class vcf.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，F(xiàn)ILTER，INFO，F(xiàn)ORMAT，也就是vcf中的一行位點信息。接下來對這些屬性一一說明：
CHROM：染色體名稱，類型為str。
POS：位點在染色體上的位置，類型為int。
ID：一般是突變的rs號，類型為str。如果是‘.’，則為None。
REF：參考基因組在該位點上的堿基，類型為str。
ALT：在該位點的測序結(jié)果。是_AltRecord類的子類實例的列表。類型為list。_AltRecord類有4個子類，代表了突變的幾種類型：如snp，indel，structual variants等。所有的實例都可以進行比較（僅限于相等的比較，沒有大于小于之說），部分子類沒有實現(xiàn)str方法，也就是說不能轉(zhuǎn)成字符串。
QUAL：該位點的測序質(zhì)量，類型為int或float。
FILTER：過濾信息。將FILTER列按分號分隔形成的字符串列表，類型為list。如果未給出參數(shù)則為None。
INFO：該位點的一些測試指標。將‘=’前的參數(shù)作為鍵，后面的參數(shù)作為值，構(gòu)建成的字典。類型為dict。
FORMAT：基因型信息。保存vcf的FORMAT列的原始形式，類型為str。
>>> 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信息作為鍵，后面對應的信息做為值，構(gòu)建成的字典（CallData對象），以及sample名稱，這兩個值組成一個Call對象，共同構(gòu)成samples的一個元素。這樣就把sample和基因型信息給關聯(lián)起來，按下標訪問每一個Call對象。samples類型為list。
start：突變開始的位置
end：突變結(jié)束的位置
alleles：該位點所有的可能情況，由REF和ALT參數(shù)組成的列表（REF類型是str，ALT參數(shù)是_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對象方法：
- 對象之間比較大小方法：根據(jù)染色體名稱和位置信息比較。Python3只實現(xiàn)了‘=’和‘<’的比較。
- 迭代方法：對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數(shù)目。
call_rate：已識別的sample數(shù)目占sample總數(shù)的比例。
num_hom_ref，num_hom_alt，num_het，num_unknown：四種基因型的數(shù)量
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文件，構(gòu)建結(jié)構(gòu)化信息
class Reader(fsock=None, filename=None, compressed=None, prepend_chr=False, strict_whitespace=False, encoding='ascii')
在讀vcf文件時，總共有六個參數(shù)可供選擇，如上圖所示。
fsock：目標文件的文件對象，可以用open(文件名)得到這個文件對象。
filename：文件名，當fsock和filename同時存在時，優(yōu)先考慮fsock。
compressed：是否要解壓，不提供參數(shù)時由程序自行判斷（以文件名是否以.gz結(jié)尾判斷是否需要解壓）。
prepend_chr：在保存染色體名稱時，是否加前綴‘chr’，默認不加，如果vcf文件的染色體名稱本來沒有前綴‘chr’，可設置為True，自動加上。
strict_whitespace：是否嚴格以制表符‘\t’分隔數(shù)據(jù)。True則表示嚴格按制表符分，F(xiàn)alse表示可以夾雜空格。
encoding：文件編碼。
>>> vcf_reader = vcf.Reader(open('vcf/test/example-4.0.vcf', 'r'))
#fsock
>>> vcf_reader = vcf.Reader(filename=r'D:\test\example.hc.vcf.gz')
#filename
頭文件信息主要保存在Reader對象的屬性中，包括alts，contigs，filters，formats，infos，metadata。
alts使用實例：
>>> vcf_reader = vcf.Reader(filename=r'D:\test\example.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對象實現(xiàn)了兩個方法：
next()：獲得下一行的數(shù)據(jù)，也就是返回下一個_Record對象。可以顯式調(diào)用next()得到下一行數(shù)據(jù)，也可以直接迭代Reader對象，它會自動調(diào)用next()函數(shù)以獲得下一行數(shù)據(jù)。
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:\test\example.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:\test\example.hc.vcf.gz'
vcf_reader = vcf.Reader(filename=filename)
# 調(diào)用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"
... ...
這個庫實現(xiàn)的所有功能，都可以自己寫代碼實現(xiàn)，而且實現(xiàn)方法比較簡單。之所以要用這個庫來處理vcf文件，是因為這個庫考慮的東西可能比我們自己了解的更多，其實現(xiàn)也可能比我們自己的代碼更加完備合理。
還有，重復造車總歸是不好的。
以上就是python PyVCF文件處理VCF文件格式實例詳解的詳細內(nèi)容，更多關于python PyVCF處理VCF格式的資料請關注腳本之家其它相關文章！
相關文章
這篇文章主要介紹了python導入不同目錄下的自定義模塊過程解析,文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值,需要的朋友可以參考下2019-11-11
這篇文章主要為大家介紹了Python強制轉(zhuǎn)換，具有一定的參考價值，感興趣的小伙伴們可以參考一下，希望能夠給你帶來幫助2021-12-12
Pandas數(shù)據(jù)形狀df.shape的實現(xiàn)
本文主要介紹了Pandas數(shù)據(jù)形狀df.shape的實現(xiàn)，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧2022-07-07
這篇文章主要介紹了Python中的迭代器與生成器使用及說明，具有很好的參考價值，希望對大家有所幫助。如有錯誤或未考慮完全的地方，望不吝賜教2022-12-12
這篇文章主要為大家介紹了python下grpc與protobuf的編寫使用，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步早日升職加薪2022-04-04
這篇文章主要介紹了讀Json文件生成pandas數(shù)據(jù)框詳情，文章圍繞主題展開詳細的內(nèi)容介紹，具有一定的參考價值，感興趣的小伙伴可以參考一下2022-08-08
pycharm 2019 最新激活方式（pycharm破解、激活）
這篇文章主要介紹了最新2019pycharm激活方式（pycharm破解、激活）,本文給大家介紹的非常詳細，具有一定的參考借鑒價值，需要的朋友可以參考下2020-01-01
tensorflow中tf.reduce_mean函數(shù)的使用
這篇文章主要介紹了tensorflow中tf.reduce_mean函數(shù)的使用，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧2020-04-04
python使用dabl幾行代碼實現(xiàn)數(shù)據(jù)處理分析及ML自動化
這篇文章主要為大家分享一個名為?dabl?的開源?Python?工具包，它可以自動化機器學習模型開發(fā)，包括數(shù)據(jù)預處理、特征可視化和分析、建模2021-11-11
Python使用Keras OCR實現(xiàn)從圖像中刪除文本
這篇文章主要為大家介紹了如何在Python中利用Keras OCR實現(xiàn)快速地從圖像中刪除文本，從而作為圖像分類器的預處理步驟，需要的可以參考一下2022-03-03