基因组与Python --PyVCF 好用的vcf文件处理器

本文链接：https://blog.youkuaiyun.com/qq_31088877/article/details/76921339

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类