企业官网建设流程全解析

pysam 坐标与索引实战指南：从 0-based 半开区间到 BAI/TBI/CSI 安全索引

【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000+ scientists worldwide. 165 ready-to-use validated skills plus 100+ scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills

本指南以 scientific-agent-skills 仓库中 pysam 坐标与索引参考文档 为核心骨架，面向所有使用 Python 处理 SAM/BAM/CRAM、VCF/BCF、FASTA 与 tabix 表格的开发者和 AI Agent。坐标写错与索引过期是"看起来合理但结果错误"的基因组学结果最常见两大根因，读完本文你将掌握 pysam 的坐标契约、各格式坐标换算表、索引类型选择（BAI/TBI 与 CSI）以及一套可直接落地的高质量边界测试方案，并学会通过仓库内置脚本与源码验证每一步结论。

适用版本说明：本文所有 API 行为与示例均以 pysam 0.24.0（内嵌 HTSlib/samtools/bcftools 1.23.1）为基准，与仓库 SKILL.md 中声明的一致性要求一致，安装时建议使用uv pip install "pysam==0.24.0"固定版本。

Pysam 坐标规则：两条必须刻进脑子的约定

pysam 的 Python API 对数值型参数与属性统一采用0-based、半开区间（half-open）：

[start, stop)

第一个碱基是0，start包含、stop不包含，区间长度为stop - start。唯一的主要例外是文本形式的 samtools 风格 region 字符串，它采用1-based、闭区间（inclusive）：

chr1:100-199

下面两种写法指向完全相同的 100 个碱基：

file.fetch("chr1", 99, 199) file.fetch(region="chr1:100-199")

该规则统一适用于以下 API，无一例外：

• AlignmentFile.fetch()、count()、count_coverage()、pileup()
• VariantFile.fetch()
• FastaFile.fetch()
• TabixFile.fetch()

最常见的坑：不要把数值形式的VariantFile.fetch()参数当作 VCF 文本坐标（VCF 的POS是 1-based）直接传入。仓库 variant_files.md 中专门强调："Do not subtract one from a numericstartpassed tofetch()"——数值参数已经是 0-based，再减一就是双重换算错误。

从源码与脚本层面也能印证这一点：inspect_hts.py 中对各类文件统一通过pysam.AlignmentFile/VariantFile/FastaFile/TabixFile读取 contig 与长度信息，其contig_records()直接使用references与lengths元组，全程以 Python 原生坐标为准；而 variant_summary.py 的--region参数则在帮助文本中明确标注为"1-based inclusive samtools region, for example chr1:1-1000000"，并且在内部通过variants.fetch(region=region)传递——这正是数值坐标与 region 字符串两条路径并存的典型工程实践。

格式坐标换算表：一张表解决跨格式对齐

不同生物信息学格式的坐标约定各不相同，参考文档给出了一张完整换算表，这是跨 BAM/VCF/BED/GFF 集成时最容易出错的地方：

对于 VCF 结构变异、符号等位基因（symbolic allele）、断点（breakend）以及带INFO/END的记录，务必使用VariantRecord.start与VariantRecord.stop，而不是用len(REF)反推区间——结构变异记录的 REF 长度与实际跨度往往不一致。这一点在 variant_files.md 中有更细致的展开：VariantRecord的pos是 1-based 位置、start是 0-based 闭区间起点、stop是 0-based 开区间终点、rlen是参考跨度，多等位基因、<DEL>/<DUP>/<INS>等符号等位基因与*跨缺失都需要单独处理，绝不能用"单碱基计数法"套用到 indel 或符号等位基因上。

单碱基位置换算

一个 1-based 的位置p，对应的 Python 单碱基区间为：

start = p - 1 stop = p

对于 VCF 记录，可直接验证并提取该位置的参考碱基：

assert record.start == record.pos - 1 base = fasta.fetch(record.contig, record.start, record.start + 1)

其中record.start是 0-based 的，而record.pos是 1-based 的，二者恒差 1。这段代码同时是"用 VCF 坐标去 FASTA 里取碱基"的标准姿势——注意FastaFile.fetch()的数值参数同样是 0-based 半开区间。

重叠查询 vs 完全包含

需要牢记：region fetch 本质上是重叠查询（overlap query）。一条比对记录或变异可以在请求区间之前就开始、之后才结束，但仍然"重叠"该区间。因此 fetch 返回的记录并不保证完全落在区间内。

如果业务逻辑要求完全包含，需要自己写过滤：

def fully_contained(read, start: int, stop: int) -> bool: return ( read.reference_start is not None and read.reference_end is not None and read.reference_start >= start and read.reference_end <= stop )

注意reference_end是 0-based 排他终点（由 CIGAR 推导，见 alignment_files.md），未比对或无 CIGAR 的记录其reference_start/reference_end可能是None或哨兵值，所以先做非空判断。

对于点突变类逻辑，还必须明确定义 deletion、reference skip（CIGARN）、符号等位基因和断点各自的"重叠"含义——例如一个 deletion 记录"覆盖"某个点，是指其起点、还是包含删掉的那段序列？参考文档提醒"Define exactly what 'overlap' means"，这是决定分析正确性的语义决策点。

pileup()还有一个额外的陷阱：如果不设置truncate=True，由于读段跨区间边界，pileup 可能输出请求区间之外的列。参考文档与 alignment_files.md 都一致强调，需要精确区间 pileup 时必须显式传入truncate=True。仓库 SKILL.md 给出了完整的精确区间 pileup 示例：

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())

同时注意 pileup 的默认参数极具误导性：默认min_base_quality为 13、默认max_depth为 8000，配对重叠检测与孤儿读段过滤默认开启——这些默认值会悄悄改变统计结果，正式分析时必须显式声明过滤语义。

解析器（Parser）坐标：三种对象三种约定

pysam 的解析器对象（用于TabixFile(parser=...)）各自做了坐标归一化：

• asBed().start/.end：0-based，半开（与 BED 文件本身一致）
• asGTF().start/.end：以 Python 坐标约定暴露（即 0-based）
• asVCF().pos：解析器特有的轻量字段，完整 VCF 记录语义请使用VariantFile

sequence_files.md 给出完整的解析器清单：asTuple()（类元组字段）、asBed()（BED 字段，0-based start/end）、asGTF()（GTF/GFF 类字段与属性）、asVCF()（轻量 tabix VCF 解析器），并强调"对于完整的 VCF 语义请使用VariantFile，而不是TabixFile(asVCF())"。

自定义 tabix 索引时，存在三套完全不同的概念，参考文档专门点名：

• seq_col、start_col、end_col是0-based 的列索引（Python 下标）
• 文件中的坐标默认按1-based编码，除非设置zerobased=True
• 之后TabixFile.fetch()的数值查询坐标仍然是 0-based

pysam.tabix_index( "custom.tsv.gz", seq_col=0, start_col=1, end_col=2, zerobased=True, )

三者分别是：Python 列下标、存储表格中的坐标编码、查询时的坐标编码。把这三者混为一谈是自定义 tabix 场景最常见的 bug 来源。

Contig 身份核对：坐标换算解决不了的问题

坐标换算只解决"位置数字"的语义，解决不了 contig 命名不一致。参考文档列出的常见不匹配包括：

• chr1与1（是否带chr前缀）
• 线粒体命名差异（chrM、MT、M）
• 交替位点（alt loci）与 decoy contig
• 组装版本差异（例如 GRCh37 与 GRCh38）
• contig 顺序与长度差异

一段可复用的核对逻辑：

alignment_contigs = dict(zip(bam.references, bam.lengths)) fasta_contigs = dict(zip(fasta.references, fasta.lengths)) shared = alignment_contigs.keys() & fasta_contigs.keys() length_mismatches = { name: (alignment_contigs[name], fasta_contigs[name]) for name in shared if alignment_contigs[name] != fasta_contigs[name] }

参考文档给出明确的操作准则：不要跨任意组装静默地添加或剥离chr前缀，必须使用经过人工审阅的显式映射。这条规则与 CRAM 场景一脉相承——cram_and_performance.md 指出 CRAM 的参考身份由@SQ头中的M5MD5 值与UR字段表示，contig 名称本身不足以确认参考一致性，需核对 contig 长度与 M5 标签。

仓库内置的 inspect_hts.py 恰好提供了现成的元数据核对入口：它对 alignment/variant/fasta/fastx/tabix 五类文件输出contigs（名称与长度）、has_index、sort_order、reference_count等 JSON 摘要，CRAM 输入强制要求--reference，是分析前核对 contig 身份与索引状态的理想第一步：

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

索引矩阵：什么数据配什么索引

参考文档给出了完整的索引矩阵：

几个关键推论：

1. 索引是特定文件的视图。数据文件一旦变化（重写、追加、修改），必须重建索引。过期索引可能报错，也可能静默返回错误或不完整的区域结果——后一种情况更难排查。
2. FASTA 的.fai是"布局索引"而非坐标排序索引，BGZF 压缩的 FASTA 还需要.gzi压缩偏移索引；普通 gzip 不适合做索引随机访问（详见 sequence_files.md）。
3. SAM、普通（未压缩的）VCF、FASTQ 在这些 API 下没有随机访问，只能顺序扫描。

BAI/TBI 与 CSI：大基因组必须选对

传统的 BAI（BAM 索引）与标准 TBI 索引最大坐标上限约为2^29（512 Mi 碱基）。这对人类基因组足够，但对某些植物、动物和合成参考基因组就不够用了。CSI 是参数化的索引格式，支持更大的坐标范围。

创建 BAM 的 CSI 索引：

import pysam pysam.index("-c", "large-reference.bam", catch_stdout=False)

创建 tabix 的 CSI 索引：

pysam.tabix_index( "large-reference.bed.gz", preset="bed", csi=True, min_shift=14, )

创建 VCF/BCF 的 CSI 索引：

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

参考文档的决策建议：当参考基因组大小未知或可能很大时优先使用 CSI，同时确认下游工具支持 CSI（部分旧工具只认 BAI/TBI）。这条建议在仓库的写规则里被提升为硬性规范：SKILL.md 的 Writing Rules 明确要求 "Use CSI rather than BAI/TBI when references or coordinates exceed legacy index limits"。

这里需要特别说明catch_stdout=False的用途：pysam 的 samtools/bcftools 命令分发器默认会捕获 stdout，对于大型或二进制输出，应使用工具自身的-o选项配合catch_stdout=False，或使用save_stdout=...，避免把完整输出塞进 Python 内存（SKILL.md Wrapped samtools and bcftools 一节）。

先排序再索引：索引不会替你排序

索引操作本身不会排序记录。BAM 必须按坐标序排好才能建索引；VCF 必须按 contig/位置序排好；tabix 表必须先排序后tabix_index()——Python 的tabix_index()函数不会校验排序状态。

BAM 的排序 + 建索引：

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

VCF 的排序 + 建索引：

import pysam.bcftools pysam.bcftools.sort( "-Oz", "-o", "sorted.vcf.gz", "input.vcf", catch_stdout=False, ) pysam.bcftools.index( "--csi", "sorted.vcf.gz", catch_stdout=False, )

注意传递给命令分发器的每个命令行 token 都要作为独立字符串传入（例如"-@", "4"），这也是 SKILL.md 强调的规范；绝不可以通过拆分不可信的 shell 命令字符串来拼装分发器参数。

安全创建 Tabix 索引：两条不可省略的约定

推荐分两步走，压缩与建索引分离：

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

这里有一个隐蔽的破坏性行为：直接调用tabix_index("regions.bed")可能会自动创建regions.bed.gz并删除原始文件。如果依赖这种一步式路径，必须设置keep_original=True。sequence_files.md 对此同样给出了警告："Lettingtabix_index()remove the uncompressed source unexpectedly" 被列在常见陷阱中。

此外，默认不要传force=True。已存在的输出应该触发人工审查，而不是被静默替换。仓库的所有脚本（如 inspect_hts.py、variant_summary.py）在写输出时都拒绝覆盖已有文件，例如 variant_summary.py 的write_report()使用destination.open("x", encoding="utf-8")独占创建模式，并在 main 中校验--output不得覆盖输入文件——"写新路径、拒绝静默覆盖"是仓库一致遵循的工程原则。

创建索引前还必须确认输入要求：按 contig 与坐标排序、使用 BGZF 而非普通 gzip、选择正确的 preset（常用bed、gff、sam、vcf），preset 决定了列含义与坐标约定。

非标准与远程索引位置

当索引文件不在标准位置、或自动推导 URL 不可靠时，通过index_filename显式传入索引路径：

with pysam.AlignmentFile( "sample.bam", "rb", index_filename="indexes/sample.csi", ) as bam: ... with pysam.VariantFile( "cohort.vcf.gz", index_filename="indexes/cohort.vcf.gz.csi", ) as variants: ...

远程随机访问（如 HTTP(S) 上的 BAM/VCF）还额外依赖：

• 支持相关网络/插件功能的 HTSlib 构建
• 可达的索引文件
• 服务器支持 range 请求
• 数据与索引 URL 的稳定性

cram_and_performance.md 给出了远程访问的更完整清单：先发一个小范围已知 region 查询验证、核对精确的索引 URL、检查内容版本化与不可变性、确认重试与超时行为、估算访问模式产生的请求数量——并提醒"大量零碎随机查询可能比一次性下载或暂存文件更慢更贵"，且绝不要将 bearer token 或密钥写进提交的 URL。

参考文档还提示：在远程或 CRAM 访问之前，先阅读 cram_and_performance.md——CRAM 解码需要匹配的参考 FASTA（pysam 0.24 默认不再联系 EBI 参考服务器），远程 CRAM 叠加了参考获取与网络两重不确定性。

索引检查：不要只检查文件是否存在

对比对文件，显式校验索引状态：

with pysam.AlignmentFile("sample.bam", "rb") as bam: if not bam.has_index(): raise ValueError("random access requires a BAM/CRAM index") bam.check_index()

check_index()在索引缺失、不可用或文件已关闭时抛出异常（详见 alignment_files.md），适合"随机访问是硬需求"的场景。

对于变异文件与 tabix 文件，构造器会自动打开发现到的索引；当索引不可用时，region fetch 会直接失败。参考文档给出的核心建议是：重新打开输出文件并对已知区域做实际查询，而不是仅仅检查索引文件是否存在——文件存在不等于内容与数据匹配、不等于索引未过期。

仓库脚本的实现可作为参考：inspect_hts.py 对 alignment 文件通过alignments.has_index()判断索引，并在有索引时调用get_index_statistics()输出 mapped/unmapped 统计；对 variant 文件则通过variants.index is not None判断。注意这些是"索引层面的检查"，get_index_statistics()返回的是索引中记录的统计信息，不是重新扫描每条记录的结果（alignment_files.md 明确说明）。

边界测试：把正确性变成可回归的检查项

对于重要 pipeline，参考文档建议至少测试以下边界：

• contig 的第一个碱基
• 区间精确的起点与终点
• 横跨查询边界的记录（overlap 语义）
• 零长度或非法区间
• contig 末端
• 当预期使用 CSI 时，超过 512 Mi 碱基的高坐标
• 缺失与别名 contig
• region 字符串与数值两种写法的等价性

其中最有价值的可执行不变量是数值坐标与 region 字符串的等价性检查：

numeric = list(file.fetch(contig, start, stop)) region = f"{contig}:{start + 1}-{stop}" textual = list(file.fetch(region=region))

对同一个已索引文件、同一个非空有效区间，这两次查询应当选出等价的记录集合。这段代码把"坐标契约"变成了可断言的测试，是任何 indexed 文件读取逻辑的黄金回归用例。结合仓库的测试组织方式（每个 skill 在 tests 目录下有对应测试目录），建议将这类断言固化进项目的自动化测试中。

实战总结：一份可复用的坐标与索引检查清单

综合参考文档与仓库源码，推荐在每次基因组分析前执行以下核对流程：

1. 确认格式、压缩、排序与索引状态：用python skills/pysam/scripts/inspect_hts.py <file>快速导出元数据 JSON，确认has_index、sort_order、contig 列表与长度。
2. 统一坐标心智模型：Python API 数值参数一律 0-based 半开；samtools region 字符串一律 1-based 闭区间；BED 原样、GFF/GTF 左端点减一、VCFPOS-1起算。
3. 结构变异不要用len(REF)反推：使用VariantRecord.start/.stop；符号等位基因与断点单独定义"重叠"语义。
4. pileup 精确区间务必truncate=True，并显式声明min_base_quality、max_depth、stepper等过滤参数。
5. 核对 contig 身份：对比references/lengths，禁止跨组装静默增删chr前缀。
6. 排序在前、建索引在后；tabix 两步走（tabix_compress+tabix_index），不传force=True；数据变更后重建索引。
7. 大基因组选 CSI：超过2^29坐标上限或参考大小未知时使用 CSI，并确认下游工具兼容。
8. 用真实 region 查询验证，配合数值/region 等价性断言做边界回归测试。

这套清单覆盖了参考文档的全部核心要点，并且每一环都能在仓库中找到对应的源码或脚本佐证：坐标契约见 SKILL.md 与 coordinates_and_indexing.md，pileup 细节见 alignment_files.md，VCF 记录语义见 variant_files.md，tabix 与 FASTA 见 sequence_files.md，远程/CRAM/线程见 cram_and_performance.md，元数据检查脚本见 inspect_hts.py 与 variant_summary.py。把坐标与索引当作一等公民对待，是杜绝"貌似合理实则错误"基因组结果的起点。

【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000+ scientists worldwide. 165 ready-to-use validated skills plus 100+ scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills