alterlab-pysam

Pysam

概述

Pysam 是一个 Python 模块，用于读取、处理和写入基因组数据集。它通过符合 Python 风格的 htslib 接口读写 SAM/BAM/CRAM 比对文件、VCF/BCF 变异文件以及 FASTA/FASTQ 序列。查询 tabix 索引文件，执行 pileup 覆盖度分析，并运行 samtools/bcftools 命令。

何时使用此技能

在以下情况应使用此技能：

• 处理测序比对文件（BAM/CRAM）
• 分析遗传变异（VCF/BCF）
• 提取参考序列或基因区域
• 处理原始测序数据（FASTQ）
• 计算覆盖度或测序深度
• 实现生物信息学分析流水线
• 对测序数据进行质量控制
• 执行变异检测和注释工作流

不适用的情况

场景	请改用
运行完整的变异检测流水线（比对 -> GATK/DeepVariant -> 注释）	alterlab-nf-core-sarek
将多个样本的变异存储为可查询数组并进行查询	alterlab-tiledbvcf
覆盖度轨迹、生成 bigWig，以及类似 deepTools 的 BAM 汇总	alterlab-deeptools
序列记录解析、翻译、比对对象（不涉及 BAM/VCF）	alterlab-biopython
从 RNA-seq FASTQ 进行转录本级定量	alterlab-rnaseq-quant

快速入门

安装

uv pip install pysam

基本示例

读取比对文件：

import pysam

# Open BAM file and fetch reads in region
samfile = pysam.AlignmentFile("example.bam", "rb")
for read in samfile.fetch("chr1", 1000, 2000):
    print(f"{read.query_name}: {read.reference_start}")
samfile.close()

读取变异文件：

# Open VCF file and iterate variants
vcf = pysam.VariantFile("variants.vcf")
for variant in vcf:
    print(f"{variant.chrom}:{variant.pos} {variant.ref}>{variant.alts}")
vcf.close()

查询参考序列：

# Open FASTA and extract sequence
fasta = pysam.FastaFile("reference.fasta")
sequence = fasta.fetch("chr1", 1000, 2000)
print(sequence)
fasta.close()

核心功能

1. 比对文件操作（SAM/BAM/CRAM）

使用 AlignmentFile 类处理已比对的测序读段。适用于分析比对结果、计算覆盖度、提取读段或进行质量控制。

常见操作：

• 打开并读取 BAM/SAM/CRAM 文件
• 获取特定基因组区域的读段
• 按比对质量、标志或其他条件筛选读段
• 写入筛选或修改后的比对结果
• 计算覆盖度统计信息
• 执行 pileup 分析（逐碱基覆盖度）
• 访问读段序列、质量分数和比对信息

**参考：**有关以下内容的详细文档，请参阅 references/alignment_files.md：

• 打开并读取比对文件
• AlignedSegment 属性和方法
• 使用 fetch() 按区域获取数据
• 用于覆盖度分析的 pileup 分析
• 创建并写入 BAM 文件
• 坐标系统和索引
• 性能优化建议

2. 变异文件操作（VCF/BCF）

使用 VariantFile 类处理变异检测流水线生成的遗传变异。适用于变异分析、筛选、注释或群体遗传学研究。

常见操作：

• 读取和写入 VCF/BCF 文件
• 查询特定区域的变异
• 访问变异信息（位置、等位基因、质量）
• 提取样本的基因型数据
• 按质量、等位基因频率或其他条件筛选变异
• 使用附加信息注释变异
• 对样本或区域进行子集筛选

**参考：**有关以下内容的详细文档，请参阅 references/variant_files.md：

• 打开并读取变异文件
• VariantRecord 属性和方法
• 访问 INFO 和 FORMAT 字段
• 使用基因型和样本
• 创建并写入 VCF 文件
• 筛选变异并生成子集
• 多样本 VCF 操作

3. 序列文件操作（FASTA/FASTQ）

使用 FastaFile 随机访问参考序列，使用 FastxFile 读取原始测序数据。适用于提取基因序列、根据参考序列验证变异，或处理原始读段。

常见操作：

• 按基因组坐标查询参考序列
• 提取感兴趣基因或区域的序列
• 读取包含质量分数的 FASTQ 文件
• 验证变异的参考等位基因
• 计算序列统计信息
• 按质量或长度筛选读段
• 在 FASTA 和 FASTQ 格式之间转换

参考： 详细文档请参阅 references/sequence_files.md，内容包括：

• FASTA 文件访问和索引
• 按区域提取序列
• 处理基因的反向互补序列
• 顺序读取 FASTQ 文件
• 质量分数转换和筛选
• 使用 tabix 索引文件（BED、GTF、GFF）
• 常见序列处理模式

4. 综合生物信息学工作流

Pysam 擅长整合多种文件类型，进行全面的基因组分析。常见工作流会结合比对文件、变异文件和参考序列。

常见工作流：

• 计算特定区域的覆盖度统计信息
• 根据比对读段验证变异
• 使用覆盖度信息注释变异
• 提取变异位置周围的序列
• 根据多个条件筛选比对结果或变异
• 生成用于可视化的覆盖度轨迹
• 对多种数据类型执行质量控制

参考： 详细示例请参阅 references/common_workflows.md，内容包括：

• 质量控制工作流（BAM 统计信息、参考序列一致性）
• 覆盖度分析（逐碱基覆盖度、低覆盖度检测）
• 变异分析（注释、按读段支持情况筛选）
• 序列提取（变异上下文、基因序列）
• 读段筛选和子集提取
• 集成模式（BAM+VCF、VCF+BED 等）
• 复杂工作流的性能优化

核心概念

坐标系统

重要： Pysam 使用 从 0 开始、左闭右开 的坐标（Python 约定）：

• 起始位置从 0 开始（第一个碱基的位置为 0）
• 结束位置不包含在范围内
• 区域 1000-2000 包含碱基 1000-1999（共 1000 个碱基）

例外： fetch() 中的区域字符串遵循 samtools 约定（从 1 开始）：

samfile.fetch("chr1", 999, 2000)      # 0-based: positions 999-1999
samfile.fetch("chr1:1000-2000")       # 1-based string: positions 1000-2000

VCF 文件： 文件格式使用从 1 开始的坐标，但 VariantRecord.start 从 0 开始。

索引要求

随机访问特定基因组区域需要索引文件：

• BAM 文件：需要 .bai 索引（使用 pysam.index() 创建）
• CRAM 文件：需要 .crai 索引
• FASTA 文件：需要 .fai 索引（使用 pysam.faidx() 创建）
• VCF.gz 文件：需要 .tbi tabix 索引（使用 pysam.tabix_index() 创建）
• BCF 文件：需要 .csi 索引

没有索引时，可使用 fetch(until_eof=True) 顺序读取。

文件模式

打开文件时指定格式：

• "rb" - 读取 BAM（二进制）
• "r" - 读取 SAM（文本）
• "rc" - 读取 CRAM
• "wb" - 写入 BAM
• "w" - 写入 SAM
• "wc" - 写入 CRAM

CRAM：pysam 0.24 / htslib 1.22 以来的两项行为变化

这些变化不易察觉，因此在归咎于数据之前，请先检查它们。

1. 默认输出 CRAM 版本现在是 3.1，而非 3.0。 使用 "wc" 写入的文件可能无法被较旧版本的 samtools 或固定使用较旧 htslib 版本的下游工具读取。如果无法控制使用方，请明确写入 3.0：

   out = pysam.AlignmentFile("out.cram", "wc", header=src.header,
                             reference_filename="ref.fa",
                             format_options=["version=3.0"])
   

   format_options 接受 str 列表（较旧版本的 pysam 接受 bytes）。

2. CRAM 参考序列不再自动从 EBI 获取。 CRAM 会相对于参考序列存储读段，因此现在读取没有匹配 FASTA 的 CRAM 会失败，而不是静默下载。打开文件时传入 reference_filename=，或将 REF_PATH / REF_CACHE 环境变量设置为本地缓存。

性能注意事项

1. 随机访问操作始终使用已建立索引的文件
2. 按列分析时使用 pileup()，而不是重复执行 fetch 操作
3. 计数时使用 count()，而不是手动迭代并计数
4. 分析相互独立的基因组区域时并行处理
5. 显式关闭文件以释放资源
6. 在没有索引的情况下进行顺序处理时使用 until_eof=True
7. 除非有必要，否则避免使用多个迭代器（如有需要，使用 multiple_iterators=True）

常见问题

1. **坐标混淆：**请记住，不同上下文中使用的坐标系统可能是从 0 或从 1 开始的
2. **缺少索引：**许多操作需要索引文件——请先创建索引
3. 部分重叠：fetch() 会返回与区域边界重叠的 reads，而不仅仅是完全包含在区域内的 reads
4. 未限定范围的 pileup：pileup(chrom, start, stop) 会为重叠 reads 覆盖的每个位置生成列，而不只是 start..stop 范围内的位置——传入 truncate=True 可限制为请求的窗口
5. 碱基质量被静默忽略：count_coverage() 默认使用 quality_threshold=15，因此会忽略低质量碱基；设置 quality_threshold=0 可统计所有碱基
6. **迭代器作用域：**请保持 pileup 迭代器引用有效，以避免出现 "PileupProxy accessed after iterator finished" 错误
7. **编辑质量分数：**更改 query_sequence 后，无法原地修改 query_qualities——请先创建副本
8. **流式处理限制：**流式处理仅支持 stdin/stdout，不支持任意 Python 文件对象
9. **线程安全：**虽然 I/O 期间会释放 GIL，但尚未全面验证线程安全性
10. **CRAM 未提供参考序列：**请参阅上面的 CRAM 说明——现在必须使用 reference_filename= 或 REF_PATH

命令行工具

Pysam 提供了对 samtools 和 bcftools 命令的访问：

# Sort BAM file
pysam.samtools.sort("-o", "sorted.bam", "input.bam")

# Index BAM
pysam.samtools.index("sorted.bam")

# View specific region
pysam.samtools.view("-b", "-o", "region.bam", "input.bam", "chr1:1000-2000")

# BCF tools
pysam.bcftools.view("-O", "z", "-o", "output.vcf.gz", "input.vcf")

错误处理：

try:
    pysam.samtools.sort("-o", "output.bam", "input.bam")
except pysam.SamtoolsError as e:
    print(f"Error: {e}")

资源

references/

各项主要功能的详细文档：

• alignment_files.md - SAM/BAM/CRAM 操作完整指南，包括 AlignmentFile 类、AlignedSegment 属性、fetch 操作、pileup 分析和比对数据写入

• variant_files.md - VCF/BCF 操作完整指南，包括 VariantFile 类、VariantRecord 属性、基因型处理、INFO/FORMAT 字段和多样本操作

• sequence_files.md - FASTA/FASTQ 操作完整指南，包括 FastaFile 和 FastxFile 类、序列提取、质量分数处理，以及基于 tabix 索引的文件访问

• common_workflows.md - 实用的生物信息学工作流示例，展示如何整合多种文件类型，包括质量控制、覆盖度分析、变异验证和序列提取

获取帮助

有关特定操作的详细信息，请参阅相应的参考文档：

• 处理 BAM 文件或计算覆盖度 → alignment_files.md
• 分析变异或基因型 → variant_files.md
• 提取序列或处理 FASTQ → sequence_files.md
• 整合多种文件类型的复杂工作流 → common_workflows.md

官方文档：https://pysam.readthedocs.io/