tooluniverse-epigenomics

ゲノミクスおよびエピゲノミクスのデータ処理

⚠️ 最優先ルール：ロング形式のメチル化 CSV では、ユニークな位置ではなく行数を数える

入力がロング形式のメチル化 CSV（(sample, CpG_position) ごとに 1 行） （例：Pos, Chromosome, MethylationPercentage の各列）の場合、「フィルタリングで除外されるサイト数」は、ほとんどの場合、除外される行数を意味し、ユニークな位置の除外数ではありません。この 2 つの答えには、およそ n_samples 倍の差があります。

質問の表現	意味
「フィルタリングで除外されるサイト数は？」	除外される行数（= サンプル数 × フィルター条件を満たさない位置数）
「フィルターを通過するユニークな CpG サイト数は？」	ユニークな位置（Pos で重複を除去してからフィルタリング）

❌ 誤り：df.drop_duplicates(["Pos"]).query("MethylationPercentage<10 or >90") の後に len(filtered) → ユニークな位置数をカウントします（通常 100～1500）

✅ 正しい：df.query("MethylationPercentage<10 or MethylationPercentage>90") の後に len(df) - len(filtered) → 行数をカウントします（通常 10k～30k）

データに 1000 以上の位置 × 20 以上のサンプルがあるのに、答えが 2000 未満の場合は、早い段階で重複を除去しています。報告する前に、質問の名詞を読み直してください。

---

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

以下の手順に従う前に、データフォルダーを調べてください。

• *_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 倍になることもよくあります）。

エピゲノミクス解析向けに、Python による計算（pandas、scipy、numpy、pysam、statsmodels）と ToolUniverse のアノテーションツールを組み合わせた、本番環境対応のスキルです。

推測せず、調べる

科学的な事実に確信が持てない場合は、まずデータベースを検索してください。

使用する場面

メチル化データ、ChIP-seq ピーク、ATAC-seq、マルチオミクス統合、ゲノム全体のエピゲノミクス統計に使用します。キーワード：メチル化、CpG、ChIP-seq、ATAC-seq、ヒストン、クロマチン、エピジェネティクス。

対象外：RNA-seq の DEG、バリアントコール、遺伝子エンリッチメント、タンパク質構造。

---

主な原則

1. データを最優先 - 解析の前に読み込み、内容を確認します
2. 質問に基づく - 具体的な数値の答えを導き出します
3. 座標系を意識する - ゲノムビルド（hg19/hg38/mm10）と chr プレフィックスを追跡します
4. 統計的厳密さ - FDR 補正と効果量によるフィルタリングを行います
5. CpG の特定 - Illumina プローブ ID とゲノム座標を解析します

PRIMARY SCRIPT — methylation_density.py（CpG 密度に関する質問では最初に使用）

ロング形式のメチル化 CSV（Pos, Chromosome, MethylationPercentage）と 染色体長 CSV を組み合わせる場合は、独自に pandas で処理する前に、必ず同梱のスクリプトを実行してください。スクリプトは一般的な指標をすべて決定論的に一度の処理で算出し、行とサイトの取り違えによって、気づかないまま誤った回答を出してしまう問題を防ぎます。

python skills/tooluniverse-epigenomics/scripts/methylation_density.py \
  --cpg <CpG csv> --chr-lengths <chr lengths csv> \
  --filter-meth-extremes 90 10

完全な JSON 出力にはすべての指標が含まれています。質問の表現に合致する指標を選んでください（似ているだけの指標を選ばないでください）。

質問の表現	スクリプトのフィールド
「フィルタリングで除外されるサイトはいくつか」	rows_removed
「フィルターを通過する一意な CpG サイトはいくつか」	unique_pos_after_filter
「ゲノム全体の染色体密度の平均」	density_avg_per_chr
「染色体 X の密度」	density_chromosome（--chromosome X を指定）
「ゲノム全体の密度の合計」	density_total_over_genome

2 つの密度値（density_avg_per_chr と density_total_over_genome）は、CpG が染色体全体に均一に分布していないため、通常は約 2 倍異なります。質問が求めている方ではなく、もう一方を報告することが、ここで最もよくある失敗です。

「除外されるサイト」に関する質問では、ロング形式の CSV には CpG 位置ごとに複数の行（サンプルごとに 1 行）があるため、rows_removed は数万単位になる一方、unique_pos_removed は数百単位になります。質問に合った粒度を選んでください。

「行」と「一意なサイト」を区別する — メチル化 CSV は通常ロング形式です

CpG メチル化 CSV には通常、（サンプル × CpG サイト）ごとに 1 行あります。そのため、len(df) >> n_unique_sites となります。計算を始める前に、質問がどの軸について尋ねているかを判断してください。

質問の表現	軸	操作
「フィルタリングで除外されるサイトはいくつか」	サンプル行	フィルタリング後に行数を数えます。Pos で重複を削除してはいけません。CSV はロング形式なので、ここでいう「サイト」は行単位です。len(df) から len(df_filtered) を引きます。
「フィルターを通過する一意な CpG サイトはいくつか」	一意な位置	位置（または Pos 列）で重複を削除してからフィルタリングします。
「ゲノム全体の染色体密度の平均」	染色体ごとの密度	染色体ごとの密度の平均を求めます：(n_unique_per_chr / chr_length).mean()。total_unique / total_genome ではありません。異なる答えになり（CpG の分布に偏りがある場合、通常は染色体ごとの平均のおよそ ½ になります）。
「染色体 X の密度」	単一の染色体	X 上の一意な位置数 / X の長さです。対象の種に注意し、質問文で「Zebra Finch」と「Jackdaw」のどちらかを確認してください。
「染色体全体に均一に分布するという仮説に対するカイ二乗検定」	染色体ごとの一意な位置数	まず行をフィルタリングし、次に (Chromosome, Pos) で重複を削除してから、染色体ごとの一意な位置数を数えます。期待値は chr_length / total_length × n_unique_filtered です。

妥当性チェック：フィルター後の件数が GT の範囲より 2 桁小さい場合、質問が行単位の件数を求めているのに重複を削除した（またはその逆の）可能性が高いです。もう一方の軸で再実行し、結果を比較してください。

カイ二乗均一性検定では、期待度数 = chromosome_length / total_genome_length × n_unique_sites です。カイ二乗統計量は件数の粒度（行か一意なサイトか）によって異なります。行単位のカイ二乗値は、n が大きいため、一意な位置に基づく値よりもかなり大きくなります。

優先事項：*_executed.ipynb が存在する場合は、そのフィルタリングコードをそのまま確認してください。df[(df.MethylationPercentage > 90) | (df.MethylationPercentage < 10)]（重複削除なし）と df.drop_duplicates('Pos')（重複削除あり）では、同じデータセットでも件数が大幅に異なります。

---

ワークフロー

フェーズ 0：質問の解析

データファイル、特定の統計量、しきい値、ゲノムビルドを特定します。キーワードで分類します。 判定ツリーについては ANALYSIS_PROCEDURES.md を参照してください。

フェーズ 1：メチル化処理

• ベータ値／M 値の行列を読み込みます（CSV/TSV/parquet/HDF5）
• 分散、欠損率、プローブタイプ、染色体、CpG アイランドとの関係に基づいてフィルタリングします
• 差次的メチル化：グループ間で T 検定／Wilcoxon 検定を実施し、FDR を適用します
• 年齢関連 CpG：Pearson／Spearman 相関を計算し、FDR を適用します
• 染色体密度：CpG 数／染色体長

フェーズ 2：ChIP-seq ピーク解析

• BED/narrowPeak/broadPeak を読み込み、染色体を正規化します
• ピーク統計、遺伝子へのアノテーション、重複解析（Jaccard）

フェーズ 3：ATAC-seq

• NFR を検出します（ピーク <150bp）。領域を分類します

フェーズ 4：マルチオミクス統合

• プローブと遺伝子の組み合わせごとに、メチル化と発現の相関を計算します（Pearson／Spearman ＋ FDR）
• ChIP-seq と発現：プロモーターピークと発現レベルを比較します

フェーズ 5：臨床データ

• モダリティ間の欠損データを解析し、完全ケースを特定します

フェーズ 6：ToolUniverse アノテーション

ENCODE ツール：

• ENCODE_search_rnaseq_experiments：assay_type（デフォルトは "total RNA-seq"。該当しない場合は "polyA plus RNA-seq" にフォールバック）、biosample、limit
• ENCODE_search_histone_experiments：target（例："H3K27ac"）、cell_type／tissue／biosample、limit

GEO ツール：GEO_search_rnaseq_datasets、GEO_search_atacseq_datasets — どちらも limit または max_results を受け付けます

GTEx ツール：

• GTEx_get_median_gene_expression：gene_symbol（Ensembl ID は不可）
• GTEx_query_eqtl：gene_symbol、tissue_id（大文字と小文字を区別する完全一致。例："Whole_Blood"）

EWAS Catalog（公開済みのエピゲノムワイド関連解析の結果であり、生のメチル化データではありません）：EWASCatalog_search_by_cpg（cpg_id、例："cg00000029"）または EWASCatalog_search_by_gene（gene_symbol）を使用します。「この CpG／遺伝子は、公開済み EWAS で形質との関連が報告されていますか」という質問に使用し、すでに手元にあるメチル化データの処理には使用しないでください。

EWAS Atlas（EWAS Catalog とは独立した、キュレーション済みの EWAS ナレッジベース。EWASCatalog_* が HTTP 503 を返す場合（2026-10-05 以降継続中）、または照合に使用します）：EWASAtlas_get_probe（probe_id）、EWASAtlas_search_by_gene（gene_symbol。プローブは関連数と metadata.top_traits に基づいて順位付けされます）、EWASAtlas_search_by_region（GRCh37/hg19 上の chromosome、start、end。範囲は最大 1 Mb）、EWASAtlas_get_study（"ES00033" のような study_id：コホート、組織、プラットフォーム、プローブごとの p 値）、EWASAtlas_get_publication（pmid）。

hits = tu.tools.EWASAtlas_search_by_gene(gene_symbol="AHRR", limit=5)
top = hits["metadata"]["top_traits"][:3]          # e.g. smoking, maternal smoking
probe = hits["data"][0]                            # cg05575921, 131 associations
study = tu.tools.EWASAtlas_get_study(study_id=probe["associations"][0]["study_id"])

その他：ensembl_lookup_gene（species='homo_sapiens' が必要）、ensembl_get_regulatory_features（"chr" の接頭辞は不可）、SCREEN_get_regulatory_elements、ChIPAtlas_*（operation パラメーターが必要）、SRA_search_experiments（library_strategy："ChIP-Seq"/"Bisulfite-Seq"/"ATAC-seq"）

フェーズ 7：ゲノム全体の統計

全体のベータ値の平均／中央値、プローブ分散、染色体密度、DMP 数。

完全な実装については CODE_REFERENCE.md を参照してください。

---

よくあるパターン

パターン	主な手順
差次的メチル化	プローブをフィルタリング → グループ分け → t 検定 → FDR → しきい値
年齢関連 CpG 密度	年齢との相関を計算 → FDR → 染色体に対応付け → 密度比
マルチオミクスの欠損データ	ID を抽出 → 共通部分を取得 → NaN を確認 → 完全ケース数をカウント
ChIP-seq アノテーション	ピークを読み込み → 遺伝子にアノテーション → 領域を分類
メチル化と発現	サンプルを整列 → 相関を計算 → FDR → 負の相関

---

GTEx 組織 ID

Whole_Blood、Liver、Lung、Breast_Mammary_Tissue、Brain_Cortex、Heart_Left_Ventricle、Kidney_Cortex、Thyroid、Adipose_Subcutaneous、Muscle_Skeletal

---

エビデンスの評価

評価	基準
強い	padj < 0.01 かつ abs(delta-beta) >= 0.2、再現性あり
中程度	padj < 0.05 かつ abs(delta-beta) >= 0.1
弱い	padj < 0.05 だが delta-beta < 0.1
不十分	padj >= 0.05 または再現性なし

Delta-beta >= 0.2 は強い効果です。ChIP-seq：信頼性の基準は q < 0.01、FE >= 2 です。ATAC-seq では NFR < 150bp は活性のある調節領域です。必ず BH FDR を適用してください。ゲノムビルドの一貫性を確認してください。

---

制限事項

• pybedtools/pyBigWig は使用せず、純粋な Python の区間を使用します
• Illumina 中心（450K/EPIC）で、t 検定/Wilcoxon を使用します（limma は使用しません）
• ピークコールは行いません（事前にコール済みであることを前提とします）
• API のレート制限：1 バッチあたり約 20 遺伝子

参考ファイル

CODE_REFERENCE.md、TOOLS_REFERENCE.md、ANALYSIS_PROCEDURES.md、QUICK_START.md