mcpField
G
genoray-api

genoray-api

github.com/d-laubデータ・データベース6128 回インストール5 週間前に更新
リポジトリ

genoray-api について

genorayは、VCF/BCF、PLINK 2 PGEN、疎なバリアントストアを扱う、NumPyを中心としたゲノムバリアント検索ツールキットです。領域検索、変換、フィルタリング、変異シグネチャ解析を行う開発者や研究者向けです。

主な機能

  1. 01VCF/BCFとPGENの読み込み
  2. 02ゲノム領域の検索
  3. 03バリアントデータをSVAR2に変換
  4. 04式によるバリアントの絞り込み
  5. 05指定フィールドと用量の抽出
  6. 06変異シグネチャの解析

活用シーン

  1. 01複数のゲノム領域で変異を検索
  2. 02コホートを疎なSVAR2ストアに変換
  3. 03バリアントフィールドを絞り込み抽出
  4. 04変異を分類してシグネチャをフィット

インストール

コマンドをコピーしてターミナルで実行します。

$npx skills add https://github.com/d-laub/genoray --skill genoray-api

インストール後、Claude Code で/genoray-apiを入力して読み込まれたか確認します。

詳細

genoray 公開 API

genoray は、VCF/BCF(cyvcf2)、PGEN(pgenlib)、およびスパース memmap 形式(SparseVar / .svar)に対応する、NumPy を優先した範囲クエリーレイヤーです。

公開インターフェース

import genoray で公開されるものは、以下に限られます。

  • genoray.PGEN — PLINK 2 PGEN リーダー
  • genoray.Reference — インデックス付き FASTA リファレンスゲノムリーダー
  • genoray.VCF — VCF/BCF リーダー
  • genoray.Filter — cyvcf2 レコード述語(record)と、それに対応する .gvi polars 式(expr)をまとめた VCF フィルター値オブジェクト
  • 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-way マージによる N 個の単一サンプル VCF/BCF → 1 つの SVAR2 ストアへの変換(reference/no_reference は from_vcf と同様にサポートし、欠落サイトは hom-ref で埋めます。regions=/merge_overlapping=/regions_overlap= をサポートしますが、samples= はサポートしません。コホートはファイルセットです)、from_vcf/from_pgen と同様に regions=/samples=/merge_overlapping=/regions_overlap= をサポートする from_svar1 による SVAR1(SparseVar)→ SVAR2 のネイティブ移行(VCF/htslib は読み込みません。二対立遺伝子の SVAR1 のみ対応)、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 を指定します)。なお、PGEN には VCF INFO/FORMAT がないため、info_fields=/format_fields= はサポートしません。from_svar1 は fields= により SVAR1 の既存フィールドを選択的に引き継ぎます(既定の None = すべて、[] = なし、または名前のサブセット)。読み戻しは fields=/with_fields/available_fields でオプトインし、decode の結果に付加できます)
  • genoray.InfoField / genoray.FormatField — SparseVar2.from_vcf の INFO/FORMAT フィールドを 1 つ設定する frozen dataclass(name、dtype=None、default=None)。単独の str 名を指定すると、推定された既定値が使われます
  • genoray.DosageField — SparseVar2.from_pgen の PGEN 用量 FORMAT フィールドを設定する frozen dataclass(name="dosage"、source="self"|Path、dtype="f16"|"f32"="f32"、default=None)
  • genoray.Tuning — 明示的なスケジューリング設定 6 項目(concurrent_chroms、reader_workers、overshard、dense_cap、merge_threads、sample_interval。すべて int | None)を持つ frozen dataclass。すべての SparseVar2.from_* ライターに tuning= として渡します。従来の環境変数による設定(「SparseVar2 — quick reference」の「Tuning」と、docs/source/svar.md の「migrating from environment variables」を参照)に代わるものです。環境変数で genoray を設定することはできません
  • genoray.exprs — .gvi インデックス用 polars フィルター式
  • genoray.cosmic_signatures — COSMIC リファレンスシグネチャーの取得/キャッシュ
  • genoray.fit_signatures — 変異カタログのシグネチャー再フィット。strategy= で Forward()(既定。スパースな前方選択)または Spa()(SigProfilerAssignment に忠実な後方削除)を選択します
  • genoray.Forward / genoray.Spa — fit_signatures の strategy dataclass(frozen)
  • genoray.Strategy / genoray.Criterion / genoray.Metric / genoray.ActivityScale — 関連する型エイリアス。genoray.SPA_CONNECTED_GROUPS には Spa の既定の共起シグネチャーグループ 4 つが格納されています

これ以外は公開されていません。genoray._vcf など、_ で始まるものは内部用です。ユーザーコードからインポートしないでください。

詳細の参照先

推測するのではなく、以下を優先して参照してください。

  • docs/source/index.md — 完全な使用例(VCF、PGEN、フィルタリング、チャンク分割)を含む解説
  • docs/source/svar.md — SparseVar の使用方法
  • genoray/__init__.py — 公開 API の確認用
  • 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) — 同じ範囲を読み取ったときに返されるバリアントの識別情報を、range(u32。行が属するクエリ範囲)、CHROM、POS(1-based)、REF、ALT の pl.DataFrame として返します。これは read/read_ranges のバリアント軸と 1:1 で対応し、リーダーの filter が適用されます。ここでの ALT は、.gvi での保存形式にかかわらず、常に list[str] です。また、CHROM はクエリに指定した contig と異なる場合でも、ファイル内の表記どおりのコンティグ名です。生成元の PGEN がなくなったあとも成果物を使う必要がある場合は、var_idxs と非公開の _index の代わりにこちらを使用してください
  • genoray/_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] は、ハードコール .pgen 自体(source="self")または別の .pgen から読み取ったサンプルごとの用量トラックを FORMAT フィールドとして保存します。from_vcf と同様に regions=/samples=/merge_overlapping=/regions_overlap= に対応)、from_vcf_list(単一サンプルの VCF/BCF を N 個受け取り、ネイティブな k-way マージで 1 つの 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= に対応しますが、範囲を事前に絞り込むのではなく、レコードごとに regions フィルターを適用します)、n_samples/available_samples/contigs/ploidy メタデータ。読み取り/クエリメソッドはミックスインにあります:genoray/_svar2_decode.py(decode — 選択したフィールドごとに 1 つの 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 — .gvi に付属するフィールドではなく、コンティグごとの Rust サイドカーを使用する、COSMIC 変異シグネチャのワークフロー。SparseVar のものと同様です)
  • genoray/_svar2_fields.py — InfoField/FormatField/DosageField データクラスと FieldDtype、および from_vcf(info_fields=, format_fields=) で使用されるヘッダー/データ型の検証;_parse_cli_field_specs(内部用 — bcftools 形式の INFO/x/FORMAT/x/FMT/x CLI フィールド文字列を解析します。genoray write vcf --fields CLI で使用されます);StoredField(固定データクラス:name、category、dtype、default、key)は、SparseVar2.available_fields が返す読み取り側マニフェストのエントリー型です。トップレベルの genoray からはエクスポートされず、この辞書を通してのみ参照できます
  • genoray/_cli/__main__.py — genoray CLI(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_array
  • genoray/exprs.py — 事前構築済みフィルター式の全セット(現在 7 種:is_snp、is_indel、is_biallelic、is_symbolic、is_breakend、is_imprecise、ILEN) 署名、キーワード引数、または形状が不明な場合は、第一原理から推測するのではなく、ソース内の docstring を読んでください。

全体に適用される規約

  • 範囲は 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 など)で使われる場合は 1 つのチャンクに直接上限を設けるため、同じ意味だと思い込む前に、下記「変換」の各項目を確認してください。
  • 染色体名は自動的に正規化されます。ファイルの表記にかかわらず、"chr1" と "1" のどちらも使用できます(ContigNormalizer)。
  • 欠損ジェノタイプは -1(int)、欠損ドサージュは np.nan(float32)です。
  • デフォルトの倍数性は 2 です。SparseVar.from_vcf/from_pgen(および genoray write-svar1)は haploid=True/--haploid を受け付けます。これはハプロタイプを OR で集約してサンプルごとに 1 つの半数体コールにし、メタデータに ploidy=1 を記録します(非フェーズ化体細胞データを想定しています)。
  • 返される配列はすべて NumPy 配列です。mode によって返される配列が選択されます。

サンプルアクセサー — 正式名と使用方法が異なる理由

available_samples(list[str])は、「ファイル内の全サンプル」を取得する正式なアクセサーで、4 つのリーダー(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= キーワード引数を受け取ります。

2 つの使用方法が異なる理由(パフォーマンス): VCF/PGEN でサンプルをサブセット化するとコストがかかります。バックエンドのリーダーが再初期化されるため、意図的に状態を変更する set_samples() 呼び出しを 1 回行い、その後の読み取りで再利用します。一方、SparseVar/SparseVar2 ではサブセット化のコストはほぼありません(メモリーマップ済みデータに対するインデックス選択にすぎません)。そのため、永続的なリーダー状態ではなく、軽量な呼び出しごとの samples= キーワード引数として提供されています。これは意図的な違いであり、不整合ではありません。一方を他方に合わせて「修正」しないでください。

モード定数 — 注意点

モードはトップレベルの名前ではなく、クラス属性です。

genoray.VCF.Genos8           # not genoray.Genos8
genoray.PGEN.GenosPhasingDosages

クラスで利用可能なモードを確認するには、_vcf.py/_pgen.py のクラス本体を読んでください(先頭付近の Genos を検索します)。

1 つのモードに複数の配列が含まれる場合、返されるタプルの順序は定数名の順序に従います。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 に指定してください。パスを 1 つだけ渡した場合、ハードコールと用量データの両方がそこから読み込まれます(ハードコールは用量のしきい値から推定されます。詳細は 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 は、ソースのフィルターを継承して適用します。SVAR をフィルタリングするには、VCF/PGEN をフィルタリングしてください。

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] を返します。

ragged の結果に含まれる各リーフ値はバリアントインデックスです。これは svar.index(少なくとも CHROM, POS, REF, ALT (list[str]), ILEN を含む polars の DataFrame)の行番号です。インデックスから 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 ベース(VCF の規約)ですが、クエリ座標は 0 ベースの半開区間です。混同しないでください。

SparseVar2(.svar2)— クイックリファレンス

SparseVar2 は次世代のスパースバリアントストアです(VariantKey スタイルのインラインエンコーディングと、バリアントごとの密/スパースのコストモデル)。2 つの部分で構成されます。変換(以下の from_vcf)でストアを書き込み、範囲クエリ(さらに下の decode / region_counts / read_ranges)で読み戻します。座標はすべて、genoray のほかの箇所と同様に、0 ベースの半開区間 [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 型で、キーワード専用の frozen dataclass に定義されています。

フィールド意味最小値
concurrent_chroms並行して変換するコンティグ数1
reader_workers並行処理する各コンティグで、独立したインデックス付きシャードリーダーを実行する数1
overshardリーダーごとの作業単位数(コンティグの正確なレコード数が不明な場合にのみ参照されます)1
dense_capリーダーと実行器の間にある dense-chunk チャネルの深さ1
merge_threadsコンティグごとの var_key マージ末尾処理に使用するスレッド数1
sample_intervalモニターのサンプリング間隔(秒)。0 にすると無効になります0

None(すべてのフィールドのデフォルト値)は「プランナーに選択させる」という意味で、「無効」という意味ではありません。Tuning を指定しない場合にプランナーが選ぶ値と同じです。設定したフィールドは反映されるか、拒否されるかのどちらかで、暗黙に縮小されることはありません。たとえば、明示的に指定した reader_workers が max_mem の制限に収まらない場合、収まる値に引き下げられるのではなく、InsufficientMemory が発生します。これがこのオブジェクトの要点です。以前の環境変数は、読み取らないバックエンドに無視されることがありましたが、Tuning が無視されることはありません。tuning= には Tuning インスタンスまたは None を指定する必要があります。たとえば dict を渡すと、フィールド名がもっともらしいキーワード引数に見えてしまう可能性があるため、パラメーター名を示す TypeError が発生します。

すべてのバックエンドがすべての設定項目を利用できるわけではありません。利用できない項目を設定すると、黙って無視されるのではなく、該当するフィールド名を示す ValueError が発生します。

フィールドfrom_vcffrom_pgenfrom_vcf_listfrom_svar1
concurrent_chromsはいはいいいえはい
reader_workersはいいいえいいえいいえ
overshardはいいいえいいえいいえ
dense_capはいはいはいはい
merge_threadsはいはいはいはい
sample_intervalはいはいはいはい

理由:

  • reader_workers/overshard はシャード化された VCF でのみ利用できます。 from_pgen は reader_workers=1 に固定されており、コンティグ内をシャード化することはありません。pgenlib はデコード中に GIL を保持するため、そこでコンティグ内をシャード化しても純粋なオーバーヘッドになります。from_vcf_list(コンティグごとに入力ファイルごとのファイル記述子をすでに 1 つ使用します)も、from_svar1(VCF/htslib リーダー自体を使用しません)も、コンティグ内をシャード化しません。
  • concurrent_chroms は from_vcf_list では利用できません。このパイプラインは設計上、コンティグを順次処理するためです。

設定したフィールドは、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= — 変換対象を、インデックス化された 1 つ以上の VCF 取得区間に 制限します。領域文字列には、既存の genoray 規約を使用します("chrom:start-end" は 1 始まりの両端を含む形式で、0 始まりの半開区間に変換されます。タプル/BED/フレーム形式の入力は、すでに 0 始まりの半開区間です)。 領域が重複している場合、merge_overlapping=True でない限りエラーになります。regions_overlap は、bcftools の --regions-overlap に対応する 3 つのモードから選択します。"pos" (デフォルト。POS が [start,end) 内)、"record"(POS が [start,end+1) 内。領域の最終塩基にある indel も保持されます)、または "variant"(アンカートリミング後のバリアント範囲が領域と重複)です。"variant" モードでは、複数アレルのレコードは、そのアレルのいずれかが実際に領域と重複していれば全体が保持されます。個々の非重複アレルは削除されません。現在、"variant" ではコンティグごとに 1 つの領域までしか指定できません。コンティグごとに複数の領域を指定するとエラーになります。"pos"/"record" を使用するか、別々に変換してください。

  • samples= — VCF サンプルを名前で選択し、並べ替えます。呼び出し元の順序を維持し、最初に出現したものだけを残して重複を除去し、未知の名前には ValueError を送出します。available_samples とデコードされるすべての列は、各サンプルの元の VCF ヘッダー内の位置にかかわらず、呼び出し元が指定した順序と完全に一致します。

  • reference(REF の検証と indel の左寄せに使用する 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 の場合、対象外の(<DEL>/<INS>/… のようなシンボリック ALT およびブレークエンド)ALT はエラーにせず削除されます。厳格なデフォルト設定では、最初に見つかった ALT でエラーになります。この層では 2 つのクラスを区別できません。「シンボリックのみ」と「ブレークエンドのみ」を切り替える個別の設定はありません。

  • 削除された対象外 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 の dosage 形式に類するデータは、代わりに format_fields=(例:DS FORMAT フィールド)で渡します。haploid= による OR 集約も、max_mem に基づくチャンク分割もありません(代わりに chunk_size を使用)。現時点では、この 2 つは SparseVar(SVAR 1.0)でのみ利用できます。

  • threads=None — スレッドの合計予算です(None の場合は自動検出)。コンティグ間の並行処理と、プランナーを介したコンティグごとのリーダー数を制御します。

  • tuning.reader_workers(tuning=Tuning(reader_workers=N) で設定。上記の「チューニング」を参照)— 並行処理中のコンティグごとの独立したインデックス付きシャードリーダーの数です。コンティグ内の読み込み並列度を設定する項目です。None のままにすると、コア予算から算出されます。使用可能なコアの 4 分の 1 はマージ処理の最後に確保され、コンティグ間の並行度は深さ(コンティグあたり約 8 リーダー)を優先して選択され、残りがリーダーに割り当てられます。明示的に指定する場合は None または >= 1 の整数である必要があります。1 未満の値を指定すると、変換開始前の Tuning 構築時に ValueError が送出されます。max_mem には収まるものの、それ以外の理由で適用できない値を指定した場合は、暗黙に減らされるのではなく InsufficientMemory が送出されます。有効な値であれば、出力は常にバイト単位で同一です。(N 個の単一サンプル VCF をマージする from_vcf_list 経路では、コンティグ内でのシャーディングは行わず、reader_workers/overshard は拒否されます。上記の適用表を参照してください。)スケーリングの数値については、docs/source/svar.md の「並列変換」を参照してください。

  • signatures=False — True の場合、書き込み中にすべての SNP/indel を SBS96/ID83 変異タイプコードに分類し、コンティグごとの mutcat サイドカーを保存します(書き込み時の密なデータ/var_key コストモデルに組み込まれます)。参照が必要です(reference=)。no_reference=True と併用すると ValueError が送出されます。SVAR2 の mutcat サイドカーを読み取るための公開 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 は int 幅を指定できません)と観測範囲の両方に照らして検証されます。オーバーフロー、または f16 の約 65504 の範囲を超える値がある場合は ValueError が発生します。f16 は唯一の損失を伴う選択肢であり、明示的に指定する必要があります。
    • default — VCF で欠損しているエントリーに書き込む値です。指定しない場合は、選択された幅の端に予約値が使われます(符号付き幅では INT*_MIN、符号なし幅では u*::MAX — 観測範囲が非負の場合、自動幅縮小では符号なしが優先されます — 浮動小数点幅では NaN)。Flag フィールドが欠損することはありません(存在しない場合 ⇒ false/0)。
    • FORMAT は独立して情報を保持するのではなく、遺伝子型に対応付けられます:FORMAT 値は遺伝子型にコールがある場合にのみ保存されます。var_key ルーティングされた変異ではキャリアのコールごとに 1 値、dense ルーティングされた変異ではサンプルごとの完全な密列(非キャリアのスロットには default/予約値を入力)として保存されます。非キャリアの FORMAT 値(たとえば、ref/ref 遺伝子型に対する補完されたドサージ)は、このバージョンでは設計上破棄されます。独立して情報を保持する FORMAT ストリームは、今後の仕様に先送りされています。
    • 読み込みパス: 以下の「INFO/FORMAT フィールドの読み込み(SVAR2)」を参照してください。SparseVar2(path, fields=…)/.with_fields(…)/.available_fields を指定すると、decode() でこれらの値をデコードして取得できます。
  • 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 1 つ分をまかなえない予算は、コホートが非常に小さい場合でも ValueError で拒否されます。この下限はバックエンド固有です。VCF の則における LP の係数はベースラインが約 457 MB、並行 contig ごとに約 111 MB ですが、contig ごとの括弧内にある kappa 項がこれら 2 つの値を完全に上回るため、非常に小さなコホートでも実際の下限は、予算が厳しいときにプランナーが実際に選ぶ 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 の条件で評価すると、1つの同時処理コンティグに必要な最小 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でチャンクMBあたりの w の計算量が kappa * (2w - 1) から w * (kappa + 2) + 8 に変更され、from_vcf では reader_workers が固定のデフォルト値ではなくなりました。None の場合は値が導出され、plan_sharded はコンティグを諦める前に w_max から 1 まで w を順に調べるため、予算が厳しい場合は拒否されるのではなく、より小さな w が選ばれます。現在のデフォルトで実際に適用される下限は cc=1, w=1 の値であり、以下の w=3 の値よりもかなり低くなっています。「旧計算式」列は、直上で引用した 2026-08-11 の改訂値です。「現在の計算式」列は、現在 reader_workers=3 を明示した場合に w * (kappa + 2) + 8 の計算量で実際に必要となる値です。この2つは同じ値ではないため、同一のものとして読まないでください。

    コホート旧計算式、w=3(2026-08-11 の改訂)現在の計算式、w=3現在の実際の下限(w=1)
    S=4,0001,380 MB1,421 MB1,016 MB
    S=128,00026,400 MB27,833 MB14,864 MB
    S=500,000101,480 MB107,069 MB56,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 を小さくしてください」という2つの対処法が示されます。ただし、差は旧計算式の大きな開きではなく、わずか(52.4 GB 対 56.4 GB)になりました。128 GB のホスト(104,858 MB)では、もはや惜しくも届かないという状況ではありません。cc=1, w=2 で計画され、必要量は 81,739 MB、つまり約23 GBの余裕があります(w=3 では 107,069 MB が必要となり、予算をわずかに超えるため、w=2 で止まります)。

  • progress=False/log_level="info"/log_filter=None — 書き込み時の進捗表示とログ出力の設定です。from_vcf/from_pgen/from_vcf_list/from_svar1/write_view で共通です。progress=True にすると、リアルタイムで進捗が表示されます。ターミナルまたは Jupyter では rich のバー(処理中のコンティグごとに1行)が表示され、それ以外では、コンティグごとに約5秒に1回の頻度に制限された簡潔なハートビート行("chr1 42% (12,345/29,000) ...")が表示されます。progress の設定にかかわらず、コンティグの処理が完了するたびに、1行の "[svar2] chrom done: N kept, M excluded (Ts)" という要約が出力されます。ただし、log_level="off" の場合は出力されません。 log_level は、構造化された書き込み時のログ行に対する最低重大度です。 大文字と小文字を区別せず、次の値を受け付けます。"off"、"critical"、"error"、"warning"、 "info"(デフォルト)、"debug"、logging モジュールの整数定数 (logging.DEBUG など。中間の整数値は、logging 自身の判定と同様に、より重大度の高い名前付きレベルに 切り上げられます)、またはその整数を10進表記した文字列("10")。CLI から整数の 経路を利用する場合は、この文字列表記を使います。標準のランクは次のとおりです。 off=0, error=1, warning=2, info=3, debug=4。"critical" は受け付けますが、 "error" のエイリアスであり、独立した CRITICAL ランクはありません。 "warn" は拒否されます(ValueError)。これは意図的な仕様です。Python 3.13 では logging.warn() が削除されたため、genoray でもこの表記は受け付けません。 "warning" を使用してください。"off" は、contig ごとのサマリーや 進捗状況の表示を含むすべてを無効にします(処理は完全な no-op となり、オーバーヘッドはゼロです)。 "info" ではさらに、スレッド予算の選択、contig ごとの開始・終了、 および参照配列と照合した contig 名の解決(参照配列とソースでの表記が異なる場合)が含まれます。 "debug" ではさらに、レコードごとの詳細(REF/FASTA の不一致で除外されたレコードや、 左寄せされる各 indel)が出力されます。正確なマッピングについては parse_log_level を 参照してください。

    log_filter は、genoray 自身の stderr 診断出力に適用される、tracing 形式の EnvFilter ディレクティブ文字列です。log_level(上記の構造化された進捗・サマリーイベントのみを 制御します)とは独立しています。None(デフォルト)では stderr に何も出力されません。 解析できないディレクティブが指定された場合は、エラーを発生させずに stderr への出力を抑制します。 log_level では表現できない、モジュールやターゲットを指定したトレースに使用します。例: log_filter="genoray::monitor=trace"。

    genoray を設定する環境変数はありません。 log_level/log_filter (および上記のすべての Tuning フィールド)は、このパイプラインを設定していた従来の 8 個の環境変数をすべて置き換えます。完全な対応表と、ログトレース用変数の genoray::monitor=trace ディレクティブを log_filter="genoray::monitor=trace" に移行する方法については、 docs/source/svar.md の「環境変数からの移行」を参照してください。

    構造化ログ行では、メッセージの後に 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 のサブコンティグのシャーディングは無効です(コンティグごとにリーダーは 1 つ)。また、threads によって出力バイト列が変わることはありません。理由(実測値、chr21c:約 100 万バリアント × 3202 サンプル):単一リーダーでの変換はすでに高速(約 33 秒)で、ボトルネックは共有エグゼキューター/ライターとリファレンス I/O にあり、pgenlib のデコードにはありません。そのためシャーディングでこの下限を上回ることはできず、実測では遅くなりました(threads=24 で 44.9 秒、シリアル実行で 32.6 秒)。GIL を解放するビルドの pgenlib(>=0.94.x、prange によるデコードの並列化)に切り替えても効果はありません。変換時間は OMP_NUM_THREADS が 1~32 の範囲で約 33 秒のままであり、デコードの並列化による利点はありません。シャーディングの仕組み自体は存在し、100 万バリアントまでバイト単位で同一になることを検証済みです。今後ボトルネックがデコードに移った場合にのみ、再び有効化できます。
  • 二倍体のみ — ploidy= キーワード引数はありません(from_vcf のデフォルト ploidy=2 が暗黙的に固定されています)。
  • chunk_size=None — from_vcf の固定デフォルト 25_000 とは異なり、ここでの None はサンプル数からバリアント数の予算を算出します(パックされた密なチャンクのコストは chunk_size * n_samples * 2 / 8 バイトです)。そのため、200 サンプルでは問題ない固定値を使っても、50 万サンプルでメモリを使い果たすことはありません。明示的な int を渡すと上書きできます。算出値が 256 バリアントを下回る場合は警告します。詳細は、以下の from_vcf_list の chunk_size の項を参照してください。ここでは dosages は n_format_fields に数えられます。
  • max_mem: int | str | None = None — 同時実行プランナーのバイト単位の予算です。コホートのベースラインメモリに、同時に処理する各コンティグの処理中チャンクバッファーを加えた量が予算内に収まるよう、同時に変換するコンティグ数を選択します(既存のコア数による上限に加え、予算にかかわらず同時実行コンティグ数は最大 8 に制限されます)。文字列形式は、モジュールレベルの max_mem の規則(上記)と同じです("4g"、"512m"、"2GB"。parse_memory で解析されます)。また、(上記の)from_vcf の max_mem と同じくプロセス全体を対象とする意味です。どちらのパイプラインにも適合済みの同時実行プランナーがあり、同じ方法で予算を同時実行数に割り当てます。ただし、RAM の計算式の係数はそれぞれ別に適合されています(PGEN のチャンクは両方のハプロタイプを一度にデコードするため、バリアントあたりのコストが高くなります)。from_vcf_list の max_mem も同じくプロセス全体の予算を意味しますが、その処理経路には割り当て先となる同時実行プランナーがありません(コンティグは厳密に順番に処理されます)。そのため、代わりに予算からチャンクごとの chunk_size を算出します。詳細は以下の項を参照してください。**デフォルトの None は、上限なしではなく、検出された予算を意味します。検出された予算は、cgroup のメモリ制限の 80%(cgroup 外では /proc/meminfo の合計値)です。**これは max_mem 導入前のプランナーから意図的にデフォルトの動作を変更したものです。検出自体に失敗した場合(cgroup の制限がなく、読み取り可能な /proc/meminfo もない場合。macOS では常に該当)、genoray は警告を出し、エラーを発生させる代わりに、従来どおりコア数のみを上限とするプランニングにフォールバックします。予算を増減するには明示的な値を渡してください。上限なしに近いプランニングにするには、非常に大きな値を渡します。**実用上の下限:**プランナーの RAM 計算式には、コホートのベースライン項として約 2.7 GB の固定値があります(PGEN で適合した係数で、from_vcf の約 457 MB より大きくなっています)。そのため、ベースラインと同時実行するコンティグ 1 つ分のチャンクバッファーをまかなえない予算は、コホートがごく小さくても ValueError になります。このベースラインはコホートサイズに応じて増加します(~0.0158 MB/sample)。したがって、小規模コホートだけの問題ではありません。約 50 万サンプルでは、この項だけで約 10.6 GB と予測されるため、より小さいホストで検出された予算を使うと変換が拒否されます。その場合は、ホストに合った max_mem を明示的に渡してください。コンティグごとのチャンクのコストは チャンクバッファー 2 個分です。VCF の処理経路では 10 個分として計上されます。8 チャンクの並べ替え待ち上限は from_vcf の仕組みであり、コンティグごとにリーダーが 1 つの PGEN パイプラインでは適用されないためです。このため、以前は ValueError になっていた 2 つのしきい値の間の予算でも、プランニングできるようになりました。
  • regions=/merge_overlapping=/regions_overlap= — .pvar のバリアントインデックス範囲を 1 つ以上指定して変換対象を制限するもので、規則、意味、3 つの重なりモード("pos"/"record"/"variant")は from_vcf と同じです。from_vcf と同様に、"variant" モードでは、領域と実際に重なるアレルが 1 つでもあれば、多アレルレコード全体を保持します。
  • samples= — 名前で .psam のサンプルを選択し、並べ替えます(規則は from_vcf と同じです)。呼び出し元が指定した順序を保持し、最初に出現したものを残して重複を除去し、未知の名前には ValueError を発生させます。各サンプルの元の .psam 内の位置にかかわらず、available_samples とデコードされたすべての列は、呼び出し元が指定した順序と正確に一致します。
  • info_fields=/format_fields= はありません — PGEN には FORMAT がなく、.pvar の INFO 抽出も実装されていません。
  • dosages=Sequence[DosageField] — サンプルごとの dosage トラックを FORMAT フィールドとして保存します。各 DosageField(name="dosage", source="self"|Path, dtype="f16"|"f32"="f32", default=None) のうち、source="self" はハードコールの .pgen(上記の source)自体から dosage を読み込みます。Path は別の .pgen から読み込みます(たとえば、ハードコールと dosage の両方が同じファイルにあると pgenlib が dosage からハードコールを導出するため、別ファイルに分けてある VAF/CCF ファイル)。そのファイルはハードコールの .psam と同じサンプルを含み、ハードコールの .pvar のバリアントと 1 対 1 で対応している必要があります。ほかの FORMAT フィールドと同様に、遺伝子型に対応付けて保存されます。var_key による振り分けでは、非キャリアの dosage は破棄されます(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 が返すアレルコードの順序に従って解決されます。これは、フェーズ化されていない GT に関して from_vcf にも当てはまる注意点です。
  • progress=False/log_level="info"/log_filter=None — from_vcf(上記)と同じです。tuning= は concurrent_chroms、dense_cap、merge_threads、sample_interval を受け付けます。reader_workers/overshard を設定するとエラーになります(from_pgen はコンティグごとに単一のリーダーを固定します。「チューニング」を参照してください)。
複数の単一サンプル 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

サイトリストが異なる N 個の単一サンプル VCF/BCF から、ネイティブの k-way マージによって 1 つの SVAR2 ストアを構築します。bcftools merge も中間のマルチサンプル VCF も不要です。

  • regions=/merge_overlapping=/regions_overlap= — from_vcf と同じ規則、セマンティクス、3 種類の重複モード("pos"/"record"/"variant")が、マージ内のすべての入力ファイルに同じように適用されます。from_vcf と同様に、"variant" モードでは、いずれかのアレルが領域と実際に重複する場合、多アレルレコード全体が保持されます。

  • samples= パラメーターはありません — from_vcf/from_pgen/from_svar1 とは異なり、from_vcf_list には名前でサブセットできるコホートがありません。各入力ファイルはすでに単一サンプルであり、コホートは sources 経由で渡されたファイルの集合そのものです。

  • 各入力ファイルは単一サンプルである必要があります — サンプル列はちょうど1つです。0個または2個以上のファイルがある場合は ValueError が発生します。そのサンプルのVCFヘッダー名がストア内でのサンプル名になります。入力ファイル間でサンプル名が重複している場合も ValueError が発生します。

  • sources — モジュールレベルの _resolve_vcf_sources によって解決される、次の3種類のいずれかです。

    • Sequence[str | Path] — 指定された順序の明示的なファイル一覧です。
    • 単一のディレクトリ Path — そのディレクトリ直下(再帰なし)にある、bgzip圧縮されたすべてのVCF(*.vcf.gz/*.vcf.bgz)に続いて、すべての *.bcf を取得します。各グループ内は natsort 順です。
    • 単一のファイル Path — .vcf.gz/.vcf.bgz/.bcf は単一ファイルとして扱い、それ以外はマニフェストとして扱います(1行につき1パス。空行と # で始まるコメント行はスキップされ、相対パスはマニフェストの親ディレクトリを基準に解決されます)。
    • 解決後のファイルが0個の場合は ValueError が発生します。
  • サイトが存在しない場合は hom-ref 0 になります。 ファイルAには呼び出されているものの、ファイルBにはまったく存在しないサイトは、Bのサンプルでは 0(hom-ref)で埋められます。

  • ファイル内の ./. はマージ後に判別できません。 SVAR2のスパースレイアウトではALTを持つエントリーのみを格納するため、欠損ハプロタイプとhom-refハプロタイプはいずれもエントリー数0となり、decode や region_counts では区別できません。-1 の欠損値センチネルは genoray.VCF/genoray.PGEN の規約であり、SVAR2のデコードには含まれません。(マージ処理中は区別されますが、遺伝子型をスパースなキャリアビットグリッドにパックする際に破棄されます。これは from_vcf と同じ動作であり、2つの経路の整合性が保たれます。)

  • マージはアトム単位で結合します。正規化された (pos, ref, alt) アトムが完全に一致する場合に限り、ファイル間で共有される1つの行になります。位置だけが一致していても結合されません。

  • 各入力ファイルのレコードは、コンティグごとに位置順にソート済みである必要があります(単一入力に対する from_vcf と同じ前提です)。ソートされていないファイルがある場合、k-wayマージを黙って破損させるのではなく、該当ファイルと位置を示す ValueError が発生します。

  • すべての入力ファイルで同じコンティグ命名方式を使用する必要があります(すべて chr1 形式、またはすべて 1 形式とし、混在させないでください)。マージではファイルごとにコンティグ名の文字列を完全一致で照合するため、命名方式が混在するコホートでは、誤った表記のコンティグ上でコホートの半数のサンプルがすべてゼロとしてデコードされるストアが黙って生成されるのではなく、事前に ValueError が発生します(競合するファイル名と表記を示します)。

  • N個すべての入力ファイルを同時に開きます(コンティグごとにファイルあたり1つのファイル記述子を使用します)。Nが大きい場合(おおよそ N > (soft RLIMIT_NOFILE - 64) / 2、既定のソフトリミット1024では多くの場合N ≈ 480)、ulimit -n による対処方法を示す ValueError が発生します。これにより、上限付近の任意のファイルについてhtslibから分かりにくい「.tbiまたは.csiファイルはありますか?」というエラーが出るのを防ぎます。非常に大規模なコホート向けのバッチ/階層型マージへのフォールバックはありません(今後の課題です)。代わりにオープンファイル数の上限を引き上げてください。

  • no_reference=True がサポートされています。from_vcf/from_pgen と同様に、REFの検証と左揃えを省略し、各アトムのREFをレコード自身のREFバイト列から再構築します。reference と no_reference のどちらか一方のみを指定するチェック、および signatures と no_reference の併用不可条件は、いずれも from_vcf と同じです。

    • このエントリーポイント固有の注意点: マージは各コンティグ内で、正規化された各アトムの (pos, ref, alt) をキーとするk-way結合を行います。左揃えを省略すると、複数ファイルで共有されるサイトが1つの出力行に結合されるのは、すべての入力で同一の表現になっている場合に限られます。つまり、アンカーベースとパディングが同じである必要があります(たとえば、すべてのファイルが同じコーラーから作成された、またはすべて同じリファレンスを使って bcftools norm を実行済みである場合です)。異なる正規化で同じインデルをエンコードしている2つのファイルは、no_reference では結合されず、1つの共有行ではなく、出力ストア内の別々の変異になります。これはエラーが発生する失敗ケースではありません。from_vcf_list で no_reference を使用する前に、上流の正規化に一貫性があることを確認してください。
  • info_fields=/format_fields= — from_vcf と同じ宣言APIです(sources の最初のファイルのヘッダーに基づいて解決されます)。サイトごとにソース列がN個あるため、マージ時のセマンティクスは単一ファイルの変換とは異なります。

    • INFOフィールドは最初のキャリア優先でマージされます。 サイトが複数のファイルに存在する場合、保存されるINFO値は、そのアトムを持つファイルのうち番号が最も小さいもの(sources の順序が最も早いもの)から取得されます。最後のファイルからの値でも、キャリアの値を集約した値(例:最大値/合計値)でもありません。
    • FORMATフィールドはサンプルごとに保持されます。 from_vcf とまったく同じく、各サンプルは自身のファイルの値を取得します。アトムをまったく持たないサンプルでは、そのフィールドの既定値(予約済みセンチネル値/NaN、または明示的な default=)が使われます。
  • chunk_size=None — from_vcf の固定値 25_000 という既定値とは異なり、ここではコホートサイズに基づく予算ベースのチャンクサイズ(_auto_chunk_size)が算出されます。これにより、パック済みの密なチャンク1つが、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 を示します)。例外は発生させず、その小さい値を維持します。たとえば、既定の約256 MiBの予算では、ploidy=2, n_format_fields=7 の場合に n_samples≈37,118 を超えると警告が発生します。解消するには max_mem を引き上げるか、要求する format_fields を減らしてください。このヘルパーを介して chunk_size を導出する3つの変換関数(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 はコンティグごとにN個のファイルをすべて一度に開くため)。このオーバーヘッドは入力ファイル数に応じて増加し、max_mem では制限されません。「プロセス全体」を、RSS合計の厳密な保証と解釈しないでください。 破壊的変更: 以前のバージョンでは、この max_mem は処理中の密なチャンク1つのバイト数上限として扱われていました。現在はプロセス全体の予算になっています。以前は「密なチャンク1つで最大512 MiBを使用する」という意味だった max_mem="512MiB" のような指定は、現在では「密なチャンクの項で最大512 MiBを使用する」という意味になり、導出される chunk_size はかなり小さくなります。引き継いだ値は調整し直してください。

    from_vcf_list には、シャーディングされた from_vcf リーダーのような、適合済みの並行処理プランナーがありません。コンティグは厳密に順次変換されるため、並行処理を計画する代わりに、チャンクごとのバイト数目標を導出します。 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。これはリーダー/エグゼキューターのチャネル容量に加えて、単一のリーダースレッドと単一のエグゼキュータースレッドがそれぞれチャネル外に保持できるチャンク1つ分です)。そのため、重複した定数を使った場合のように、オーケストレーターの実際のアーキテクチャとの不整合がひそかに生じることはありません。ただし、どちらの場合も適合済みのメモリ法則ではありません(シャーディングパスとは異なり、このパイプラインにはそのような法則はありません)。_DENSE_CHUNK_TARGET_BYTES は上限のままなので、大きな予算を指定してもチャンクは現在のサイズを超えて大きくなりません。厳しい予算の場合にのみ、チャンクがそれより小さくなります。さらに、導出された目標値も最悪ケースの上限です。見積もりではすべてのバリアントが密に振り分けられると仮定するため、バリアントが疎に振り分けられるコホート(例:プライベートな体細胞変異のコール)では、使用量は大幅に少なくなります。chunk_size が明示的に渡された場合は無視されます。これには、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" の場合、不正なレコードは、そのレコードを含むファイル内でのみ除外され、マージ後のサイト全体が除外されることはありません。コンティグごとのログには、すべての入力ファイルで除外された合計数が報告されます。

  • progress=False/log_level="info"/log_filter=None —(上記の)from_vcf と同じです。tuning= で指定できるのは dense_cap、merge_threads、sample_interval のみです。concurrent_chroms、reader_workers、overshard を設定するとエラーになります(このパイプラインはコンティグを順次処理し、コンティグ内でシャーディングすることはありません。「Tuning」を参照してください)。

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 のみ対応 — ソースストアに多対立遺伝子のバリアントがある場合、ValueError が発生します(SVAR1 の geno==1 モデル)。まず SVAR1 ストアを二対立遺伝子として再作成してください。
  • regions=/merge_overlapping=/regions_overlap= — 規約、セマンティクス、3 つのオーバーラップモード("pos"/"record"/"variant")は from_vcf/from_pgen と同じです。"variant" モードでは、いずれかのアレルが実際に領域と重なる場合、レコード全体を保持します(ただし SVAR1 は二対立遺伝子のみのため、判定対象となる ALT は常に 1 つです)。from_pgen と異なり、SVAR1 には事前に検索範囲を絞り込むためのディスク上のカバー範囲インデックスがありません。選択した contig のローカルバリアントは引き続きすべてスキャンされます。出力を実際に制限するのはレコードごとのフィルターであるため、範囲を限定したスキャンではなく contig 全体のスキャンが必要です。
  • samples= — SVAR1 のサンプルを名前で選択し、並べ替えます(規約は from_vcf/from_pgen と同じです)。呼び出し側の順序を保持し、最初に現れた重複だけを残し、名前が不明な場合は ValueError が発生します。各サンプルの SVAR1 での元の位置にかかわらず、available_samples とデコードされた各列は、呼び出し側が指定した順序と完全に一致します。
  • fields=Sequence[str] | None — SVAR1 の FORMAT フィールド(例:dosages)から引き継ぐものを、SVAR1 内の名前で選択します。None(デフォルト)ではすべてを引き継ぎ(従来どおり、デフォルトで損失なく変換)、[] では何も引き継がず、名前のサブセットを指定するとそのフィールドだけを引き継ぎます。不明な名前を指定すると、利用可能なフィールドを列挙した ValueError が発生します。mutcat はこの方法では選択できず、常に破棄されます。SVAR1 のものを引き継ぐ代わりに参照配列からシグネチャを再計算するには、signatures=True を渡してください。
  • フィールドの同等性に関する注意: SVAR1 には非キャリアの FORMAT 値が保存されないため、フィールド出力が from_vcf とバイト単位で一致するのは、var_key(キャリアのみ)にルーティングされたバリアントに限られます。dense にルーティングされたバリアントでは、非キャリアのセルにはソース 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 を読み込まず、contig 内でシャーディングもしません。「Tuning」を参照してください)。
範囲クエリ

完成したストアを開き、コンティグごとにクエリします。構築時に meta.json を読み込み、コンティグごとにネイティブリーダーを開きます。.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 の不透明な文字列、バリアント軸のオフセットオブジェクトは共通です)。純粋な欠失の場合、ALT は空です(参照塩基は再出力されません)。seqpro が必要です。
  • region_counts(contig, regions) は、デコードを伴わないカウント(オフセットの差分+密なマスクの popcount)を返します。これは SparseVar.var_ranges を簡略化した代替機能です(SVAR2 には統一されたバリアントテーブルがないため、バリアントの インデックス は存在しなくなりました)。
  • クエリは コンティグごとに実行します。コンティグをまたぐバッチ処理は呼び出し側で行います。領域は (start, end) ペアのイテラブルです。
  • decode/region_counts/read_ranges の contig 引数では、別の命名方式も使用できます。chr 付き/なし(chr1 ↔ 1)およびミトコンドリアの別名 {M, MT, chrM, chrMT} は、ContigNormalizer によりストア内の表記に解決されます。解決できないコンティグを指定すると ValueError が発生します。

ユーザー向けの SVAR2 クエリ API は、上記の decode/region_counts/read_ranges です。read_ranges(contig, starts, ends, samples=None) は検索と取得を統合した処理です。starts/ends は並列する 1D 配列で(SparseVar.read_ranges と同様)、samples は名前によってサブセットを選択し、並べ替えます。生の 2 チャンネル BatchResult → numpy dict を返します。これはフィールドが固定された 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-dict ワイヤーコントラクトです。公開 API の一部ではなく、セマンティックバージョニングの対象外です。予告なく変更または削除される場合があるため、ユーザーコードから呼び出さないでください。

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 は正規キー(下記の available_fields を参照)の Sequence[str] です。 with_fields は同じストアを参照する新しい SparseVar2 を返します。元のインスタンスは変更されません。fields=None(コンストラクターのデフォルト)では何も選択されません。フィールドはオプトイン方式です。
  • available_fields -> dict[str, StoredField] — ストアの meta.json で宣言されたすべてのフィールドを、正規キーで返します。INFO と FORMAT を通じて一意な場合はフィールド名のみ、両方のカテゴリで使われている場合は bcftools 形式の INFO/DP / FORMAT/DP になります。StoredField(トップレベルの genoray ではなく genoray._svar2_fields で定義)は凍結されたデータクラスで、name、category("info"/"format")、dtype(np.dtype)、default(float | None)、key を持ちます。
  • decode(contig, regions) は、選択したフィールドごとに 1 つの Ragged を、返されるレコード Ragged に pos/ilen/allele とともに追加します。すべてが同じバリアント軸オフセットオブジェクトを共有し、形状は (R, S, P, None) です。 フィールドのデータには rag["KEY"](Ragged.__getitem__)でアクセスします。rag.fields["KEY"] ではありません。Ragged.fields はレコード上のフィールド名を表す list[str] にすぎません。
  • データ型は保存時のまま保持されます。 SVAR2 は書き込み時に整数フィールドの幅を自動的にロスレスで縮小するため、たとえば AC フィールドが int8 として返される場合があります。読み込み時に幅が拡張されることはありません。
  • 欠損値には、書き込み時に設定した場合はそのフィールドの default が、それ以外の場合は予約済みのセンチネル値(浮動小数点数では NaN、整数では iinfo.min/iinfo.max)が使われます。値はそのまま返され、変換されません。
  • FORMAT フィールドは遺伝子型に対応付けられます(上記の from_vcf の format_fields= を参照)。decode() が出力するのはキャリアレコードのみなので、書き込み時に「非キャリアの値は保存されない」という注意点は、この読み込みでは表面化しません。
変異シグネチャ(SBS96 / DBS78 / ID83)

.gvi に付加されたフィールドの代わりに contig ごとの Rust サイドカーを使用しますが、COSMIC のワークフローは SparseVar と同じです(上記の「変異カタログ」を参照)。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 付きと chr なし、およびミトコンドリアの別名 {M, MT, chrM, chrMT} — これらは ContigNormalizer を使ってストア内の表記に解決されます。要求された contig がすべて解決できない場合は ValueError が発生します。SparseVar.annotate_mutations と異なり、write_back= の切り替えオプションはありません — SVAR2 は常にサイドカーをディスクに保存します。任意の gtf= では GTF/GFF 遺伝子モデルのパスを指定できます。指定すると各 SNV は転写鎖クラス(feature == "gene" のフットプリントから判定)でも分類され、 strand.bin サイドカーに保存されます。これにより、以下の "SBS192"/"SBS384" カタログが利用可能になります。各サイドカーと meta.json はアトミックに書き込まれるため(同じディレクトリ内の一時ファイルを作成してからリネーム)、クラッシュしてもファイルが途中で切れた状態にはなりません。再実行は常に行われるため、サイドカーが見つからない場合や破損が疑われる場合の復旧手段にもなります。gtf= なしで再実行すると、以前に書き込まれた strand.bin は削除されます(新たに GTF なしで注釈を付けた際に、古い鎖クラスが残ることはありません)。
  • mutation_matrix(kind, *, count="allele"|"sample", contigs=None) -> pl.DataFrame — 固定された COSMIC コードブック順の MutationType 列と、サンプルごとの列を返します。kind ∈ {"SBS96", "DBS78", "ID83", "SBS192", "SBS384"}。 count="allele" は非参照アレルのコピーをすべてカウントします。count="sample" は、contig をまたいで論理和を取り、各サンプルで各カテゴリを最大 1 回カウントします。 contigs= では合計対象をサブセットに限定できます(別名も受け付け、重複は統合されます。ストアに存在しない名前を指定すると ValueError が発生します)。注釈が必要なのは指定した contig のみです。None(デフォルト)ではすべての contig が対象になります。対象範囲内の contig に注釈がない(ディスク上にサイドカーがない)場合は ValueError が発生します — 先に annotate_mutations または from_vcf(..., signatures=True) で注釈を付けてください。"SBS192"/"SBS384" にはさらに鎖の注釈(annotate_mutations(..., gtf=...))が必要で、対象範囲内の contig に鎖の注釈がない場合は ValueError が発生します。サイドカーが途中で切れている場合、または古くなっている場合(長さがストアのバリアント数と一致しない場合)は、annotate_mutations の再実行を促す OSError が発生します。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 に渡されます。SPA のアルゴリズムを使うには strategy=Spa() を指定するか、None のままにして max_delta/min_activity/criterion による前方選択の簡略指定を使用します。 fit_signatures は 1 回だけ呼び出されるため、strategy= と簡略指定用の 3 つのキーワード引数のいずれかを同時に渡すと、 fit_signatures と同様に ValueError が発生します。
  • v1 と同じ分類ルールを使用します(Rust 分類器を共有):DBS78 は、孤立した隣接する同一ハプロタイプの SNV ペアからのみ生じます — 3 個以上連続する隣接 SNV は個別の SBS96 エントリーのままとなり、2 bp を超えるネイティブ MNV は SBS96 に分解され、孤立した各ダブレットは(構成する SNV ごとではなく)1 回だけカウントされます。
  • 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 — 変異間距離のカットオフ(塩基単位)です。変異は、最も近い隣接変異までの距離が、そのサンプルのカットオフ <= の場合にクラスター化されます。すべてのサンプルに適用する単一の float、またはサンプル名からカットオフへのマッピングを指定します。マッピングには、すべての available_samples エントリを含め、他のエントリを含めてはいけません(それ以外の場合は ValueError)。カットオフは > 0 である必要があります。カットオフは呼び出し元が指定するもので、ストアにはありません。上流ではシミュレーションによってサンプルごとに導出し、imds.pickle に保存します。同じ値を指定すると、上流で分類されるすべてのイベントに同じラベルが付与されます(SigProfilerClusters v1.2.2)。
  • vaf_field — 呼び出しごとの VAF/CCF を格納する、任意の FORMAT フィールドです(正規キーまたは available_fields にあるフィールド名のみを指定)。指定すると、VAF の整合性が決定木に反映され、失敗したイベントは貪欲法で再分割されます。2 バイトまたは 4 バイトの float dtype(f16/f32)を持つ FORMAT フィールドに解決される必要があります。それ以外の場合は ValueError になります。欠損値は、上流の -1.5 センチネルとして扱われます。VAF を使用しない場合は、省略します(デフォルトは None)。
  • vaf_cut — 隣接する変異が同一イベントにとどまるための |delta VAF| の最大値です(上流のデフォルトは 0.1、CCF 設定では 0.25)。> 0 である必要があります。
  • contigs — 注釈をサブセットに限定します(chr 付きと chr なしの名前、およびミトコンドリアのエイリアスは、annotate_mutations と同様に解決されます)。ストアに存在しない名前を指定すると ValueError になります。対象外の contig はすべて 255 で埋められ、その領域でもフィールドをデコードできる状態に保ちます。None(デフォルト)を指定すると、すべての contig に注釈を付けます。再実行時は常に再処理されます。
  • meta.json に cluster_version(現在は 1)、cluster_contigs、cluster_cutoff、cluster_vaf_field、cluster_vaf_cut を記録し、cluster_class FORMAT エントリ(u8)を追加します。

ラベルのコードブック — SNP キャリアの呼び出しごとに 1 つのラベルを付与します(1/2 遺伝子型には、記録された代替アレルごとにラベルを付与します。indel はサブ分類されません)。

codelabelupstream
0nonclusteredクラスター化ファイルに含まれない
1doubletClassIA
2mbsClassIB
3omikliClassIC
4kataegisClassII
5otherClassIII
255not_annotated対象外 contig / 非 SNP 呼び出し

既存のフィールド機構を使って読み戻します — with_fields(["cluster_class"]) の後に decode(contig, ranges) を実行します。このフィールドは pos/ilen/allele と同じバリアント軸に沿って格納されるため、rag["cluster_class"] はストア内の呼び出しごとのフラットな順序で並び、rag["pos"] と 1:1 で対応します。ディスク上では、var_key_* ストリームは呼び出しごとに 1 つのラベルを格納します。dense_* ストリームは (dense_row, sample) ごとに 1 つのラベルを、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 付きと 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 ごとに 1 つの単一 contig ストアへ分割し、out_dir/{contig}.svar2 に出力します。出力パスを self.contigs の順序で返します。contig ごとに subset_contigs を 1 回呼び出すことで実装されています。
  • SparseVar2.concat(output, sources, *, mode="copy", overwrite=False) -> None(クラスメソッド)— contig セットが互いに重複しないストアを 1 つに結合します。sources はパス(または SparseVar2 インスタンス)のシーケンスです。すべてのソースで samples、ploidy、format_version、fields が一致している必要があります。いずれかが一致しない場合、または複数のソースに同じ contig 名がある場合は、ValueError が発生します。結合後の contig リストは natsorted で並べ替えられ、sources の指定順序には依存しません。
  • mode(3 つすべてのメソッドで使用)は共通の 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 内部のバグ(ワーカースレッドのパニック)です。ご報告ください。

CLI

genoray write には 3 つのサブコマンド 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

3 つの write サブコマンドはすべて、--regions/-r、--regions-file/-R、--samples/-s、--samples-file/-S、--merge-overlapping、--regions-overlap(pos/record/variant)、--reference と --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 は拒否されます。Python 3.13 で logging.warn() が削除されたため、これは慣例ではなく非推奨の表記です。パラメーターには固定の選択肢リストではなく、検証済みの 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」を参照してください)。さらに、write vcf、write pgen、write svar1 では --concurrent-chroms N も使用できます。vcf-list 形式では使用できません(from_vcf_list は設計上、コンティグを順番に処理します。「Tuning」を参照してください)。Python の tuning= キーワード引数と同様に、明示的に指定したフラグ値はそのまま適用されるか、拒否されます。暗黙のうちに小さく調整されることはありません。また、バックエンドで使用できないフラグを渡すと、無視されずにエラーになります。

write vcf では、さらに --reader-workers N と --overshard N も指定できます(単一ファイル入力のみ。ディレクトリ/マニフェストと併せていずれかを渡すとエラーになります。vcf-list 形式ではコンティグ内のシャーディングは行いません)。単一ファイル形式では、どちらもそのまま 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 形式ではエラーになります(各入力ファイルがすでにサンプルを1つだけ提供しているため、サブセット化するコホートがありません)。--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 です(PGEN は二倍体専用のため、--ploidy はありません)。--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(ここでは独立したフラグです。SVAR2 の write サブコマンドでは単一の --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)にはさらに、write と同じ6種類の名前を指定する --log-level LEVEL と、 write と同じセマンティクスの --log-filter DIRECTIVE があります。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 ごとに1つの region(SVAR2 には contig 長のメタデータがないため、SparseVar2.contigs を使用)で、範囲は [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 への切り替えでは、キャリアのコールごとに1つの値が保存され、キャリアではないサンプルの FORMAT 値を格納する場所がありません。そのため、FORMAT フィールドを引き継ぐビューでソースの dense バリアントを再ルーティングすると、その値が気付かれないまま失われます。"auto" は FORMAT を扱う場合は忠実性を優先し、それ以外の場合(ジェノタイプのみ/INFO のみのビューでは、失われる可能性のあるサンプルごとの格納領域がありません)はサイズ最適化のために再ルーティングします。

    --reference(サブセット上で mutcat をゼロから再計算)と -@/--threads(同時にスライスする contig 数を制限。省略時は自動検出)は、--reroute と --no-reroute のどちらでも実際に機能します。どちらの経路にも「受け付けるが無視される/使用されない」という但し書きはなくなりました。--progress、--log-level、--log-filter もすべてここで実際に機能します。contig ごとに1行を表示する簡易表示とログレベルのセマンティクスについては、以下の「write_view 進捗バー」を参照してください。write_view の基盤にある reroute= キーワード引数が受け付けるのは "auto"、True、False のみです。それ以外の値(例:reroute=1)を指定すると、reroute=False のスライサーに暗黙的にフォールスルーするのではなく、ValueError が発生します。

  • 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 の docstring を参照)と --overwrite を受け付けます。

フィルタリング

VCF:filter= に genoray.Filter(record=, expr=) の値オブジェクトを渡します。 record は、ジェノタイプのスキャン中に適用される Callable[[cyvcf2.Variant], bool] です。expr は .gvi インデックスに適用される、対応する polars の pl.Expr です。VCF では両方の要素が必要で、常に一緒に扱われるため、互いに食い違うことはありません。

構築後に VCF のフィルターを変更するには、Filter(または解除する場合は None)を vcf.filter セッターに代入します。メモリー内のインデックスは無効化されます。 ゲッターは現在適用されている 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

従来の 2 引数コンストラクター(filter= とは別に polars 式をキーワード引数で渡す形式)と、タプル値を返す vcf.filter ゲッター/セッターは、3.0.0 で削除されます。レコード述語と polars 式を別々に渡しているコードは、上記の単一の Filter(record=, expr=) オブジェクトを渡す形式に移行してください。

PGEN:.gvi インデックスの列を対象に、ブールマスクを返す polars の pl.Expr を渡します。組み込み式は genoray.exprs にあります(すべての一覧)。

  • is_snp(すべての ALT アレルの ILEN が 0 の場合は True。null の ILEN がある行は False)
  • is_indel(すべての ALT アレルの ILEN が 0 以外の場合は True。null の ILEN がある行は False)
  • is_biallelic
  • is_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 アレルごとに 1 つの値)

記号 SV における ILEN の意味。 精密な <DEL>/<INS>/<DUP> の場合、ILEN はインデックス構築時に INFO フィールドから計算されます。<DEL> では -|SVLEN|、<INS>/<DUP> では +|SVLEN| となり、SVLEN がない場合は |END - POS| にフォールバックします。VCF では INFO フィールドはヘッダーで宣言された列から読み込まれ(oxbow 経由)、PGEN では PVAR の INFO 文字列から解析されます。記号以外の ALT には、len(ALT) - len(REF) が使われます。

サイズを特定できない記号アレルの ILEN は null です。 次の場合、アレルのサイズを特定できません。INFO フラグ IMPRECISE が設定されている場合、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 を参照してください。2 つの 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 では、呼び出し元が対象とするコンティグをフィルタリングします(デフォルトでは FASTA 内のすべてのコンティグ)。
  • fetch(contig, start, end) — 0 ベースの半開区間 [start, end) です。コンティグの範囲外の位置は N で埋められます。NDArray[np.uint8] を返します。
  • contig_array(contig) — キャッシュされたコンティグ全体の配列を NDArray[np.uint8] で返します。fetch と同じ、メモリ上の 1 コンティグキャッシュを共有します。chr 付き、または chr なしの名前を受け付けます。
  • コンティグ名に依存しません。ContigNormalizer により、"chr1" と "1" のどちらも正しく解決されます。
  • 一度にメモリにキャッシュされるコンティグは 1 つです。コンティグごとに順次アクセスすると効率的です。

変異カタログ(SBS-96/DBS-78/ID-83)

write_view の進捗バー

SparseVar.write_view(..., progress=False) は、オプトインの progress キーワードを受け付けます。 True にすると、ビューの書き込み中にフェーズ単位の rich 進捗バーが表示されます(主要な各ステップで 1 ティック:カウント、遺伝子型、各保持フィールド、インデックスの構築、および reference= が指定されている場合の変異アノテーション)。デフォルトは False で、バーもオーバーヘッドも発生しないため、ライブラリやパイプラインの呼び出し元への影響はありません。genoray view svar1 CLI では、同じオプションを --progress として利用できます(こちらもデフォルトではオフです)。

genoray view svar1 in.svar out.svar -r chr1:1-1000 -s A,B --progress

バーは見た目だけのものです。有効にしても無効にしても、出力バイト列、スキーマ、dtype は同一です。

SparseVar2.write_view(..., progress=False, log_level="info", log_filter=None) は、from_* ライターと同様に、書き込みの進捗をリアルタイムで表示します(上記の「変換」を参照)。ただし、from_* ライターとは異なり、write_view にはサンプリング対象となるレコードストリームがないため、進捗は大まかで、コンティグごとに 1 行表示され、コンティグ内でのバーの移動はありません。ターミナルまたは Jupyter では、progress=True にすると、処理中および完了したコンティグのリストがリアルタイムで更新されます。それ以外の環境では、各コンティグの処理が完了するたびに、簡潔な「chrom done」行が表示されます。progress の値にかかわらず、コンティグの処理完了時には、log_level="off" でない限り、1 行の "[svar2] chrom done: N kept, 0 excluded (Ts)" サマリーが表示されます(スライスで変異が除外されることはないため、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 スタイルの変異カタログをサポートしています。ワークフローは次のとおりです。

  1. svar.annotate_mutations(reference) を 1 回呼び出して、すべての変異を分類し、.svar ディレクトリに mutcat.npy を書き込みます。
  2. 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 — 指定した場合(コンティグ名のリスト)、そのコンティグ上の変異のみが分類されます。 それ以外のすべてのコンティグ上のエントリーは NOT_ANNOTATED(センチネル -4)として記録され、 それらのコンティグがリファレンスから取得されることはありません。名前の照合には ContigNormalizer が使用されるため、chr1 と 1 のどちらも使用できます。 .svar インデックスに存在しない指定コンティグは警告を出してスキップされます。インデックスには存在するものの リファレンスには存在しないコンティグをリストに含めるとエラーになります(除外するにはリストから 省いてください)。None(デフォルト)ではすべてのコンティグを分類します。write_back=True の場合、 正規化された対象範囲は metadata.json に mutcat_contigs として記録されます(None = すべて)。
  • write_back=True — mutcat.npy を保存し、metadata.json を更新するため、以降の SparseVar(dir, fields=["mutcat"]) のオープン時にこのフィールドが読み込まれます。 注: write_view は出力に mutcat を位置ベースでコピーすることはありません (以下を参照)。サブセット上で再計算するには write_view に reference= を渡すか、 出力ビューで 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"] が設定されます。

分類内容:

変異タイプチャネル
孤立した SNVSBS-96(トリヌクレオチドコンテキスト)
同じハプロタイプ上で隣接する SNV ペアDBS-78(5' エントリー = DBS コード、3' エントリー = DBS_PARTNER センチネル)
3 個以上の隣接する SNV の連続SBS(それぞれ独立したまま。DBS への統合なし)
VCF 内のネイティブな 2 bp MNVDBS-78
2 bp を超える MNV、シンボリック変異、非 ACGTUNCLASSIFIED
挿入/欠失ID-83(サイズ、リピートコンテキストによる分類)
contigs= の範囲外のコンティグ上にある変異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" — 各サンプルにつき、各カテゴリーを最大 1 回数えます(有無)。
  • MutationType という文字列カラムに続き、サンプルごとに 1 つの Int64 カラムを持つ Polars の DataFrame を返します。行は固定された 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
-1DBS_PARTNER — 隣接する SNV ペアの 3' 側。カウント対象外です
-2UNCLASSIFIED — シンボリック変異/複合変異/2 bp を超える MNV/非 ACGT
-3MISSING — 予約済みセンチネル(コード空間では定義されていますが、annotate_mutations v1 では出力されません。SparseVar は ALT を持つエントリーのみを保存するため、コールなしのスロットは不定長フィールドに現れません)
-4NOT_ANNOTATED — contigs= による注釈範囲外のコンティグ上にあるエントリー。カウント対象外です

注釈済みのファイルを読み込むには:

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 の統合は、同じハプロタイプ上の孤立した隣接ペアにのみ適用されます。隣接する SNV が 3 個以上連続する場合は、個別の 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 列(正規のコードブック順)と、シグネチャごとに 1 列を返します。ID83 では genome は無視されます。
  • fit_signatures(catalogue, reference, *, strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame — catalogue を reference に対して再フィッティングします。MutationType をキーに結合して行をそろえます(カタログにリファレンスにないタイプがある場合は ValueError を発生させます)。サンプルごとに 1 行を返します。列は Sample、シグネチャごとに 1 つの Float 列(活性。選択されていない場合は 0.0)、および cosine_similarity です。 n_jobs=1(デフォルト)は逐次実行、n_jobs=-1 はすべてのコアを使用します。結果は n_jobs/backend にかかわらず同一です。 strategy= はアルゴリズムを選択します。max_delta/min_activity/criterion とは併用できません(両方が指定されると、デフォルト値であっても ValueError が発生します)。これら3つは strategy=Forward(...) の恒久的な簡略指定であり、非推奨ではありません。
    • Forward(max_delta=0.01, min_activity=0.005, criterion="cosine")(strategy を省略した場合のデフォルト) — genoray 独自のスパースな貪欲法による前向き選択です。NNLS + 条件に基づく追加 + 最小活性値による枝刈りを行います。 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)をゼロから再実装したものです。SigProfiler への依存はなく、numpy/scipy/polars のみを使用します。すべての参照シグネチャに対して飽和させ、相対 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" を指定すると、Forward と同様に、スケーリングされていない NNLS の重みを返します。上流の SPA とは異なり、決定論的で乱数生成器には触れません。測定したすべての変異負荷で Forward より多くの真のシグネチャを復元しますが、Forward の "cosine" 条件と同様にスケール不変であるため、一貫した推定量ではありません。 サンプルごとの処理は Forward よりおよそ1桁遅くなります(実際の 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 自体と同様にそのまま渡されます。SigProfilerAssignment のアルゴリズムを使うには strategy=Spa() を指定します。strategy を None(デフォルト)のままにして、max_delta/min_activity/criterion で Forward の簡略指定を設定することもできます。strategy= と max_delta/min_activity/criterion のいずれかを併せて渡すと、fit_signatures と同様に ValueError が発生します(省略されたキーワード引数にはデフォルト値が使用されます)。 SparseVar2.assign_signatures のシグネチャは同一ですが、さらに contigs= を受け付け、ストア内のコンティグの一部にカタログを絞り込めます。
リフィット戦略

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 の 4 つの SBS グループ(genoray.SPA_CONNECTED_GROUPS:SBS2/13、SBS7a-d、SBS10a/b、SBS17a/b)が使用されます。独自のグループをシーケンスのシーケンスとして渡すか、False を指定して無効にできます。

Forward(criterion="cosine") も Spa も、変異負荷を考慮しません。一貫性が重要な全ゲノムカタログには、Forward(criterion="bic") を使用してください。

出力スキーマは戦略によって変わりません。Sample、参照シグネチャごとに 1 つの Float 列、末尾に cosine_similarity が含まれます。

対象外(v1):デノボ抽出、機会正規化、ブートストラップ信頼区間、プロット。

よくある間違い

間違い修正方法
genoray.Genos8genoray.VCF.Genos8(クラス属性)
vcf.read(..., phasing=True)VCF() コンストラクターで phasing=True を設定します
dosage_field= を指定せずに VCF から用量を読み取るコンストラクターで 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 レコード)を返します
mutcat フィールドなしで mutation_matrix を呼び出す先に 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 では、reference= を指定すると reroute=True と reroute=False のどちらでも mutcat が一から再計算されます。
FORMAT フィールドを fields= に渡して SparseVar2.write_view(..., reroute=True) を呼び出し、フィールドが除外または拒否されると想定する現在は reroute=True と reroute=False のどちらでも fields が引き継がれます(以前は reroute=True で ValueError が発生していました)。代わりに、デフォルトの reroute="auto" に注意してください。FORMAT フィールドが引き継がれる場合は常に reroute=False に解決されます。これは、dense から var_key への切り替えでは、非キャリアサンプルの FORMAT 値を格納する場所がないためです。
write_view の output にソースデータセットのディレクトリを渡す(overwrite=True を指定した場合も同様)ValueError が発生します。インプレースで書き込むと、ビューが書き込まれる前にソースが削除されるためです。別の output パスを指定してください。
annotate_mutations に FASTA パスを直接渡すサポートされています。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 を参照してください。

コメント