遍历变异位点：基础循环与字段提取

为什么要把迭代放在首位？

拿到一个VCF文件，第一件事肯定是看看里面到底有什么变异。不管是过滤、统计还是输出，都得先把每条记录读出来。pysam 在这块做得特别轻量，你甚至不用写什么 fetch 语句——直接 for 循环就完事了。

我刚开始用 pysam 时，还想着按区间读是不是更快。后来发现，在不知道具体位置的情况下，直接迭代整文件往往是最直观的选择。说白了，它就是 Python 的迭代器，内存占用很低。

for 循环的基本写法

我们先用 pysam.VariantFile 打开 VCF，然后把它当成一个可迭代对象。

import pysam

vcf = pysam.VariantFile("sample.vcf.gz")
for record in vcf:
    # 每一条 record 就是一个变异
    pass

注意这里的 vcf 对象本身就是一个迭代器，不需要再调用 .fetch()。如果你想遍历整个文件，直接 for 就行。这也符合 Python 的“打开即迭代”习惯。

从 record 里拿到核心字段

每条 VariantRecord 里都藏着 contig、位置、REF、ALT 这些关键信息。我列了一个常见字段的获取方式：

这里有个小细节：record.start 是 0-based 的，而 record.stop 是 1-based 的。如果你要用 BED 格式或者某些下游工具，记得统一坐标系。我就在项目里踩过这个坑，算出来的位置总是偏一位，折腾了半天才意识到。

一个完整的遍历例子

下面这段代码打印出每个变异的基本信息，顺便统计一下总条数：

vcf = pysam.VariantFile("input.vcf.gz")
count = 0
for rec in vcf:
    chrom = rec.contig
    pos = rec.start + 1          # 换成 1-based 更直观
    ref = rec.ref
    alts = rec.alts if rec.alts else ("",)
    alt_str = ",".join(alts)
    print(f"{chrom}\t{pos}\t{ref}\t{alt_str}")
    count += 1
print(f"Total variants: {count}")
vcf.close()

你可能会想，为什么不直接用 rec.pos？其实 pysam 里没有直接叫 pos 的属性，只有 start 和 stop。我习惯用时加上 1，这样打印出来的坐标和 VCF 文件里一致。

性能注意事项——别踩的坑

遍历上千万条记录的 VCF 时，性能往往会成为瓶颈。我分成几点来讲：

• 文件格式一定要压缩并索引：别用原始 .vcf，尽量用 bgzip 压缩成 .vcf.gz，然后用 bcftools index 或 tabix 建立索引。没有索引时 pysam 会从头读到尾，速度差一个数量级。我曾在一个 30 万条的小文件上忘了索引，跑了一分多钟，加上索引后不到 5 秒。
• 避免在循环内反复解析复杂类型：比如获取基因型时，rec.samples[0].alleles 返回的是元组列表，如果你只要某个样本的 GT，直接用 rec.samples[0].get("GT") 更快。
• 谨慎使用 rec.format 或 rec.info 的默认访问：当你写 rec.info.DP 时，pysam 会做字段查找——如果这个 INFO 在 header 里不存在，它会抛 KeyError。我习惯先判断一下 if "DP" in rec.info，或者直接用 rec.info.get("DP")，省下不必要的异常处理开销。
• 能用数组尽量别用列表：例如获取所有样本的 GT 时，pysam 提供了 rec.get_format("GT")，返回的是 numpy 数组，处理大量样本时比逐行拆列表快很多。

SVG 流程图：遍历变异记录的典型路径

下图总结了从打开文件到获取字段再到性能检查的整体流程。你可以把它当作思维导图来记：

性能实测对比

我在一个 2.4G 的 WGS VCF 文件上（约 800 万条记录）对比了几种常见的遍历方式，结果如下：

你可以看到，索引比什么都重要。第二条和第三条的差别几乎就是是否用 bgzip 的区别。至于正则解析——我想不出来有什么理由要那么做，除非你在处理纯文本且没有 pysam 的环境下。

一个小技巧：用上下文管理器自动关闭文件

每次手动 close 容易漏掉，尤其是在脚本中途异常退出时。pysam 的 VariantFile 支持 with 语句，我强烈建议这样写：

with pysam.VariantFile("input.vcf.gz") as vcf:
    for rec in vcf:
        # 你的处理逻辑
        pass  # 退出 with 后自动关闭文件

这样就算循环里出了异常，资源也能正确释放。我在生产流水线上就是这么写的，省心很多。

避坑指南——我曾经遇到的一个“怪事”

有一次我写了个脚本，循环里每次都调用 len(rec.alts)，结果跑了一整天没跑完。后来才意识到，rec.alts 是个 tuple，但在这之前我为了获取 ALT 的逗号分隔字符串，写了 ",".join(rec.alts)，也在循环里重复做了。其实完全可以提到循环外面，或者只做一次。如果你只需要判断 ALT 数量，直接 len(rec.alts) 没问题，但要小心 rec.alts 可能是 None（多等位基因行？不，None 表示没有 ALT，比如参考位点）。

嗯，更准确的说法是：VCF 定义每个变异至少有一个 ALT，但有时候你会遇到空值。所以我一般先判断 if rec.alts is None: continue，避开边界情况。

怎么优雅地获取多个样本的基因型？

如果你拿到的是群体 VCF，样本数量可能是几百甚至上千。逐样本循环效率不高，pysam 提供了 rec.get_format_array() 返回 numpy 数组。举个例子：

with pysam.VariantFile("cohort.vcf.gz") as vcf:
    for rec in vcf:
        gt_arr = rec.get_format_array("GT")   # 形状 (n_samples, 2) 或 (n_samples,)
        # 后续可以用 numpy 矩阵操作
        ...

这个方法比 for sample in rec.samples: ... 快很多，因为避免了每次构建 Python 对象。我习惯在处理千人基因组级别的 VCF 时用这个。