mcpField
T
tooluniverse-variant-analysis

tooluniverse-variant-analysis

github.com/mims-harvardデータ・データベース1.8k372 回インストール2 日前に更新
リポジトリ

tooluniverse-variant-analysis について

バリアント解析とアノテーションを行うToolUniverseスキルです。再解析の前に既存の結果や指定スクリプトを確認する手順を示し、VCF集計、コーディングバリアントのフィルタリング、証拠検索を扱います。

主な機能

  1. 01事前計算済みの結果を先に確認
  2. 02集計とフィルタリング用スクリプト
  3. 03VCFの多アレルレコードに対応
  4. 04臨床情報と集団データの検索

活用シーン

  1. 01VCFやHaplotypeCallerのSNP・indel集計
  2. 02複数サンプルのコーディングバリアント平均算出
  3. 03バリアント注釈とClinVar分類割合の評価

MCP 対応

VCFの統計と正規化に加え、バリアント、集団頻度、遺伝子量感受性を調べるToolUniverseツールを案内します。

インストール

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

$npx skills add https://github.com/mims-harvard/tooluniverse --skill tooluniverse-variant-analysis

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

同梱ファイル

SKILL.md と一緒にインストールされる参考資料やスクリプト。

.env.template
python_implementation.py
QUICK_START.md
test_variant_analysis.py
references/
annotation_guide.md
mutation_classification_guide.md
sv_cnv_analysis.md
variant_notation_conversion.md
vcf_filtering.md
scripts/
annotate_variants.py
coding_variant_filter.py
filter_variants.py
gatk_haplotypecaller_pipeline.py
parse_vcf.py
variant_fraction.py
test_data/
cohort2.vcf
sample.vcf

詳細

バリアント解析とアノテーション

ルールゼロ — まず事前計算済みの結果を確認する

以下の指示に従う前に、データフォルダーを調べて、次のファイルがないか確認します。

  • *_executed.ipynb → tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' で読み込み、セルの出力を信頼できる回答として引用します
  • 事前計算済みの結果ファイル(*results*、*deseq*、*enrich*、*stats*、*_simplified.csv などの名前の CSV/TSV)→ 直接読み込み、求められている値を報告します
  • 正規の解析スクリプト(analysis.R、run_*.py、find_*.R、*.Rmd)→ そのまま実行し、出力を確認します

上記のいずれも存在しない場合に限り、このスキルの以下の再解析手順に従います。生データから再実行すると、公開済みの回答とは異なる数値になり、処理時間も大幅に長くなります(多くの場合、ターン数が5~10倍になります)。


主要スクリプト — まずこれらを使用する

同梱のこれらのスクリプトには、質問固有の注意点(分母の選択、倍数性の既定値、多アレルの分割、複数行の Excel ヘッダー、非コード領域の許可リスト)が反映されています。解析しやすいラベル付きの KEY=VALUE 行を出力するため、場当たり的な pandas/awk の出力を使うよりも優先してください。

スクリプト使用する場合
gatk_haplotypecaller_pipeline.py「BAM から HaplotypeCaller がコールした SNP/indel の数」を尋ねられた場合。BWA インデックス作成 → アラインメント → ソート → インデックス作成 → HaplotypeCaller の処理に対応し、既存の BAM から開始(アラインメントをスキップ)することも、既存の VCF のみをカウントすることもできます。既定値の --ploidy 2 は GATK 自体の既定値に一致します(「HaplotypeCaller がコールした」ものの GT の多くはこの設定で生成されています)。原核生物を明示的に半数体としてコールする場合は --ploidy 1 を指定します。多アレルの分割と bcftools 形式の SNP/indel 判定が組み込まれています。
coding_variant_filter.py「イントロン、遺伝子間領域、UTR バリアントを除外した後のサンプルあたりの CHIP/コーディングバリアントの平均数」を尋ねられた場合。2段階の正規フィルターを適用します。(1)Zygosity == Reference の行がある場合は除外します(この行はカウントを約10倍に増やします)。(2)イントロン、遺伝子間領域、UTR、上流、下流の SO 用語を除外します。2行の VarSeq Excel ヘッダーと、サンプルごとのフォルダーまたは統合 CSV に対応しています。
variant_fraction.py「VAF < X のバリアントのうち、Y とアノテーションされた割合」を尋ねられた場合。分母はコーディング領域のサブセットのみです(synonymous/missense/splice_region/stop_gained/lost/start_lost/frameshift/inframe indel)。すべてのレコードではありません。

既存の VCF/BCF をスクリプトを書かずにカウントする場合(特に MCP サーバー内で、シェルから bcftools を連結して実行するのが難しい場合)、VCFStatsTool ツールでは、以下の正規手順を1回の構造化された呼び出しで実行できます。VCF_summary_stats(レコード/SNP/indel/MNP/ts-tv/サンプルごとの統計)、VCF_count_variants(PASS/QUAL/領域/式フィルター適用後のカウント)、VCF_normalize(多アレルの分割と任意の左寄せを行い、処理前後のカウントを報告)。これらのツールは内部で bcftools を実行するため、数値は以下に記載されたシェルコマンドと一致します。CLI 出力を解析する代わりに、決定的な JSON 結果が必要な場合に使用してください。

ワークスペースの分離(重要)

gatk_haplotypecaller_pipeline.py と coding_variant_filter.py は、 the input data folder ディレクトリ内への書き込みを拒否します。 これらのディレクトリは慣例上、読み取り専用です。HaplotypeCaller の中間 BAM/VCF および スクリプト内部の一時ファイルには、必ず --workdir /tmp/<run_dir>(またはデータフォルダー外の書き込み可能なパス)を指定してください。

入力データフォルダー内の参照 FASTA、FASTQ、既存の BAM/VCF ファイルは読み取り専用です。 データフォルダー内の BAM に .bai インデックスを追加する必要がある場合、スクリプトはその BAM を workdir にコピーします。

具体的な実行例

サンプルのソート済み BAM で HaplotypeCaller を再実行します(「HaplotypeCaller は BAM 内で SNP/indel をいくつ特定しましたか」という質問に対する標準的な方法です。デフォルト以外のフラグや、質問と一致しないフィルタリング後の結果で生成されている可能性がある、同梱の *_raw_variants.vcf のカウントより優先してください)。

python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --reference <data-folder>/REF.fna \
  --bam <data-folder>/SAMPLE_sorted.bam \
  --workdir /tmp/hc_run --sample-name SAMPLE

FASTQ からの完全なパイプライン(BWA + ソート + HaplotypeCaller、約 5~10 分):

python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --reference <data-folder>/REF.fna \
  --fastq-r1 <data-folder>/SAMPLE_1.fastq.gz --fastq-r2 <data-folder>/SAMPLE_2.fastq.gz \
  --workdir /tmp/hc_run --sample-name SAMPLE

既存の VCF のみをカウントします(「variants.vcf にはレコードがいくつありますか」のように、質問でそのファイルが明示されている場合に限ります。「HaplotypeCaller は SNP をいくつ特定しましたか」には使用しないでください。同梱ファイルの倍数性/フィルタリングが質問と一致しない可能性があります)。

python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --vcf <data-folder>/SAMPLE_variants.vcf

イントロン/遺伝子間領域/UTR のフィルタリング後の、サンプルあたりの CHIP バリアント数の平均(サンプルごとに 2 行のヘッダーがある VarSeq Excel ファイルのフォルダー):

python skills/tooluniverse-variant-analysis/scripts/coding_variant_filter.py \
  --dir <data-folder>/CHIP_DP10_GQ20_PASS --pattern '*.xlsx' --header-rows 2

単一の結合済み CSV に同じフィルタを適用します:

python skills/tooluniverse-variant-analysis/scripts/coding_variant_filter.py \
  --file all_samples.csv --sample-col sample --header-rows 1
grep する出力キー

gatk_haplotypecaller_pipeline.py:SNP_COUNT_ALLELES、INDEL_COUNT_ALLELES、 TOTAL_RECORDS、PLOIDY、VCF_PATH。

coding_variant_filter.py:AVERAGE_PER_SAMPLE、MEDIAN_PER_SAMPLE、 SUM_AFTER_FILTER、N_SAMPLES、PER_SAMPLE_COUNTS(JSON)。

質問が「サンプルあたりの平均」の場合は、AVERAGE_PER_SAMPLE を報告してください(SUM_AFTER_FILTER ではありません)。コホート全体の合計は、N_SAMPLES × サンプルあたりの平均です。平均を尋ねられているのに合計を報告すると、一般的な CHIP コホートでは約 80 倍の誤差になります。スクリプトは常に両方を出力するので、質問の表現に合った方を選んでください。

倍数性:生物種ではなく、質問のパイプラインに合わせる

GATK HaplotypeCaller のデフォルトは --sample-ploidy 2 です。公表されている「HaplotypeCaller は SNP/indel をいくつ特定しましたか」という回答の多くは、原核生物であっても、このデフォルト設定で HC を実行して得られています。そのため、スクリプトもデフォルトでは倍数性を 2 に設定しています。倍数体 1 のコールを明示的に求められた場合(例:「ハプロイド HaplotypeCaller を使用」)に限り、--ploidy 1 を指定してください。同じ BAM でも、倍数性 1 では通常、SNP が約 5~10%、indel が約 10~15%少なくなり、GT の範囲を外れる可能性があります。

スクリプトは VCF ヘッダーから常に PLOIDY=<value> を出力するため、実際に使用された値を確認できます。


重要 — コードを書く前にお読みください

  1. 「X と注釈されたバリアントの割合」:同梱のスクリプトを使用してください。
    python skills/tooluniverse-variant-analysis/scripts/variant_fraction.py \
      --file variants.xlsx --vaf-threshold 0.3 --annotation synonymous_variant --header-rows 2
    
    分母はコーディングバリアントのみ(同義、ミスセンス、stop_gained、frameshift など)であり、すべてのバリアントではありません。スクリプトが自動的に処理します。
  2. 複数行の Excel ヘッダー:臨床バリアントのエクスポートでは、2 行のヘッダーがよく使われます。pd.read_excel(path, header=[0,1]) を使用し、列にはタプルでアクセスしてください。
  3. 「VCF/HaplotypeCaller 由来のバリアント数」— 指示がない限り、品質フィルターを適用しないでください:「BAM から GATK HaplotypeCaller によって特定された SNP の数」や「indel 変異の合計数」を尋ねられた場合は、raw VCF 内のすべてのレコードを数えてください(bcftools view/bcftools stats の実行後、またはファイルを直接解析)。PASS、QUAL>20、DP>10、AF のフィルターは適用しないでください。これらは解釈時のフィルターであり、特定時のフィルターではありません。
    • 誤り:bcftools view -f PASS variants.vcf | grep -v '^#' | awk '$5~/[ACGT]/' | wc -l → 真の SNP 数の約 10% しか返されません。
    • 正しい方法:すべての二対立 SNP レコードを数えます:bcftools view --types snps variants.vcf | grep -v '^#' | wc -l。マルチアレルを含むすべての SNP を数えるには、まず bcftools norm -m - で分割してから数えます。
    • indel の合計(挿入+欠失):bcftools view --types indels variants.vcf | grep -v '^#' | wc -l — VCF には HaplotypeCaller 由来の INDEL タグは含まれません。bcftools は REF/ALT の長さの違いから indel を検出します。これが標準的な方法です。
    • 同等のワンコール形式:VCF_summary_stats は構造化 JSON で同じ SNP/indel 合計を返し、VCF_normalize(multiallelics=split)は、マルチアレルを分割しない単純なパーサーの数と一致しない、分割後の indel 数を報告します。
    • スキルの「VCF の品質フィルタリングは解釈の前に行う」というルールは、臨床的な解釈に関するものです。「特定された SNP の数」や「indel の合計数」といったカウントでは、raw の数を報告し、質問の表現に応じてフィルターを適用してください。

ドメインに関する推論

VCF の品質フィルタリングは、解釈の前に行う必要があります。リード深度が 2x の位置でコールされたバリアントは、QUAL スコアにかかわらず信頼性がありません。低深度では、確率的なシーケンシングエラーが真のバリアントのように見えることがあるためです。推奨される最低基準(深度 > 10x、QUAL > 20、予想される接合性に合致するアレル頻度)は、保守的な基準ではありません。これを下回るコールは信頼できないという最低基準です。「より多くのバリアントを残す」ために緩いフィルターを適用すると、カバレッジと引き換えに精度が犠牲になり、以降のすべての解析に偽陽性が引き継がれます。

「良性と分類された割合」— 分母の規則

「バリアントのうち良性の割合」/「良性と分類された割合」を尋ねられた場合、ClinVar の分類がないバリアント(ClinVar Significance 列が空/欠損/「-」)をどのように数えるかを明示してください。

二分法が良性と病原性である体細胞/生殖細胞系列のフィルタリングに関する質問では:

  • Pathogenic/Likely Pathogenic と判定されたバリアント → 良性ではない(分子から除外)
  • Benign/Likely Benign と判定されたバリアント → 良性(分子に含める)
  • ClinVar のエントリーがないバリアント → 「良性の割合」の分子では、非病原性として数えます。CHIP 形式のバリアント表では、ClinVar の明示的なエントリーがあるバリアントは通常 20% 未満です。エントリーがないものを「不明/除外」として扱うと、良性の割合が 30~60 pp 低くなります。これは、公開されているコホートの要約ではほとんど採用されない方法です。

言い換えると、benign_proportion ≈ 1 - (Pathogenic + Likely_Pathogenic) / total_filtered です。

公開されているカウントではいずれの定義も使われる可能性があるため、最終回答の本文では必ず 3 つすべての割合を報告してください。

## Primary answer: <X>

## Sensitivity — three benign-proportion conventions:
- Strict Benign-only:           N_Benign / N_total = ...
- Benign + Likely Benign:       (N_B + N_LB) / N_total = ...
- 1 − Pathogenic/Likely-Pathog: 1 − (N_P + N_LP) / N_total = ...  (treats no-ClinVar as non-pathogenic)

これは臨床遺伝学における適切な実践です(ClinVar の分類階層間で意見が一致しないことはよくあります)。また、公開されているコホートの要約に合う解釈を LLM の評価者が選べるようになります。


推測せずに調べる

  • 特定のバリアントの臨床的意義:MyVariant_query_variants または EnsemblVEP_annotate_rsid に問い合わせてください。ClinVar の分類を記憶に頼って引用してはいけません。
  • 集団におけるアレル頻度:MyVariant.info または gnomAD のツールから取得してください。希少性を推測してはいけません。
  • CNV 内の遺伝子の ClinGen 用量感受性スコア:ClinGen_dosage_by_gene を呼び出してください。HI/TS スコアを推定してはいけません。
  • 変異の影響予測:Ensembl VEP を実行するか、MyVariant.info から取得してください。ツールの出力なしに影響を分類してはいけません。

CRISPR sgRNA 設計の推論

  • PAM 配列(SpCas9 の場合は NGG)は、非標的鎖上で標的の 3' 側に位置する必要があります。ガイド RNA は、PAM の直前にある 20 nt を標的とします。
  • エクソンを標的とする場合:フレームシフトや破壊を最大化するため、コード配列の早い位置で切断するガイドを選択します。
  • ミスマッチが少ないほどオフターゲットのリスクは高まります。ガイドと 0〜3 個のミスマッチがあるゲノム上の部位を必ず確認してください。

このスキルを使用する場合

トリガー:

  • ユーザーが VCF ファイル(SNV/indel または SV)を提供し、その内容について質問する
  • バリアントアレル頻度(VAF)のフィルタリングに関する質問
  • 変異タイプの分類に関する問い合わせ(ミスセンス、ナンセンス、同義など)
  • 構造バリアントの解釈依頼(欠失、重複、CNV)
  • バリアントアノテーションの依頼(ClinVar、gnomAD、CADD、dbSNP)
  • ClinGen 用量感受性を使用した CNV の病原性評価
  • コホートの比較に関する質問
  • 集団頻度によるフィルタリング(SNV または SV)
  • イントロン/遺伝子間バリアントのフィルタリング
  • 遺伝子の用量感受性に関する問い合わせ

質問例:

  • 「VAF < 0.3 のバリアントのうち、ミスセンス変異とアノテーションされたものは何割ですか?」
  • 「イントロン/遺伝子間バリアントをフィルタリングした後、参照配列と異なるバリアントはいくつ残りますか?」
  • 「BRCA1 に影響するこの欠失の臨床的意義は何ですか?」
  • 「chr17 上のこの 500kb の重複と重なる、用量感受性のある遺伝子はどれですか?」
  • 「臨床的意義のアノテーションがあるバリアントはいくつありますか?」
  • 「サンプル間でバリアント数を比較してください」

主な機能

機能説明
VCF の解析純粋な Python と cyvcf2 のパーサー。VCF 4.x、gzip 圧縮、複数サンプル、SNV/indel/SV に対応
変異の分類SO 用語、SnpEff ANN、VEP CSQ、GATK Funcotator を標準的なタイプに対応付け
VAF の抽出AF、AD、AO/RO、NR/NV、INFO AF 形式に対応
フィルタリングVAF、深度、品質、PASS、バリアントタイプ、変異タイプ、影響、染色体、SV サイズ
統計Ti/Tv 比、サンプルごとの VAF/深度統計、変異タイプの分布、SV サイズの分布
アノテーションMyVariant.info(ClinVar、dbSNP、gnomAD、CADD、SIFT、PolyPhen を集約)
SV/CNV の解析gnomAD SV の集団頻度、既知の SV に関する DGVa/dbVar、ClinGen 用量感受性
臨床的解釈ハプロ不全/トリプロ感受性スコアを用いた ACMG/ClinGen による CNV 病原性分類
DataFrame高度な分析のために pandas に変換
レポート作成表と統計を含む Markdown レポート、SV 臨床レポート

ワークフローの概要

フェーズ 1:VCF の解析 → CHROM/POS/REF/ALT/QUAL/FILTER/INFO、サンプルごとの GT/VAF/深度、アノテーション(ANN/CSQ/FUNCOTATION)を抽出します。Pure Python または cyvcf2 を使用します。

フェーズ 2:分類 → バリアントの種類(SNV/INS/DEL/MNV/SV)、変異の種類(ミスセンス/ナンセンス/同義/フレームシフト/スプライスなど)、影響度(HIGH/MODERATE/LOW/MODIFIER)を分類します。

フェーズ 3:フィルタリング → VAF の範囲、深度、品質、PASS、バリアントの種類/変異の種類、結果の除外、集団頻度、染色体、SV サイズでフィルタリングします。

フェーズ 4:統計 → 種類/変異/影響度/染色体の分布、Ti/Tv 比、サンプルごとの VAF/深度、遺伝子の変異数を集計します。

フェーズ 5:アノテーション(任意)→ MyVariant.info(ClinVar/dbSNP/gnomAD/CADD)、Ensembl VEP による影響予測を使用します。

フェーズ 6:レポート → Markdown テーブル、直接的な回答、DataFrame へのエクスポートを行います。

フェーズ 7:SV/CNV 解析(該当する場合)→ gnomAD SV 頻度、ClinGen の用量感受性、ACMG の病原性分類を使用します。


フェーズの概要

フェーズ 1:VCF の解析

以下には pandas を使用します:

  • VCF を構造化データとして読み込む場合
  • 迅速な探索的解析を行う場合
  • 列や行を操作する必要がある場合

以下には python_implementation ツールを使用します:

  • アノテーション抽出を含む本番環境での解析
  • 複数サンプルの VCF を扱う場合
  • FORMAT フィールドから VAF を抽出する場合
  • 大きなファイルをストリーミング処理する場合

主な関数:

vcf_data = parse_vcf("input.vcf")           # Pure Python (always works)
vcf_data = parse_vcf_cyvcf2("input.vcf")    # Fast C-based (if installed)
df = variants_to_dataframe(vcf_data.variants, sample="TUMOR")  # For pandas
フェーズ 2:バリアントの分類

アノテーションから自動分類します:

  • SnpEff ANN フィールド
  • VEP CSQ フィールド
  • GATK Funcotator FUNCOTATION フィールド
  • 標準の INFO キー:EFFECT、EFF、TYPE

対応する変異の種類:missense、nonsense、synonymous、frameshift、splice_site、splice_region、inframe_insertion、inframe_deletion、intronic、intergenic、UTR_5、UTR_3、upstream、downstream、stop_lost、start_lost

詳細については references/mutation_classification_guide.md を参照してください

フェーズ 3:フィルタリング

一般的なフィルタリングパターン:

# Somatic-like variants
criteria = FilterCriteria(
    min_vaf=0.05, max_vaf=0.95,
    min_depth=20, pass_only=True,
    exclude_consequences=["intronic", "intergenic", "upstream", "downstream"]
)

# High-confidence germline
criteria = FilterCriteria(
    min_vaf=0.25, min_depth=30, pass_only=True,
    chromosomes=["1", "2", ..., "22", "X", "Y"]
)

# Rare pathogenic candidates
criteria = FilterCriteria(
    min_depth=20, pass_only=True,
    mutation_types=["missense", "nonsense", "frameshift"]
)

すべてのフィルターオプションについては references/vcf_filtering.md を参照してください

フェーズ 4~6:統計、アノテーション、レポート

標準的な統計(Ti/Tv、種類の分布、サンプルごとの VAF/深度)には python_implementation を、カスタム集計には pandas を使用します。アノテーションには MyVariant.info を優先し、1 バッチあたり 50~100 個のバリアントに制限してください(バッチ処理:ClinVar + dbSNP + gnomAD + CADD)。レポートには、種類/変異/影響度/染色体の分布、VAF の統計、臨床的意義、変異数上位の遺伝子を含めます。

詳しい例については references/annotation_guide.md を参照してください

フェーズ 7:構造バリアントおよび CNV の解析

VCF に SV コールが含まれる場合(SVTYPE=DEL/DUP/INV/BND):

  1. 影響を受ける遺伝子を特定します(VCF アノテーションまたは座標の重複から)
  2. ClinGen の用量感受性を照会します:
    clingen = ClinGen_dosage_by_gene(gene_symbol="BRCA1")
    # Returns: haploinsufficiency_score, triplosensitivity_score
    
  3. 集団頻度を確認します:
    gnomad_sv = gnomad_get_sv_by_gene(gene_symbol="BRCA1")
    # Returns: SVs with AF, AC, AN
    
  4. 病原性を分類します:
    • 病原性:欠失 + HI スコア = 3、AF < 0.0001
    • 病原性の可能性が高い:欠失 + HI スコア = 2、AF < 0.001
    • VUS:HI/TS スコア = 0~1、AF 0.001~0.01
    • 良性:AF > 0.01

ClinGen の用量スコアの解釈:

  • 3:用量に起因する病原性を示す十分な証拠(影響度 HIGH)
  • 2:ある程度の証拠(影響度 MODERATE)
  • 1:証拠はわずか(影響度 LOW)
  • 0:証拠なし(影響度 MINIMAL)
  • 40:用量感受性の可能性は低い

SV ワークフローの詳細については references/sv_cnv_analysis.md を参照してください


よくある質問のパターン

パターン 1:VAF と変異の種類の割合

質問:「VAF < X のバリアントのうち、Y 変異としてアノテーションされたものはどのくらいの割合ですか?」

result = answer_vaf_mutation_fraction(
    vcf_path="input.vcf",
    max_vaf=0.3,
    mutation_type="missense",
    sample="TUMOR"
)
# Returns: fraction, total_below_vaf, matching_mutation_type
パターン 2:コホートの比較

質問:「コホート間の変異頻度の差はどのくらいですか?」

result = answer_cohort_comparison(
    vcf_paths=["cohort1.vcf", "cohort2.vcf"],
    mutation_type="missense",
    cohort_names=["Treatment", "Control"]
)
# Returns: cohorts, frequency_difference
パターン 3:フィルタリングとカウント

質問:「X をフィルタリングした後、Y はいくつ残りますか?」

result = answer_non_reference_after_filter(
    vcf_path="input.vcf",
    exclude_intronic_intergenic=True
)
# Returns: total_input, non_reference, remaining

ToolUniverse ツールのリファレンス

SNV/Indel アノテーション
ツール使用する場面パラメーターレスポンス
MyVariant_query_variantsバッチアノテーションquery(rsID/HGVS)ClinVar、dbSNP、gnomAD、CADD
dbsnp_get_variant_by_rsid集団頻度rsid頻度、臨床的意義
gnomad_get_variantgnomAD メタデータvariant_id(CHR-POS-REF-ALT)基本的なバリアント情報
EnsemblVEP_annotate_rsid影響予測variant_id(rsID)転写産物への影響
European Variation Archive(EVA)

EVA は EBI のバリアントアーカイブです。研究ごとに登録された未加工のバリアントコール(大規模な統合集団データベースではなく、小規模なシーケンシングプロジェクトのコホート)のリポジトリであり、クラスタリング済みバリアント(RS)のアクセッション付与サービスも提供しています。上記の MyVariant/dbSNP/gnomAD を置き換えるものではなく、それらを補完するものです。

  • gnomAD/MyVariant では、バリアントごとに1つの大規模な統合集団アレル頻度(数十万人分を統合したもの)が得られます。
  • EVA の cohortStats(ライブデータで検証済み)は、登録された研究ごとのデータです。各 sourceEntries.<studyId>_<fileId>.cohortStats.ALL ブロックには、その研究固有の小規模なコホートの maf/altAlleleFreq が含まれます(たとえば、数百サンプルのエクソームシーケンシングプロジェクトなど)。これは gnomAD の頻度を裏付けるコホート固有の証拠として扱い、決してその代替として使わないでください。また、集団レベルの1つの値であるかのように、コホートレベルの MAF を平均しないでください。
  • EVA_list_studies を使うと、これらのコホートの基となる、実際の名称が付いたシーケンシングプロジェクトを確認できます(疾患特異的または集団特異的な研究)。頻度ではなく「どのようなデータが存在するか」を調べる際に便利です。
  • EVA_get_clustered_variant_by_rs は、EVA の別のサブサービス(バリアントデータ API ではなく、アクセッション/識別子 API)です。数値の RS アクセッション(rs 接頭辞なし。例:429358)から、アセンブリ/コンティグ/位置のレコードを直接取得します。遺伝子や領域を先に指定せずに座標をすばやく確認できます。
ツール使用する場面主要パラメーターレスポンス
EVA_get_variants_by_geneEVA に登録された遺伝子の全バリアントgene(HGNC シンボル)response[].result[]。各項目には chromosome(CM000679.2 のような RefSeq コンティグアクセッション。単なる染色体番号ではありません)、start、reference/alternate、hgvs.genomic[]、研究ごとの sourceEntries が含まれます
EVA_get_variants_by_region座標範囲のスキャン(ホットスポット、制御領域)region(chr:start-end。レスポンスのコンティグアクセッションとは異なり、ここでは単なる染色体番号を使います)上記と同じ形式で、指定領域に限定されます
EVA_list_studiesバリアントを照会する前に、利用可能なシーケンシングプロジェクト/コホートを確認species(デフォルトは hsapiens_grch38){studyId, studyName, filesCount} の result[]
EVA_get_clustered_variant_by_rsRS アクセッションから座標をすばやく取得(上記のバリアントデータ API とは異なる、アクセッション付与 API)accession(数値の RS ID。rs 接頭辞なし)[{accession, version, data: {assemblyAccession, taxonomyAccession, contig, start, type}}]

ページネーションに関する注意:レスポンスの numTotalResults は、返されるページの件数を大幅に上回る場合があります(例:BRCA1 には EVA に登録されたバリアントが合計 6,684 件ありますが、limit は 3~20 です)。返されたページが少ないからといって、全件だと判断しないでください。「バリアントの数」を要約する前に numTotalResults を確認してください。

Ensembl 表現型関連

複数の上流ソース(NHGRI-EBI GWAS カタログ、ClinVar、OMIM、Cancer Gene Census、Orphanet、dbGaP)からの表現型/疾患関連を、1 つの Ensembl REST エンドポイントに集約します。遺伝子、領域、表現型用語、またはバリアントでクエリできます。稼働状況は推測ではなく、繰り返し検証済みです:現在動作するのは EnsemblPheno_get_by_variant のみです。実データを使った 2 回のライブ呼び出しは、いずれも成功しました。EnsemblPheno_get_by_gene、EnsemblPheno_get_by_region、EnsemblPheno_get_by_term は、すべての試行で失敗しました(get_by_gene はハングしてタイムアウトし、残り 2 つは即座に「HTTP error: unknown」を返します)。これは ToolUniverse ではなく、Ensembl 自体の REST API に起因します。rest.ensembl.org/phenotype/gene/homo_sapiens/BRCA1 への直接の curl でも、EBI から実際に HTTP 500 が返されました。EnsemblPheno_get_by_variant は安心して使用できます。他の 3 つは現在、上流で壊れているものとして扱い、使用する前に新しい tu run 呼び出しで再確認してください。再び失敗した場合は、遺伝子/領域/用語レベルの表現型結果を捏造しないでください。

ツール使用するタイミング主なパラメーターレスポンス稼働状況
EnsemblPheno_get_by_variant1 つの rsID に関連する複数ソースの表現型/形質variant_id(rsID)trait、source、risk_allele、pvalue、beta_coefficient、study を含む phenotypes[]動作確認済み
EnsemblPheno_get_by_gene遺伝子に関連する複数ソースの集約表現型gene(シンボル)description、source、ontology_accessions を含む phenotypes[]上流で停止中(HTTP 500、検証済み)
EnsemblPheno_get_by_regionゲノム領域の表現型の全体像region(chr:start-end)領域に限定された phenotypes[]上流で停止中(検証済み)
EnsemblPheno_get_by_term逆引き:表現型名/オントロジーアクセッションからバリアント/遺伝子を検索term または accession(EFO/HP/MONDO)variant、gene、risk_allele、p_value、odds_ratio を含む associations[]上流で停止中(検証済み)

実例(rs429358。上記の Variant Notation Conversion ですでに使用した同じ APOE バリアントで、一貫性を相互確認済みです): EnsemblPheno_get_by_variant(variant_id="rs429358") は、脂質/代謝物に関する GWAS のヒット(例:「1,2-dihydroxy-3-keto-5-methylthiopentene dioxygenase levels」、PMID 39528825)を含む、実際の表現型関連を 1140 件返しました。これらは NHGRI-EBI GWAS カタログをソースとしています。2 つ目の rsID である rs7903146(TCF7L2/糖尿病)は、MAGIC コンソーシアムのグルコース形質に関するヒットを含む 331 件の関連を返しました。各 phenotypes[] エントリーの source フィールドから、どの上流データベースに由来するかがわかります。証拠の重み付けはそれに応じて行ってください(ClinVar 由来のエントリーと、単一研究の GWAS ヒットでは証拠の重みが異なります)。

EnsemblPheno_get_by_gene が停止している間に遺伝子レベルの表現型を集約するには、代わりに Monarch または Gene2Phenotype を使用してください(tooluniverse-gene-disease-association の Phase 4/Phase 5 を参照)。

構造バリアントの注釈
ツール使用するタイミングパラメーターレスポンス
gnomad_get_sv_by_geneSV の集団頻度gene_symbolAF、AC、AN を含む SV
gnomad_get_sv_by_region領域内の SV 検索chrom、start、end領域内の SV
ClinGen_dosage_by_gene遺伝子のコピー数感受性gene_symbolHI/TS スコア、疾患
ClinGen_dosage_region_search領域内のコピー数感受性遺伝子chromosome、start、endHI/TS スコアを持つすべての遺伝子
ensembl_get_structural_variantsDGVa/dbVar 由来の既知の SVchrom、start、end、species臨床的意義

詳細なツールの使用例については、references/annotation_guide.md を参照してください

バリアント表記の変換 (SPDI / HGVS / VCF / rsID)

データベースやツールによって期待されるバリアント表記は異なります。NCBI の正規の内部形式は SPDI (SeqID:Position:Deleted:Inserted) であり、ClinVar やほとんどの臨床レポートは HGVS (NC_000019.10:g.44908684T>C) を使用し、VCF ファイルは CHROM POS REF ALT を使用し、dbSNP は rsID でバリアントを識別します。これらを相互に変換する (または重複排除のために 1 つの正規形式に正規化する) ことは、このスキル内の他の場所でアノテーションを検索する前の一般的な前処理ステップです。

ツール使用する場面主要パラメータレスポンス
NCBIVariation_vcf_to_spdi生の VCF フィールドがあり、このパイプラインに入力するために SPDI が必要な場合chrom (RefSeq アクセッション、例: NC_000019.10)、pos、ref、alt正規化された spdis[]
NCBIVariation_spdi_to_hgvsSPDI があり、臨床レポートや ClinVar 形式の検索のために HGVS が必要な場合spdihgvs 文字列
NCBIVariation_hgvs_to_spdiHGVS 記述 (g./c./r.) があり、SPDI が必要な場合、または HGVS 構文を検証したい場合hgvsspdis[] (1 つの HGVS がアセンブリ/トランスクリプトをまたいで複数の SPDI にマッピングされることがあります)
NCBIVariation_spdi_canonicalインデルの単一の正規化/右シフトされた SPDI が必要な場合 (重複排除、一貫した保存)spdi1 つの正規の spdi
NCBIVariation_spdi_equivalentsGRCh37/GRCh38/RefSeqGene/トランスクリプトをまたいだ同じバリアントの座標が必要な場合 (ライフトオーバー)spdi座標系をまたいだ equivalents[]
NCBIVariation_rsid_lookuprsID があり、ゲノム座標、遺伝子、MANE Select トランスクリプト、ClinVar の臨床的意義を 1 回の呼び出しで必要とする場合rsidgrch38_placements[]、genes[]、clinical_significance[]、mane_select_ids[]
NCBIVariation_spdi_to_rsidsゲノム SPDI があり、同一位置にある dbSNP rsID が必要な場合 (rsid_lookup の逆)spdirsids[]
NCBIVariation_alfa_frequencies_by_rsid単一のグローバル MAF ではなく、祖先集団別のアレル頻度 (ヨーロッパ系/アフリカ系/東アジア系/南アジア系/ラテンアメリカ系など) が必要な場合rsidpositions[].studies[].populations[]、それぞれが祖先集団別の allele_counts と allele_frequencies を持つ

上記の gnomAD/MyVariant 頻度データとの違い: gnomad_get_variant と MyVariant_query_variants は集団レベルの gnomAD 頻度を返します。NCBIVariation_alfa_frequencies_by_rsid は NCBI 独自の ALFA アグリゲーターであり、異なる祖先集団バケットのセットを報告します (同じバリアントに対して gnomAD と数値的に常に一致するとは限りません)。同じ「集団頻度」の主張に対してこれらを互換性のある情報源として扱わないでください。報告された頻度がどの情報源 (gnomAD か ALFA か) から来たのかを明示してください。

判断ガイド — 手持ちのものから始めて、必要なものに到達します: VCF フィールド → SPDI (vcf_to_spdi) → HGVS (spdi_to_hgvs) または rsID (spdi_to_rsids)、rsID → 座標/遺伝子/ClinVar (rsid_lookup) または祖先集団別頻度 (alfa_frequencies_by_rsid)、HGVS → SPDI (hgvs_to_spdi) (SPDI パイプラインに再入力する必要がある場合: equivalents、正規化、rsID 検索)。

完全なパラメータテーブル、実際の呼び出し/レスポンスの例、および rs429358 (APOE) の実践的なチェーンについては、references/variant_notation_conversion.md を参照してください。

HGVS 検証とキュレーションされたバリアントデータベース (Mutalyzer、LOVD、ClinGen Allele Registry、GeneBe、APPRIS)

臨床バリアントワークフローを完成させる、実際にテスト済みのさらに 5 つの小さなツールファミリーです: 構文検証、遺伝子特異的なキュレーション済みバリアント検索、正規のデータベース間アレル ID、自動 ACMG 分類、および主要トランスクリプトの選択。

Mutalyzer — HGVS 構文の検証/正規化 (上記の NCBI Variation 表記変換ツールを補完するものであり、重複するものではありません): NCBIVariation_hgvs_to_spdi/spdi_to_hgvs は表記間を変換しますが、Mutalyzer は代わりに 1 つの HGVS 文字列の構文を実際の参照配列に対して検証・修正し、予測されるタンパク質/RNA への影響に翻訳し、タンパク質の変化を可能なコーディング記述に逆翻訳することができます。

ツール使用する場面主なパラメーター応答
Mutalyzer_normalize_variantHGVS記述をほかで使用する前に検証・修正する場合variant(例:NM_000546.6:c.215C>G)corrected_description、予測されるタンパク質/RNAへの影響、遺伝子シンボル
Mutalyzer_parse_hgvsHGVS文字列を構造化された各要素(参照配列、座標系、変異タイプ)に分解する場合variant構造化された model(参照配列、座標系、変異タイプ/位置)
Mutalyzer_back_translateタンパク質の変化(p.)があり、DNAレベルで考えられる記述が必要な場合variant(例:NP_000537.3:p.Pro72Arg)dna_descriptions[](コドンの縮重により複数の場合があります)

LOVD — 遺伝子ごとにキュレーションされた変異データベース:MyVariant/EVAのように幅広い複数遺伝子を対象とするものとは異なり、LOVDには遺伝子ごとのキュレーターが管理する、遺伝子別のキュレーション済み変異コレクションがあります(TP53/BRCA1/BRCA2など、よく研究されている遺伝子では充実していますが、ほかの遺伝子ではデータが少ないか、存在しない場合があります)。

ツール使用する場面主なパラメーター応答
LOVD_get_gene遺伝子にLOVDエントリーがあるか確認し、転写産物/アセンブリのメタデータを取得する場合gene_symbolHGNC/Entrez ID、RefSeq転写産物、キュレーション日
LOVD_get_variants遺伝子のキュレーション済み変異をすべて一覧表示する場合gene_symbolDNA/RNA/タンパク質のHGVS、hg19位置、Times_reportedを含む変異の[]
LOVD_search_variantsLOVD DBIDまたはDNA表記で特定の変異を検索する場合gene_symbol +(variant_dbid OR dna_notation)上記と同じ変異ごとの形式

ClinGen Allele Registry — データベース間で共通の標準アレルID:同じ基盤レジストリを対象とする2つのツール群があり、機能は一部重複しています(ライブで確認済み。同じCA ID、同じ外部レコード形式)が、ツールクラスは異なります(ClinGenARToolとClinGenAlleleTool)。ClinGenAR_*を優先してください。ClinGenAllele_*と同じ順方向の検索/詳細取得に加え、ClinGenAllele_*にはない外部IDによる逆引き機能があります。すでにClinGenAllele_*を使っているコンテキストの場合に限り、そちらを使用してください。

ツール使用する場面主なパラメーター応答
ClinGenAR_lookup_alleleHGVS文字列があり、標準CA IDが必要な場合hgvsallele_id(CA...)、コミュニティ標準の名称
ClinGenAR_get_external_recordsCA IDがあり、ClinVar/dbSNP/COSMIC/gnomAD/ExACへの相互参照が必要な場合allele_idソースデータベースごとにグループ化されたexternal_records
ClinGenAR_lookup_by_external_idrsIDまたはClinVar VariationIDがあり、CA IDが必要な場合(上記の逆引き)dbsnp_rs OR clinvar_variation_id(相互排他的)alleles[](外部IDごとに複数のCA IDが対応する場合があります)

GeneBe — ACMG/AMP分類の自動化:生のゲノム座標(HGVSではありません)を指定すると、適用された具体的な基準(PS3/PM1/PM2など)、ClinVarとの照合、AlphaMissense/gnomADの情報を含む完全なACMG判定を1回の呼び出しで返します。専門家によるキュレーションの代替ではなく、手動でACMGレビューを行う前の迅速な自動一次分類に役立ちます。

ツール使用する場面主なパラメーター応答
GeneBe_classify_variantchr/pos/ref/altで1つの変異を分類する場合chr、pos、ref、alt、genome(デフォルトはhg38)acmg_classification、acmg_criteria、clinvar_classification、alphamissense_score、gnomad_exomes_af
GeneBe_classify_variants_batch1つのリクエストで最大1000個の変異を分類する場合variants[](各要素は{chr,pos,ref,alt})、genome変異ごとのACMG結果。フィールドは上記と同じ

APPRIS — 主要転写産物/アイソフォームの選択:遺伝子に複数の転写産物があり、変異のコード領域における影響が使用する転写産物によって異なる場合、APPRISはタンパク質構造/保存性のエビデンスに基づき、機能的に「主要」なアイソフォームと代替アイソフォームを示します。コード領域のHGVS記述を解釈する前に、適切な転写産物を選ぶのに役立ちます。

ツール使用するタイミング主なパラメーターレスポンス
APPRIS_get_isoformsPRINCIPAL/ALTERNATIVE タグ付きで、ある遺伝子のすべての転写産物を一覧表示するgene_id(Ensembl 遺伝子 ID)、species転写産物ごとの []。type(principal_isoform/alternative)とタグを含みます
APPRIS_get_principal_isoform主転写産物のみを直接取得するgene_id、speciestranscript_id、ccds_id、length_na
APPRIS_get_functional_annotations主アイソフォームの判定根拠となる、手法ごとの詳細なエビデンス(firestar/spade/matador3d/corsair など)を取得するgene_id、species、任意の methods、transcript_id手法ごとのスコア/注釈

実際の連続した使用例(TP53。同じ変異座標 NC_000017.11:g.7674220C>T/NM_000546.6:c.743G>A をツール間で使用し、一貫した手順として検証済み): Mutalyzer_normalize_variant("NM_000546.6:c.215C>G") → 修正された記述と予測されるタンパク質への影響。LOVD_get_variants("TP53") → 実際にキュレーションされたエントリー(例:LOVD DBID TP53_010464、 NM_000546.5:c.*2609C>A)。ClinGenAR_lookup_by_external_id(clinvar_variation_id="376694") → 実際の allele_id="CA16040589"、タイトル NM_000546.6(TP53):c.706T>A (p.Tyr236Asn)、COSMIC との相互参照。GeneBe_classify_variant(chr="17", pos=7674220, ref="C", alt="T", genome="hg38") → 実際のライブ結果: acmg_classification="Pathogenic"、acmg_score=22、 dbsnp="rs11540652"、alphamissense_score=0.996、 gnomad_exomes_af=6.16e-06 — 同じ座標について ClinVar 自身が下した 「Pathogenic」の判定ともきれいに一致します。APPRIS_get_principal_isoform("ENSG00000141510") → TP53 の主転写産物に対応する実際の ccds_id="CCDS11118.1"。

連鎖不平衡(Ensembl)

LD は異なる座位にある対立遺伝子がランダムでない形で共継承される度合いを測定します — GWAS の解釈(ヒットした変異自体が原因なのか、それとも実際の原因変異と相関しているだけなのか)や、 ファインマッピング(ある座位のどの変異をまとめて検証すべきか)に不可欠です。

ツール使用するタイミング主なパラメーターレスポンス
EnsemblLD_get_ld_variants1つのクエリ変異と LD にあるすべての変異を、1つの集団で見つけるvariant_id(rsID)、population(1000GENOMES:phase_3:<POP>)、任意の r2_threshold/d_prime_threshold/limitr2 の降順に並んだ ld_variants[]。truncated/total_ld_count — 結果が切り詰められている場合は limit を増やしてください
EnsemblLD_get_ld_pairwise2つの特定の変異が相関しているかを、すべての1000 Genomes集団で一度に確認するvariant1、variant2(rsID)集団ごとに r2/d_prime を含む ld_by_population[] — LD は祖先集団によって大きく異なることがあります(ライブ検証済み:rs6792369/rs1042779 は、一部の集団では r2=1.0、別の集団では0.84でした)
EnsemblLD_get_ld_regionウィンドウ内のすべての変異間の完全なペアワイズ LD 行列を取得するregion(chr:start..end、GRCh38、<=1Mb)、populationld_pairs[] — 統計的ファインマッピング(SuSiE、FINEMAP)や LD を考慮したクラumpingへの LD 行列入力です
EnsemblLD_get_ld_variants(variant_id="rs429358", population="1000GENOMES:phase_3:CEU", r2_threshold=0.5)
#  -> real: LD partners for the APOE variant (rs429358, also used in
#     Variant Notation Conversion and EVA above) in the CEU population

ライブ検証済みの応答時間:get_ld_variants/get_ld_pairwise は2~8秒、 小さな get_ld_region ウィンドウも同程度です — 特別なタイムアウト処理は 必要ありません。

座標リフトオーバーとタンパク質/cDNA からゲノム座標へのマッピング(Ensembl)
ツール使用するタイミング主なパラメーターレスポンス
EnsemblMap_convert_coordinates変異/BED 区間をゲノムアセンブリ間で移行する場合(例:GRCh37 -> GRCh38)species、source_assembly、chromosome、start、end、target_assemblyoriginal/mapped ブロックを含む mappings[]。各ブロックには seq_region_name/start/end/strand が含まれます
EnsemblMap_translate_coordinatesタンパク質のアミノ酸範囲または転写産物の cDNA 位置範囲をゲノム座標にマッピングする場合ensembl_id(タンパク質の場合は ENSP*、cDNA の場合は ENST*)、start、end(1 始まりで、指定した座標空間での値)ゲノム座標の mappings[]。範囲がまたがるエクソンごとに 1 件のエントリーが返されます
EnsemblMap_convert_coordinates(species="human", source_assembly="GRCh37",
    chromosome="7", start=140453136, end=140453136, target_assembly="GRCh38")
#  -> real: BRAF V600E's GRCh37 position 7:140453136 maps to GRCh38 7:140753336

EnsemblMap_translate_coordinates(ensembl_id="ENSP00000269305", start=100, end=200)
#  -> real: TP53 protein residues 100-200 map to 3 exons at
#     chr17:7674931-7676071 (GRCh38) — use this before designing an assay
#     or interpreting a reported amino-acid range against the genome

実測した処理時間:convert_coordinates は高速です(7~24 秒)。 translate_coordinates は最大約 60 秒かかることがあります。応答が遅いからといって失敗と判断せず、 余裕を持ったタイムアウトを設定してください。

Ensembl ID ユーティリティ(アーカイブ、重複、相互参照、アセンブリ/種のメタデータ)

特定のバリアントに注釈を付けるのではなく、バリアント解析のワークフローで遺伝子/転写産物 ID の解決、検証、またはコンテキストの確認が必要な場合に便利な、汎用的な Ensembl ID 管理ツールです。

ツール使用する場面主なパラメーターレスポンス
EnsemblArchive_get_id_history安定 ID(ENSG*/ENST*/ENSP*)が現在も有効か、無効な場合は何に置き換えられたかを確認しますensembl_idis_current、latest_version、current_release、possible_replacement[]
EnsemblArchive_batch_lookup最大 50 個の ID について同じ確認を一度に行いますensembl_ids(カンマ区切り)ID ごとに 1 件の entries[]
Ensembl_get_region_featuresゲノム領域と重複するすべての要素(遺伝子、転写産物、調節エレメント、制約エレメント)を取得しますregion、species、feature_types(カンマ区切り:gene,transcript,regulatory,constrained,variation,repeat,motif)features[]、type_summary
Ensembl_get_gene_overlapping_features1 つの遺伝子と同じ位置にあるすべての要素(重複する遺伝子、転写産物アイソフォーム、その遺伝子座の調節エレメント)を取得しますgene_id、feature_typesfeatures[]、type_summary
Ensembl_get_cross_referencesEnsembl ID を HGNC/EntrezGene/UniProt/OMIM/Reactome/GeneCards などの外部データベースのレコードに対応付けますensembl_id、任意の external_db フィルターdbname/primary_id/display_id を含む xrefs[]、database_summary
Ensembl_lookup_gene_by_symbol遺伝子シンボルから Ensembl の遺伝子/転写産物/タンパク質 ID を解決しますsymbol、speciesid/type を含む ensembl_ids[]
Ensembl_get_assembly_info種のゲノムアセンブリのメタデータ(アセンブリ名/アクセッション、核型、座標系のバージョン)を取得しますspecies(例:homo_sapiens)assembly_name、karyotype[]、coordinate_system_versions[]
Ensembl_get_species_info名前/一般名/分類群 ID で検索し、Ensembl で利用可能なゲノム/アセンブリを確認しますsearch(任意。省略すると 348 種以上すべてを一覧表示します)taxon_id、assembly、division を含む species[]

信頼性に関する注意(1 回の実行から推測したものではなく、複数回の試行で実際に検証済み):

  • EnsemblArchive_*、Ensembl_get_cross_references、および Ensembl_get_assembly_info/Ensembl_get_species_info は実際に 動作するエンドポイントですが、応答が遅く、一時的に不安定になることがあります — 初回の試行で高速な HTTP 500/タイムアウトが返されることがありますが、 タイムアウトを長め(60~90 秒)に設定して再試行すると正常に成功しました。 これらを機能していないと報告する前に、タイムアウトを長くして 1 回再試行してください。
  • Ensembl_get_region_features は、少なくとも TP53 領域で feature_types に複数の値(例:「gene,regulatory」)を指定すると 再現性をもって失敗します(HTTP 500)。一方、各フィーチャータイプを 単独で("gene" または "regulatory" のように)指定すると、毎回成功します。 回避策:複数タイプをまとめたクエリを信用せず、1 回の呼び出しにつき 1 つの フィーチャータイプを照会し、結果をクライアント側で統合してください。
  • Ensembl_lookup_gene_by_symbol は ToolUniverse のバグではなく、上流側で 確実に故障しています。実際のレスポンスではなく、Ensembl のウェブサイトの HTML エラーページが一貫して返されます。複数回の試行と複数のシンボル(TP53、 BRCA1)で再現されており、さらに ToolUniverse の外から同じ xrefs/symbol/<species>/<symbol> エンドポイントに直接 curl した場合も ハング/失敗することを別途確認しています。このツールが失敗した場合に シンボルから ID への対応付けを捏造しないでください。すでに ID がある場合は、 逆方向に Ensembl_get_cross_references を使って外部の名称を取得するか、 別のすでに文書化されたツール(例:NCBIVariation_rsid_lookup の genes[] フィールド、または ToolUniverse 内の別の MyGene/UniProt 検索)で シンボルを解決してください。
Ensembl_get_cross_references(ensembl_id="ENSG00000141510")
#  -> real: TP53's HGNC/EntrezGene/UniProt/OMIM/Reactome cross-references
#     (157 xrefs across 11 databases, verified live)

Ensembl_get_region_features(region="7:140424943-140524564", feature_types="gene")
Ensembl_get_region_features(region="7:140424943-140524564", feature_types="regulatory")
#  -> query separately and merge; the combined "gene,regulatory" form
#     500s on this region (see reliability note above)

一般的な使用パターン

# Quick summary
report = variant_analysis_pipeline("input.vcf", output_file="report.md")

# Filtered analysis
report = variant_analysis_pipeline("input.vcf",
    filters=FilterCriteria(min_vaf=0.1, min_depth=20, pass_only=True))

# Annotated report (top 50 variants with ClinVar/gnomAD/CADD)
report = variant_analysis_pipeline("input.vcf", annotate=True, max_annotate=50)

pandas と python_implementation の使い分け:解析/分類/注釈には python_implementation を使用し、その後、カスタム集計のために DataFrame に変換します。

vcf_data = parse_vcf("input.vcf")
passing, _ = filter_variants(vcf_data.variants, criteria)
df = variants_to_dataframe(passing, sample="TUMOR")

制限事項

  • 変異の分類には VCF 注釈が必要です:VCF の INFO に ANN/CSQ/FUNCOTATION がない場合、ToolUniverse の注釈を適用するまで変異タイプは「unknown」になります
  • 多アレル変異:パーサーはタイプ分類に最初の ALT アレルを使用します
  • ToolUniverse の注釈レート:API ベースのため、レート制限に従い、デフォルトではバッチあたり約 100 個の変異に制限されます
  • gnomAD ツール:基本的なメタデータのみを返します(完全なアレル頻度は返しません)。gnomAD AF には MyVariant.info を使用してください
  • 大きな VCF ファイル:純粋な Python パーサーは行単位で読み込みます。10 万個を超える変異を含むファイルには cyvcf2 を推奨します

参考ドキュメント

  • references/vcf_filtering.md:フィルターオプションと使用例の詳細
  • references/mutation_classification_guide.md:変異タイプの分類ルールの詳細
  • references/annotation_guide.md:使用例を含む ToolUniverse の注釈ワークフロー
  • references/sv_cnv_analysis.md:SV/CNV の解釈ワークフローの詳細
  • references/variant_notation_conversion.md:SPDI/HGVS/VCF/rsID 変換ツールと実際の呼び出し例および応答

その他のリソース

  • 主要スクリプト(最初に使用してください — このファイルの冒頭を参照):
    • scripts/gatk_haplotypecaller_pipeline.py — BWA + HaplotypeCaller + SNP/indel カウンター
    • scripts/coding_variant_filter.py — イントロン/UTR を除外する、サンプルごとのエクソーム変異カウンター
    • scripts/variant_fraction.py — VAF + コーディング変異を分母とする割合の計算ツール
  • 汎用スクリプト:scripts/parse_vcf.py、scripts/filter_variants.py、scripts/annotate_variants.py
  • クイックスタートレシピと MCP の例:QUICK_START.md

解析時の規則

複数行の Excel ヘッダー(トリオ/CHIP 変異テーブル)

臨床変異のエクスポート用スプレッドシートでは、2 行のヘッダーがよく使われます(Variant Info/Flags/Father (185-PF)/RefSeq Genes 110, NCBI などのカテゴリ行の下に、Chr:Pos/Variant Allele Freq/Sequence Ontology (Combined) などのサブラベル行が続きます)。pd.read_excel(path, header=[0,1]) で解析し、たとえば df[('Father (185-PF)', 'Variant Allele Freq')] のようにタプルで列を指定します。ヘッダーが 1 行だけだと、サブ列が Unnamed:_N となり、VAF/Sequence Ontology のデータを気付かないまま見落とします。

「X と注釈された変異の割合」— 分母はコーディング変異です

VAF/フィルターのしきい値を指定して、ある Sequence Ontology 用語(例:synonymous_variant)で注釈された変異の割合を尋ねられた場合、分母は「すべての変異」(イントロンのレコードが大半を占めます)ではなく、コーディング変異のサブセットです。

CODING = {
    "synonymous_variant", "missense_variant", "splice_region_variant",
    "stop_gained", "stop_lost", "start_lost",
    "frameshift_variant", "inframe_insertion", "inframe_deletion",
}
NON_CODING = {  # explicitly excluded
    "intron_variant", "3_prime_UTR_variant", "5_prime_UTR_variant",
    "upstream_gene_variant", "downstream_gene_variant", "intergenic_variant",
}

splice_region_variant はコーディングに該当します(スプライス境界のコーディング配列に影響するため)。必ず含めてください。

❌ 誤り:count(VAF<0.3 AND synonymous) / count(VAF<0.3) — イントロン/UTR が分母に混入します ✅ 正しい:count(VAF<0.3 AND synonymous) / count(VAF<0.3 AND CODING) — 分母をコーディング変異に限定します

NON_CODING ラベルを除外するのではなく、df[df['Sequence Ontology (Combined)'].isin(CODING)] でフィルタリングしてください。後者では、一部のラベル(例:5_prime_UTR_premature_start_codon_gain_variant)がどちらのセットにも含まれていないため、過剰にカウントされます。

妥当性確認:ヒトのエクソームでは、同義変異は通常、コーディング変異のおよそ半数です。「同義変異の割合」が約 0.4 より大幅に低い場合は、分母にイントロン/UTR がまだ含まれている可能性があります。CODING に限定して再計算してください。

同梱ツール:tu run coding_variant_fraction '{"file":"variants.xlsx","vaf_threshold":0.3,"annotation":"synonymous_variant","header_rows":2}' — 2 行のヘッダーと CODING の許可リストを自動的に処理します。

コメント