基因組與Python --PyVCF 好用的vcf文件處理器

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

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

通常處理vcf文件時，在讀取，處理階段總是會寫很多重復代碼，核心的任務代碼很少。當然，如果僅僅是找位點的CHROM，POS，ID，REF，ALT，QUAL這幾個參數時，這樣做也可以。因為vcf格式規范，這幾個參數的結構相對簡單。但是如果處理頭文件信息，或者處理INFO，FORMAT參數時，要寫比較復雜的正則表達式，這樣做不僅繁瑣，而且容易出錯。

Python的PyVCF庫解決了這個問題，它通過正則表達式把vcf文件信息轉換成結構化的信息，簡化了vcf文件的處理過程，方便后續提取相關參數及處理。

PyVCF庫的安裝，cmd界面：

 

pip install PyVCF

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

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 record
Record(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對象------位點信息的儲存形式

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，FILTER，INFO，FORMAT，也就是vcf中的一行位點信息。接下來對這些屬性一一說明：

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

POS：位點在染色體上的位置，類型為int。

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

REF：參考基因組在該位點上的鹼基，類型為str。

ALT：在該位點的測序結果。是_AltRecord類的子類實例的列表。類型為list。_AltRecord類有4個子類，代表了突變的幾種類型：如snp，indel，structual variants等。所有的實例都可以進行比較（僅限於相等的比較，沒有大於小於之說），部分子類沒有實現str方法，也就是說不能轉成字符串。

QUAL：該位點的測序質量，類型為int或float。

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

INFO：該位點的一些測試指標。將‘=’前的參數作為鍵，后面的參數作為值，構建成的字典。類型為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信息作為鍵，后面對應的信息做為值，構建成的字典（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.1
0/1
13115 13116 13116
T [G] ['T', G]

_Record對象方法：

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

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

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

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

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

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

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

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

9. 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 record
Record(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_called
1
>>> record.call_rate
1.0
>>> record.num_hom_ref
0
>>> record.aaf
[0.5]
>>> record.num_het
1
>>> record.heterozygosity
0.5
>>> record.var_type
'snp'
>>> record.var_subtype
'ts'
>>> record.is_snp
True
>>> record.is_indel
False

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:\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.alts
OrderedDict([('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:\test\example.hc.vcf.gz')
>>> vcf_reader.next()
<vcf.model._Record object at 0x0000000003ED8780>
>>> record = vcf_reader.next()
>>> print record
Record(CHROM=chr1, POS=10347, REF=AACCCT, ALT=[A])
>>> for record in vcf_reader:
    print record
Record(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) # 調用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文件，是因為這個庫考慮的東西可能比我們自己了解的更多，其實現也可能比我們自己的代碼更加完備合理。

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