第3章 pysam 读取 VCF 文件

VCF 解析的起点就一个函数：pysam.VariantFile。这句话虽短，但背后有不少门道。我最初接触 pysam 时，以为它跟普通文件打开差不多——指定路径、选个模式就行。结果第一次跑就报了个「Not a VCF file」的错误，折腾了半天才发现是压缩格式没处理好。所以这一章咱们把打开、读头、遍历这些基本操作彻底聊透。

本章核心目标：学会用 pysam 打开 VCF 文件，理解三种访问模式，能正确读取头信息并高效遍历所有记录。

pysam.VariantFile 读取 VCF 流程 VariantFile() 模式选择 r / w / a 读取 header for variant in vcf vcf.close()

3.1 用 VariantFile 打开 VCF

方法签名很简单：

import pysam

vcf = pysam.VariantFile("input.vcf.gz", "r")
# 或直接给路径，pysam 会自动判断格式
vcf = pysam.VariantFile("input.vcf")

嗯，这里要特别说一下：pysam 支持 .vcf、.vcf.gz、.bcf 甚至是 ftp:// 和 s3:// 路径（当然要装对应库）。我遇到过一位同事，用 pysam 读本地的 .vcf 文件，结果忘了加模式参数，默认行为跟预期不一致——代码没报错，但写操作把原文件覆盖了。所以，显式传模式参数是个好习惯。

3.2 模式选择：r / w / a

模式参数跟普通文件打开类似，但含义略有区别：

模式 含义 典型场景
"r" 只读 分析已有 VCF / BCF
"w" 写入（覆盖） 创建新 VCF 或覆盖旧文件
"a" 追加 向已有 VCF 末尾添加记录（少用）
"rb" / "wb" 二进制读写（BCF） 处理 BCF 格式时需指定

我个人习惯：写工具时总是先用 "r" 打开做验证，确认文件可读且头内容正确后，再用 "w" 创建输出文件。这样能避免中途发现文件损坏而白跑。

3.3 读取头信息

打开后第一件事就是检查头。我经常在头里找 ##INFO 和 ##FORMAT 的定义，看看文件是不是我预期的版本。头信息存在 vcf.header 对象里：

vcf = pysam.VariantFile("sample.vcf.gz", "r")
header = vcf.header

# 查看头记录的原始文本
print(header.records[:3])

# 获取所有样本名
samples = list(header.samples)
print(f"样本数: {len(samples)}, 样本名: {samples[:5]}")

# 查看特定字段定义
info_def = header.info.get("DP")
print(f"DP 定义: {info_def}")  # 输出类型、描述等

# 查看 contig 列表
for contig in header.contigs.values():
    print(f"{contig.name}  length={contig.length}")

vcf.close()

你想想看，有时候同一个 VCF 文件被不同工具处理过，头信息里可能会残留下游工具的注释。我记得有一次跑流程，发现某条 INFO 字段被两个工具重复定义了，类型还不一致。后来养成习惯，先打印 header.records 过滤出 ##INFO 看看有没有重复，能省掉不少 debug 时间。

注意：头信息在写入模式下会自动生成基础行（##fileformat=VCFv4.2、#CHROM 等），但如果你想自定义 ##INFO / ##FORMAT，必须显式通过 header.add_line() 添加。我在生产中遇到过因为忘记添加而写出不规范的 VCF，被下游工具拒收。

头信息的核心结构可以用这个表格快速回忆：

属性 类型 说明
header.records list 所有原始头行（包括 ## 行和 #CHROM 行）
header.samples list 样本名列表（按文件中的列顺序）
header.contigs OrderedDict contig 信息，键为 contig 名，值为 Contig 对象
header.info dict-like INFO 字段定义，键为字段名
header.formats dict-like FORMAT 字段定义
header.filters dict-like FILTER 字段定义

3.4 遍历记录

遍历才是 VCF 处理的主战场。pysam 把 VariantFile 对象设计成可迭代的，每条记录是一个 VariantRecord 对象。基本写法：

vcf = pysam.VariantFile("input.vcf.gz", "r")
for record in vcf:
    # record 类型: pysam.VariantRecord
    chrom = record.chrom
    pos = record.pos       # 1-based
    rid = record.id
    ref = record.ref
    alt = record.alts      # tuple of str
    qual = record.qual
    filter = record.filter

    # 获取 INFO 字段
    dp = record.info.get("DP")
    af = record.info.get("AF")

    # 获取 FORMAT / 样本数据
    for sample in vcf.header.samples:
        gt = record.samples[sample]["GT"]
        ad = record.samples[sample]["AD"]

    # 只处理前几条做测试
    if record.pos > 1000000:
        break
vcf.close()

这个 for record in vcf 看起来很直观，但如果你要处理超大 VCF（比如全基因组 80GB 的 WGS 文件），内存并不会爆掉——pysam 底层是流式读取，一次只缓存一个记录。我原来用 pandas 读 VCF，几百兆就崩了；换成 pysam 流式遍历，几十 GB 的文件安稳跑完。

性能提醒：如果只关注某个区域，比如 chr1:100000-200000，可以用 vcf.fetch("chr1", 100000, 200000)。pysam 会自动使用索引（.tbi / .csi）跳转，比全量遍历快几十倍。我曾经在做基因面板靶向分析时，用 fetch 把 5 分钟的处理降到了 3 秒。

关于 record.id 还有个小坑：有些 VCF 文件不包含 ID 列（全是 .），此时 record.id 是 None。保险做法：

rid = record.id if record.id is not None else f"{record.chrom}:{record.pos}"

3.5 实战示例：快速扫描 VCF 摘要

下面这段代码，能在一分钟内给你一个百万行 VCF 的全景：总变异数、每条染色体的分布、变异类型统计。

import pysam
from collections import Counter

def scan_vcf(vcf_path):
    vcf = pysam.VariantFile(vcf_path, "r")
    print(f"样本数: {len(vcf.header.samples)}")
    print(f"Contig 数: {len(vcf.header.contigs)}")

    total = 0
    chr_counter = Counter()
    type_counter = Counter()

    for record in vcf:
        total += 1
        chr_counter[record.chrom] += 1
        if len(record.ref) == 1 and record.alts and len(record.alts[0]) == 1:
            type_counter["SNP"] += 1
        elif len(record.ref) > 1 or (record.alts and any(len(a) > 1 for a in record.alts)):
            type_counter["Indel"] += 1
        else:
            type_counter["Other"] += 1

    vcf.close()

    print(f"\n总变异数: {total:,}")
    print("\n染色体分布 (Top 5):")
    for chrom, cnt in chr_counter.most_common(5):
        print(f"  {chrom}: {cnt:,}")
    print(f"\n变异类型:")
    for t, cnt in type_counter.most_common():
        print(f"  {t}: {cnt:,}")

    return {"total": total, "by_chrom": chr_counter, "by_type": type_counter}

# 运行
stats = scan_vcf("my_data.vcf.gz")

这里用了 collections.Counter 做计数，简单高效。你可能会问：为什么不用 pandas ？因为在流式读取中维护一个 Pandas DataFrame 会累积所有行到内存，违反了轻量原则。等你想分析统计结果时，再从本地的 SQLite 或 Parquet 里读出来就行。

教训分享：我曾经写过一个 VCF 统计脚本，一开始图省事把所有记录的 CHROM、POS、TYPE 存进列表再统计，结果一个 5M 行的 VCF 直接吃掉 12GB 内存。后来改成流式计数+定期落盘的策略，内存降到 200MB。记住：处理生物数据时，能流式就不要缓存。

最后提一个容易被忽略的操作：用 with 语句管理上下文。pysam 的 VariantFile 支持上下文管理器：

with pysam.VariantFile("data.vcf.gz", "r") as vcf:
    for record in vcf:
        # 处理记录
        pass
# 自动关闭，无需显式 close()

这个写法更 Pythonic，也避免了你忘记关文件导致的 fd 泄漏。我在多线程处理时吃过亏——开了 100 个线程各连一个 VCF，忘了关闭，结果操作系统把进程 kill 了。从那以后，with 就成了我的默认选择。

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