📚 课程目录第2章 VCF解析

2. VCF 解析：Python 读取 VCF

做 VCF 分析，第一步就是把它读进内存里。这部分我踩过不少坑，今天把两种最主流的 Python 工具——pysam 和 cyvcf2——给你讲透。

2.1 两种主流工具的定位

你可能会问：为什么有了 pysam 还要用 cyvcf2？说白了，这俩工具的设计目标不一样。

特性	pysam	cyvcf2
核心依赖	htslib (C 库)	htslib + cython
处理速度	中等	快（约 3-5x）
字段访问	通过列索引或字典	直接 Python 属性访问
Bcf 支持	良好	原生支持
学习曲线	稍陡	平缓

我个人习惯：如果是快速原型验证，直接用 cyvcf2，代码写得快；如果是需要精细控制过滤逻辑，或者和已有的 SAMtools/BCFtools 脚本配合，pysam 更合适。

2.2 安装与环境配置

安装其实很简单，但有个坑我得先说。我曾经在旧版 conda 环境里直接 pip install cyvcf2，结果编译报错了一下午，后来发现是 htslib 版本冲突。

推荐用 conda，它会把依赖一起管好：

conda create -n vcf_analysis python=3.9
conda activate vcf_analysis

# 安装 pysam
conda install -c bioconda pysam

# 安装 cyvcf2
conda install -c bioconda cyvcf2

建议优先用 conda 安装。pip 版有时会缺 htslib 的动态库，conda 版的预编译包更稳定。如果遇到“libhts.so 找不到”的错误，八成是环境变量 LD_LIBRARY_PATH 没设好。

2.3 pysam 基本读取操作

pysam 封装了 htslib，所以底层操作和 BCFtools 非常像。你想想看，它的 VariantFile 其实就是 C 语言中 htslib 的 bcf_sr_t 的 Python 包装。

打开一个 VCF 文件：

import pysam

# 打开文件，'r' 是读模式
vcf = pysam.VariantFile("sample.vcf.gz", "r")

# 查看头部信息
print("样本列表:", vcf.header.samples)
print("第1条染色体的长度:", vcf.header.contigs.get("chr1").length)

遍历记录是核心操作。这里有一个重点：记录是按染色体顺序排序的，但如果你没建索引，它不会报错，只是遍历时会非常慢。

for record in vcf:
    print(f"位置: {record.chrom}:{record.pos}")
    print(f"参考等位基因: {record.ref}")
    print(f"替代等位基因: {record.alts}")
    print(f"质量值: {record.qual}")
    print(f"过滤器状态: {record.filter.keys()}")
    break  # 只演示第一条

访问基因型信息时，我早期犯过一个错误——直接拿 .gt_types 来比较，结果发现数字含义和我想的不一样。正确方式如下：

# 访问第一个样本的基因型
sample_name = vcf.header.samples[0]
genotype = record.samples[sample_name]

# gt_type: 0=纯合参考, 1=杂合, 2=纯合替代, 3=缺失
print(f"基因型类型: {genotype.gt_type}")
print(f"等位基因深度: {genotype.ad}")
print(f"基因型质量: {genotype.gq}")

record.filter.keys() 返回的是通过/未通过集合的键，不是所有标签。如果一个位点的 FILTER 列为 PASS，那 keys() 返回 {'PASS'}；如果是 . 或空，就返回空集合。别拿它来判断“是否有过滤标签”。

2.4 cyvcf2 快速上手

cyvcf2 的设计哲学是“让 VCF 操作像操作普通 Python 对象一样自然”。我用它写过一个小工具，40 行的脚本处理了 1.2 亿行 VCF，性能非常惊艳。

基本读取和 pysam 类似：

from cyvcf2 import VCF

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

# 查看样本
print("样本数量:", len(vcf.samples))

# 遍历前两条记录
for i, variant in enumerate(vcf):
    if i >= 2:
        break
    print(f"CHROM: {variant.CHROM}, POS: {variant.POS}")
    print(f"REF: {variant.REF}, ALT: {variant.ALT}")
    print(f"QUAL: {variant.QUAL}, FILTER: {variant.FILTER}")

访问基因型是 cyvcf2 的强项。它的 num_called 和 gt_types 是属性，不用调用方法：

# 获取第一个样本的基因型
sample_gt = variant.genotypes[0]
print(f"等位基因: {sample_gt[:2]}")  # [0,1] 表示杂合
print(f"是否定相: {sample_gt[2]}")   # True/False

# 更常用的方式
gt_type = variant.gt_types[0]  # 0/1/2/3
depth = variant.gt_depths[0]
qual = variant.gt_quals[0]

如果要写高性能的生产级脚本，cyvcf2 是首选。它在底层做了类型缓存和批量内存分配，同样大小的文件，cyvcf2 能比 pysam 快 2-3 倍，而且语法更简洁。

2.5 两种工具的选择策略

我结合项目经验，总结了一个简单的决策树：

代码量优先（< 50 行） → 选 cyvcf2。写起来像在写 JSON 解析，省时间。
需要精细过滤（复杂过滤器） → 选 pysam。它和 BCFtools 的过滤表达式兼容，可以直接传递正则或过滤串。
处理大量样本（1000+） → 必须用 cyvcf2。大规模样本时 pysam 的内存碎片化问题会更严重。
和 R/Julia 做交叉验证 → 两者都可以，但我推荐 pysam，因为它输出的是 htslib 的标准字节流，其他语言工具兼容性更好。

2.6 避坑指南：索引与压缩

VCF 文件在处理前有两个硬性要求：必须是 .vcf.gz 格式（用 bgzip 压缩），必须有 .tbi 或 .csi 索引。缺少任何一个，代码都会在遍历时报错或无限卡死。

我曾经接手过一个项目，同事把 VCF 文件用 gzip 压缩了就直接用，结果 pysam 始终报“truncated file”错误。后来才发现，VCF 必须用 bgzip（块压缩），不能用普通 gzip：

# 正确压缩方式
bgzip -c input.vcf > input.vcf.gz
tabix -p vcf input.vcf.gz

# 错误方式（会导致解析失败）
gzip -c input.vcf > input.vcf.gz

可以用 file 命令快速检查：file sample.vcf.gz，如果输出包含 "BGZF" 或 "blocksorted"，就是对的；如果只显示 "gzip compressed data"，那就要重新压缩。

建立索引不只是为了速度——没有索引，pysam 的 fetch() 方法会直接抛异常。而 cyvcf2 虽然没有强制要求索引，但如果你文件超过 100 MB 还没有索引，遍历效率会低到让人崩溃。

2.7 核心流程总结

这个流程我实际跑过上百次。记住一句话：索引是命。没有索引，再好的工具也快不起来。