VCF解析实践：使用Python读取VCF

做生物信息的人，谁也绕不开VCF文件。无论是看单样本的突变位点，还是处理大队列WGS数据，第一步都是把VCF读进来、提字段、做分析。我这些年接触过的工具不少，但要说灵活性和可定制性，还得是Python。今天我们就聊聊怎么用Python把VCF里的变异位点信息干净地提出来。

为什么用Python解析VCF

你可能会问，VCF是文本文件，直接写awk不也能处理吗？嗯，简单场景确实可以。但你想想看，VCF的INFO字段结构复杂，样本格式字段更是五花八门，光靠文本切割很容易踩坑。Python有专用的操作库，能直接把VCF解析成结构化对象，省心不少。

说白了，专业的事交给专业的库。Python社区有两个库最常用：vcfpy 和 pysam。它们各有侧重点，下面我会分别演示。

两个主流库：vcfpy vs pysam

先做个对比，让你有个整体印象：

我个人的建议很简单：做原型开发和数据探索时用vcfpy，代码写起来顺手；到了生产环境，尤其是几百G的WGS数据，还是pysam更靠谱，速度差距肉眼可见。

环境准备

安装不复杂，pip一行搞定：

pip install vcfpy pysam

记得pysam在Linux/macOS下需要编译C扩展，Windows用户建议用conda：

conda install -c bioconda pysam vcfpy

嗯，这里要注意版本兼容性。我遇到过pysam新版对某些老VCF格式支持有变化的情况，建议先在测试数据上跑一遍再上全量数据。

使用vcfpy读取VCF

vcfpy的API很直观。打开文件、遍历记录、访问字段，三步走：

import vcfpy

# 打开VCF文件
reader = vcfpy.Reader.from_path('input.vcf')

# 遍历每一条变异记录
for record in reader:
    chrom = record.CHROM
    pos = record.POS
    ref = record.REF
    alts = [a.value for a in record.ALT]
    qual = record.QUAL
    print(chrom, pos, ref, alts, qual)

reader.close()

你看，record.CHROM、record.POS 这些属性名和VCF列名几乎一一对应，上手很快。不过要注意，ALT字段返回的是 Allele 对象列表，需要用 .value 拿到字符串。我在项目里第一次用的时候忘了这步，debug了好一会儿。

INFO字段是个字典结构，访问起来也很方便：

# 提取INFO字段
dp = record.INFO.get('DP')
af = record.INFO.get('AF')
print(f"DP={dp}, AF={af}")

样本级别的数据通过 record.calls 访问：

for call in record.calls:
    sample_name = call.sample
    gt = call.data.get('GT')
    print(sample_name, gt)

使用pysam读取VCF

pysam的用法略有不同，它用上下文管理器打开文件，字段名也小写：

import pysam

with pysam.VariantFile('input.vcf') as vcf:
    for record in vcf:
        chrom = record.chrom
        pos = record.pos
        ref = record.ref
        alts = record.alts
        qual = record.qual
        print(chrom, pos, ref, alts, qual)

注意这里 record.alts 直接返回元组，不需要再转换。pysam的INFO访问也更像字典：

    dp = record.info.get('DP')
    af = record.info.get('AF')

样本数据通过 record.samples 访问：

    for sample_name, sample_data in record.samples.items():
        gt = sample_data['GT']
        print(sample_name, gt)

我个人习惯在大队列数据处理时用pysam，速度优势很明显。之前跑一个10x WGS的VCF，pysam比vcfpy快了将近3倍，内存也更省。