关于 genoray-api
genoray 是以 NumPy 为核心的基因组变异查询工具包,支持 VCF/BCF、PLINK 2 PGEN 和稀疏变异存储,面向需要查询区间、转换数据、筛选变异或分析突变特征的开发者与研究人员。
主要功能
- 01读取 VCF/BCF 和 PGEN
- 02查询基因组区间
- 03转换变异数据为 SVAR2
- 04使用表达式筛选变异
- 05提取指定字段和剂量
- 06分析突变特征
适用场景
- 01查询多个基因组区域的变异
- 02将队列转换为稀疏 SVAR2 存储
- 03筛选并提取变异字段
- 04分类突变并拟合特征
安装
复制下面的命令,在终端中运行。
npx skills add https:/ / github.com/ d-laub/ genoray --skill genoray-api安装后在 Claude Code 中输入/genoray-api确认已加载。
详细介绍
genoray 公共 API
genoray 是一个以 NumPy 为先的范围查询层,支持 VCF/BCF(cyvcf2)、PGEN(pgenlib)和稀疏 memmap 格式(SparseVar / .svar)。
公共接口
import genoray 恰好会公开以下内容:
genoray.PGEN— PLINK 2 PGEN 读取器genoray.Reference— 带索引的 FASTA 参考基因组读取器genoray.VCF— VCF/BCF 读取器genoray.Filter— VCF 筛选值对象,将 cyvcf2 记录谓词(record)与匹配的.gvipolars 表达式(expr)打包在一起genoray.SparseVar— 稀疏.svar读取器/写入器genoray.SparseVar2— 新一代稀疏变异存储(通过from_vcf将 VCF/BCF 转换为 SVAR2(支持regions=/samples=/merge_overlapping=/regions_overlap=),通过from_pgen将 PLINK2 PGEN 转换为 SVAR2;通过from_vcf_list中的原生 k 路合并,将 N 个单样本 VCF/BCF 合并为一个 SVAR2 存储(支持reference/no_reference,与from_vcf相同;缺失位点填充为纯合参考型;支持regions=/merge_overlapping=/regions_overlap=,但不支持samples=——队列由文件集构成);通过from_svar1原生迁移 SVAR1(SparseVar)→ SVAR2(不读取 VCF/htslib;仅支持双等位基因 SVAR1;支持regions=/samples=/merge_overlapping=/regions_overlap=,与from_vcf/from_pgen相同);通过decode/region_counts/read_ranges进行范围查询;通过annotate_mutations/mutation_matrix/assign_signatures支持突变特征(SBS96/DBS78/ID83),也可以在写入时通过from_vcf(signatures=True)/from_pgen(signatures=True)/from_svar1(signatures=True)进行分类;写入时通过from_vcf(info_fields=, format_fields=)/from_vcf_list(info_fields=, format_fields=)提取标量数值型 INFO/FORMAT 字段(from_vcf_list的 INFO 字段采用首个携带者优先合并,FORMAT 字段按样本合并);from_pgen则通过dosages=Sequence[DosageField]将每个样本的剂量轨迹存储为 FORMAT 字段(可从硬调用.pgen本身读取,使用source="self",也可从单独的.pgen读取)——它仍不支持info_fields=/format_fields=(PGEN 没有 VCF INFO/FORMAT 字段);from_svar1可通过fields=选择性传递 SVAR1 中已有的字段(默认值None= 全部,[]= 无,或指定名称子集)——读取时通过fields=/with_fields/available_fields选择启用,并附加到decode的结果中)genoray.InfoField/genoray.FormatField— 冻结数据类(name、dtype=None、default=None),用于配置SparseVar2.from_vcf的单个 INFO/FORMAT 字段;直接使用str名称时,则采用推断出的默认值genoray.DosageField— 冻结数据类(name="dosage"、source="self"|Path、dtype="f16"|"f32"="f32"、default=None),用于配置SparseVar2.from_pgen的 PGEN 剂量 FORMAT 字段genoray.Tuning— 包含六个显式调度参数(concurrent_chroms、reader_workers、overshard、dense_cap、merge_threads、sample_interval,均为int | None)的冻结数据类,通过tuning=传递给每个SparseVar2.from_*写入器;用于替代已移除的旧环境变量配置(参见下文“SparseVar2 — quick reference”中的“Tuning”,以及docs/source/svar.md中的“migrating from environment variables”)——没有任何环境变量可以配置 genoraygenoray.exprs— 用于.gvi索引的 polars 筛选表达式genoray.cosmic_signatures— 获取/缓存 COSMIC 参考特征genoray.fit_signatures— 突变目录特征重拟合;strategy=选择Forward()(默认,稀疏前向选择)或Spa()(忠实于 SigProfilerAssignment 的反向消除)genoray.Forward/genoray.Spa—fit_signatures策略数据类(冻结)genoray.Strategy/genoray.Criterion/genoray.Metric/genoray.ActivityScale— 相关类型别名;genoray.SPA_CONNECTED_GROUPS包含Spa默认的四组共现特征组
其他内容均不属于公共接口。任何以下划线 _ 开头的内容(例如 genoray._vcf)都是内部实现——请勿在用户代码中导入。
详细信息的查阅位置
请优先查阅以下内容,不要靠猜测:
docs/source/index.md— 叙述性导览,包含完整示例(VCF、PGEN、筛选、分块)docs/source/svar.md— SparseVar 用法genoray/__init__.py— 确认公共接口genoray/_vcf.py—VCF类:构造函数、read、chunk、类顶部附近的模式常量;get_record_info(contig=None, start=None, end=None, fields=None, info=None, lazy=False)— 获取指定范围或整个文件中的非 FORMAT 记录级字段(包括 INFO),返回pl.DataFrame(lazy=True时返回pl.LazyFrame)。INFO 字段会作为顶层列返回,而不是嵌套的INFO结构体,列名使用表头中各 ID 的拼写形式。info=名称与表头声明的 INFO ID 不区分大小写地匹配;None表示全部,[]表示不选,表头未声明的名称会引发ValueError(此前会静默丢弃该名称以及所有其他 INFO 字段,#139)。若 INFO ID 与非 INFO 列名冲突(例如INFO/QUAL),也会引发错误,而不会产生含义不明的扁平化结果。genoray/_pgen.py—PGEN类:构造函数、read、chunk、read_ranges、chunk_ranges、类顶部附近的模式常量;var_records(contig, starts=0, ends=POS_MAX)— 返回与对相同范围执行读取时相同的变体标识,结果为pl.DataFrame,包含range(u32,表示该行所属的查询范围)、CHROM、POS(1 起始)、REF、ALT,与read/read_ranges的变体轴严格 1:1 对齐,并遵循读取器的filter。此处的ALT始终为list[str],无论.gvi中的存储形式如何;CHROM为文件中实际使用的 contig 名称(可能与查询的contig不同)。若产物必须比生成它的 PGEN 存活得更久,请使用此方法,而不是var_idxs加私有的_indexgenoray/_svar.py—SparseVar:__init__、from_vcf、from_pgen、read_ranges、read_ranges_with_length(contig, starts=0, ends=POS_MAX, samples=None)(保证长度的范围读取;返回类型与read_ranges相同,即Ragged或添加字段后的记录)、with_fields、annotate_mutations、mutation_matrix、assign_signatures、annotate_with_gtf(gtf, level_filter=1, write_back=True, *, strand_encoding=None, codon_null_token=None)(GTF CDS 注释入口,返回包含varID/gene_id/strand/codon_pos的pl.DataFrame)、cache_afs()(计算AF列并持久化到.gvi索引;返回None)genoray/_svar2.py—SparseVar2:__init__(path, *, fields=None)、with_fields(fields)(使用所选字段创建同一存储上的新读取器)、available_fields(dict[str, StoredField],在__init__中设置)、from_vcf(VCF/BCF → SVAR2 转换入口;signatures=在写入期间进行分类,info_fields=/format_fields=在写入期间提取数值标量字段;支持regions=/samples=/merge_overlapping=/regions_overlap=)、from_pgen(PLINK2 PGEN → SVAR2 转换入口;仅支持二倍体,不支持ploidy=/info_fields=/format_fields=;dosages=Sequence[DosageField]将每个样本的剂量轨迹存为 FORMAT 字段,数据读取自硬调用.pgen本身(source="self")或另一个.pgen;与from_vcf一样支持regions=/samples=/merge_overlapping=/regions_overlap=)、from_vcf_list(N 个单样本 VCF/BCF → 通过原生 k 路合并写入一个 SVAR2 存储;sources接受Sequence/目录/清单,由模块级_resolve_vcf_sources解析;支持reference/no_reference(no_reference会跳过左对齐,因此跨文件合并要求输入预先标准化);支持info_fields=/format_fields=— INFO 按最先出现的携带者优先合并,FORMAT 保留为逐样本数据;与from_vcf一样支持regions=/merge_overlapping=/regions_overlap=,但不支持samples=— 队列由文件集合决定)、from_svar1(SVAR1(SparseVar)→ SVAR2 原生迁移入口;不读取 VCF/htslib,ploidy来自 SVAR1 元数据,仅支持双等位 SVAR1,不支持info_fields=/format_fields=(这些仅适用于 VCF)— 使用fields=Sequence[str] | None选择要迁移的 SVAR1 字段(默认None= 全部,[]= 不选,子集仅迁移指定名称,未知名称会引发ValueError);无法通过此方式选择mutcat,且它始终会被丢弃;与from_vcf/from_pgen一样支持regions=/samples=/merge_overlapping=/regions_overlap=,但区域筛选逐条记录进行,而不是预先缩小覆盖范围);n_samples/available_samples/contigs/ploidy元数据。读取/查询方法位于 mixin 中:genoray/_svar2_decode.py(decode— 为每个选定字段附加一个Ragged,region_counts)、genoray/_svar2_batch.py(公开的read_ranges;仅供 gvl 内部使用的_overlap_batch/_find_ranges/_gather_ranges),以及genoray/_svar2_mutcat.py(annotate_mutations、mutation_matrix、assign_signatures— COSMIC 突变特征工作流,仿照SparseVar,但由按 contig 划分的 Rust 侧车文件提供支持,而不是附加到.gvi的字段)genoray/_svar2_fields.py—InfoField/FormatField/DosageField数据类,以及FieldDtype和from_vcf(info_fields=, format_fields=)使用的表头/数据类型验证;_parse_cli_field_specs(内部使用 — 解析 bcftools 风格的INFO/x/FORMAT/x/FMT/xCLI 字段字符串,供genoray write vcf --fieldsCLI 使用);StoredField(冻结的数据类:name、category、dtype、default、key)是读取端清单条目类型,由SparseVar2.available_fields返回 — 不会在顶层genoray导出,只能通过该字典访问genoray/_cli/__main__.py—genorayCLI(index、write vcf/write pgen/write svar1(全部 → SVAR2)、顶层write-svar1(旧版 VCF/PGEN → SVAR1)、view/view svar1、concat、split)genoray/_signatures/—cosmic_signatures、fit_signatures、Forward、Spa、Strategy、Metric、ActivityScale、SPA_CONNECTED_GROUPS。fit_signatures(..., strategy=Forward()|Spa())用于选择重新拟合算法;Forward是 genoray 自有的贪心前向选择(默认算法),Spa则重新实现 SigProfilerAssignment 的cosmic_fit(从饱和特征集合开始进行后向剔除,然后进行添加-移除优化)genoray/_reference.py—Reference:from_path、fetch、contig_arraygenoray/exprs.py— 预构建筛选表达式的完整集合(目前有 7 个:is_snp、is_indel、is_biallelic、is_symbolic、is_breakend、is_imprecise、ILEN) 当签名、关键字参数或形状不明确时,请阅读源码中的文档字符串,而不要从基本原理推断。
通用约定
- 范围采用从 0 开始、左闭右开格式
[start, end)。 max_mem接受"4g"、"512m"、"2GB"等字符串——但SparseVar2.from_vcf、SparseVar2.from_pgen和SparseVar2.from_vcf_list的max_mem除外:它们都是整个进程级别的规划预算,而非单个分块的上限;在假设它们与此名称在其他地方的含义相同之前,请先查看下方「转换」中的对应条目(VCF.chunk/chunk_ranges、PGEN.chunk/chunk_ranges等情况下,该参数会直接限制单个分块)。- 染色体名称会自动规范化:无论文件采用哪种命名约定,
"chr1"和"1"都可用(ContigNormalizer)。 - 缺失基因型 =
-1(int)。缺失剂量 =np.nan(float32)。 - 默认倍性为 2;
SparseVar.from_vcf/from_pgen(以及genoray write-svar1)接受haploid=True/--haploid,会将单倍型 OR 合并为每个样本的单个单倍体调用,并在元数据中记录ploidy=1(适用于无相位的体细胞数据)。 - 所有返回数组均为 NumPy 数组;
mode用于选择要返回的数组。
样本访问器——规范名称及两种写法不同的原因
available_samples(list[str])是规范的「文件中所有样本」访问器——四种读取器(VCF、PGEN、SparseVar、SparseVar2)都提供此属性。
VCF 和 PGEN 还提供:
current_samples——当前选中的子集(只读属性)。set_samples(samples) -> Self——有状态调用,会就地修改读取器以选择子集(或传入None恢复所有样本),然后返回self。
SparseVar 和 SparseVar2 没有 current_samples/set_samples。
相反,每个读取方法(read_ranges、read_ranges_with_length 等)
都接受作为单次调用关键字参数的 samples=。
两种写法不同的原因(性能): 对 VCF/PGEN 进行样本子集选择的成本较高——因为这会重新初始化后端读取器——因此有意采用有状态的 set_samples() 调用,只调用一次并在多次读取中重复使用。对 SparseVar/SparseVar2 而言,子集选择的成本几乎为零(只是在已内存映射的数据上选择索引),因此将其作为轻量的单次调用关键字参数 samples= 提供,而非持久化的读取器状态。这种差异是有意设计,并非不一致——不要为了统一而「修正」其中一种写法。
模式常量——注意事项
模式是类属性,不是顶层名称:
genoray.VCF.Genos8 # not genoray.Genos8
genoray.PGEN.GenosPhasingDosages
要查看某个类支持的模式,请阅读 _vcf.py / _pgen.py 中的类定义(在文件顶部附近搜索 Genos)。
当某个模式包含多个数组时,返回元组的顺序与常量名称中的顺序一致。PGEN.GenosPhasingDosages 返回 (genos, phasing, dosages);VCF.Genos8Dosages 返回 (genos, dosages)。
VCF——快速参考
vcf = genoray.VCF(
"file.vcf.gz",
phasing=True, # constructor-time, not per-read
dosage_field="DS", # required to read dosages; FORMAT field with Number=A
filter=genoray.Filter(
record=lambda v: ..., # cyvcf2.Variant -> bool
expr=~genoray.exprs.is_symbolic, # matching .gvi index predicate
),
)
# Single range
arr = vcf.read("chr1", start=0, end=1_000_000, mode=genoray.VCF.Genos8)
# Chunked
for chunk in vcf.chunk("chr1", start=0, end=1_000_000,
max_mem="2g", mode=genoray.VCF.Genos8Dosages):
...
phasing=False时的形状:(samples, ploidy=2, variants)。phasing=True时的形状:(samples, ploidy+1=3, variants)——沿倍性轴的第 3 行为0(无相位)/1(有相位),与 cyvcf2 一致。- 剂量数组会移除倍性轴:
(samples, variants),dtype 为float32。 - VCF 特意不提供
read_ranges——基准测试表明它不会提升吞吐量。 read(out=...)仅适用于 VCF——传入预分配数组即可就地填充。PGEN 随机访问读取会分配新数组,不提供out=缓冲区。
PGEN — 快速参考
pgen = genoray.PGEN(
"hardcalls.pgen", # hardcalls live in the main path
dosage_path="dosages.pgen", # optional; defaults to the main path
filter=genoray.exprs.is_snp & genoray.exprs.is_biallelic,
)
重要提示:如果你有一个仅含剂量的 PGEN 和一个单独的硬调用 PGEN,硬调用文件应放在主路径中,剂量文件应放在 dosage_path 中。如果只传入一个路径,硬调用和剂量都会从该文件中读取(硬调用根据剂量阈值推断,详见 PLINK 2 文档)。
首次构建时,会在 PGEN 旁创建一个 .gvi 索引文件。不要删除它。
# Single range
genos = pgen.read("chr2", start=0, end=1000)
# Multiple ranges in one call (PGEN-only optimization)
data, offsets = pgen.read_ranges(
"chr2",
starts=[0, 1000, 2000],
ends=[1000, 2000, 3000],
mode=genoray.PGEN.GenosPhasingDosages,
)
# `data` matches the mode (tuple when mode bundles multiple arrays)
# `offsets` shape: (n_ranges + 1,). Slice range i with: arr[..., offsets[i]:offsets[i+1]]
# Chunked variants of both
for chunk in pgen.chunk("chr2", 0, 1000, max_mem="4g"): ...
for range_iter in pgen.chunk_ranges("chr2", starts, ends, max_mem="4g"):
for chunk in range_iter: ...
基因型 dtype:int32。剂量 dtype:float32。定相信息是一个形状为 (samples, variants) 的独立 bool 数组——不是基因型数组中的额外一行(与 phasing=True 时的 VCF 不同)。
SparseVar(.svar)——快速参考
构建:
# From a configured VCF reader
vcf = genoray.VCF("file.vcf.gz", dosage_field="DS")
genoray.SparseVar.from_vcf("out.svar", vcf, max_mem="4g",
with_dosages=True, overwrite=True)
# Or from a PGEN
genoray.SparseVar.from_pgen("out.svar", "file.pgen", max_mem="4g")
# Unphased somatic data: collapse to a single haploid call per sample (ploidy=1)
genoray.SparseVar.from_vcf("out.svar", vcf, max_mem="4g", haploid=True)
SparseVar.from_vcf / from_pgen 会继承并应用源文件的过滤条件——对 VCF/PGEN 进行过滤,即可过滤 SVAR。
SparseVar.from_vcf / from_pgen 接受 regions=、samples=、merge_overlapping=、regions_overlap=,可在转换期间按区域和/或样本筛选(语义与 SparseVar.write_view 相同);样本子集会从输出中移除 MAC=0 的变体。
读取:
# Plain ragged: data is just variant indices
svar = genoray.SparseVar("out.svar")
ragged = svar.read_ranges("chr1", starts=[0, 50_000], ends=[10_000, 60_000],
samples=["S1", "S2"])
# shape: (ranges, samples, ploidy, ~variants) — last axis is ragged
# With extra fields attached
svar = genoray.SparseVar("out.svar", fields={"dosages": np.float32})
# or, on an existing instance:
svar_with = svar.with_fields({"dosages": np.float32})
result = svar_with.read_ranges("chr1", [0], [10_000])
result.genos # Ragged of variant indices (uint32)
result.dosages # Ragged of dosages (float32)
with_fields(False) 会移除所有额外字段,之后的读取会再次返回普通的 Ragged[V_IDX_TYPE]。
不规则数组结果中的每个叶值都是一个变体索引——即 svar.index 中的行号。svar.index 是一个 polars DataFrame,至少包含 CHROM, POS, REF, ALT (list[str]), ILEN。要将索引映射回 chrom/pos/ref/alt,请按行索引该 DataFrame。
v_idxs = ragged[0, 0, 0].to_numpy()
rows = svar.index[v_idxs.tolist()].select("CHROM", "POS", "REF", "ALT")
svar.index.POS 是 1-based(VCF 约定),而查询坐标采用 0-based 半开区间。不要混淆两者。
SparseVar2(.svar2)——快速参考
SparseVar2 是新一代稀疏变体存储(采用 VariantKey 风格的内联编码和逐变体稠密/稀疏成本模型)。它分为两部分:转换(下文的 from_vcf)用于写入存储;范围查询(下文的 decode / region_counts / read_ranges)用于读取数据。与 genoray 中其他地方一样,所有坐标均采用 0-based 半开区间 [start, end)。
调优(genoray.Tuning)
每个 SparseVar2.from_* 写入方法(from_vcf、from_pgen、
from_vcf_list、from_svar1)都接受 tuning: Tuning | None = None。
Tuning 是显式的调度参数对象,用于替代旧的环境变量配置——没有任何环境变量可以配置
genoray;现在每个参数都是构造函数参数(或 CLI 标志,见下文“CLI”)。
from genoray import SparseVar2, Tuning
SparseVar2.from_vcf(
"out.svar2", "file.vcf.gz", "ref.fa",
tuning=Tuning(reader_workers=4, dense_cap=64),
)
包含 6 个字段,均为 int | None,仅限关键字参数,属于冻结的数据类:
| 字段 | 含义 | 最小值 |
|---|---|---|
concurrent_chroms | 并发转换的 contig | 1 |
reader_workers | 每个并发 contig 使用的独立索引分片读取器数量 | 1 |
overshard | 每个读取器的工作单元数(仅在 contig 没有精确记录数时使用) | 1 |
dense_cap | 读取器与执行器之间密集块通道的深度 | 1 |
merge_threads | 用于每个 contig 的 var_key 合并尾段的收集线程数 | 1 |
sample_interval | 监控采样间隔(秒);0 表示禁用 | 0 |
None(每个字段的默认值)表示**“由规划器选择”,而不是“关闭”——与完全不传入
Tuning 时规划器会选择的值相同。你设置的字段会按指定值执行,或直接拒绝,绝不会静默缩减**:如果显式设置的
reader_workers 无法满足 max_mem,就会引发 InsufficientMemory,而不会降低该值以适应内存限制。
这正是该对象的意义所在——旧环境变量可能会被从未读取它们的后端静默忽略;Tuning 不会如此。
tuning= 必须是 Tuning 实例或 None——例如,传入 dict 会引发 TypeError 并指出参数名称,因为其中的字段名称可能会被误认为是有效的关键字参数。
并非每个后端都能使用所有参数。设置了后端无法使用的参数时,会引发 ValueError 并指出有问题的字段,而不会静默忽略:
| 字段 | from_vcf | from_pgen | from_vcf_list | from_svar1 |
|---|---|---|---|---|
concurrent_chroms | 是 | 是 | 否 | 是 |
reader_workers | 是 | 否 | 否 | 否 |
overshard | 是 | 否 | 否 | 否 |
dense_cap | 是 | 是 | 是 | 是 |
merge_threads | 是 | 是 | 是 | 是 |
sample_interval | 是 | 是 | 是 | 是 |
原因:
reader_workers/overshard仅适用于分片 VCF。from_pgen将reader_workers=1固定为常量,且绝不会在 contig 内分片——pgenlib在解码过程中会持有 GIL,因此在该处进行亚 contig 分片纯属额外开销。from_vcf_list(每个 contig 的每个输入文件已各占用一个文件描述符)和from_svar1(完全不使用 VCF/htslib 读取器)也都不会在 contig 内分片。from_vcf_list不提供concurrent_chroms,因为该处理流水线按设计会依次遍历各个 contig。
你设置的字段会在 pipeline config 日志行中标记为 explicit;保持为 None 的字段则会标记为
planner——因此,单看该日志行就能知道哪些参数来自你,哪些来自规划器。
在适用的所有场景中,Tuning 纯粹用于控制调度:无论取值如何,输出的字节内容都完全相同(由存储哈希预言机验证)。
转换
from genoray import SparseVar2
dropped = SparseVar2.from_vcf(
"out.svar2", "file.vcf.gz", "ref.fa", # reference: validates REF + left-aligns indels
overwrite=True,
)
# Pre-normalized input (e.g. `bcftools norm`'d): skip REF validation/left-align
dropped = SparseVar2.from_vcf("out.svar2", "file.vcf.gz", no_reference=True)
签名:from_vcf(out, source, reference=None, *, regions=None, samples=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, ploidy=2, chunk_size=25_000, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, info_fields=None, format_fields=None, check_ref="e", progress=False, tuning=None, log_level="info", log_filter=None, max_mem=None) -> int
-
source— 经过 bgzip 压缩的 VCF(.vcf.gz,或等效写法.vcf.bgz) 或 BCF(.bcf)。如果未找到.csi/.tbi,则自动创建索引(.csi)。对于 PLINK2 PGEN 源,请改用from_pgen(见下文)。 -
regions=/merge_overlapping=/regions_overlap=— 将转换范围限制为一个或多个带索引的 VCF 获取区间。区域字符串使用现有的 genoray 约定("chrom:start-end"为从 1 开始且包含两端,转换为从 0 开始且右端不包含;元组/BED/帧输入已是从 0 开始且右端不包含)。 区域重叠时会报错,除非merge_overlapping=True。regions_overlap选择三种模式之一,与 bcftools--regions-overlap一致:"pos"(默认;POS 位于[start,end))、"record"(POS 位于[start,end+1),因此会保留位于区域最后一个碱基的插入缺失变异),或"variant"(锚点裁剪后的变异范围与区域重叠)。在"variant"模式下,如果多等位记录中任意一个等位基因确实与区域重叠,就会保留整条记录;不会丢弃单个不重叠的等位基因。"variant"模式目前要求每条染色体最多指定一个区域;每条染色体有多个区域时会报错——请使用"pos"/"record",或分别转换。 -
samples=— 按名称选择并重新排序 VCF 样本:保留调用方指定的顺序,对重复项仅保留首次出现的项,遇到未知名称时抛出ValueError。无论样本在原始 VCF 标头中的位置如何,available_samples和每个解码后的列都会严格匹配调用方指定的顺序。 -
reference(用于验证 REF 并对插入缺失变异进行左对齐的 FASTA 路径)和no_reference=True(信任已标准化的输入,跳过验证和左对齐)必须且只能指定其中一个——两者都传或都不传都会抛出ValueError。 -
reference=FASTA 使用的染色体命名方式可以与变异源不同(例如,源中为chr1,FASTA 中为1;或任一方的线粒体染色体命名为M/MT/chrM/chrMT);genoray 会在验证 REF/执行左对齐前,根据 FASTA 自身的命名方式解析源中的染色体名称(包括chr前缀和线粒体别名)。无论 FASTA 使用何种命名,输出存储都会保留源中的染色体命名方式。 -
skip_out_of_scope=False— 设为True时,会丢弃超出范围的 ALT(符号型<DEL>/<INS>/… 和断点),而不是报错;严格模式下,默认会在遇到第一个此类 ALT 时报错。在这一层无法区分这两类——没有单独控制“仅符号型”或“仅断点”的选项。 -
返回被丢弃的超出范围 ALT 数量,类型为
int(除非skip_out_of_scope=True,否则始终为0)。 -
check_ref: Literal["e", "x"] = "e"— 当记录的 REF 与参考 FASTA 不一致时采用的策略(no_reference=True时忽略)。"e"(默认)会抛出异常并中止构建,与bcftools norm --check-ref e一致。"x"会丢弃有问题的记录(包括越过染色体末端的 REF),并继续处理,同时记录每条染色体的数量。 比对不区分大小写(软掩蔽的小写参考碱基也视为匹配)。任何其他值都会在转换开始前抛出ValueError。 -
此处不接受
dosages=关键字参数(与下文的from_pgen不同)——VCF 中类似剂量的数据应通过format_fields=传入(例如DSFORMAT 字段)。不支持haploid=的 OR 合并,也不支持基于max_mem的分块(请改用chunk_size)——这两项目前仍仅适用于SparseVar(SVAR 1.0)。 -
threads=None— 线程总预算(若为None,则自动检测)。该参数决定染色体并发数,并由规划器决定每条染色体的读取器数量。 -
tuning.reader_workers(通过tuning=Tuning(reader_workers=N)设置,见上文“Tuning”)— 每条并发处理的染色体使用的独立索引分片读取器数量,用于设置子染色体读取并行度的参数。若保留为None,则根据核心预算推导:预留四分之一可用核心用于合并收尾,优先选择较高的并发深度(每条染色体约 8 个读取器),其余核心用于读取器。显式值必须为None或大于等于 1 的整数;小于 1 的值会在Tuning构造时、转换开始前抛出ValueError;若某个值符合max_mem限制但因其他原因无法满足,则会抛出InsufficientMemory,而不是静默降低该值。任何有效值生成的输出均逐字节相同。(from_vcf_list,即合并 N 个单样本 VCF 的路径,不会对单条染色体进行分片,并且不接受reader_workers/overshard——见上文的适用性表。)有关扩展数据,请参阅docs/source/svar.md中的“Parallel conversion”。 -
signatures=False— 设为True时,会在写入过程中将每个 SNP/插入缺失变异分类为 SBS96/ID83 突变类型代码,并为每条染色体存储一个mutcat侧文件(计入写入过程的密集数据/var_key 成本模型)。需要参考序列(reference=);与no_reference=True同时使用时会抛出ValueError。目前尚无用于读取 SVAR2mutcat侧文件的公开 API(与下文的SparseVar.annotate_mutations/mutation_matrix不同)——此标志仅控制是否写入该侧文件。 -
info_fields=/format_fields=— 默认值为None,类型分别为Sequence[str | InfoField]/Sequence[str | FormatField]。在写入过程中,将 INFO/FORMAT 中的标量数值字段提取到存储中:from genoray import SparseVar2, InfoField, FormatField SparseVar2.from_vcf( "out.svar2", "file.vcf.gz", "ref.fa", info_fields=["AC", InfoField("AF", dtype="f16")], format_fields=[FormatField("DS", default=0.0)], )- 范围:仅限标量数值类型。 标头中的
Type必须为Integer、Float或Flag;Number必须为1、双等位拆分的A,或0(仅限Flag的 INFO)。其他任何值(Number=R/G/.、String/Character字段)都会在配置时、转换开始前引发ValueError。裸str名称使用推断的 默认值(dtype=None,无default);传入InfoField/FormatField可覆盖。 dtype(FieldDtype = Literal["bool","i8","u8","i16","u16","i32","u32","f16","f32"]):None(默认值)会自动解析——Integer/Flag会无损自动缩窄为可容纳观测到的全局范围 的最小位宽(另加一个保留的缺失值哨兵);Float始终解析为f32(不会静默降精度)。 显式指定的dtype会在转换时根据标头类型(例如,Float不能使用整数位宽) 和观测范围进行校验——溢出或超出f16约 65504 的范围会引发ValueError。f16是唯一有损选项,必须显式请求。default——写入 VCF 缺失条目的值;否则使用所选位宽极值处的保留哨兵 (有符号位宽使用INT*_MIN,无符号位宽使用u*::MAX——自动缩窄会在观测范围非负时 优先选择无符号类型;浮点位宽使用NaN)。Flag字段永不缺失(不存在 ⇒false/0)。- FORMAT 与基因型对齐,而非独立无损:FORMAT 值仅在基因型有调用的位置存储——对于
按 var_key 路由的变体,每个携带者调用对应一个值;对于密集路由的变体,则存储每个样本的
完整密集列(非携带者位置填入
default/哨兵)。此版本会有意丢弃非携带者的 FORMAT 值(例如,ref/ref 基因型处插补得到的剂量);独立无损的 FORMAT 数据流将推迟到未来规范。 - **读取路径:**请参阅下方“读取 INFO/FORMAT 字段(SVAR2)”——通过
SparseVar2(path, fields=…)/.with_fields(…)/.available_fields可选择解码并读回这些字段。
- 范围:仅限标量数值类型。 标头中的
-
max_mem: int | str | None = None——并发规划器的字节预算:一次转换多少个 contig, 由此决定,使队列基线内存加上每个并发 contig 的在途分块缓冲区不超出预算 (此外还受现有核心数上限约束)。字符串格式与上文模块级max_mem约定相同 ("4g"、"512m"、"2GB",由parse_memory解析),且与from_vcf_list的max_mem具有相同的整个进程含义(见下文)——from_vcf_list没有可供其使用的已拟合并发规划器(其 contig 严格按顺序运行),因此会改为根据此预算推导自身的 每分块chunk_size。None(默认值)表示使用检测到的预算——cgroup 内存限制的 80% (或 cgroup 外的/proc/meminfo总量)——而非不设上限。这是相较于max_mem之前的规划器 有意做出的默认行为变更。如果检测本身失败(无 cgroup 限制且/proc/meminfo不可读取—— 在 macOS 上始终如此),genoray 会发出警告,并回退到旧的仅受核心数限制的规划方式,而非引发错误。 传入显式值可提高或降低预算;也可以传入非常大的值,以近似不设上限的规划。实际下限: 规划器的 RAM 公式包含固定的队列基线项和每个并发 contig 的项,因此任何无法覆盖基线加一个并发 contig 的预算都会引发ValueError,即使队列规模很小也不例外。该下限因后端而异:VCF 公式的原始 LP 系数约为 457 MB 基线,加上每个并发 contig 约 111 MB;但每个 contig 项中的kappa项完全压过 这两个数值,因此,即使队列规模很小,实际下限——按规划器在预算紧张时实际选取的cc=1, w=1点计算——也约为 1.0 GB(S=4,000 时为 1,016 MB;见下表),而非约 600 MB——实践中远低于此值的预算都会被拒绝。from_pgen的下限约为 2.7 GB,另加每个并发 contig 约 210 MB,因此其下限接近 3 GB。 2026-08-11 的内存预算模型调整使from_vcf的实际内存下限大约变为调整前的 4 倍,在固定的示例参数点chunk_size=25_000, reader_workers=3, cc=1下评估:单个并发 contig 所需的最小max_mem,在 S=4,000 时从约 1.15 GB 升至约 1.38 GB,在 S=128,000 时从约 7.8 GB 升至约 26.4 GB,在 S=500,000 时从约 27.9 GB 升至约 101.5 GB。随后,任务 3 将每个w的开销从每 chunk-MB 的kappa * (2w - 1)改为w * (kappa + 2) + 8,而且from_vcf不再将reader_workers固定为默认值——None会推导该值,并且plan_sharded会从w_max开始向下扫描w,直至1,之后才放弃某个 contig,因此预算紧张时现在会选择较小的w,而不是拒绝执行。当前默认值下,已发布的规划器实际强制要求的下限是cc=1, w=1,远低于下面任一w=3数值。“旧模型”列是上文引用的 2026-08-11 调整后数值;“当前模型”列则是如今显式指定reader_workers=3时,按w * (kappa + 2) + 8开销计算出的实际成本——两者数值不同,不要将它们视为同一个数值:群组 旧模型, w=3(2026-08-11 调整后)当前模型, w=3当前实际下限( w=1)S=4,000 1,380 MB 1,421 MB 1,016 MB S=128,000 26,400 MB 27,833 MB 14,864 MB S=500,000 101,480 MB 107,069 MB 56,408 MB 显式指定
reader_workers会使内存下限高于这个w=1最低值(根据上面的w * (kappa + 2) + 8开销,较大的w成本更高),因为显式指定的值会被遵守,或被拒绝,而不会悄悄降级为适合当前预算的w——传入reader_workers=3就要求满足上面的当前模型w=3数值(例如 S=500,000 时为 107,069 MB),否则会引发PlanError::InsufficientMemory,绝不会悄悄降级为w=1。整体方向仍然安全(要求更高只会导致分配过多内存或直接拒绝规划,不会导致 OOM)。在 S=500,000 时,64 GB 主机(max_mem默认为检测到的 RAM 的 80%,即52,429 MB)仍低于w=1的下限56,408 MB,因此仍会引发PlanError::InsufficientMemory——错误信息会指出两种解决方法:“提高max_mem或降低chunk_size”,不过现在差距已很小(52.4 GB 对比 56.4 GB),而不像旧模型那样相差悬殊。128 GB 主机(104,858 MB)则不再是差一点就够:它会按cc=1, w=2进行规划,所需内存为81,739 MB——约有 23 GB 的余量。(它会停在w=2,因为w=3需要107,069 MB,略超预算。) -
progress=False/log_level="info"/log_filter=None— 写入时的进度和日志记录,由from_vcf/from_pgen/from_vcf_list/from_svar1/write_view共用。progress=True会显示实时进度:在终端或 Jupyter 中显示rich进度条(每个正在处理的 contig 占一行);其他环境中则输出精简的心跳行,每个 contig 大约每 5 秒限流一次("chr1 42% (12,345/29,000) ...")。无论progress的值是什么,每个 contig 完成后都会输出一行"[svar2] chrom done: N kept, M excluded (Ts)"摘要,除非log_level="off"。log_level是结构化写入时日志行的最低严重级别。 可接受的值(不区分大小写):"off"、"critical"、"error"、"warning"、"info"(默认值)、"debug"、logging模块的整数常量 (logging.DEBUG及其他类似常量——介于两者之间的整数会向上取整到 严重级别更高的命名级别,与logging自身的过滤规则一致),或该整数的 十进制写法("10"),这是 CLI 使用整数路径的方式。规范级别:off=0, error=1, warning=2, info=3, debug=4。接受"critical",但它 是"error"的别名——不存在单独的 CRITICAL 级别。"warn"不受支持(ValueError),这是有意为之:logging.warn()已在 Python 3.13 中移除,因此 genoray 也不接受这种写法——请使用"warning"。"off"会禁用所有内容,包括每个 contig 的摘要和进度渲染(纯空操作,零开销);"info"还会显示线程预算选择、每个 contig 的开始/结束信息,以及当 contig 名称 与来源文件中的写法不同时,根据参考文件解析 contig 名称的信息;"debug"还会 显示每条记录的详细信息(因 REF/FASTA 不匹配而排除的记录、每个经过 左对齐处理的 indel)。具体映射请参阅parse_log_level。log_filter是可选的tracing风格EnvFilter指令字符串,用于筛选 genoray 自身输出到 stderr 的诊断信息,与log_level相互独立(后者仅控制上述结构化 进度/摘要事件)。None(默认值)会让 stderr 保持静默;无法解析的指令会让 stderr 静默,而不会引发异常。可用它指定log_level无法表达的模块/目标级 tracing,例如log_filter="genoray::monitor=trace"。没有环境变量可以配置 genoray。
log_level/log_filter(以及上文中的每个Tuning字段)取代了此前用于配置此流水线的全部 8 个环境变量—— 完整映射请参阅docs/source/svar.md中的“从环境变量迁移”,其中包括如何将日志 tracing 变量的genoray::monitor=trace指令迁移为log_filter="genoray::monitor=trace"。结构化日志行会在消息后以
key=value的形式内联显示字段 (例如:pipeline config concurrent_chroms=8 reader_workers=4 exact_counts=true planned_units=32)。
从 PGEN 转换
from genoray import SparseVar2
dropped = SparseVar2.from_pgen(
"out.svar2", "file.pgen", "ref.fa", # reference: validates REF + left-aligns indels
overwrite=True,
)
签名:from_pgen(out, source, reference=None, *, regions=None, samples=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, chunk_size=None, max_mem=None, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, dosages=None, check_ref="e", progress=False, tuning=None, log_level="info", log_filter=None) -> int
source— 一个.pgen文件。变异元数据从同级的.pvar/.pvar.zst读取,样本名称从同级的.psam读取。reference/no_reference、skip_out_of_scope、overwrite、long_allele_capacity、signatures和check_ref的含义均与from_vcf(见上文)相同,并返回相同的int(被丢弃的范围外 ALT)。- 与
from_vcf不同,PGEN 子 contig 分片已禁用(每个 contig 仅有一个读取器),且threads永远不会改变任何输出字节。原因 (实测数据:chr21c,约 1M 个变异 x 3202 个样本):单读取器转换 已经很快(约 33s),瓶颈在共享执行器/写入器和参考数据 I/O, 而非 pgenlib 解码;因此分片无法突破这一耗时下限,且实测 更慢(threads=24时为 44.9s,串行时为 32.6s)。将pgenlib升级为可释放 GIL 的构建版本(>=0.94.x,该版本通过prange并行化解码) 也没有帮助:在OMP_NUM_THREADS1..32 范围内,转换耗时始终约为 33s,因此解码并行化毫无收益。分片机制仍然存在, 且输出字节完全一致(已通过 1M 个变异验证);只有未来变更使瓶颈转移到解码时才会重新启用。 - 仅支持二倍体 — 没有
ploidy=关键字参数(from_vcf的默认值ploidy=2在此处 隐式固定)。 chunk_size=None— 与from_vcf固定的默认值25_000不同,此处的None会根据样本数推导变异数量预算(一个打包的密集分块占用chunk_size * n_samples * 2 / 8字节),避免适用于 200 个样本的固定常量在 500k 个样本时耗尽内存。传入显式int可覆盖此值。 如果推导值低于 256 个变异,则会发出警告 — 详情见下方from_vcf_list的chunk_size条目。此处的dosages计入n_format_fields。max_mem: int | str | None = None— 并发规划器的字节预算:同时转换多少个 contig。规划器会选择合适数量,使队列基线内存加上每个并发 contig 的在途分块缓冲区都不超出预算 (除现有的核心数限制外,并发 contig 数量还会限制为最多 8 个)。字符串形式与上文模块级max_mem约定相同("4g"、"512m"、"2GB",由parse_memory解析),且含义与from_vcf的max_mem(见上文)相同,均指整个进程的预算 — 两条流水线都有拟合好的并发规划器,且以相同方式将预算用于并发,只是采用各自拟合的 RAM 规律系数(PGEN 分块会同时解码两个单倍型,因此每个变异的成本更高)。from_vcf_list的max_mem也表示整个进程的预算,但该路径没有并发规划器可供分配预算(其 contig 严格按顺序运行),因此会根据预算推导每个分块的chunk_size— 见下方对应条目。None(默认值)表示使用检测到的预算 — cgroup 内存限制的 80%(或 cgroup 外/proc/meminfo中的总内存)— 并非不设上限。 这是相对于max_mem规划器引入前的有意默认行为变更。如果检测本身失败(没有 cgroup 限制,且无法读取/proc/meminfo— 在 macOS 上始终如此),genoray 会发出警告,并退回到仅受核心数限制的旧规划方式,而不会报错。传入显式值可提高或降低预算,也可传入很大的值以近似不设上限的规划。实际下限: 规划器的 RAM 规律包含一个固定的队列基线项,约为 2.7 GB (PGEN 自身拟合的系数,高于from_vcf的约 457 MB),因此,任何无法覆盖基线加上一个并发 contig 的分块缓冲区的预算都会被拒绝并抛出ValueError,即使是很小的队列也不例外。该基线会随队列大小增长(~0.0158 MB/sample),因此这并非仅与小队列有关:在约 500k 个样本时,基线本身就预测约 10.6 GB,因此较小主机上的检测到的预算会导致转换被拒绝 — 此时请传入符合主机内存的显式max_mem。每个 contig 的分块开销为2 个分块缓冲区,而非 VCF 路径计入的 10 个:8 个分块的重排积压上限是from_vcf的机制,PGEN 流水线(每个 contig 一个读取器)从不执行此限制,因此过去会抛出ValueError的两个预算界限之间,现在可以用于规划。regions=/merge_overlapping=/regions_overlap=— 约定、语义和三种重叠模式("pos"/"record"/"variant")均与from_vcf相同,可将转换限制在一个或多个.pvar变异索引范围内。与from_vcf一样,"variant"模式下,只要多等位记录的任意一个等位基因与区域确实重叠,就会保留整条记录。samples=— 按名称选择并重新排序.psam样本(约定与from_vcf相同):保留调用方指定的顺序,去重时保留首次出现项,遇到未知名称则抛出ValueError。无论各样本在原始.psam中的位置如何,available_samples和每个解码后的列都严格按照调用方指定的顺序排列。- 不支持
info_fields=/format_fields=— PGEN 不包含 FORMAT,且尚未实现.pvarINFO 提取。 dosages=Sequence[DosageField]— 将每个样本的剂量轨道存储为 FORMAT 字段。每个DosageField(name="dosage", source="self"|Path, dtype="f16"|"f32"="f32", default=None):source="self"从硬调用.pgen(上文的source)本身读取剂量;Path则从单独的.pgen读取(例如,由于 pgenlib 在硬调用与剂量同存于一个文件时会根据剂量推导硬调用,因此将 VAF/CCF 文件单独保存)— 该文件必须与硬调用.psam的样本相同,并与硬调用.pvar的变异逐一对应。 存储方式与任何 FORMAT 字段一样,按基因型对齐:使用 var_key 路由时,非携带者的剂量会被丢弃(对 VAF/CCF 类字段无影响;非携带者的值约为 0)。读取方式与其他 FORMAT 字段相同 — 见下文“读取 INFO/FORMAT 字段(SVAR2)”。from genoray import SparseVar2, DosageField SparseVar2.from_pgen("out.svar2", "cohort.pgen", "ref.fa", dosages=[DosageField(name="DS", source="self")]) # separate dosage file (e.g. VAF stored as dosage): SparseVar2.from_pgen("out.svar2", "hardcalls.pgen", "ref.fa", dosages=[DosageField(name="VAF", source="vaf.pgen")])- 未定相杂合子会按
pgenlib返回的等位基因编码顺序解析单倍型——这与from_vcf对未定相GT的处理存在同样的注意事项。 progress=False/log_level="info"/log_filter=None——与上文的from_vcf相同。tuning=接受concurrent_chroms、dense_cap、merge_threads、sample_interval;设置reader_workers/overshard会引发异常(from_pgen固定为每条 contig 使用单个读取器——参见上文“Tuning”)。
从单样本 VCF 列表转换
from genoray import SparseVar2
# Explicit list
dropped = SparseVar2.from_vcf_list("out.svar2", ["s1.vcf.gz", "s2.bcf"], "ref.fa")
# A directory of single-sample files
# (non-recursive: all *.vcf.gz/*.vcf.bgz, then all *.bcf)
dropped = SparseVar2.from_vcf_list("out.svar2", "vcfs/", "ref.fa")
# A manifest file (one path per line; blank/`#`-comment lines skipped;
# relative entries resolved against the manifest's directory)
dropped = SparseVar2.from_vcf_list("out.svar2", "manifest.txt", "ref.fa")
签名:from_vcf_list(out, sources, reference=None, *, regions=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, ploidy=2, chunk_size=None, max_mem=None, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, info_fields=None, format_fields=None, check_ref="e", progress=False, tuning=None, log_level="info", log_filter=None) -> int
通过原生 k 路归并,从 N 个样本各异且位点列表不同的单样本 VCF/BCF 构建 一个 SVAR2 存储库——无需使用 bcftools merge,也无需生成中间的多样本 VCF。
-
regions=/merge_overlapping=/regions_overlap=——与from_vcf使用相同的约定、语义和三种重叠模式("pos"/"record"/"variant"),并以相同方式应用于归并中的每个输入文件。与from_vcf一样,"variant"模式下,只要多等位基因记录中的任一等位基因与区域真正重叠,就会保留整条记录。 -
没有
samples=参数——与from_vcf/from_pgen/from_svar1不同,from_vcf_list无法按名称筛选队列样本:每个输入文件本身就是单样本,队列正是通过sources传入的文件集合。 -
每个输入文件必须是单样本文件——恰好包含一列样本;如果任一文件的样本列为零列或多于一列,则引发
ValueError。该样本的 VCF 头部名称将成为存储中的样本名称;不同输入文件间的样本名称重复时,引发ValueError。 -
sources——由模块级函数_resolve_vcf_sources解析,支持以下三种形式:Sequence[str | Path]——显式指定的文件,按给定顺序排列。- 单个目录
Path——目录内所有经过 bgzip 压缩的 VCF(*.vcf.gz/*.vcf.bgz),然后是所有*.bcf文件(不递归查找);每组均按natsort排序。 - 单个文件
Path——.vcf.gz/.vcf.bgz/.bcf文件按单个文件处理;其他文件均视为清单(每行一个路径,跳过空行和以#开头的注释行,相对路径以清单所在目录为基准解析)。 - 若解析结果中没有文件,则引发
ValueError。
-
缺失的位点 → 纯合参考基因型
0。 如果位点在文件 A 中有调用,但文件 B 中完全不存在,则 B 的样本在该位点填入0(纯合参考基因型)。 -
合并后无法观察到文件内的
./.。 SVAR2 的稀疏布局只存储携带 ALT 的条目,因此缺失单倍型和纯合参考单倍型都会产生零条目,无法通过decode或region_counts区分。-1缺失值哨兵是稠密genoray.VCF/genoray.PGEN的约定,不属于 SVAR2 的解码。(这种区别在合并过程中确实存在,但将基因型打包到稀疏载体位网格时会被丢弃——这与from_vcf一致,确保两条路径保持一致。) -
合并按原子变异执行连接——只有标准化后的
(pos, ref, alt)原子完全匹配时,变异才会作为所有文件共享的一行;仅位置相同并不足够。 -
每个输入文件的记录都必须已按每条 contig 的位置排序(与
from_vcf对其单个输入文件的假设相同)——未排序的文件会引发ValueError,并指出出错文件及相关位置,而不会静默破坏 k 路归并。 -
所有输入文件必须采用相同的 contig 命名方案(全部使用
chr1风格或全部使用1风格,不能混用)——合并会按每个文件中的确切字符串匹配 contig,因此如果一个队列混用不同方案,会立即引发ValueError并指出冲突文件和名称,而不会静默生成一个存储,使队列中一半样本在“名称不匹配”的 contig 上解码为全零。 -
会同时打开全部 N 个输入文件(每个 contig 的每个文件占用一个文件描述符)——当 N 较大时(大致为
N > (soft RLIMIT_NOFILE - 64) / 2;在默认软限制 1024 下通常约为 N ≈ 480),会引发ValueError并提示使用ulimit -n解决,而不是遇到 htslib 更难理解的“is there a .tbi or .csi file?”错误,该错误会指向接近限制的某个随机文件。没有可供超大型队列回退使用的分批/分层合并方式(后续工作)——请改为提高打开文件限制。 -
支持
no_reference=True,与from_vcf/from_pgen一致:跳过 REF 校验和左对齐,并从记录自身的 REF 字节重建每个原子变异。reference/no_reference必须且只能指定一个,以及signatures与no_reference不兼容的检查,也与from_vcf完全相同。- 此入口特有的注意事项:由于合并是按每条 contig 执行的 k 路连接,并以每个原子变异标准化后的
(pos, ref, alt)为键,跳过左对齐意味着只有当所有输入对同一位点的表示完全一致时,该位点才会合并为一个输出行——锚定碱基和填充方式必须相同(例如,所有文件都来自同一个调用器,或都已使用同一参考基因组运行bcftools norm)。如果两个文件用不同的标准化方式编码同一个 indel,在no_reference下它们不会合并:输出存储中会静默保留为两个单独的变异,而不是一个共享行。这不会引发错误——使用from_vcf_list并依赖no_reference前,请先确认上游标准化方式一致。
- 此入口特有的注意事项:由于合并是按每条 contig 执行的 k 路连接,并以每个原子变异标准化后的
-
info_fields=/format_fields=——声明 API 与from_vcf相同(根据sources中第一个文件的头部解析)。由于每个位点现在有 N 个来源列,合并语义与单文件转换有所不同:- INFO 字段采用首个携带文件优先。 位点由多个文件共享时,存储的 INFO 值取自按
sources顺序最靠前且包含该原子变异的文件——不是最后一个文件,也不是对各携带文件的值进行聚合(例如取最大值或求和)。 - FORMAT 字段按样本保存,与
from_vcf完全相同:每个样本使用其所在文件的值;完全不携带该原子变异的样本则使用该字段的默认值(保留哨兵值/NaN,或显式指定的default=)。
- INFO 字段采用首个携带文件优先。 位点由多个文件共享时,存储的 INFO 值取自按
-
chunk_size=None——不同于from_vcf固定的25_000默认值,这里的None会根据队列大小,通过_auto_chunk_size推导基于预算的块大小,使每个打包后的稠密块保持在由max_mem推导的单块目标上限内(见下文max_mem),最高约为 256 MiB;这与from_pgen/from_svar1已采用的默认值相同。预算同时考虑打包后的基因型网格和任何暂存的format_fields(每个变异n_format_fields * n_samples * 4字节——这一项可能是网格的32 * F / ploidy倍,例如在 F=7、ploidy=2 时为 112 倍),因此在大型队列中请求 FORMAT 字段会相应缩小自动块大小。适用范围:它只限制稠密块占用的内存;在典型队列大小下,这只是峰值 RAM 的一小部分——这是大型队列的防护措施,并不能解决输入文件数量导致的整体 RAM 扩展问题。传入整数可覆盖自动值并固定数量。如果推导出的chunk_size小于 256 个变异,_auto_chunk_size会发出警告(指出预算、推导出的chunk_size以及n_samples/ploidy/n_format_fields),并保留这个较小值,而不是引发错误——例如,在ploidy=2, n_format_fields=7时,若默认预算约为 256 MiB,n_samples≈37,118以上就会触发警告。提高max_mem或减少请求的format_fields可消除警告。所有通过此辅助函数推导chunk_size的三个转换器(from_pgen、from_svar1、from_vcf_list)均采用此行为。 -
max_mem: int | str | None = None——整个进程可使用的字节预算(字符串格式与模块级max_mem约定相同,例如"4g";由parse_memory解析)——含义与from_vcf的max_mem相同(见上文),并非每个块的上限。此预算只用于确定稠密块的大小——不会限制读取器每个输入文件的额外开销(from_vcf_list会对每条 contig 一次性打开全部N个文件);这部分开销随输入文件数量增长,完全不受max_mem限制。不要将“整个进程”理解为对总 RSS 的硬性保证。 破坏性变更: 早期版本将此max_mem视为单个在途稠密分块的字节数上限。现在它是整个进程的预算。以前,max_mem="512MiB"表示“允许单个稠密分块最多使用 512 MiB”;现在则表示“稠密分块这一项最多使用 512 MiB”,由此计算出的chunk_size会小得多。请重新调整任何沿用的值。from_vcf_list没有像分片式from_vcf读取器那样经过拟合的并发规划器——它严格按顺序转换 contig——因此,它不会规划并发,而是计算每个分块的字节目标:chunk_target = min(_DENSE_CHUNK_TARGET_BYTES, max_mem // (concurrent_jobs * in_flight_chunks_per_job))。concurrent_jobs和in_flight_chunks_per_job并非硬编码的 Python 字面量——它们会实时从 Rust 扩展模块读取(_core.VCF_LIST_CONCURRENT_CHROMS,当前为1;以及_core.VCF_LIST_DENSE_CHANNEL_CAP + 2,当前为8:读取器/执行器通道的容量,加上单个读取器线程和单个执行器线程各自在通道之外可能持有的一个分块),因此不会像重复定义的常量那样,与编排器的实际架构悄然脱节——无论如何,这都不是拟合出的内存规律(这条流水线不存在这样的规律,不同于分片路径的 RAM 规律)。_DENSE_CHUNK_TARGET_BYTES仍是一个上限,因此较大的预算不会让分块超过当前大小——只有预算紧张时,分块才会缩小。除此之外,计算出的目标仍是一个最坏情况上限:估算时假设每个变异都走稠密路径,因此变异走稀疏路径的队列(例如仅私有的体细胞调用)占用会少得多。显式传入chunk_size时,此设置将被忽略——包括完全跳过内存预算检测(不尝试,也不警告)。None(默认值)表示使用检测到的预算,检测方式与from_vcf相同;如果检测失败,则会退回历史上的固定约 256 MiB 稠密分块目标(并发出警告),而不是抛出错误。 -
ploidy、skip_out_of_scope、threads、overwrite、long_allele_capacity、signatures、check_ref的含义均与from_vcf相同,返回值也是相同的int(被丢弃的范围外 ALT 数量)。合并期间会对每个输入文件分别应用check_ref(no_reference=True时忽略):在"x"模式下,错误记录只会从其所属文件中排除(不会排除整个合并位点),每个 contig 的日志会报告所有输入文件中被排除记录的总数。 -
progress=False/log_level="info"/log_filter=None——与from_vcf相同(见上文)。tuning=只接受dense_cap、merge_threads、sample_interval——设置concurrent_chroms、reader_workers或overshard会引发错误(此流水线按顺序遍历 contig,且不会在 contig 内分片——见上文“调优”)。
从 SVAR1 转换
from genoray import SparseVar2
dropped = SparseVar2.from_svar1(
"out.svar2", "old.svar", "ref.fa", # reference: validates REF + left-aligns indels
overwrite=True,
)
签名:from_svar1(out, source, reference=None, *, regions=None, samples=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, chunk_size=None, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, fields=None, check_ref="e", progress=False, tuning=None, log_level="info", log_filter=None) -> int
将现有的 SVAR 1.0(SparseVar)存储原生迁移到 SVAR2——不读取
VCF,也不使用 htslib;SVAR1 本身已是稀疏格式,因此会从 SVAR1 的数组重建变异
记录,并复用与 from_vcf 相同的转换流程。
source—SparseVar存储目录(SVAR1)。reference/no_reference、skip_out_of_scope、overwrite、long_allele_capacity、signatures和check_ref的含义均与from_vcf(上文)相同,并返回相同的int(被丢弃的范围外 ALT 数量)。ploidy从 SVAR1 的元数据中读取——没有ploidy=关键字参数。chunk_size=None会根据队列规模推导变异数量预算,方式与from_pgen/from_vcf_list相同(_auto_chunk_size),并在低于 256 个变异的情况下发出相同 的警告——详情请参阅上文from_vcf_list中的chunk_size条目。预算会计入fields=携带的 FORMAT 字段(默认全部计入),因此宽字段的 SVAR1 存储会推导出更小的块, 与from_pgen处理dosages=时相同。- 仅支持双等位基因 SVAR1 — 如果源存储包含多等位基因变异(SVAR1 的
geno==1模型),则会引发ValueError;请先以双等位基因方式重新创建 SVAR1 存储。 regions=/merge_overlapping=/regions_overlap=— 约定、 语义以及三种重叠模式("pos"/"record"/"variant")均与from_vcf/from_pgen相同。"variant"模式下,只要记录中的任意等位基因与区域真正 重叠,就会保留整条记录(不过 SVAR1 本身仅支持双等位基因,因此这里只会判断单个 ALT)。 与from_pgen不同,SVAR1 没有磁盘上的覆盖范围索引可供预先缩小搜索范围——仍会完整扫描 所选染色体上的局部变异;真正限制输出的是逐条记录过滤,因此需要完整扫描整条染色体, 而不是仅扫描受范围限制的区域。samples=— 按名称选择并重新排序 SVAR1 样本(约定与from_vcf/from_pgen相同):保留调用方指定的顺序,对重复项仅保留首次出现的项,遇到未知名称时引发ValueError。available_samples和所有解码后的列都严格遵循调用方指定的顺序, 与各样本在原始 SVAR1 中的位置无关。fields=Sequence[str] | None— 选择要传递的 SVAR1 FORMAT 字段(例如dosages),按其 SVAR1 名称指定:None(默认值)传递全部字段(此前默认无损的行为),[]不传递任何字段,名称子集只传递对应字段——遇到未知名称时会引发ValueError,并列出可用字段。无法通过此方式选择mutcat,它始终会被丢弃——请传入signatures=True,以从参考中重新计算签名,而不是传递 SVAR1 中的签名。- **字段一致性注意事项:**由于 SVAR1 从不存储非携带者的 FORMAT 值,因此字段输出仅在 var_key
(仅携带者)路由的变异中与
from_vcf逐字节一致——对于密集路由的变异,非携带者单元格会填入 字段的默认/缺失哨兵值,而不是源 VCF 中的实际值。基因型流本身(而非字段)在归一化方式一致时, 无论采用何种路由方式,都与from_vcf逐字节一致。 - 没有
info_fields=/format_fields=关键字参数 — 这些名称仅适用于 VCF; 请改用fields=(上文)来选择要传递的 SVAR1 字段。 progress=False/log_level="info"/log_filter=None— 与from_vcf(上文)相同。tuning=接受concurrent_chroms、dense_cap、merge_threads、sample_interval;设置reader_workers/overshard会引发错误(from_svar1不读取 VCF/htslib,且从不会在染色体内分片 ——请参阅上文“Tuning”)。
区间查询
打开一个已完成的存储,然后按 contig 查询。构造时会读取 meta.json,并为每个 contig 打开一个原生读取器,提供 .available_samples(列表;与 VCF/PGEN/SparseVar 共享的规范样本名称访问器)、.n_samples、.contigs、.ploidy、.format_version。
from genoray import SparseVar2
sv = SparseVar2("out.svar2")
regions = [(0, 40), (1_000, 2_000)] # 0-based half-open [start, end)
# Analysis path — decode to a seqpro Ragged record (one call per contig)
rag = sv.decode("chr1", regions) # fields pos (i32), ilen (i32), allele (ALT bytes)
# + one per selected field (see "Reading
# INFO/FORMAT fields" below); shape
# (R, S, P, None); pure-DEL ALT is empty
# Decode-free per-(region, sample, ploid) variant count — replaces SVAR 1.0's var_ranges
counts = sv.region_counts("chr1", regions) # np.ndarray, shape (R, S, P)
decode(contig, regions)返回一个seqpro.rag.Ragged,其布局与 gvl 的RaggedVariants逐字节相同(pos/ilen为数值型,allele为不透明字符串 ALT,共享一个变异轴 offsets 对象)。纯缺失变异的 ALT 为空(不会重新输出参考碱基)。需要seqpro。region_counts(contig, regions)是无需解码的计数方式(offset 差值 + 稠密掩码 popcount)——可简化替代SparseVar.var_ranges(SVAR2 没有统一变异表,因此变异索引已不复存在)。- 查询按 contig 分别进行——跨 contig 批处理由调用方负责。regions 是由
(start, end)对组成的可迭代对象。 decode/region_counts/read_ranges的contig参数接受替代命名方案——带chr前缀与不带前缀(chr1↔1),以及线粒体别名{M, MT, chrM, chrMT}——并通过ContigNormalizer解析为存储自身使用的名称。无法解析的 contig 会引发ValueError。
面向用户的 SVAR2 查询 API 是 decode / region_counts / read_ranges(如上所述)。read_ranges(contig, starts, ends, samples=None) 是融合的搜索与收集操作;starts/ends 是平行的一维数组(与 SparseVar.read_ranges 相同),samples 按名称选择/重新排序子集。它返回原始的双通道 BatchResult → numpy 字典,这是一个 TypedDict,包含固定字段集:vk_pos/vk_key/vk_off、dense_pos/dense_key/dense_range/dense_present/dense_present_off、lut_bytes/lut_off,以及标量 n_regions/n_samples/ploidy。
SparseVar2 还提供 _overlap_batch/_find_ranges/_gather_ranges(以下划线开头)——这是仅供 gvl 内部使用的 numpy 字典线协议,用于写入时重叠缓存所采用的搜索/收集拆分。它们不属于公共 API,不受 semver 约束,并且可能随时更改或消失,恕不另行通知;请勿在用户代码中调用它们。
读取 INFO/FORMAT 字段(SVAR2)
使用 from_vcf(info_fields=…, format_fields=…) 写入的字段(见上文)需要选择启用后才能读回——默认不会解码(每个字段都会产生额外 I/O 开销)。
sv = SparseVar2("out.svar2")
sv.available_fields # {"AF": StoredField(...), "DS": StoredField(...)}
sv = sv.with_fields(["AF", "DS"]) # or SparseVar2("out.svar2", fields=["AF", "DS"])
rag = sv.decode("chr1", [(0, 10_000)])
rag["AF"] # Ragged, sharing offsets with pos/ilen/allele
SparseVar2(path, *, fields=None)/.with_fields(fields)—fields是由规范键组成的Sequence[str](见下文的available_fields)。with_fields会基于同一存储返回一个新的SparseVar2,不会就地修改原对象。fields=None(构造函数的默认值)表示不选择任何字段——字段必须选择启用。available_fields -> dict[str, StoredField]— 存储的meta.json中声明的所有字段,键为规范形式:如果某个字段名在 INFO 和 FORMAT 中唯一,则使用裸字段名;否则使用 bcftools 风格的INFO/DP/FORMAT/DP。StoredField(定义于genoray._svar2_fields,未从genoray顶层导出)是一个冻结的数据类,包含:name、category("info"/"format")、dtype(np.dtype)、default(float | None)、key。decode(contig, regions)会为每个选定字段附加一个Ragged,与pos/ilen/allele一同放入返回的记录Ragged中——它们都共享同一个变异轴偏移对象,形状为(R, S, P, None)。 使用rag["KEY"](Ragged.__getitem__)访问字段数据,不要使用rag.fields["KEY"]——Ragged.fields只是记录上字段名称组成的list[str]。- 保留存储时的 dtype。 SVAR2 会在写入时无损地自动缩窄整数型字段,因此例如
AC字段读回时可能是int8;读取时不会拓宽类型。 - 缺失值为写入时设置的字段
default;若未设置,则使用保留哨兵值(浮点数为NaN,整数为iinfo.min/iinfo.max)——按原样返回,不会转换。 - FORMAT 字段与基因型对齐(见上文
from_vcf的format_fields=)——decode()只会输出携带变异的记录,因此写入路径中“非携带者的值不会存储”这一注意事项在此读取接口中不会造成影响。
突变特征(SBS96 / DBS78 / ID83)
与 SparseVar 使用相同的 COSMIC 工作流(见上文“突变目录”),但由每条染色体的 Rust 侧车文件提供支持,而不是附加在 .gvi 上的字段。在调用 mutation_matrix 前必须进行注释——可以事后注释,也可以在调用 from_vcf 时传入 signatures=True(见上文):
sv = SparseVar2("out.svar2")
ref = genoray.Reference.from_path("hg38.fa")
sv.annotate_mutations(ref) # post-hoc; writes the mutcat sidecar
sv.annotate_mutations(ref, contigs=["chr1"]) # restrict to a subset of contigs
df = sv.mutation_matrix("SBS96") # count="allele" (default)
df = sv.mutation_matrix("DBS78", count="sample")
act = sv.assign_signatures("SBS96") # mutation_matrix + fit_signatures
annotate_mutations(reference, *, gtf=None, contigs=None) -> None—reference是genoray.Reference或 FASTA 路径;contigs=None(默认值)会注释每个 contig。contigs=接受替代命名方式 — 带chr前缀或不带前缀,以及线粒体别名{M, MT, chrM, chrMT}— 通过ContigNormalizer解析为存储区自身的 拼写;如果所有请求的 contig 都无法解析,则会引发ValueError。 与SparseVar.annotate_mutations不同,此处没有write_back=开关 — SVAR2 始终将 sidecar 持久化到磁盘。gtf=可选地提供 GTF/GFF 基因模型路径;提供后,每个 SNV 还会按转录链类别进行分类(基于feature == "gene"足迹),并持久化到strand.binsidecar,从而解锁下文的"SBS192"/"SBS384"目录。每个 sidecar 和meta.json都会以原子方式写入(同目录临时文件 + 重命名),因此崩溃 不会留下截断文件;重新运行时始终会无条件重写,因此也可用于恢复缺失或可疑的 sidecar。 不带gtf=重新运行时,会移除先前写入的strand.bin(不会将过时的链类别 保留在不带 gtf 的新注释中)。mutation_matrix(kind, *, count="allele"|"sample", contigs=None) -> pl.DataFrame— 包含MutationType列(固定的 COSMIC 编码顺序)以及每个样本各一列。kind ∈ {"SBS96", "DBS78", "ID83", "SBS192", "SBS384"}。count="allele"会计数每个非参考等位基因拷贝;count="sample"会对每个样本中的每个类别最多计数一次,并在 contig 间进行 OR 合并。contigs=将求和范围限制为一个子集(接受替代命名方式;重复项会合并;存储区中不存在的名称会引发ValueError),且只要求这些 contig 已注释;None(默认值)使用每个 contig。如果范围内的 contig 尚未注释(磁盘上没有 sidecar),则会引发ValueError— 请先通过annotate_mutations或from_vcf(..., signatures=True)进行注释。"SBS192"/"SBS384"还需要链注释(annotate_mutations(..., gtf=...));如果范围内的 contig 缺少链注释,则会引发ValueError。如果 sidecar 被截断或已过时(长度与存储区的变异计数不匹配),则会引发OSError,提示你重新运行annotate_mutations。assign_signatures不接受"SBS192"/"SBS384"— 见下文。assign_signatures(kind, *, reference=None, count="allele", contigs=None, strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame— 先调用mutation_matrix(kind, count=..., contigs=...),再调用genoray.fit_signatures(...)。reference接受pl.DataFrame、TSV 路径或None(默认使用genoray.cosmic_signatures(kind))。contigs=(contig 子集,接受替代命名方式)会传递给mutation_matrix。strategy=/max_delta=/min_activity=/criterion=会原样传递给fit_signatures,与顶层函数相同 — 传入strategy=Spa()可使用 SPA 算法;也可保留为None,以使用max_delta/min_activity/criterion前向选择简写形式。只会调用一次fit_signatures,因此同时传入strategy=和任意一个简写参数都会引发ValueError,与fit_signatures完全相同。- 分类规则与 v1 相同(共用 Rust 分类器):DBS78 仅由孤立的相邻同单倍型 SNV 对产生 — ≥3 个 相邻 SNV 的连续片段仍作为单独的 SBS96 条目,原生长度 > 2bp 的 MNV 会 原子化为 SBS96,且每个孤立双联体只计数一次(而非按其组成的每个 SNV 各计数一次)。
- SVAR2 不提供公开的读取接口来访问每个基因型的原始
mutcat代码(不同于 v1 的fields=["mutcat"])— 仅公开聚合后的mutation_matrix输出。
链特异性目录(SBS192 / SBS384)
SparseVar2 还支持转录链偏倚目录,这些目录在注释时需要基因模型(GTF):
sv2.annotate_mutations(reference, gtf="gencode.v45.annotation.gtf.gz")
sbs384 = sv2.mutation_matrix("SBS384") # 384 rows: [T, U, N, B] x 96
sbs192 = sv2.mutation_matrix("SBS192") # 192 rows: the {T, U} sub-view = SBS384[:192]
- SBS384 = 96 个三核苷酸通道 x 4 个链类别,采用 SigProfiler 顺序
[T, U, N, B]:Transcribed(转录链)、Untranscribed(非转录链)、Nontranscribed(非转录区,基因间区)、Bidirectional(双向,位置被两条链上的基因覆盖)。 - SBS192 是
{T, U}子视图(SBS384[:192])。 - 链判定规则(嘧啶折叠):基因区内的 SNV 仅当其 ref/alt 碱基对中的嘧啶位于基因的编码链上时,才归为非转录链;否则归为转录链。基因覆盖范围取自
feature == "gene"行(完整基因体);可预先筛选 GTF 以限定生物类型。 - 如果未指定
gtf=,mutation_matrix("SBS192"/"SBS384")会报错。写入时的from_vcf(..., signatures=True)不处理链信息;如需获取链特异性目录,可事后调用annotate_mutations(reference, gtf=...)。 assign_signatures("SBS192"/"SBS384")会引发NotImplementedError:COSMIC 未发布链特异性参考集。请使用mutation_matrix进行链偏倚分析。
突变簇(SigProfilerClusters)
annotate_clusters 将 SigProfilerClusters 的逐突变子类分类(doublet / MBS / omikli / kataegis / other)写入已完成的存储,作为 cluster_class FORMAT 字段。此操作仅用于事后处理(没有写入时标志),且无需参考基因组——分类器不会使用序列上下文。
sv = SparseVar2("out.svar2")
sv.annotate_clusters(imd_cutoff=1000.0) # no-VAF mode
sv.annotate_clusters(imd_cutoff=1000.0, vaf_field="VAF") # VAF mode
sv.annotate_clusters(imd_cutoff={"S1": 500.0, "S2": 1_000.0}, vaf_field="VAF")
sv.annotate_clusters(imd_cutoff=1000.0, vaf_field="VAF", contigs=["chr1"])
# The original instance is stale for the new field (`available_fields` is read
# at construction); re-open the store and opt in before decoding.
sv = SparseVar2("out.svar2").with_fields(["cluster_class"])
rag = sv.decode("chr1", [(0, 1_000_000)])
labels = rag["cluster_class"] # 1:1 with rag["pos"], the per-call flat order
- 签名:
annotate_clusters(*, imd_cutoff: float | Mapping[str, float], vaf_field: str | None = None, vaf_cut: float = 0.1, contigs: Sequence[str] | None = None) -> None。 imd_cutoff— 突变间距离阈值(单位:碱基):当突变与最近邻的距离<=其样本的阈值时,该突变会被归入突变簇。可以是应用于所有样本的单个浮点数,也可以是从样本名称映射到阈值的映射;映射必须且只能包含available_samples中的每个条目(否则会引发ValueError);阈值必须> 0。阈值由调用方提供,而非存储提供:上游通过模拟为每个样本推导阈值,并将其保存在imds.pickle中——输入相同的值,上游分类的每个事件都会得到相同的标签(SigProfilerClusters v1.2.2)。vaf_field— 可选的 FORMAT 字段(规范键或裸名称,参见available_fields),用于保存每个变异的 VAF/CCF。指定后,VAF 一致性会参与决策树判断,未通过的事件会被贪心地重新拆分。必须解析为数据类型为 2 或 4 字节浮点数(f16/f32)的 FORMAT 字段;否则会引发ValueError。缺失值按上游的-1.5哨兵值处理。不指定(默认值None)则使用无 VAF 路径。vaf_cut— 相邻突变仍归入同一事件所允许的最大|delta VAF|(上游默认值为0.1;0.25是其 CCF 设置)。必须> 0。contigs— 将注释限制到指定子集(名称别名按annotate_mutations中的规则解析:带chr前缀与不带前缀的名称,以及线粒体别名);存储中不存在的名称会引发ValueError。所有范围外的 contig 都会填充 255,以确保该字段在这些位置仍可解码。None(默认值)表示注释每个 contig。重新运行时始终会执行注释。- 在
meta.json中写入cluster_version(当前为1)、cluster_contigs、cluster_cutoff、cluster_vaf_field、cluster_vaf_cut,并添加cluster_classFORMAT 条目(u8)。
标签编码表——每个 SNP 携带者调用对应一个标签(1/2 基因型中的每个记录的 alt 都各有一个标签;indel 不会被分类):
| 编码 | 标签 | 上游 |
|---|---|---|
| 0 | nonclustered | 不在聚类文件中 |
| 1 | doublet | ClassIA |
| 2 | mbs | ClassIB |
| 3 | omikli | ClassIC |
| 4 | kataegis | ClassII |
| 5 | other | ClassIII |
| 255 | not_annotated | 范围外 contig / 非 SNP 调用 |
通过现有字段机制读回——
with_fields(["cluster_class"]),然后执行 decode(contig, ranges)。该字段与 pos/ilen/allele 共用同一变异轴,因此
rag["cluster_class"] 按存储中的逐调用扁平顺序排列,与 rag["pos"] 一一对应。在磁盘上,var_key_* 流为每次调用存储一个标签;dense_* 流按 (dense_row, sample) 存储一个标签,顺序为 dense_row * n_samples + sample,非携带者为 255;indel 流的值始终为 255。范围内的 SNP 调用都会获得标签:对于上游丢弃的事件(聚类单例、无 VAF 的失败事件),genoray 会写入 other;非聚类的 SNP 调用则标记为 nonclustered。
按 contig 合并和拆分
SVAR2 的 contig 在磁盘上完全独立,因此重新合并或选取整个 contig 只需进行低成本的元数据重写和文件复制操作。与 write_view(见下方 CLI 部分)不同,这些方法都不会重新运行转换或 var_key/dense 成本模型。
from genoray import SparseVar2
sv = SparseVar2("out.svar2")
sv.subset_contigs("chr1.svar2", "chr1") # single contig
sv.subset_contigs("subset.svar2", ["chr1", "chr2"]) # multiple, source order preserved
paths = sv.split_by_contig("by_contig/") # one store per contig, out_dir/{contig}.svar2
SparseVar2.concat("merged.svar2", ["chr1.svar2", "chr2.svar2"]) # disjoint-contig merge
subset_contigs(output, contigs, *, mode="copy", overwrite=False) -> None— 写入一个仅包含contigs的新存储(单个 contig 名称或名称序列)。contigs接受其他命名方式——带chr前缀或不带前缀,以及线粒体别名{M, MT, chrM, chrMT}——并通过ContigNormalizer转换为存储自身使用的名称。仅重写元数据并复制保留 contig 的目录,同时保留源存储的 contig 顺序。如果任何名称无法与self.contigs匹配,或output指向此存储自身的路径(拒绝原地选取,类似于write_view的原地操作防护),则引发ValueError。如果output已存在且overwrite=False,则引发FileExistsError。split_by_contig(out_dir, *, mode="copy", overwrite=False) -> list[Path]— 拆分为每个 contig 一个单 contig 存储,路径为out_dir/{contig}.svar2;按self.contigs的顺序返回输出路径。实现方式是针对每个 contig 调用一次subset_contigs。SparseVar2.concat(output, sources, *, mode="copy", overwrite=False) -> None(类方法)— 将 contig 集合互不相交的存储合并为一个。sources是路径序列(或SparseVar2实例序列);所有源必须在samples、ploidy、format_version和fields上保持一致。上述任一项不一致,或某个 contig 名称出现在多个源中,都会引发ValueError。合并后的 contig 列表会按natsorted排序,不受sources的传入顺序影响。mode(适用于全部三种方法)是共用的Mode字面量 —"copy"|"hardlink"|"symlink"|"move"— 用于控制将每个 contig 目录迁移到输出存储的方式。
错误
genoray 按类别引发标准 Python 内置异常:
ValueError— 输入内容有误:未找到 contig/sample、REF 与参考 FASTA 不一致,或skip_out_of_scope=False时出现符号型 ALT 或断点 ALT。FileNotFoundError— 缺少必需的输入文件。OSError— 存储的边车文件损坏或被截断,或底层磁盘 I/O 失败。RuntimeError— genoray 内部错误(工作线程发生 panic);请报告此问题。
CLI
genoray write 有三个子命令——write vcf、write pgen、write svar1——且这三个命令的目标格式都是 SVAR2。不再支持不指定类型、自动检测输入的 genoray write SOURCE OUT;必须指定源类型。此前写入 SVAR 1.0(SparseVar)的路径位于顶层 genoray write-svar1 命令(使用连字符,不是 write 子命令)中——和以前一样,它接受 VCF 或 PGEN 源。genoray view 仍默认使用 SVAR2;此前 SVAR 1.0 的行为可通过 view svar1 使用。genoray concat/genoray split 仅支持 SVAR2(没有对应的 SVAR1 版本)。
genoray write vcf / genoray write pgen / genoray write svar1
# write vcf — VCF/BCF (or a directory/manifest of single-sample VCFs/BCFs) → SVAR2
genoray write vcf file.vcf.gz out.svar2 --reference ref.fa
genoray write vcf file.vcf.gz out.svar2 --no-reference
genoray write vcf file.vcf.gz out.svar2 --reference ref.fa --skip-symbolics-and-breakends --threads 4
genoray write vcf file.vcf.gz out.svar2 --reference ref.fa --fields INFO/AF --fields FORMAT/DP
genoray write vcf vcf_dir/ out.svar2 --no-reference --regions chr1:1-1000 # vcf-list form
# write pgen — PLINK2 PGEN → SVAR2 (no --ploidy; PGEN is diploid-only)
genoray write pgen file.pgen out.svar2 --reference ref.fa
genoray write pgen file.pgen out.svar2 --no-reference --regions chr1:1-1000 --samples A,B
genoray write pgen file.pgen out.svar2 --reference ref.fa --dosages DS=self
genoray write pgen hardcalls.pgen out.svar2 --reference ref.fa --dosages VAF=vaf.pgen
# write svar1 — SVAR1 (SparseVar) → SVAR2
genoray write svar1 store.svar out.svar2 --no-reference --samples A,B
genoray write svar1 store.svar out.svar2 --no-reference --fields dosages
genoray write svar1 store.svar out.svar2 --no-reference --empty-fields
# write-svar1 (legacy, top-level) — VCF or PGEN → SVAR 1.0, dosages, --haploid, --max-mem
genoray write-svar1 file.vcf.gz out.svar --max-mem 4g --haploid
这三个 write 子命令都支持 --regions/-r、--regions-file/-R、
--samples/-s、--samples-file/-S、--merge-overlapping、
--regions-overlap(pos/record/variant)、--reference XOR
--no-reference(必填)、--chunk-size、--threads/-@、--overwrite、
--long-allele-capacity(高级选项)、单个 --skip-symbolics-and-breakends
标志(映射到 skip_out_of_scope=;SVAR2 核心无法将符号 ALT(<DEL>、<INS>、…)
或断点扩展为核苷酸,因此会一并丢弃,并在设置该标志时输出一行
Dropped {n} out-of-scope (symbolic/breakend) ALT alleles.)、--check-ref {e,x}(默认值为 e,与 --no-reference 一起使用时会被忽略;e 会在首次遇到 REF/FASTA 不一致时中止,x 会丢弃有问题的记录并继续——与 bcftools norm --check-ref 一致)、
--progress/--no-progress + --log-level LEVEL(映射到
progress=/log_level=;行为见上文“Conversion”——默认值为 --no-progress --log-level info。LEVEL 可以是 off、critical、error、warning、info、debug,不区分大小写——即 parse_log_level 接受的所有名称。warn 不会被接受:logging.warn() 已在 Python 3.13 中移除,因此它是已弃用的拼写,而非惯例。该参数经过验证,类型为 str,而不是固定选项列表,正是为了确保它与作为唯一事实来源的 parse_log_level 保持一致;该函数也接受的整数 logging 常量同样可以通过此标志传入,写法如下:--log-level 10 == --log-level debug),以及
--log-filter DIRECTIVE(映射到
log_filter=;采用 tracing 风格的 EnvFilter 字符串,例如 --log-filter "genoray::monitor=trace",且独立于 --log-level;默认未设置/None)。write vcf 的 vcf-list 形式会将这些参数传递给 from_vcf_list;单文件形式则传递给 from_vcf。
每个 write 子命令还会提供 Tuning 中始终适用的调节项,形式为
--dense-cap N、--merge-threads N 和 --sample-interval N(映射到
tuning=Tuning(dense_cap=, merge_threads=, sample_interval=),见上文“Tuning”)。--concurrent-chroms N 还可用于 write vcf、
write pgen 和 write svar1,但不适用于 vcf-list 形式(from_vcf_list
按设计会按顺序遍历 contig;见上文“Tuning”)。与 Python 的 tuning= 关键字参数一样,显式指定的标志值会被遵循或拒绝——绝不会静默降低;传入后端无法使用的标志会引发错误,而不会被忽略。
write vcf 还接受 --reader-workers N 和 --overshard N(仅适用于单文件输入;如果与目录/清单一起传入其中任何一个,都会引发错误——vcf-list 形式绝不会在 contig 内分片)。单文件形式会将两者直接传递给 tuning=Tuning(reader_workers=N, overshard=N),因此在此处传入 N < 1 会引发相同的 ValueError——CLI 不会单独检查下限。
genoray write vcf(SparseVar2.from_vcf/from_vcf_list):source是单个.vcf.gz/.vcf.bgz/.bcf文件时使用from_vcf;其他情况(目录,或非.vcf.gz/.vcf.bgz/.bcf文件)使用 vcf-list 形式(包含单样本 VCF/BCF 文件的目录,或列出这些文件的清单),并调用from_vcf_list——.svar(SVAR1)源应使用write svar1,而不是此命令。--samples/--samples-file仅适用于单文件形式——用于 vcf-list 形式时会报错(每个输入文件本身恰好包含一个样本,因此不存在可供筛选的队列)。--fields(-f,可重复)接受 bcftools 风格的INFO/x/FORMAT/x/FMT/x规范,由_parse_cli_field_specs解析并作为info_fields=/format_fields=传递;默认未设置(不携带任何字段,仅保留基因型)。--chunk-size默认值为25000。此处接受--ploidy(默认值为 2)。genoray write pgen(SparseVar2.from_pgen):source是.pgen文件。不支持--ploidy(PGEN 仅支持二倍体)。--dosages(可重复)接受NAME=self(从source本身读取剂量)或NAME=/path/to/vaf.pgen(从单独的 PGEN 读取),每项都会转换为DosageField(name=NAME, source=...),并作为dosages=传递。--chunk-size默认使用根据内存推导的值(None)。--max-mem(默认值为None= 检测到的预算,而非不设上限)与from_pgen(max_mem=)使用相同的整个进程并发规划器预算(见上文)。genoray write svar1(SparseVar2.from_svar1):source是*.svar(SVAR1)目录。--fields(可重复)用于选择要保留的 SVAR1 FORMAT 字段(默认:全部);--empty-fields会覆盖--fields,不保留任何字段。--chunk-size默认使用根据内存推导的值(None)。genoray write-svar1(顶层旧版命令,SparseVar.from_vcf/from_pgen→ SVAR 1.0):沿用之前的行为,不变——VCF 或 PGEN 源(自动检测)、--dosages(VCF 的 FORMAT 字段名,或 PGEN 的剂量.pgen路径)、--max-mem(默认值为"1g")、--haploid、--no-symbolic/--no-breakend(此处为独立标志,与 SVAR2write子命令的单个--skip-symbolics-and-breakends不同)、--threads/-@、--overwrite。不支持--regions/--samples/--fields/--reference/--check-ref/--progress/--log-level/--log-filter,也不支持任何基于Tuning的标志(--dense-cap/--merge-threads/--sample-interval/--concurrent-chroms/--reader-workers/--overshard)——这些都仅适用于 SVAR2 的write(旧版SparseVar.from_vcf/from_pgen后端不接受progress=/log_level=/log_filter=/tuning=)。
genoray view
# SVAR2 (default) — thin CLI over SparseVar2.write_view
genoray view in.svar2 out.svar2 -r chr1:1-1000 -s A,B
genoray view in.svar2 out.svar2 -r chr1:1-1000 # all samples
genoray view in.svar2 out.svar2 -s A,B # all variants (one region per contig)
genoray view in.svar2 out.svar2 -r chr1:1-1000 --no-reroute # representation-preserving, low-memory view
genoray view in.svar2 out.svar2 -r chr1:1-1000 --reroute # force the size-optimal re-route
# SVAR 1.0 (previous default) — unchanged SparseVar.write_view CLI
genoray view svar1 in.svar out.svar -r chr1:1-1000 -s A,B --progress
两个子命令共用 -r/--regions、-R/--regions-file、
-s/--samples、-S/--samples-file、-f/--fields、--merge-overlapping、
--regions-overlap、--overwrite、-@/--threads、--progress/
--no-progress 选项,也共用相同的空操作保护(必须至少指定 regions 或
samples 之一)和互斥检查(--regions/--regions-file
和 --samples/--samples-file 分别互斥)。genoray view
(SVAR2)还提供 --log-level LEVEL(与上文 write 中的六个名称相同)和
--log-filter DIRECTIVE(语义与上文 write 相同);genoray view svar1 则没有这些选项,因为其 SparseVar.write_view 后端没有
log_level=/log_filter= 参数。两个 view 子命令都没有基于 Tuning 的标志——
write_view 不是转换流水线,因此 Tuning 的调节项均不适用。
-
genoray view(SVAR2,是对SparseVar2.write_view的轻量封装):省略--regions/--regions-file时,“所有变异”默认按每条 contig 一个区域处理 (SparseVar2.contigs,因为 SVAR2 没有 contig 长度元数据),范围为[0, 2**31 - 1)——所有真实 POS 都小于该值。--fields默认值为None,表示不保留任何字段(仅基因型)——即使存储中有 INFO/FORMAT 字段, 也始终可以成功。--reroute和--no-reroute都使用相同的切片后端, 并以相同方式传递--fields/--reference——两者不再分别对应 保留字段与仅保留基因型:--reroute会在子集上重新运行 var_key/dense 路由成本模型—— 大小最优(根据子集的样本数/携带者数,将每个变异重新路由到较小的表示形式)。--no-reroute(reroute=False)直接切片每个变异现有的 磁盘表示(不运行成本模型,按字节级切片)——无论子集的样本数/携带者数如何, 都保留原表示形式。建议用于体细胞/全稀有变异队列(几乎每个变异都已经通过 var_key 路由)或内存受限的运行。- 两个标志都省略(默认)时为
"auto":只要携带任何 FORMAT 字段,就采用--no-reroute的行为;否则采用--reroute的行为。原因:dense→var_key 转换每个携带者调用只存储一个值,且没有位置存放非携带者样本的 FORMAT 值,因此在 携带 FORMAT 字段的视图中重新路由源 dense 变异会导致该值被静默丢弃——只要涉及 FORMAT,"auto"就优先保证保真度;否则(仅基因型/仅 INFO 视图)选择大小最优的 重新路由,因为没有会丢失的逐样本位置。
--reference(在子集上从头重新计算mutcat)和-@/--threads(限制并发切片的 contig 数;省略时自动检测)在--reroute和--no-reroute下都确实生效——两条路径都不再存在“接受但忽略/未使用” 的情况。--progress、--log-level和--log-filter在此处也都确实生效—— 有关逐 contig 一行的粗粒度显示和日志级别语义,请参阅下文的“write_view进度条”。底层write_view的reroute=参数只接受"auto"、True或False——任何其他值(例如reroute=1)都会引发ValueError,而不会静默地转由reroute=False切片器处理。 -
genoray view svar1:SVAR 1.0 行为不变——“所有变异”默认根据SparseVar的_contig_stats确定(每条 contig 的范围为[0, pos_max + 1));--fields默认包含所有可用字段(使用显式的空选择可不携带任何字段);没有--reference/--reroute/--log-level选项;--progress显示真正的阶段级进度条(见下文)。
genoray concat / genoray split
genoray concat merged.svar2 part1.svar2 part2.svar2 # disjoint-contig merge
genoray split in.svar2 out_dir/ # explode into out_dir/{contig}.svar2
genoray split in.svar2 subset.svar2 --contigs chr1,chr2 # subset into one store
两者都接受 --mode(Literal["copy", "hardlink", "symlink", "move"],默认值为
"copy"——参阅 SparseVar2.concat/split_by_contig/subset_contigs
文档字符串)和 --overwrite。
过滤
VCF:将 genoray.Filter(record=, expr=) 值对象传给 filter=。
record 是一个 Callable[[cyvcf2.Variant], bool],在基因型扫描期间应用;expr 是匹配的 polars pl.Expr,应用于 .gvi 索引 — VCF 要求两部分都提供,并将它们打包在一起,以免两者不一致。
要在构造后更改 VCF 的过滤器,请将 Filter(或用于清除过滤器的 None)赋给 vcf.filter setter;内存中的索引会失效。
getter 返回当前生效的 Filter | None,因此 vcf.filter = vcf.filter 可以往返赋值。
from genoray import VCF, Filter
vcf = VCF("file.vcf", filter=Filter(
record=lambda v: not v.INFO.get("SVTYPE"), # cyvcf2 record predicate
expr=~genoray.exprs.is_symbolic, # matching .gvi index predicate
))
vcf.filter = None # clear
f = vcf.filter # -> Filter | None
原先的双参数构造方式(在 filter= 旁单独传入 polars 表达式的关键字参数)以及返回元组的 vcf.filter getter/setter 已在 3.0.0 中移除 — 请将任何分别传入记录谓词和 polars 表达式的代码迁移为上方所示的单个 Filter(record=, expr=) 对象。
PGEN:传入返回布尔掩码并作用于 .gvi 索引列的 polars pl.Expr。genoray.exprs 中的内置表达式(完整列表):
is_snp(当所有 ALT 等位基因的 ILEN == 0 时为 True;含有任何nullILEN 的行 → False)is_indel(当所有 ALT 等位基因的 ILEN != 0 时为 True;含有任何nullILEN 的行 → False)is_biallelicis_symbolic(当任意 ALT 是 VCF 4.x 符号等位基因,即以<开头时为 True)is_breakend(当任意 ALT 是 VCF 4.x 断点等位基因,采用配对断点 / 单断点表示法时为 True,例如G[chr2:321[、]chr2:321]G、.TGCA、TGCA.。它是与符号等位基因不同的一类 ALT —is_symbolic不会标记断点)is_imprecise(当任意 ALT 的 ILEN 为null时为 True — 表示无法确定大小的符号等位基因或断点)ILEN(一个List[Int32]表达式 — 每个 ALT 等位基因对应一个值,不是布尔值)
符号 SV 的 ILEN 语义。 对于精确的 <DEL> / <INS> / <DUP>,ILEN 在构建索引时根据 INFO 字段计算:<DEL> 使用 -|SVLEN|,<INS> / <DUP> 使用 +|SVLEN|(若缺少 SVLEN,则回退为 |END - POS|)。对于 VCF,INFO 字段通过 oxbow 从 header 声明的列中读取;对于 PGEN,则从 PVAR INFO 字符串中解析。非符号 ALT 使用字面值 len(ALT) - len(REF)。
无法确定大小的符号等位基因的 ILEN 为 null。 以下情况会导致等位基因无法确定大小:设置了 IMPRECISE INFO 标记;SVLEN / END 均缺失;符号类型不受支持(<BND>、<CNV>、<INV>、<*> / <NON_REF>);或者 ALT 是采用配对断点 / 单断点表示法的断点(例如 G[chr2:321[)。在 NumPy 物化时,null ILEN 会强制转换为 0(视为点变异)。
过滤建议(使用 filter= — PGEN 使用裸 pl.Expr,VCF 使用 genoray.Filter):
-
~genoray.exprs.is_symbolic— 丢弃所有符号等位基因(无论是否精确)。 适用于无法将任何符号 ALT 展开为字面序列的单倍型消费者(例如genvarloader):# PGEN pgen = genoray.PGEN("file.pgen", filter=~genoray.exprs.is_symbolic) # VCF (both halves required, bundled in a Filter) vcf = genoray.VCF( "file.vcf.gz", filter=genoray.Filter( record=lambda rec: not any(a.startswith("<") for a in rec.ALT), expr=~genoray.exprs.is_symbolic, ), ) -
~genoray.exprs.is_imprecise— 保留精确的符号 SV(大小 / 跨度准确),仅丢弃无法确定大小的等位基因(包括始终无法确定大小的断点)。适用于可查询精确 SV 的范围 / 重叠查询:pgen = genoray.PGEN("file.pgen", filter=~genoray.exprs.is_imprecise) -
对于单倍型消费者,丢弃所有无法展开的 ALT(符号等位基因和断点)—
~is_symbolic无法捕获断点:hap_safe = ~genoray.exprs.is_symbolic & ~genoray.exprs.is_breakend pgen = genoray.PGEN("file.pgen", filter=hap_safe)
其他情况下,请针对 .gvi 架构编写 pl.col(...) — 请阅读 genoray/exprs.py 了解可用列。使用 & / | 组合两个 exprs 表达式时无需导入 polars;只有在构建自定义谓词时才需要 import polars as pl。
参考 — 快速参考
genoray.Reference 是一个基于 pysam 的索引 FASTA 读取器,用于为突变目录分类提供侧翼上下文。
ref = genoray.Reference.from_path("hg38.fa") # auto-creates .fai if absent
ref = genoray.Reference.from_path("hg38.fa", contigs=["chr1", "chr2"])
seq: np.ndarray = ref.fetch("chr1", start=1_000_000, end=1_000_010)
# returns uint8 NDArray, 0-based half-open [start, end)
# bytes(seq) gives the ASCII sequence
主要属性:
from_path(fasta, contigs=None)—fasta是一个str | Path;如果缺少.fai索引,则自动调用pysam.faidx。contigs用于筛选调用方关注的 contig(默认为 FASTA 中的所有 contig)。fetch(contig, start, end)— 采用从 0 开始、左闭右开的区间[start, end)。超出 contig 范围的位置会以 N 填充。返回NDArray[np.uint8]。contig_array(contig)— 返回完整 contig 序列,作为缓存的NDArray[np.uint8]。与fetch共用单个 contig 的内存缓存。接受带chr前缀或不带前缀的名称。- 不依赖 contig 名称格式:
"chr1"和"1"均可正确解析(底层使用ContigNormalizer)。 - 一次在内存中缓存一个 contig;按 contig 顺序访问时效率很高。
突变目录(SBS-96 / DBS-78 / ID-83)
write_view 进度条
SparseVar.write_view(..., progress=False) 支持选择启用的 progress 关键字参数。
当设为 True 时,写入视图期间会显示阶段级 rich 进度条
(每个主要步骤更新一次:计数、基因型、每个保留字段、索引
构建,以及指定 reference= 时的突变注释)。默认值为
False — 不显示进度条,也没有额外开销 — 因此不会影响库和流水线调用方。
genoray view svar1 CLI 提供相同选项 --progress
(同样默认为关闭):
genoray view svar1 in.svar out.svar -r chr1:1-1000 -s A,B --progress
进度条仅用于展示:无论是否启用,输出字节、架构和数据类型都相同。
SparseVar2.write_view(..., progress=False, log_level="info", log_filter=None)
以与 from_* 写入器相同的方式显示实时写入进度(见上文
“转换”),但有一点不同:write_view 不像 from_* 写入器那样有可供采样的逐记录数据流,因此其进度是粗粒度的 — 每个 contig 输出一行,不会在 contig 内更新进度条。在终端或
Jupyter 中,progress=True 会显示实时更新的进行中/已完成
contig 列表;在其他环境中,每个 contig 完成时会输出一条简洁的“chrom done”行。无论 progress 为何值,每个 contig 完成后都会输出一行 "[svar2] chrom done: N kept, 0 excluded (Ts)" 摘要,除非 log_level="off"(切片不会排除变异,因此
excluded 始终为 0)。log_level(包括 "critical" 别名和不接受 "warn")与 log_filter 的行为和上文的 from_* 写入器相同;请注意,write_view 没有 tuning= 参数 — 它不是转换流水线,因此 Tuning 的所有参数均不适用。genoray view CLI 提供 --progress/--no-progress、--log-level 和 --log-filter(见上文
genoray view CLI 章节);genoray view svar1(SVAR 1.0)仅提供
--progress — 其 SparseVar.write_view 后端(见上文)没有
log_level/log_filter 参数。
原子且崩溃安全的写入
写入具备崩溃安全性和原子性。from_vcf、from_pgen 和 write_view 会在隐藏的同级暂存目录(输出旁的 .<name>.tmp…)中构建
.svar 目录,并且仅在写入完全成功后才原子地重命名到目标位置;.gvi 索引文件也以相同方式写入。写入过程中发生崩溃不会留下不完整或损坏的输出;覆盖现有输出时,在替换内容完整写入前会保留原输出。输出字节不变 — 这仅是一项持久性保证。
概述
SparseVar 支持 COSMIC 风格的突变目录。工作流如下:
- 调用一次
svar.annotate_mutations(reference),对每个变异进行分类,并将mutcat.npy写入.svar目录。 - 调用
svar.mutation_matrix(kind),获取每个样本的计数矩阵。
SparseVar.annotate_mutations
svar = genoray.SparseVar("out.svar")
ref = genoray.Reference.from_path("hg38.fa")
svar.annotate_mutations(ref) # write_back=True (default)
svar.annotate_mutations(ref, write_back=False) # in-memory only; not persisted
svar.annotate_mutations("hg38.fa") # path accepted directly
签名:annotate_mutations(reference, *, contigs=None, write_back=True) -> None
reference—genoray.Reference实例或 FASTA 文件路径 (通过Reference.from_path自动包装)。contigs=None— 如果指定(contig 名称列表),则只分类这些 contig 上的变异;所有其他 contig 上的条目都标记为NOT_ANNOTATED(哨兵值-4),且不会从参考序列中获取这些 contig。 名称通过ContigNormalizer匹配(chr1/1均可)。.svar索引中不存在的指定 contig 会跳过并发出警告;索引中存在但 参考序列中不存在的已列出 contig 会引发异常(若要干净地排除它, 请将其从列表中省略)。None(默认值)会分类所有 contig。当write_back=True时,规范化后的范围会记录在metadata.json的mutcat_contigs中(None= 全部)。write_back=True— 持久化mutcat.npy并更新metadata.json,以便 后续通过SparseVar(dir, fields=["mutcat"])打开时能读取该字段。 注意:write_view绝不会按位置将mutcat复制到输出文件 (见下文);请将reference=传给write_view,以便在子集上重新计算, 或者自行对输出视图调用annotate_mutations。(这是SparseVar/v1 的行为;对于SparseVar2.write_view,无论reroute=True还是reroute=False,传入reference=都会从头开始在 子集上重新计算mutcat。)write_back=False—mutcat字段仅驻留在内存中 (svar.fields["mutcat"]);重新打开文件时将无法找到它。- 调用后,无论
write_back的值为何,都会填充svar.fields["mutcat"]。
分类内容:
| 变异类型 | 通道 |
|---|---|
| 单独的 SNV | SBS-96(三核苷酸上下文) |
| 同一单倍型上的相邻 SNV 对 | DBS-78(5' 条目 = DBS 代码,3' 条目 = DBS_PARTNER 哨兵值) |
| 连续 ≥ 3 个相邻 SNV | SBS(每个仍单独处理;不合并为 DBS) |
| VCF 中的原生 2 bp MNV | DBS-78 |
| > 2 bp 的 MNV、符号型变异、非 ACGT 变异 | UNCLASSIFIED |
| 插入/缺失 | ID-83(按大小、重复上下文分箱) |
位于 contigs= 之外的 contig 上的变异 | NOT_ANNOTATED(从所有矩阵中排除) |
SparseVar.mutation_matrix
svar = genoray.SparseVar("out.svar", fields=["mutcat"]) # pre-load field
df = svar.mutation_matrix("SBS96") # default count="allele"
df = svar.mutation_matrix("DBS78", count="sample")
df = svar.mutation_matrix("ID83", count="allele")
签名:mutation_matrix(kind, *, count="allele") -> pl.DataFrame
kind—"SBS96"、"DBS78"、"ID83"之一。count="allele"— 统计每个非参考等位基因拷贝(双倍体纯合 = 2)。count="sample"— 每个样本中每个类别最多统计一次(存在/不存在)。- 返回一个 Polars
DataFrame,包含一个MutationType字符串列,其后是 每个样本对应的一个Int64列。行按固定的 COSMIC 编码表顺序排列 (分别为 96 / 78 / 83 行)。 - 需要
mutcat字段可用:在打开时通过fields=["mutcat"]加载, 或者该字段已因之前调用annotate_mutations而存在于内存中, 或者之前调用annotate_mutations(write_back=True)后已写入磁盘。 若以上条件均不满足,则引发ValueError。
mutcat 字段
mutcat 是一个按基因型条目存储的 int16 字段(布局与
genos 相同,采用不规则数组)。int16 代码空间如下:
| 范围 | 通道 |
|---|---|
[0, 96) | SBS-96 |
[96, 174) | DBS-78 |
[174, 257) | ID-83 |
-1 | DBS_PARTNER — 相邻 SNV 对的 3' 半部分;永不计数 |
-2 | UNCLASSIFIED — 符号型 / 复杂变异 / > 2 bp 的 MNV / 非 ACGT 变异 |
-3 | MISSING — 保留的哨兵值(在代码空间中定义,但 v1 的 annotate_mutations 不会生成;SparseVar 只存储携带 ALT 的条目,因此无调用位点不会出现在不规则字段中) |
-4 | NOT_ANNOTATED — 位于 contigs= 注释范围之外的 contig 上的条目;永不计数 |
读取先前已注释的文件:
svar = genoray.SparseVar("out.svar", fields=["mutcat"])
# svar.fields["mutcat"] is a Ragged[int16] mirroring svar.genos
v1 范围限制(无链偏倚;根据 PCAWG/SigProfiler 规则校准)
- 不区分链偏倚(SBS-192 / SBS-384)——v1 仅支持
SparseVar。对于解析到转录链的目录,请使用SparseVar2.annotate_mutations(reference, gtf=...)。 - DBS 合并仅适用于同一单倍型上的孤立相邻碱基对。≥ 3 个连续相邻 SNV 的片段仍作为单独的 SBS 条目保留。
- Indel 通道(ID-83)分箱遵循 PCAWG/SigProfiler 发布的规则,并由
tests/test_mutcat.py中的单元测试固定。与 SigProfilerMatrixGenerator 的交叉验证暂不进行(它不是已声明的依赖项)。
签名拟合(COSMIC)
将目录分解为每个样本的 COSMIC 签名活性。
import genoray
from genoray import Forward, Spa
ref = genoray.cosmic_signatures("SBS96") # pooch-fetched + cached
cat = svar.mutation_matrix("SBS96") # MutationType + sample cols
act = genoray.fit_signatures(cat, ref) # activities + cosine_similarity
# strategy=: pick the refit algorithm and its parameters explicitly
act = genoray.fit_signatures(cat, ref, strategy=Forward(max_delta=0.02))
act = genoray.fit_signatures(cat, ref, strategy=Spa()) # SPA's cosmic_fit, faithfully
# convenience: mutation_matrix -> fit_signatures in one call
act = svar.assign_signatures("SBS96") # default COSMIC ref
act = svar.assign_signatures("SBS96", reference=ref, min_activity=0.01)
act = svar.assign_signatures("SBS96", reference="my_sigs.txt") # TSV path
act = svar.assign_signatures("SBS96", strategy=Spa()) # same strategy= as fit_signatures
签名:
cosmic_signatures(kind, *, version="3.4", genome="GRCh38") -> pl.DataFrame— 获取并缓存kind ∈ {"SBS96","DBS78","ID83"}对应的 COSMIC 参考集。 返回一个MutationType列(按规范码本顺序排列)以及每个签名对应的一列。genome对ID83无效。fit_signatures(catalogue, reference, *, strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame— 根据reference对catalogue重新拟合。通过联结MutationType对齐行(如果目录包含参考中缺失的类型,则引发ValueError)。每个样本返回一行:Sample、每个签名对应的一个 Float 列(活性;未选中的签名值为0.0),以及cosine_similarity。n_jobs=1(默认值)为串行;n_jobs=-1使用所有核心。无论n_jobs/backend如何设置,结果都相同。strategy=用于选择算法,不能与max_delta/min_activity/criterion组合使用(如果两者都给出,即使这些参数取默认值,也会引发ValueError)——这三个参数始终是strategy=Forward(...)的简写形式,并未弃用:Forward(max_delta=0.01, min_activity=0.005, criterion="cosine")(省略strategy时的默认值)——genoray 自有的稀疏贪心前向选择:NNLS + 由 criterion 指导的添加 + 最小活性剪枝。criterion="cosine"在余弦相似度进入平台期时停止(与尺度无关,因此无法反映突变负担);criterion="bic"则依据 Poisson-BIC 检验停止(会考虑负担,在突变数高于约 1,000 时更一致,但在负担极低时会过度选择)。返回的活性是原始 NNLS 权重——不一定总和为某个特定值。Spa(metric="l2", initial_remove_penalty=0.05, add_penalty=0.05, remove_penalty=0.01, background_sigs=("SBS1","SBS5"), connected_sigs=True, activity_scale="burden")——从头重新实现 SigProfilerAssignment 的cosmic_fit(solver="nnls"、pcawg_rule=False),纯 numpy/scipy/polars 实现,不依赖 SigProfiler:遍历每个参考特征签名进行饱和拟合,依据相对 L2 误差进行反向剪枝,然后通过增删层进行优化。background_sigs中的名称会在每个优化层开始时强制加入工作集,但不会受到绝对保护而免于移除扫描:SPA 自身的保护会在扫描过程中逐渐减弱,genoray 也复现了这一行为,因此样本不需要的背景特征签名仍会被移除(reference中不存在的名称会被忽略,因此在 DBS78/ID83 参考集上,默认的 SBS1/SBS5 不会产生影响);connected_sigs=True会在选中任一成员后强制添加共同出现的 SBS 配对特征签名(例如 SBS2/SBS13),False会禁用此功能,也可以传入自定义的Sequence[Sequence[str]],其中包含特征签名名称组。默认的activity_scale="burden"会重新缩放并取整活性,使每个样本的特征签名列之和等于其突变负担;activity_scale="raw"则返回未经缩放的 NNLS 权重,与Forward的行为相同。结果具有确定性(不涉及任何 RNG),这一点不同于上游 SPA。在所有测试的突变负担下,恢复的真实特征签名数量都多于Forward,但与Forward的"cosine"criterion 一样,它与尺度无关,因此也不是一致估计量。每个样本的运行速度大约比Forward慢一个数量级(在真实 COSMIC SBS96、86 个特征签名的数据上测得约慢 7–27 倍,具体取决于目录)。
SparseVar.assign_signatures(kind, *, reference=None, count="allele", strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame——先调用mutation_matrix(kind, count=...),再调用fit_signatures(...)。reference可接受pl.DataFrame、TSV 路径或None(默认使用cosmic_signatures(kind))。 将n_jobs/backend转发给fit_signatures,以实现按样本并行处理(n_jobs=1(默认值)表示串行;n_jobs=-1表示使用所有核心)。strategy=和criterion=的转发方式与fit_signatures本身完全相同——传入strategy=Spa()可使用 SigProfilerAssignment 的算法;也可以让strategy保持为None(默认值),并通过max_delta/min_activity/criterion配置Forward简写形式。将strategy=与max_delta/min_activity/criterion中任意参数同时传入时,会引发ValueError,与fit_signatures完全相同(未传入的关键字参数使用默认值)。SparseVar2.assign_signatures的签名完全相同,但还接受contigs=,用于将目录限制为存储中的部分 contig。
重拟合策略
fit_signatures(strategy=...) 用于选择算法。两种策略都是冻结的数据类,只包含各自的参数。
Forward(max_delta=0.01, min_activity=0.005, criterion="cosine")— 从空集合开始进行贪心正向选择。默认策略。活性是原始 NNLS 权重,不必与样本的突变负荷相加相等。Spa(metric="l2", initial_remove_penalty=0.05, add_penalty=0.05, remove_penalty=0.01, background_sigs=("SBS1","SBS5"), connected_sigs=True, activity_scale="burden")— SigProfilerAssignment 的cosmic_fit。所有默认值均来自 SPA。在默认设置activity_scale="burden"下,活性是整数值浮点数,总和与样本的突变负荷完全相等;activity_scale="raw"则会返回未缩放的 NNLS 权重。
参考集中不存在的 background_sigs 名称会被忽略,因此默认设置对 DBS78 和 ID83 不起作用。connected_sigs=True 使用 SPA 的四个 SBS 组(genoray.SPA_CONNECTED_GROUPS:SBS2/13、SBS7a-d、SBS10a/b、SBS17a/b);可传入由序列组成的自定义分组,或传入 False 以禁用。
Forward(criterion="cosine") 和 Spa 都不会考虑突变负荷。对于需要一致性的全基因组目录,请使用 Forward(criterion="bic")。
输出模式不受策略影响:Sample、每个参考特征一个 Float 列,以及末尾的 cosine_similarity。
范围之外的功能(v1):从头提取、机会归一化、Bootstrap 置信区间、绘图。
常见错误
| 错误 | 修正 |
|---|---|
genoray.Genos8 | genoray.VCF.Genos8(类属性) |
vcf.read(..., phasing=True) | 在 VCF() 构造函数中设置 phasing=True |
从 VCF 读取剂量,却未设置 dosage_field= | 在构造函数中传入 dosage_field="DS"(或适当的 Number=A 字段) |
| 同时有硬调用时,将仅含剂量的 PGEN 放在主路径中 | 主路径放硬调用,剂量放在 dosage_path= 中 |
导入 from genoray._vcf import VCF | 使用 from genoray import VCF |
期望 VCF 具有 read_ranges | VCF 不支持此功能;循环调用单范围 read,或使用 PGEN/SparseVar |
将 svar.index["POS"] 当作从 0 开始的坐标 | 它从 1 开始;与查询坐标比较时减去 1 |
调用 read_ranges 后假设返回扁平数组 | PGEN 返回 (data, offsets);SparseVar 返回 Ragged(或带有 fields 的 awkward 记录) |
调用 mutation_matrix 时没有 mutcat 字段 | 先运行 annotate_mutations,或使用 fields=["mutcat"] 打开 |
期望 mutation_matrix 自动运行注释 | 它不会自动运行;请单独调用 annotate_mutations |
重新打开 SparseVar 后丢失 mutcat 字段 | 在 annotate_mutations 中使用 write_back=True(默认值);然后使用 SparseVar(dir, fields=["mutcat"]) 打开 |
调用 write_view 并期望输出中包含 mutcat | write_view 绝不会按位置复制 mutcat(子集化会使 DBS 邻接码失效)。向 write_view 传入 reference=,以便在子集上重新计算;或者自行对输出视图调用 annotate_mutations。在未传入 reference= 时,显式在 fields= 中包含 "mutcat" 会引发 ValueError。对于 SparseVar2.write_view,无论 reroute=True 还是 reroute=False,reference= 都会从头重新计算 mutcat。 |
将 FORMAT 字段传入 SparseVar2.write_view(..., reroute=True) 的 fields=,并期望该字段被丢弃或拒绝 | 现在 reroute=True 和 reroute=False 都会保留 fields(此前 reroute=True 会引发 ValueError)。应留意默认值 reroute="auto":只要保留了任何 FORMAT 字段,它就会解析为 reroute=False,因为 dense→var_key 转换时,没有位置存放非载体样本的 FORMAT 值。 |
将源数据集目录作为 output 传给 write_view(即使设置了 overwrite=True) | 会引发 ValueError——原地写入会在视图写入前删除源数据。请传入不同的 output 路径。 |
将 FASTA 路径直接传给 annotate_mutations | 支持——它会通过 Reference.from_path 自动封装 |
在 SparseVar2.decode() 的结果上使用 rag.fields["AF"] | Ragged.fields 是包含名称的 list[str],不是映射;请使用 rag["AF"] 索引字段本身 |
期望 SparseVar2.decode() 包含 INFO/FORMAT 字段 | 字段需显式启用——向 SparseVar2(...) 传入 fields=[...],或先调用 .with_fields([...]) |
何时需要更新此技能
任何新增、删除、重命名或更改公共名称(即无需下划线即可通过
import genoray 访问的名称)语义的 PR,都必须同时更新此技能
和代码更改。请参阅项目的 CLAUDE.md。