基因组与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
或者从github/jamescasbon/PyVCF⽹站上下载安装包，⾃⾏安装。
PyVCF库的导⼊：
import vcf
PyVCF库的名字为vcf，导⼊之后可以使⽤其⽅法对vcf⽂件做处理。
PyVCF库详细介绍：
使⽤实例：
>>> import vcf
>>> vcf_reader = vcf.Reader(filename=r'D:\test\example.')
>>> 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对象------位点信息的储存形式
del._Record(CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, sample_indexes, samples=None)
_Record是del中的⼀个对象，除了它还有_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' <type 'float'> -2.743
<type 'str'> GT:AD:DP:GQ:PL
除了这些基本属性之外，_Record对象还有⼀些其他属性：
samples：把FORMAT信息作为键，后⾯对应的信息做为值，构建成的字典（CallData对象），以及sample名称，这两个值组成⼀个Call
对象，共同构成samples的⼀个元素。这样就把sample和基因型信息给关联起来，按下标访问每⼀个C
all对象。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, d
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])
>>> pe('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_calledpython怎么读取py文件
1
>>> record.call_rate
1.0
>>> record.num_hom_ref
>>> 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.') #filename
头⽂件信息主要保存在Reader对象的属性中，包括alts，contigs，filters，formats，infos，metadata。
alts使⽤实例：
>>> vcf_reader = vcf.Reader(filename=r'D:\test\example.')
>>> vcf_reader.alts
OrderedDict([('NON_REF', Alt(id='NON_REF', desc='Represents any possible alternative allele at thi
s 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.')
>>> ()
&del._Record object at 0x0000000003ED8780>
>>> record = ()
>>> 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.'
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⽂件，是因为这个库考虑的东西可能⽐我们⾃⼰了解的更多，其实现也可能⽐我们⾃⼰的代码更加完备合理。
还有，重复造车总归是不好的。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。