pysam

pysam

概览

使用 pysam 以低层级、流式方式访问 HTSlib 支持的基因组格式：

• 使用 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 调用封装的命令分发器

当前上游基准版本：pysam 0.24.0（2026 年 4 月 27 日），封装了 HTSlib/samtools/bcftools 1.23.1。更新特定版本的指南前，请先阅读 references/sources.md。

安装

为确保工作可复现，请使用固定版本：

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 构建依赖；请阅读 references/sources.md 中链接的官方安装指南。

首先确定

编写代码前：

1. 确认实际格式、压缩方式、排序顺序和可用索引。
2. 决定使用 Python 数值坐标还是区域字符串，不要混用。
3. 对于 CRAM，确认准确的参考基因组组装版本和 FASTA。
4. 优先使用带索引的区域访问；仅在有意遍历全部记录时才顺序迭代。
5. 写入时保留头信息，并默认写入新路径。
6. 明确筛选语义：比对质量/碱基质量、标志位、重叠处理、重复项处理以及 pileup 深度上限。

对于不熟悉的文件，请先使用随附的只读检查工具：

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

随附脚本

脚本	用途	常用调用方式
scripts/inspect_hts.py	仅检查比对、变异、FASTA、FASTQ 和 tabix 文件的元数据	python scripts/inspect_hts.py sample.cram --reference ref.fa
scripts/alignment_qc.py	以流式方式统计读取和 QC 汇总数据，并输出为 JSON	python scripts/alignment_qc.py sample.bam --max-records 100000
scripts/variant_summary.py	以流式方式汇总变异、FILTER 和基因型，并输出为 JSON	python scripts/variant_summary.py cohort.vcf.gz --region chr1:1-1000000
scripts/filter_alignments.py	筛选 SAM/BAM/CRAM，同时不改变记录顺序	python scripts/filter_alignments.py input.bam output.bam --exclude-secondary

所有脚本都会拒绝覆盖已有输出。使用 --help 运行各脚本，查看坐标、索引和隐私相关说明。

坐标约定

**pysam API 接受的数值坐标采用从 0 开始、左闭右开的形式。**这包括 AlignmentFile.fetch()、VariantFile.fetch()、FastaFile.fetch()、TabixFile.fetch() 和 pileup() 的数值参数。

区域字符串采用 samtools 风格：从 1 开始且两端均包含。

# The same 100 bases:
bam.fetch("chr1", 99, 199)          # [99, 199)
bam.fetch(region="chr1:100-199")    # 1-based inclusive

VCF 文本使用从 1 开始的 POS，而记录属性会同时提供两种坐标系统：

record.pos    # 1-based
record.start  # 0-based inclusive
record.stop   # 0-based exclusive

有关格式转换、重叠语义、索引选择和 contig 名称检查，请阅读 references/coordinates_and_indexing.md。

比对文件

使用上下文管理器和显式模式：

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) 按文件顺序流式读取所有记录，包括未放置的未比对读取片段，且无需索引：

with pysam.AlignmentFile("sample.bam", "rb") as bam:
    for read in bam.fetch(until_eof=True):
        ...

重要区别：

• fetch() 返回与某个区域重叠的比对记录。
• count() 统计记录数，默认使用 read_callback="nofilter"。
• count_coverage() 返回 A/C/G/T 碱基计数，默认要求碱基质量为 15，并使用 read_callback="all"。
• pileup() 提供按列排列的读取片段，并有自己的筛选、碱基质量、重叠、孤立读取片段和 max_depth=8000 默认设置。

对于精确区域的 pileup，请设置 truncate=True 并显式指定筛选条件：

with pysam.FastaFile("reference.fa") as fasta, pysam.AlignmentFile(
    "sample.bam", "rb"
) as bam:
    for column in bam.pileup(
        "chr1",
        1_000,
        2_000,
        truncate=True,
        stepper="samtools",
        fastafile=fasta,
        min_mapping_quality=20,
        min_base_quality=20,
        max_depth=100_000,
    ):
        print(column.reference_pos, column.get_num_aligned())

有关标志位、CIGAR 操作、标签、修饰碱基、记录写入、pileup 详情和迭代器生命周期，请阅读 references/alignment_files.md。

变异文件

输入格式会自动检测。数字型获取坐标仍采用从 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"))

在检索记录之前筛选样本：

with pysam.VariantFile("cohort.bcf") as variants:
    variants.subset_samples(["sample_A", "sample_B"])
    for record in variants:
        ...

更改 header 时，先复制每条记录并将其转换为目标 header，再添加新声明的 INFO/FORMAT/FILTER 字段。不要手动清空并重建 header.samples。

阅读 references/variant_files.md，了解安全使用 header、写入、样本筛选、缺失基因型、符号等位基因、过滤、转换和索引。

FASTA、FASTQ 和 Tabix

索引 FASTA 使用从 0 开始的数字型坐标：

with pysam.FastaFile("reference.fa") as fasta:
    sequence = fasta.fetch("chr1", 999, 1_099)

FastxFile 按顺序读取。persist=False 速度更快，但迭代继续后，已生成的记录会失效：

with pysam.FastxFile("reads.fastq.gz", persist=False) as reads:
    for read in reads:
        qualities = read.get_quality_array()
        ...

Tabix 输入必须按坐标排序并采用 BGZF 压缩，不能使用普通 gzip。请采用无损的两步工作流：

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

with pysam.TabixFile("regions.bed.gz", parser=pysam.asBed()) as tbx:
    for interval in tbx.fetch("chr1", 1_000, 2_000):
        print(interval.contig, interval.start, interval.end)

阅读 references/sequence_files.md，了解 FASTA/FASTQ 记录和安全创建 tabix 索引。

CRAM、远程 I/O 和线程

pysam 0.24 更改了继承自 HTSlib 的行为：

• 新写入的 CRAM 默认采用 CRAM 3.1，而非 3.0。
• HTSlib 默认不再访问 EBI 参考序列服务。
• 为确保本地读写结果稳定，优先使用 reference_filename="reference.fa"。

with pysam.AlignmentFile(
    "sample.cram",
    "rc",
    reference_filename="reference.fa",
    threads=4,
) as cram:
    for read in cram.fetch("chr1", 1_000, 2_000):
        ...

仅在有意通过 MD5 查找参考序列时配置 REF_PATH/REF_CACHE。不要假定 CRAM 文件包含所有必需信息。threads= 可加速压缩/解压缩，但不会并行执行 Python 分析。

在转换 CRAM、远程访问或并发迭代之前，请阅读 references/cram_and_performance.md。

封装的 samtools 和 bcftools

显式导入命令模块。将每个命令行标记作为单独的字符串传入：

import pysam.samtools
import pysam.bcftools

pysam.samtools.sort(
    "-@", "4", "-o", "sorted.bam", "input.bam", catch_stdout=False
)
pysam.samtools.index("-@", "4", "sorted.bam", catch_stdout=False)

pysam.bcftools.index("--csi", "variants.vcf.gz", catch_stdout=False)

分派器默认会捕获标准输出。对于大型或二进制输出，请使用工具的 -o 选项并设置 catch_stdout=False，或使用 save_stdout=...，而不要将完整输出返回到内存中。

try:
    pysam.samtools.quickcheck("-v", "sample.bam")
except pysam.SamtoolsError as error:
    messages = pysam.samtools.quickcheck.get_messages()
    raise RuntimeError(messages or str(error)) from error

对记录级逻辑使用 Python API；对 sort、index、merge、view 和 normalization 等成熟的批量操作使用分派器。切勿通过拆分不可信的 shell 命令来组合分派器参数。

写入规则

• 打开输出文件之前，先复制或构建有效的 header。
• 写入新路径；除非明确要求替换，否则不要使用 force=True。
• 如果输出将建立索引，请保留排序顺序。
• 先设置 query_sequence，再设置 query_qualities。
• 优先使用 pysam.CIGAR_OPS 枚举成员；pysam.CMATCH 等顶层常量是兼容性别名，未来将移除。
• 对比对文件使用 pysam.samtools.quickcheck() 验证输出；在下游使用变异/序列输出之前，重新打开并检查。
• 如果参考序列或坐标超出旧版索引的限制，请使用 CSI，而非 BAI/TBI。

参考资料索引

需求	阅读
比对 API、标志、CIGAR、pileup、修饰碱基	references/alignment_files.md
VCF/BCF header、记录、样本、写入	references/variant_files.md
FASTA/FASTQ 和 tabix 索引表格	references/sequence_files.md
坐标转换和索引选择	references/coordinates_and_indexing.md
CRAM 参考序列、远程 I/O、线程、性能	references/cram_and_performance.md
正确的集成分析模式	references/common_workflows.md
精简的当前 API 签名和默认值	references/api_reference.md
现有环境的升级说明	references/migration_to_0_24.md
官方文档、规范和发布来源	references/sources.md

常见故障模式

• 将数值型 VariantFile.fetch() 坐标视为 1 基坐标
• 在需要 BGZF 加 tabix/CSI 时使用普通 gzip
• 未使用索引就按区域调用 fetch
• 假设 fetch() 会包含未定位的未比对序列
• 忘记为精确的 pileup 区间设置 truncate=True
• 忽略 pileup 的默认设置，例如碱基质量 13 和深度上限 8000
• 在活动迭代器或线程之间共用同一个文件句柄
• 未使用完全匹配的参考序列解码 CRAM
• 在输出标头中声明新的 VCF 字段之前就为其赋值
• 将大量 samtools/bcftools 输出捕获到内存中
• 对 indel 或符号等位基因使用 SNP 碱基计数方法

相关技能

• polars-bio：用于基因组区间操作和在 Polars DataFrame 上进行生物信息学文件 I/O 的 Python 库 polars-bio，基于 Arrow 和 DataFusion 构建。
• tiledbvcf：使用 tiledbvcf Python API 和 CLI 在 TileDB-VCF 数据集中存储和查询基因组变异数据（create、store、export、list、stat）。
• gtars：检查并规划使用 Gtars 的工作。Gtars 是用于基因组区间的 Rust/Python/CLI 工具包：BED RegionSet 集合代数（reduce、setdiff、intersec…