pysam

用途

pysam 用于对 HTSlib 支持的基因组格式进行底层、流式访问，主要面向需要在 Python 中读取、查询、过滤或写入 SAM/BAM/CRAM、VCF/BCF、FASTA/FASTQ 以及 tabix 索引表格的场景。它封装了 HTSlib/samtools/bcftools，当前上游基线为 pysam 0.24.0（2026 年 4 月 27 日），对应 HTSlib/samtools/bcftools 1.23.1。

主要能力

AlignmentFile 与 AlignedSegment：处理 SAM/BAM/CRAM 比对文件。
VariantFile、VariantHeader、VariantRecord：处理 VCF/BCF 变异文件。
FastaFile：访问带索引的 FASTA；FastxFile：顺序读取 FASTA/FASTQ。
TabixFile：访问 BGZF 压缩、tabix 索引的 BED/GFF/GTF 或自定义表格。
pysam.samtools 与 pysam.bcftools：封装命令行调度器。

安装

使用固定版本以保证可复现：

uv pip install "pysam==0.24.0"

确认运行时：

import pysam
print(pysam.__version__)           # 0.24.0
print(pysam.__samtools_version__)  # 1.23.1

支持的 macOS 和 Linux 平台提供预编译 wheel。源码构建需要 C 编译器和 HTSlib 构建依赖。

使用方式

编码前应先确认真实格式、压缩方式、排序顺序和可用索引；区分数值 Python 坐标与区域字符串，不要混用；CRAM 需明确参考组装和 FASTA；优先使用索引区域访问，仅在有意为之的情况下使用顺序迭代；写入时保留头部并默认写入新路径；明确过滤语义（比对/碱基质量、标志、重叠处理、重复处理、pileup 深度上限）。

对于不熟悉的文件，可使用内置只读检查脚本：

python scripts/inspect_hts.py sample.bam
python scripts/inspect_hts.py cohort.vcf.gz
python scripts/inspect_hts.py reference.fa

内置脚本包括：inspect_hts.py（仅元数据检查）、alignment_qc.py（流式聚合读取/QC 计数为 JSON）、variant_summary.py（流式变异、FILTER 和基因型摘要为 JSON）、filter_alignments.py（过滤 SAM/BAM/CRAM 且不改变记录顺序）。所有脚本拒绝覆盖已有输出。

坐标约定：pysam API 接受的数值坐标为 0 基半开区间，包括 AlignmentFile.fetch()、VariantFile.fetch()、FastaFile.fetch()、TabixFile.fetch() 和 pileup() 的数值参数。区域字符串为 samtools 风格，1 基闭区间。VCF 文本使用 1 基 POS，而记录属性同时暴露两种体系：record.pos（1 基）、record.start（0 基闭）、record.stop（0 基开）。

比对文件使用上下文管理器和显式模式，例如：

import pysam

with pysam.AlignmentFile("sample.bam", "rb", threads=4) as bam:
    for read in bam.fetch("chr1", 1_000, 2_000):
        if (
            not read.is_unmapped
            and not read.is_secondary
            and not read.is_supplementary
            and read.mapping_quality >= 30
        ):
            print(read.query_name, read.reference_start, read.cigarstring)

使用 fetch(until_eof=True) 可按文件顺序流式读取所有记录，包括未定位的未比对读段，且不需要索引。注意 fetch() 返回与区域重叠的比对记录；count() 计数且默认 read_callback="nofilter"；count_coverage() 返回 A/C/G/T 碱基计数且默认碱基质量 15 加 read_callback="all"；pileup() 暴露每列读段并有自己的过滤、碱基质量、重叠、孤儿和 max_depth=8000 默认值。

变异文件输入格式自动检测，数值 fetch 坐标仍为 0 基：

import pysam

with pysam.VariantFile("cohort.vcf.gz", threads=4) as variants:
    for record in variants.fetch("chr1", 999_999, 2_000_000):
        print(record.contig, record.pos, record.ref, record.alts)
        for sample_name, call in record.samples.items():
            print(sample_name, call.get("GT"))

可在检索记录前子集样本。更改头部时，应复制每条记录并翻译到目标头部，再分配新声明的 INFO/FORMAT/FILTER 字段。

FASTA、FASTQ 与 Tabix：带索引 FASTA 使用数值 0 基坐标；FastxFile 为顺序读取，persist=False 更快但记录在迭代推进后失效；Tabix 输入必须坐标排序且为 BGZF 压缩（非普通 gzip），使用非破坏性两步工作流：

pysam.tabix_compress("regions.bed", "regions.bed.gz")
pysam.tabix_index("regions.bed.gz", preset="bed")

CRAM、远程 I/O 与线程：pysam 0.24 改变了继承的 HTSlib 行为——新写入的 CRAM 默认使用 CRAM 3.1 而非 3.0；HTSlib 默认不再联系 EBI 参考服务器；建议使用 reference_filename="reference.fa" 进行确定性本地读写。仅在有意进行基于 MD5 的参考查找时配置 REF_PATH/REF_CACHE。不要假设 CRAM 是自包含的。threads= 加速压缩/解压缩，但不并行化 Python 分析。

封装的 samtools 和 bcftools：显式导入命令模块，每个命令行 token 作为单独字符串传入。调度器默认捕获 stdout；对于大型或二进制输出，使用工具的 -o 选项配合 catch_stdout=False 或 save_stdout=...，而不是在内存中返回完整输出。记录级逻辑使用 Python API，成熟批量操作（如 sort、index、merge、view、normalization）使用调度器。切勿通过拆分不可信 shell 命令来组合调度器参数。

写入规则：打开输出前复制或构造有效头部；写入新路径，除非明确替换否则不要使用 force=True；若输出将被索引则保留排序顺序；先设置 query_sequence 再设置 query_qualities；优先使用 pysam.CIGAR_OPS 枚举成员；使用 pysam.samtools.quickcheck() 验证比对输出，并在下游使用前重新打开变异/序列输出；当参考或坐标超出旧索引限制时使用 CSI 而非 BAI/TBI。

适用场景

适用于基因组数据分析、变异检测与汇总、序列读取与处理、区域查询与过滤、格式转换与索引等生物信息学工作流。相关技能包括 polars-bio（基于 Polars DataFrame 的基因组区间操作与文件 I/O）、tiledbvcf（TileDB-VCF 数据集的存储与查询）、gtars（基因组区间工具包）。

限制

常见失败模式包括：将数值 VariantFile.fetch() 坐标当作 1 基；在需要 BGZF 加 tabix/CSI 的地方使用普通 gzip；无索引调用区域 fetch；假设 fetch() 包含未定位的未比对读段；精确 pileup 区间忘记 truncate=True；忽略 pileup 默认值（如碱基质量 13 和深度上限 8000）；在活动迭代器或线程间共享一个文件句柄；解码 CRAM 时缺少精确参考；在输出头部声明前分配新 VCF 字段；在内存中捕获大型 samtools/bcftools 输出；对 indel 或符号等位基因使用 SNP 碱基计数方法。

技能介绍

用途

主要能力

安装

使用方式

适用场景

限制