第3章 VCF文件解析：Python库与自定义解析

VCF文件解析，说白了就是把那些用制表符分隔的文本行，变成程序能直接操作的数据结构。我刚开始做变异数据处理时，第一反应就是找个现成的库——谁愿意自己造轮子呢？但跑过几个项目后发现，不同场景下的解析需求真不一样。有的要速度，有的要兼容，有的要能解析畸形数据。这一章我就把常用的解析手段摊开来讲，帮你选对工具。

3.1 主流Python库横向对比

Python生态里解析VCF的库不少，但真正经得起生产环境考验的，也就那么几个。我个人习惯在正式项目里首选cyvcf2，原因很简单——它基于htslib的C语言底层，解析速度是纯Python实现的10倍以上。你想想看，一个全基因组VCF文件动辄几十GB，解析速度直接决定了你的开发迭代效率。

下面这张表我整理了几个主流库的核心差异：

库名解析速度功能完整性维护状态推荐场景
PyVCF较慢完整维护中小型数据、教学验证
cyvcf2极快完整活跃大规模生产环境
vcfpy中等较完整维护中兼容性优先的场景

选库的时候有一个容易被忽略的点：PyVCF对VCF 4.2的INFO字段解析更细致，但遇到超大文件会明显卡顿。cyvcf2则相反，它在读取时做了惰性加载，内存占用控制得很好。我在项目中遇到过两次因为用了PyVCF导致内存溢出的情况，后来全切到了cyvcf2。

3.2 cyvcf2实战：高效解析VCF

cyvcf2的安装很简单，一行命令就行。但要注意，它依赖htslib，在macOS上可能需要先装libomp：

pip install cyvcf2

基础用法几乎是直觉式的。打开一个VCF文件，遍历每个变异记录，直接访问属性：

from cyvcf2 import VCF

vcf = VCF("input.vcf.gz")

for variant in vcf:
    chrom = variant.CHROM
    pos = variant.POS
    ref = variant.REF
    alt = ",".join(variant.ALT)
    qual = variant.QUAL
    print(f"{chrom}\t{pos}\t{ref}\t{alt}\t{qual}")

小提示： cyvcf2 的 variant.ALT 返回的是一个元组，多个ALT等位基因时需要手动拼接。我一般直接用 ",".join(variant.ALT) 转成字符串，写入DataFrame时方便很多。

cyvcf2还支持按区域查询，这个功能在大文件里特别实用。比如你只想看1号染色体前1Mb的变异：

for variant in vcf("1:1-1000000"):
    print(variant.CHROM, variant.POS, variant.REF, variant.ALT)

注意： 按区域查询要求VCF文件必须经过 bgzip 压缩并用 tabix 建立索引。缺少索引时cyvcf2会抛出 IndexError，而且提示信息不太明显。我第一次用的时候就卡在这里了，后来才发现是忘了跑 tabix -p vcf input.vcf.gz。

3.3 自定义解析器：何时需要，如何实现

现成库大部分场景够用，但总有些边界情况让你想自己动手。举几个例子：VCF文件里的FILTER字段有自定义标签，库的解析逻辑会直接丢弃未知标签；或者你需要把某个INFO子字段拆成多列，库的抽象层反而碍事。

我曾经在项目中遇到过一个坑：上游工具生成的VCF里，QUAL字段在某些行是缺失的，cyvcf2默认返回None，但下游分析流程期望的是字符串"NA"。为了一个字段改整个下游不现实，这时候自己写一个解析器反而是最干净的方案。

核心逻辑其实就几行：

def parse_vcf_line(line):
    if line.startswith("#"):
        return None
    fields = line.strip().split("\t")
    if len(fields) < 8:
        return None

    info_dict = {}
    info_fields = fields[7].split(";")
    for item in info_fields:
        if "=" in item:
            key, val = item.split("=", 1)
            info_dict[key] = val
        else:
            info_dict[item] = True

    return {
        "chrom": fields[0],
        "pos": int(fields[1]),
        "id": fields[2],
        "ref": fields[3],
        "alt": fields[4],
        "qual": fields[5] if fields[5] != "." else "NA",
        "filter": fields[6],
        "info": info_dict,
    }

这个函数保留了原始字段的完整性，QUAL为缺失值时返回"NA"，INFO字段也拆成了字典。关键是你可以在里面自由加逻辑，遇到不合规的行直接返回None，不会让整个解析崩溃。

个人经验： 自定义解析器的单元测试一定要写全。我一般会准备三组测试数据：一行标准格式、一行缺失QUAL、一行INFO为空。这三行覆盖了大部分线上异常场景。

3.4 解析结果结构化：从VCF到DataFrame

解析完VCF，下一步就是把它变成结构化数据。我个人习惯统一输出为Pandas DataFrame，因为后续的过滤、统计、知识图谱构建都能直接复用。用cyvcf2配合DataFrame，代码非常紧凑：

import pandas as pd
from cyvcf2 import VCF

vcf = VCF("input.vcf.gz")
rows = []

for variant in vcf:
    rows.append({
        "chrom": variant.CHROM,
        "pos": variant.POS,
        "ref": variant.REF,
        "alt": ",".join(variant.ALT),
        "qual": variant.QUAL,
        "filter": variant.FILTER,
    })

df = pd.DataFrame(rows)
print(df.head())

如果你需要把结果导出成JSON供下游API使用，直接调用 df.to_json() 即可。但我建议按变异记录拆成列表，方便逐条写入知识图谱：

records = df.to_dict(orient="records")

import json
with open("variants.json", "w") as f:
    json.dump(records, f, indent=2)

核心思路： 解析阶段的输出格式决定了后续整合的复杂度。我强烈建议统一使用"记录列表"结构——每条变异是一个字典，字段名与VCF标准保持一致。这样不管你是要构建知识图谱还是做统计分析，接口都是统一的。

3.5 本章知识体系

下面这张图把VCF解析的整体脉络串了起来。从源文件开始，到选择解析策略，再到最终的结构化输出和知识图谱构建，每一步都是承上启下的。你可以在自己的项目里把它当作路线图，遇到瓶颈时回头看看，就知道该在哪一层做调整。

VCF文件解析知识体系 VCF文件 策略选择 Python库解析 自定义解析 结构化存储 知识图谱构建 解析策略的选择 取决于数据规模与场景

从SVG图里你能看到，解析阶段其实是整个VCF数据整合的"腰部"——上面连着原始文件，下面连着结构化和知识图谱。选库还是自写，没有绝对的好坏，关键看你的数据量、格式规范程度和下游要求。我自己的经验是：先拿cyvcf2快速出原型，遇到定制需求再局部替换成自定义解析，两者互补着用最舒服。

← 上一章 📖 返回目录 下一章 →