引言

vcf文件的全稱是variant call file，即突變識(shí)別文件，它是基因組工作流程中產(chǎn)生的一種文件，保存的是基因組上的突變信息。通過對(duì)vcf文件進(jìn)行分析，可以得到個(gè)體的變異信息。嗯，總之，這是很重要的文件，所以怎么處理它也顯得十分重要。它的文件信息如下：

文件的開頭是一堆以“##”開始的注釋行，包含了文件的基本信息。然后是以“#”開頭的一行，共9+n個(gè)部分，前九部分標(biāo)注的是后面行每部分代表的信息，相當(dāng)于表頭。后面部分是樣本名稱，可以有多個(gè)。注釋行結(jié)束后是具體的突變信息，每一行分為9+n個(gè)部分，每部分之間用制表符（‘\t’）分隔。

通常處理vcf文件時(shí)，在讀取，處理階段總是會(huì)寫很多重復(fù)代碼，核心的任務(wù)代碼很少。當(dāng)然，如果僅僅是找位點(diǎn)的CHROM，POS，ID，REF，ALT，QUAL這幾個(gè)參數(shù)時(shí)，這樣做也可以。因?yàn)関cf格式規(guī)范，這幾個(gè)參數(shù)的結(jié)構(gòu)相對(duì)簡(jiǎn)單。但是如果處理頭文件信息，或者處理INFO，F(xiàn)ORMAT參數(shù)時(shí)，要寫比較復(fù)雜的正則表達(dá)式，這樣做不僅繁瑣，而且容易出錯(cuò)。

Python的PyVCF庫(kù)解決了這個(gè)問題，它通過正則表達(dá)式把vcf文件信息轉(zhuǎn)換成結(jié)構(gòu)化的信息，簡(jiǎn)化了vcf文件的處理過程，方便后續(xù)提取相關(guān)參數(shù)及處理。

PyVCF庫(kù)的安裝

cmd界面

pip install PyVCF

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

PyVCF庫(kù)的導(dǎo)入

import vcf

PyVCF庫(kù)的名字為vcf，導(dǎo)入之后可以使用其方法對(duì)vcf文件做處理。

PyVCF庫(kù)詳細(xì)介紹

使用實(shí)例：

>>> 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中了。它是一個(gè)可迭代對(duì)象，它的迭代元素都是一個(gè)_Record對(duì)象的實(shí)例，保存著非注釋行的一行信息，即變異位點(diǎn)的具體信息。通過它，我們可以很輕易地得到位點(diǎn)的詳細(xì)信息。

_Record對(duì)象------位點(diǎn)信息的儲(chǔ)存形式

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

_Record是vcf.model中的一個(gè)對(duì)象，除了它還有_Call，_AltRecord等對(duì)象。它的基本屬性為CHROM，POS，ID，REF，ALT，QUAL，F(xiàn)ILTER，INFO，F(xiàn)ORMAT，也就是vcf中的一行位點(diǎn)信息。接下來(lái)對(duì)這些屬性一一說(shuō)明：

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

POS：位點(diǎn)在染色體上的位置，類型為int。

ID：一般是突變的rs號(hào)，類型為str。如果是‘.’，則為None。

REF：參考基因組在該位點(diǎn)上的堿基，類型為str。

ALT：在該位點(diǎn)的測(cè)序結(jié)果。是_AltRecord類的子類實(shí)例的列表。類型為list。_AltRecord類有4個(gè)子類，代表了突變的幾種類型：如snp，indel，structual variants等。所有的實(shí)例都可以進(jìn)行比較（僅限于相等的比較，沒有大于小于之說(shuō)），部分子類沒有實(shí)現(xiàn)str方法，也就是說(shuō)不能轉(zhuǎn)成字符串。

QUAL：該位點(diǎn)的測(cè)序質(zhì)量，類型為int或float。

FILTER：過濾信息。將FILTER列按分號(hào)分隔形成的字符串列表，類型為list。如果未給出參數(shù)則為None。

INFO：該位點(diǎn)的一些測(cè)試指標(biāo)。將‘=’前的參數(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對(duì)象還有一些其他屬性：

samples：把FORMAT信息作為鍵，后面對(duì)應(yīng)的信息做為值，構(gòu)建成的字典（CallData對(duì)象），以及sample名稱，這兩個(gè)值組成一個(gè)Call對(duì)象，共同構(gòu)成samples的一個(gè)元素。這樣就把sample和基因型信息給關(guān)聯(lián)起來(lái)，按下標(biāo)訪問每一個(gè)Call對(duì)象。samples類型為list。

start：突變開始的位置

end：突變結(jié)束的位置

alleles：該位點(diǎn)所有的可能情況，由REF和ALT參數(shù)組成的列表（REF類型是str，ALT參數(shù)是_AltRecord對(duì)象的子類實(shí)例），類型是list。

>>> for record in vcf_reader:       
print record.samples, 
'\n', record.samples[0].sample, 
'\n', record.samples[0]['GT'] 
#按下標(biāo)訪問Call，按.sample訪問sample，按鍵訪問FORMAT對(duì)應(yīng)信息      
print record.start, record.POS, record.end      
print record.REF, record.ALT, record.alleles 
#注意G沒有引號(hào)，它是_AltRecord對(duì)象 
[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對(duì)象方法：

• 對(duì)象之間比較大小方法：根據(jù)染色體名稱和位置信息比較。Python3只實(shí)現(xiàn)了‘=’和‘<’的比較。
• 迭代方法：對(duì)samples里的元素進(jìn)行迭代。
• 字符串方法：只返回CHROM，POS，REF，ALT四列信息。
• genotype(name)方法，和samples按下標(biāo)訪問不同，這個(gè)方法提供按sample名稱進(jìn)行訪問的功能。
• add_format(fmt)， add_filter(flt)， add_info(info, value=True)：給相應(yīng)的屬性添加元素。
• get_hom_refs()：拿到samples中該位點(diǎn)未突變的所有sample，返回列表。
• get_hom_alts()：拿到samples中該位點(diǎn)100%突變的所有sample，返回列表。
• get_hets()：拿到samples中該位點(diǎn)基因型為雜合的所有sample，返回列表。
• get_unknown()：拿到samples中該位點(diǎn)基因型未知的所有sample，返回列表。

>>> record = next(vcf_reader)
>>> record2 = next(vcf_reader)
>>> print record > record2 
#按染色體名稱和位置進(jìn)行比較False
>>> for i in record: 
#按samples列表進(jìn)行迭代        
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名字進(jìn)行訪問
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對(duì)象還有很多有用的方法屬性：

num_called：該位點(diǎn)已識(shí)別的sample數(shù)目。

call_rate：已識(shí)別的sample數(shù)目占sample總數(shù)的比例。

num_hom_ref，num_hom_alt，num_het，num_unknown：四種基因型的數(shù)量

aaf：所有sample等位基因的頻率（即除開REF），返回列表。

heterozygosity：該位點(diǎn)的雜合度，0.5為雜合突變，0為純合突變。

var_type：突變類型，包括‘snp’，‘indel’，‘sv’（structural variant），‘unknown’。

var_subtype：更加細(xì)化的突變類型，如‘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 
#只有一個(gè)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對(duì)象------處理vcf文件，構(gòu)建結(jié)構(gòu)化信息

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

在讀vcf文件時(shí)，總共有六個(gè)參數(shù)可供選擇，如上圖所示。

fsock：目標(biāo)文件的文件對(duì)象，可以用open(文件名)得到這個(gè)文件對(duì)象。

filename：文件名，當(dāng)fsock和filename同時(shí)存在時(shí)，優(yōu)先考慮fsock。

compressed：是否要解壓，不提供參數(shù)時(shí)由程序自行判斷（以文件名是否以.gz結(jié)尾判斷是否需要解壓）。
prepend_chr：在保存染色體名稱時(shí)，是否加前綴‘chr’，默認(rèn)不加，如果vcf文件的染色體名稱本來(lái)沒有前綴‘chr’，可設(shè)置為True，自動(dòng)加上。

strict_whitespace：是否嚴(yán)格以制表符‘\t’分隔數(shù)據(jù)。True則表示嚴(yán)格按制表符分，F(xiàn)alse表示可以?shī)A雜空格。

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對(duì)象的屬性中，包括alts，contigs，filters，formats，infos，metadata。

alts使用實(shí)例：

>>> 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對(duì)象實(shí)現(xiàn)了兩個(gè)方法：

next()：獲得下一行的數(shù)據(jù)，也就是返回下一個(gè)_Record對(duì)象?？梢燥@式調(diào)用next()得到下一行數(shù)據(jù)，也可以直接迭代Reader對(duì)象，它會(huì)自動(dòng)調(diào)用next()函數(shù)以獲得下一行數(shù)據(jù)。

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

start和end都不給，就返回chrom染色體所有的突變位點(diǎn)。這個(gè)方法需要用另一個(gè)第三方Python模塊pysam來(lái)建立文件索引，如果沒有安裝這個(gè)模塊，會(huì)導(dǎo)致錯(cuò)誤。

另外，使用這個(gè)方法之后，它會(huì)將對(duì)象的可迭代范圍改成fetch()得到的突變位點(diǎn)，所以用這個(gè)方法，原來(lái)的迭代進(jìn)度就失效了。

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

這個(gè)庫(kù)還有一個(gè)Writer對(duì)象，在此就不詳細(xì)介紹了，因?yàn)榇蟛糠謱?duì)vcf文件的處理都可以用上面兩個(gè)對(duì)象的知識(shí)搞定。

綜合使用：

#!/usr/bin/env python
# -*- coding: utf-8 -*- import vcf  
# 導(dǎo)入PyVCF庫(kù) 
filename = r'D:\test\example.hc.vcf.gz' 
vcf_reader = vcf.Reader(filename=filename) 
# 調(diào)用Reader對(duì)象處理vcf文件 
for record in vcf_reader: 
# 迭代Reader對(duì)象,返回的是_Record對(duì)象   
# record是_Record對(duì)象    
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:等價(jià)        
     print "I'm not a sv"    
     if record.heterozygosity == 0.5: 
         # 判斷是否為雜合突變        
         print "I'm a heterozygous mutation"    
         ...    ...

這個(gè)庫(kù)實(shí)現(xiàn)的所有功能，都可以自己寫代碼實(shí)現(xiàn)，而且實(shí)現(xiàn)方法比較簡(jiǎn)單。之所以要用這個(gè)庫(kù)來(lái)處理vcf文件，是因?yàn)檫@個(gè)庫(kù)考慮的東西可能比我們自己了解的更多，其實(shí)現(xiàn)也可能比我們自己的代碼更加完備合理。

還有，重復(fù)造車總歸是不好的。

以上就是python PyVCF文件處理VCF文件格式實(shí)例詳解的詳細(xì)內(nèi)容，更多關(guān)于python PyVCF處理VCF格式的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！