python PyVCF檔案處理VCF檔案格式範例詳解

2022-07-01 14:02:34

引言

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:testexample.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])

呼叫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.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檔案，是因為這個庫考慮的東西可能比我們自己瞭解的更多，其實現也可能比我們自己的程式碼更加完備合理。

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

以上就是python PyVCF檔案處理VCF檔案格式範例詳解的詳細內容，更多關於python PyVCF處理VCF格式的資料請關注it145.com其它相關文章！