genoray-api について
genorayは、VCF/BCF、PLINK 2 PGEN、疎なバリアントストアを扱う、NumPyを中心としたゲノムバリアント検索ツールキットです。領域検索、変換、フィルタリング、変異シグネチャ解析を行う開発者や研究者向けです。
主な機能
- 01VCF/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 は、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)と、それに対応する.gvipolars 式(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/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) 署名、キーワード引数、または形状が不明な場合は、第一原理から推測するのではなく、ソース内の 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_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に固定されており、コンティグ内をシャード化することはありません。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=(例:DSFORMAT フィールド)で渡します。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,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を小さくしてください」という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=)が使われます。
- INFOフィールドは最初のキャリア優先でマージされます。 サイトが複数のファイルに存在する場合、保存されるINFO値は、そのアトムを持つファイルのうち番号が最も小さいもの(
-
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_classFORMAT エントリ(u8)を追加します。
ラベルのコードブック — SNP キャリアの呼び出しごとに 1 つのラベルを付与します(1/2 遺伝子型には、記録された代替アレルごとにラベルを付与します。indel はサブ分類されません)。
| code | label | upstream |
|---|---|---|
| 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"] と 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_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 アレルごとに 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 スタイルの変異カタログをサポートしています。ワークフローは次のとおりです。
svar.annotate_mutations(reference)を 1 回呼び出して、すべての変異を分類し、.svarディレクトリにmutcat.npyを書き込みます。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"]が設定されます。
分類内容:
| 変異タイプ | チャネル |
|---|---|
| 孤立した 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= の範囲外のコンティグ上にある変異 | 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 |
-1 | DBS_PARTNER — 隣接する SNV ペアの 3' 側。カウント対象外です |
-2 | UNCLASSIFIED — シンボリック変異/複合変異/2 bp を超える MNV/非 ACGT |
-3 | MISSING — 予約済みセンチネル(コード空間では定義されていますが、annotate_mutations v1 では出力されません。SparseVar は ALT を持つエントリーのみを保存するため、コールなしのスロットは不定長フィールドに現れません) |
-4 | NOT_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.Genos8 | genoray.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 を参照してください。