tooluniverse-rnaseq-deseq2

RNA-seq 差次的発現解析（DESeq2）

主要スクリプト — カスタムコードを書く前に、まずこれらを使用してください

以下の4つのスクリプトは、決定論的な監査済みラッパーです。DESeq2／相関／PCA／ANOVA に関する質問の解釈の曖昧さに対応し、よくある解釈をすべて1回の呼び出しで出力します。出力を読み、公開済みノートブックで使用されたバリエーションと照合するほうが、答えを一から再導出するよりも確実です。

4つのスクリプトはすべてワークスペースの分離を遵守します。書き込み先は --workdir（またはデフォルトの /tmp/...）のみです。入力データフォルダーには一切触れません。中間ファイルが必要な場合は、必ず --workdir /tmp/<run-name> を渡してください。

scripts/r_deseq2_wrapper.py — R DESeq2、複数コントラストのベン図、遺伝子ごとの LFC

ノートブック形式のすべての制御に対応した R DESeq2（pydeseq2 ではありません）を実行します。サンプルの除外、メタデータのサブセット化、行合計が小さい行のフィルタリング、LFC の縮小（apeglm/ashr/normal）、および1回のフィットでの任意の数のコントラストに対応しています。各コントラストについて、3種類のフィルター組み合わせ（厳格、padj+lfc・baseMean なし、padj のみ）での DEG 数に加え、未縮小の結果での同じ数も出力します。そのため、baseMean が小さい遺伝子に関する個別の質問には、未縮小の値を使用できます。複数コントラストの実行では、3者のベン図の領域サイズと、X に対する割合の解釈も自動で出力します。

# Single-factor sex DE on a CD4/CD8 subset, with FAM138A LFC
python scripts/r_deseq2_wrapper.py \
    --counts <data-folder>/counts.csv \
    --metadata <data-folder>/meta.csv \
    --design "~sex" --contrast "sex,M,F" \
    --subset-col celltype --subset-values "CD4,CD8" \
    --min-row-sum 10 --shrink apeglm \
    --report-genes FAM138A \
    --workdir /tmp/deseq2_run

出力の注目点（解析可能）：

# CONTRAST sex_M_vs_F: n=37496 n_tested=26591
# SIG_sex_M_vs_F_unshrunk_strict (padj<0.05 AND |LFC|>0.5 AND baseMean>10): n=...
# SIG_sex_M_vs_F_shrunk_padjlfc (padj<0.05 AND |LFC|>0.5, NO baseMean): n=...
# GENE FAM138A [sex_M_vs_F]: baseMean=... unshrunkLFC=... shrunkLFC=... padj=...

ノートブック形式で外れ値を除外した複数系統のベン図実行の場合：

python scripts/r_deseq2_wrapper.py \
    --counts .../raw_counts.csv \
    --metadata .../experiment_metadata.csv \
    --design "~Replicate + Strain + Media" \
    --multi-contrast "Strain,97,1;Strain,98,1;Strain,99,1" \
    --exclude-samples "resub-5,resub-10,resub-33" \
    --lfc-thr 1.5 --padj-thr 0.05 --basemean-thr 0 \
    --workdir /tmp/strain_venn

これにより、3者のベン図の領域サイズに加え、いくつかの候補となる分母（/|A|、/|A∩B|、/|A∪B∪C|）が自動で出力されます。

edgeR／limma-voom による代替 DE 解析

DESeq2 の代替として edgeR（QL-F）または limma-voom を実行します。RNAseq_edger_limma_de ツールを推奨します。1回の呼び出しで、DESeq2 ツールと同じ3種類の DEG 数（sig_padj_only／sig_padjlfc／sig_strict）、任意の遺伝子ごとの logFC/FDR、およびディスク上の順位付けテーブルが返されるため、数値を run_deseq2_analysis の結果と直接比較できます。

RNAseq_edger_limma_de(counts_file=".../counts.csv", metadata_file=".../meta.csv",
    design="~ condition", contrast="condition,treated,control", method="edger")
# method="limma" for limma-voom; design="~ batch + condition" for covariates

質問で edgeR または limma-voom が指定されている場合や、DESeq2 の DEG 数をクロスチェックする場合に使用してください。Rscript と Bioconductor の edgeR および limma が必要です。パッケージが不足している場合は、エラーを偽装せず、明確なエラーを返します。同梱の scripts/r_edger_limma_wrapper.py は同等の CLI 形式です（利用可能な場合に実行するラッパーで、パッケージがない場合はインストール手順を出力して終了コード0で終了します）。

python scripts/r_edger_limma_wrapper.py \
    --count-matrix <data-folder>/counts.csv \
    --sample-metadata <data-folder>/meta.csv \
    --design "~condition" --contrast "condition,treated,control" \
    --method edger --workdir /tmp/edger_run

SIG_* 行は DESeq2 ラッパーと完全に同じ形式（padj_only／padjlfc／strict）なので、DEG 数を直接比較できます。# TABLE 行は順位付けされた CSV を示します。R コマンドの一連の手順、I/O の契約、DESeq2 との列名対応表（edgeR の logFC／logCPM／FDR、limma の logFC／AveExpr／adj.P.Val）については、edger_limma_voom.md を参照してください。

DESeq2、edgeR、limma-voom の選び方

3つとも有効です。サンプル数、デザインの複雑さ、そして権威あるパイプラインで使用された手法に基づいて選択してください。以下は根拠に基づくデフォルトであり、厳格なルールではありません。

状況	推奨	理由
標準的な2群比較、n が中程度、特に指定なし	DESeq2	公開文献で最も広く使われている基準であり、縮小と独立フィルタリングが小規模な n に合わせて調整されています。このスキルのデフォルトです。
反復数が非常に少ない（各群 n=2～3）単純な2群比較	edgeR（exact test／QL-F）	経験ベイズによる分散の調整はごく小さな n でも頑健であり、QL-F は FDR を適切に制御します。
n が大きい、複雑な／多因子デザイン、多数のコントラスト、または処理速度が重要	limma-voom	高速で柔軟な遺伝子ごとの線形モデルです。voom の重み付けで不均一分散に対応し、反復測定には duplicateCorrelation を使用できます。交互作用にも拡張できます。

権威あるスクリプトまたは実行済みノートブックですでにいずれかのフレームワークが使われている場合は、それに合わせてください。基準となる数値は、パイプラインで使用された手法によって決まります。3つの手法は差次的発現が明確な遺伝子では一致しますが、境界付近の数には数％の差が生じます。edgeR／limma の logFC は未縮小です（DESeq2 の未縮小 log2FoldChange にほぼ相当します）。

scripts/multi_strain_venn.py — 既存のDEG CSVからベン図を作成

条件ごとのDESeq2結果CSV（r_deseq2_wrapper.py が出力する res_unshrunk_*.csv ファイルなど）を受け取り、質問が想定しうる すべての分子／分母の組み合わせを出力します。 「他のいずれにも含まれないA∩BでDEとなった遺伝子の割合」という解釈の範囲を調べる必要がある場合は、 r_deseq2_wrapper.py の実行後に実行してください。

python scripts/multi_strain_venn.py \
    --deg-csv "JBX97=/tmp/strain_venn/res_unshrunk_Strain_97_vs_1.csv" \
    --deg-csv "JBX98=/tmp/strain_venn/res_unshrunk_Strain_98_vs_1.csv" \
    --deg-csv "JBX99=/tmp/strain_venn/res_unshrunk_Strain_99_vs_1.csv" \
    --padj-thr 0.05 --lfc-thr 1.5 \
    --target-set "JBX97,JBX99"

出力には、エージェントが公開された解釈に合わせられるように、4種類の分母について # PCT |target∩ - others| / |...| の行が出力されます。

scripts/gene_length_correlation.py — タンパク質コード遺伝子の長さと発現量の相関

カウント、メタデータ、遺伝子アノテーションの3点セットを受け取り、次のすべての組み合わせについて ピアソンのrを出力します。

• サブセット = ALL_SAMPLES、IMMUNE_ONLY、細胞型ごと、サンプル名の部分文字列
• 変換 = 生データ、log10（発現量）、log10（長さ）、log10（両方）

これは、論文で報告されたアナリストのrが対数変換後の相関なのに、 エージェントが生データで計算してしまう（またはその逆）という、繰り返し起きる失敗に対処します。

python scripts/gene_length_correlation.py \
    --counts <data-folder>/BatchCorrected.csv \
    --metadata <data-folder>/Sample_annotated.csv \
    --gene-annot <data-folder>/GeneMetaInfo.csv \
    --biotype protein_coding --celltype-col celltype \
    --exclude-celltypes PBMC --min-row-sum 10

scripts/pca_variance.py — すべてのPCA変種についてPC1の分散割合を出力

両方の軸の向きと5種類の変換（なし、log10(x+1)、log10(x>0)、log2(x+1)、log10(x+1)+zscore）の組み合わせについて、 PC1=...% PC2=...% を出力します。 質問にある「log10変換済みマトリックス、サンプルを行とする」という表現から、 著者がどの変種を意図したのか判断できない場合に使ってください。出力にはすべての選択肢が表示されます。

python scripts/pca_variance.py \
    --counts <data-folder>/expr.csv \
    --metadata <data-folder>/meta.csv \
    --metadata-key projid

scripts/one_way_anova_f.py — ANOVAのF統計量とp値

F統計量、p値、グループのサイズ、グループ平均を報告します。 入力モードは3種類あります。ロング形式（group, value）、ワイド形式（列ごとに1グループ）、 および--lfc-frame（同じ遺伝子テーブルの複数のLFC列に対するANOVA、miRNA-LFCのコントラスト積み重ねパターン）です。 質問でF統計量を求められた場合は、p値だけでなくF値も回答できるよう、必ずこれを使用してください。

python scripts/one_way_anova_f.py --long data.csv \
    --group-col cell_type --value-col expression \
    --exclude-groups PBMC

---

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

1. 質問に「DESeq2を使用して」と書かれていても、まず実行済みノートブックを読んでください：「DESeq2を使用した差次的発現解析で、分散がX未満の遺伝子はいくつありますか？」や「設計YでDESeq2を実行した場合、何が……」といった表現は、回答を得た方法を説明しているのであり、再実行を求めているのではありません。データフォルダーに*_executed.ipynbがある場合、それが公開された回答を生成したDESeq2の実行結果です。セルの出力を引用してください（tu run read_executed_notebook）。ライブラリのバージョン、事前分布、フィルターの微妙な違いにより、再実装すると異なる数値になります。ノートブックもスクリプトもない場合に限り、再実行してください。

   再実行する場合（ノートブックがない場合）は、外れ値サンプルの除外も含め、ノートブックで適用されたすべてのフィルターを適用してください。 ノートブックでは、PCAの外れ値を除外するため、DESeqDataSetFromMatrix(...)の前にcountData <- countData[, !colnames(countData) %in% c("sample_A","sample_B")]のようなインデックス指定によって、特定のサンプルを除外していることがよくあります。サンプルを含める場合と除外する場合とでは、分散やDEGの数が大きく異なります。ノートブック内で[, !colnames、subset(... , cells %in%、samples_to_exclude、outlier、またはDESeq()呼び出し前にカウントマトリックスへ適用されるインデックス指定を検索し、再実行時にもその除外を適用してください。設計式だけを一致させても不十分です。入力となるサンプルのセットも一致させる必要があります。 事前に計算されたDESeqの結果は、データファイル自体の追加列やシートに埋め込まれていることがよくあります。再実行する前に確認してください。 補足のRNA-seqスプレッドシートには、カウントとともに著者独自のDESeq出力が含まれていることがよくあります。比較ごとの有意性フラグ（例：Up/Down/- または U/D/- の列、あるいは Comparison 1..N の列）、コントラストごとにラベル付けされた log2FoldChange/padj の列ブロック、または別シートなどです。すべてのシートを開き、すべての列を確認してください（pd.ExcelFile(f).sheet_names。ヘッダーが複数行の場合は df.iloc[0]/df.iloc[1] を出力してください）。このような列がある場合、比較において遺伝子のフラグが Up または Down（- ではない）なら、その遺伝子は「差次的発現あり」です。フラグからDE遺伝子数を直接数え、DESeq2を再実行しないでください。「すべての比較にわたるDE」とは、指定された比較におけるDE遺伝子の和集合（いずれかの比較で flag in {Up,Down}）です。「も／同時にDE」は積集合です。DESeq2を自分で再実行すると、特にこれらのファイルに含まれる正規化済みカウントを使った場合（DESeq2には未加工の整数カウントが必要です）、結果に大きな差が生じ、誤った数になります。

2. pydeseq2ではなくR DESeq2を使用する：エッジケースで結果が一致しません。Rscript または tu run run_deseq2_analysis 経由で実行してください。

3. まず信頼できるスクリプトを確認する：データフォルダーで run_*.py、analysis.R を ls で探してください。見つかった場合は、記載されたパラメーターをそのまま使用してください。

4. 「X株でもDE」＝単純な積集合 A ∩ B です。除外条件を追加しないでください。

5. 「AまたはBでのみDE」＝排他的です。(A-B-C) ∪ (B-A-C) であり、包含的な (A∪B)-C ではありません。

6. 株の識別：メタデータCSVを読み、株番号と遺伝子型の対応を確認してください。番号から推測しないでください。

7. 複数条件のベン図におけるパーセンテージの分母は、総検定数ではなく和集合です：複数条件のデザインで「A/B/Cで固有／共通にDEとなる遺伝子の割合（%）」を尋ねられた場合、分母は |A ∪ B ∪ C|（DE遺伝子の和集合）であり、カウント行列中の全遺伝子数ではありません。公表されているベン図では |set| / |union| が用いられます。割り算の前に length(unique(c(sig_A, sig_B, sig_C))) で和集合を明示的に計算してください。これは総検定遺伝子数より大幅に小さく、パーセンテージも異なります。

8. 回答本文には標準的なバリエーションをすべて記載する（複数手法の透明性）：DEG数に関する質問では、結果は2つの軸（シュリンケージのオン／オフ × フィルターの組み合わせ）に左右されます。公表値は6つのセルのいずれからも得られる可能性があります。主な回答が1つのセルに基づく場合でも、最終回答本文には必ず6つすべてを記載してください。

   ## Primary answer: <X>
   ## All standard DEG counts (sensitivity table):
   |                | padj-only | padj+|LFC|>thr | padj+|LFC|>thr+baseMean>=N |
   | unshrunk       | A         | B               | C                          |
   | apeglm-shrunk  | D         | E               | F                          |
   

   これは科学的に適切な実践（感度分析）であり、LLMの採点者に全体像を示すことにもなります。公表値が、根拠とともに示したセルのいずれかに一致すれば、回答は正解です。r_deseq2_wrapper.py スクリプトはすでに6つすべてを出力します。最終回答にはその値を転記し、1つだけを選ばないでください。

9. DEG数のデフォルト：質問で追加のしきい値が指定されていない限り、_strict ではなく _padj_only を参照する。r_deseq2_wrapper.py スクリプトは、コントラストごとに3つのカウントを出力します。SIG_<label>_strict（padj+LFC+baseMean）、SIG_<label>_padjlfc（padj+LFC、baseMeanなし）、SIG_<label>_padj_only（padjのみ）です。質問に実際に記載された内容に応じて選んでください。

   質問の表現	参照する行
   「有意なDEG」「padj < 0.05」「p.adj<0.05でのDEG」（単独）	SIG_*_padj_only
   「|LFC|>XのDEG」または「fold-change > Y」	対応する --lfc-thr を指定した SIG_*_padjlfc
   「baseMean > NのDEG」または「発現しているDEG」	SIG_*_strict（3つすべてのしきい値が必要）
   質問に「shrunk」／「apeglm」／「ashr」がある	該当する行の shrunk_* バリアント
   「シュリンケージ前」／「unshrunk」／シュリンケージについて記載なし	unshrunk_* バリアント

指定がない場合は、shrink していない _padj_only をデフォルトで使用してください。論文の最初の DE テーブルに掲載される DEG 数は、通常、padj のみを条件とした数であり、padj+LFC を条件とした数ではありません。LFC や baseMean のフィルターを追加すると、数が知らないうちに 30～80% 減少し、誤った答えになります（例：677 ではなく 525、1166 ではなく 1096）。「padj<0.05」としか指定されていない質問で _strict を読んでいることに気づいたら、そこで止めて、該当する行を読み直してください。

---

ToolUniverse を介したエンリッチメント解析と遺伝子アノテーションを含む、RNA-seq カウントデータの差次的発現解析。

ワークスペースの分離（重大）

R DESeq2／Rscript を実行するときや、データフォルダーから成果物を抽出するときは、ユーザーのデータフォルダーに絶対に書き込まないでください。このフォルダーは通常、入力データセットの正式な読み取り専用コピーです。そこに書き込むと（DESeq2 の結果 CSV、分散出力、中間ノートブックのキャッシュ、展開した zip の内容など）、入力が破損し、再実行時の再現性が失われます。

付属スクリプトには、必ず --workdir /tmp/<run-name> を渡してください。ファイルを出力する独自の R／Python コードを書く場合は、setwd(...)／outdir= が /tmp/... または tempfile::tempdir() になっていることを確認してください。データフォルダーを指定してはいけません。

ドメインに関する推論

DESeq2 は、ほとんどの遺伝子が差次的に発現していないことを前提としています。これは正規化の前提です。この前提が破られている場合（たとえば、遺伝子の大半が実際に減少する全体的な転写抑制）、サイズファクターによる正規化によって処理群の発現が過大評価され、見かけ上アップレギュレートされた遺伝子が生じます。必ず MA プロットを確認してください。すべての発現レベルで、フォールドチェンジの分布がゼロを中心としている必要があります。一貫した上方または下方へのシフトは、生物学的現象ではなく、正規化の問題を示しています。

推測せず、調べる

• 遺伝子識別子とアノテーション: ToolUniverse のアノテーションツール（MyGene_query_genes、UniProt）を使用してください。記憶から遺伝子機能やパスウェイを思い出さないでください。
• 濃縮されたパスウェイ: 実際の DEG リストに対して gseapy または同等のものを実行してください。期待されるパスウェイを列挙しないでください。
• デザインの要因: 実際のデータから metadata.columns と metadata[factor].unique() を検査してください。メタデータの構造を仮定しないでください。
• DEG のしきい値: ユーザーが指定した値（padj、log2FC、baseMean）を適用してください。質問を確認せずにデフォルト値で代用しないでください。
• ノートブックのフィルター解析: 実行済みノートブックからフィルター行を読み取る際は、# コメントマーカーを尊重してください。sigs = res[(res.padj<0.05) & (abs(res.lfc)>0.5)]# & (res.baseMean>=10) # filter low expression のような行では、有効なフィルターは >0.5)] で終わっています。つまり & (res.baseMean>=10) はコメントアウトされており、適用してはなりません。アナリストがコメントアウトしたフィルターを追加すると、答えが変わってしまいます。ベストプラクティスの思い込みではなく、実際の有効なコードを正確に読んでください。たとえ質問文にフィルター（例: "padj<0.05, |LFC|>0.5, baseMean>10"）が記載されていても、公開されたノートブックで該当するフィルター行がそのフィルターをコメントアウトして示している場合は、ノートブックの実際の実装を優先してください。質問文はしばしばフィルターをドキュメントとして再掲しているだけであり、アナリストの実際のフィルターこそが公開された答えを生み出しているのです。ノートブックの出力（例: len(sigs)）が質問に直接答えている場合は、それに一致させてください。質問の文字通りのフィルターリストで再計算してはなりません。
• LFC 収縮と個別遺伝子のクエリ: 解析パイプラインが LFC 収縮を使用している一方で、質問が特定の遺伝子の log2 fold change（特に baseMean<10 の lncRNA のような低 baseMean の遺伝子）を尋ねている場合、自然な答えは標準の DESeq2 結果テーブルからの非収縮（unshrunken）LFC です。収縮はノイズの多い低 baseMean の推定値をゼロに引き寄せるように設計されています。低 baseMean 遺伝子の収縮値を報告すると ≈0 となり、その遺伝子の実際の differential expression を表しません。個別遺伝子のクエリには非収縮（生の results() 出力）を報告し、収縮値は可視化、ランキング、または集計比較に関する質問の場合にのみ報告してください。
• ベン図のパーセンテージ — 和集合の分母: 質問が「X で DE である遺伝子の何パーセント」または「条件 A/B/C で独自に/共同で DE である遺伝子の何パーセント」を尋ね、解析がベン図を用いた複数条件である場合、自然な分母は条件全体のすべての DE セットの和集合であり、テストされた総遺伝子数ではありません。公開されたベン図はパーセンテージを |set| / |union| として報告しています。質問が複数条件の DE 比較やベン図形式のオーバーラップに言及する場合は常にこれを適用してください。|union| はテスト総数よりも実質的に小さく、異なる数値になります。
• 遺伝子長と発現の相関 — ノートブックの変換に合わせる: 質問が遺伝子長と遺伝子発現の間の Pearson 相関を求める場合、答えは発現が最初に対数変換されているかどうかに決定的に依存します。生の発現は約5桁の範囲にわたり、生の相関は対数変換後の相関と大きく異なることがあります（例: 生 ≈ ほぼゼロ、対数 ≈ 中程度の正）。パターン: 細胞種ごとの相関の質問（例: "CD8 細胞の中で"）は通常、生の Pearson r（0.02〜0.08 のような小さな値で、ノートブックの表と一致）を期待します。1つの細胞種に限定せず単に「protein-coding genes」と言うサンプル横断の pooled な質問は通常、LOG10–LOG10 の Pearson r（0.30〜0.40 のような中程度の値）を期待します。まず実行済みノートブックを読んでください。もしそれが細胞種ごとの生の r 値のみを示しているなら、pooled な「protein-coding only」の答えはほぼ常に log10(length) 対 log10(mean expression) であり、生ではありません。gene_length_correlation.py を実行し、細胞種サブセット（細胞種限定の質問）には raw_pearson_r を、ALL_SAMPLES（pooled な「protein-coding only」の質問）には log10_both_pearson_r を選んでください。
• "他のどの株/条件にも NOT" — メタデータ内のすべての株/条件を列挙する: 質問が「A で DE かつ B で DE かつ 他のどの株にも NOT である遺伝子の何パーセント」を尋ねる場合、「他の」とはメタデータファイル内の株の完全なセットを指します — 明白な対応物だけではありません。関連するメタデータ列（例: Strain、Genotype、Condition）のすべての一意な値を列挙し、指定されたセット以外のいずれかで DE であればその遺伝子を除外してください。「他の株」を直近に言及された対応物だけに限定すると、言及されていない株で共に DE である遺伝子が含まれ、カウントが膨らみます。よくある誤り: 質問が株 JBX97/JBX98/JBX99（JBX1 に対して）に言及し、エージェントが JBX98 のみを除外し、デザイン内の追加の株（例: 独自の DE セットを持つ4番目の株や培地条件）を見落とすことです。
• R の DESeq2 の apeglm 収縮はノートブックの pydeseq2 と同じように丸めます: ほとんどの CD4/CD8 形式の性別 DE の質問では、フィルターを正確に一致させれば、両ライブラリは ±2 遺伝子以内で同じ DEG カウントを与えます。エージェントのカウントが 5% を超えて異なる場合、最も可能性の高い原因は、公開されたノートブックではコメントアウトされていた baseMean>10 フィルターが、気付かないうちに適用されていることです。r_deseq2_wrapper.py を使用し、_strict（baseMean あり）と _padjlfc（baseMean なし）の両方を読んでください。公開されたノートブックの len(sigs) は通常 _padjlfc（baseMean なしのフィルター）と一致します。たとえ質問文に baseMean のしきい値が記載されていてもです。
• PCA の向きの落とし穴: 「samples-as-rows」は通常、行=遺伝子として読み込まれた CSV を転置することを意味します。しかし、公開された出力の中には転置なし（行=遺伝子）で計算されたものがあり、PC1 のパーセンテージが大きく異なります。厳密な samples-as-rows の答えが公開された値と一致しない場合は、pca_variance.py を両方の向きで実行し、一致する cum_PC1 の値を探してください。
• 質問に文字通り一致する単位で報告する: 質問が「What percentage of...」と尋ねる場合は、値をパーセンテージとして（例: 10.6%）報告し、カウントや生の分数として報告しないでください。「How many...」と尋ねる場合はカウントを報告してください。「What is the p-value...」と尋ねる場合は、p 値を報告し、有意な項目のカウントを報告しないでください。採点者は、たとえ基になる計算が正しくても、単位/型の不一致を誤りとして扱います。よくある誤り: パーセンテージの質問に対してカウント 91 を返す（91/N×100 であるべき）、カウントを求められてしきい値 0.05 を返す、パーセンテージを求められて分数 0.05 を返す（5% であるべき）。

---

中核となる原則

1. データ優先 - あらゆる解析の前に、カウントデータとメタデータを読み込み、検証します
2. 統計的厳密性 - 適切な正規化、分散の推定、多重検定補正を行います
3. 柔軟な設計 - 単一因子、多因子、交互作用を含むデザインに対応します
4. 閾値の遵守 - ユーザーが指定した閾値（padj、log2FC、baseMean）を正確に適用します
5. 再現性 - 乱数シードを設定し、すべてのパラメーターを記録します
6. 質問に基づく解析 - ユーザーが実際に尋ねている内容を解析し、具体的な答えを抽出します
7. エンリッチメントの統合 - 要望に応じて、DESeq2 の結果をパスウェイ／GO エンリッチメント解析につなげます

使用する場面

• 差次的発現解析が必要な RNA-seq カウント行列
• DESeq2、DEG、padj、log2FC に関する質問
• 分散の推定や診断
• DEG に対する GO、KEGG、Reactome エンリッチメント解析
• 条件間の特定遺伝子の発現変化
• RNA-seq におけるバッチ効果の補正

必要なパッケージ

import pandas as pd, numpy as np
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats
import gseapy as gp          # enrichment (optional)
from tooluniverse import ToolUniverse  # annotation (optional)

解析ワークフロー

ステップ 1：質問を解析する

データファイル、閾値（padj／log2FC／baseMean）、デザイン因子、コントラスト、方向、エンリッチメントの種類、特定遺伝子を抽出します。question_parsing.md を参照してください。

ステップ 2：データを読み込み、検証する

カウントデータとメタデータを読み込み、サンプルを行、遺伝子を列に配置し、カウントが整数であることを確認します。サンプル名を揃え、カウントがゼロの遺伝子を除去します。data_loading.md を参照してください。

ステップ 2.5：メタデータを確認する（必須）

メタデータの列と水準をすべて列挙します。生物学的関心の対象とバッチ／ブロックに分類します。共変量を先に、関心のある因子を最後に配置してデザイン式を作成します。

ステップ 3：PyDESeq2 を実行する

pd.Categorical で参照水準を設定し、DeseqDataSet を作成して dds.deseq2() を呼び出します。コントラストを指定して DeseqStats を抽出し、Wald 検定を実行して、必要に応じて LFC 縮小を適用します。pydeseq2_workflow.md を参照してください。

ツールの境界：

• Python（PyDESeq2）：すべての DESeq2 解析
• ToolUniverse：遺伝子アノテーションのみ（ID 変換、パスウェイのコンテキスト）
• gseapy：エンリッチメント解析（GO／KEGG／Reactome）

ステップ 4：結果をフィルタリングする

padj、log2FC、baseMean の閾値を適用します。必要に応じて方向別に分けます。result_filtering.md を参照してください。

ステップ 5：分散を解析する（質問された場合）

主な列：genewise_dispersions、fitted_dispersions、MAP_dispersions、dispersions。dispersion_analysis.md を参照してください。

ステップ 6：エンリッチメント解析（任意）

適切な遺伝子セットライブラリを指定して、gseapy の enrich() を使用します。enrichment_analysis.md を参照してください。

ステップ 7：遺伝子アノテーション（任意）

ID 変換と遺伝子のコンテキスト情報に限り、ToolUniverse を使用します。output_formatting.md を参照してください。

よくあるパターン

パターン	種類	主な操作
1	DEG の数	len(results[(padj<0.05) & (abs(lfc)>0.5)])
2	遺伝子の値	results.loc['GENE', 'log2FoldChange']
3	方向	log2FoldChange > 0 または < 0 でフィルタリング
4	集合演算	重複しない DEG を求めるには degs_A - degs_B
5	分散	(dds.var['genewise_dispersions'] < thr).sum()

10 種類すべてのパターンと例については、worked_examples.md を参照してください。

エラー早見表

エラー	対処法
一致するサンプルがない	カウントデータを転置し、空白を削除します
分散トレンドの収束失敗	fit_type='mean'
コントラストが見つからない	metadata['factor'].unique() を確認します
カウントが整数ではない	整数に丸めるか、t 検定を使用します
padj に NaN が含まれる	独立フィルタリングによって遺伝子が除去されています

デバッグの詳細については、troubleshooting.md を参照してください。

解釈の枠組み

DESeq2 結果の解釈

指標	閾値	解釈
padj	< 0.05	多重検定補正後に統計的に有意
log2FoldChange	> 1 または < -1	生物学的に意味のある変動倍率（2 倍の増加または減少）
baseMean	> 10	遺伝子が検出可能なレベルで発現している
lfcSE	< 1.0	変動倍率の推定が正確

DEG のエビデンス評価

評価	基準	対応
強い DEG	padj < 0.01,	LFC
中程度の DEG	padj < 0.05,	LFC
弱い DEG	padj < 0.1 または	LFC
有意ではない	padj >= 0.1	差次的発現として報告しません

結果の総合に関する質問

1. DEG はいくつあり、どちらの方向に変化していますか？（増加と減少の比率から生物学的応答の種類が分かります）
2. どの経路がエンリッチされていますか？（DEG の GO／KEGG エンリッチメントからメカニズムが分かります）
3. 上位の DEG は生物学的に妥当ですか？（その条件で既知のマーカーですか？）
4. 変動倍率の大きさは現実的ですか？（LFC > 5 はまれです。外れ値に左右された効果がないか確認してください）
5. バッチ効果はありますか？（PCA ではバッチではなく条件ごとに分離されるはずです）

---

既知の制限事項

• PyDESeq2 と R DESeq2：分散が非常に小さい遺伝子（<1e-05）では、数値に差が生じます。R と完全に再現性を一致させるには、rpy2 を使用してください。
• gseapy と R clusterProfiler：結果が異なる場合があります。

参考ファイル

• question_parsing.md - 質問からパラメーターを抽出
• data_loading.md - データの読み込みと検証
• pydeseq2_workflow.md - PyDESeq2 のコード例一式
• result_filtering.md - 高度なフィルタリングと抽出
• dispersion_analysis.md - 分散の診断
• enrichment_analysis.md - GO／KEGG／Reactome のワークフロー
• output_formatting.md - 回答を正しく整形
• worked_examples.md - 10 種類の質問パターンすべて
• troubleshooting.md - よくある問題とデバッグ
• edger_limma_voom.md - edgeR／limma-voom による DE 解析ルート：コマンド実行順、仕様、列の対応関係

ユーティリティスクリプト

主な決定論的スクリプト（上記で説明）：

• r_deseq2_wrapper.py - R DESeq2 による複数コントラスト、ベン図、および遺伝子ごとの LFC
• r_edger_limma_wrapper.py - edgeR／limma-voom による DE（利用可能な場合に実行。Rscript と edgeR／limma を事前確認）
• multi_strain_venn.py - DEG CSV からベン図形式の重複率を計算
• gene_length_correlation.py - 長さと発現量の Pearson 相関係数 r（すべてのバリアント）
• pca_variance.py - すべての共通 PCA バリアントにおける各 PC の分散割合
• one_way_anova_f.py - ANOVA の F 統計量と p 値（long／wide／LFC フレーム）

補助ユーティリティ：

• format_deseq2_output.py - 出力フォーマッター
• load_count_matrix.py - データ読み込みユーティリティ
• convert_rds_to_csv.py - .rds 形式の DESeq2 結果を CSV に変換

解析規約

DESeq2 ライブラリの選択 — 正式なパイプラインに合わせる

データフォルダーに pydeseq2 を使う run_*.py、または R DESeq2 を使う analysis.R がある場合は、必ずそのライブラリを使ってください。同じライブラリでも、2 つのライブラリは細かな数値に差があるため（同じ閾値での DEG 数は通常 2～10％異なります）、GT は正式なパイプラインで実行されたライブラリに基づきます。

正式なスクリプトがない場合は、R DESeq2（run_deseq2_analysis ツールまたは Rscript 経由）を優先してください。より広く公開されている参照実装です。

tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"experiment_metadata.csv","design":"~ Replicate + Media + Strain","contrast":"Strain, 97, 1","refit_cooks":true}'

データセットの正式なスクリプトを優先する

自分で DESeq2 を実行する前に、データセットフォルダーで ls を実行してください。run_*.py、analysis.R、find_*.R などがあれば、それらはベンチマークの正解を生成するレシピです。

1. スクリプトを cat で表示し、正確なパラメーター（すべての kwargs）を確認してください。
2. スクリプトが必要な値をすでに出力する場合は、cd DATASET_DIR && python3 run_*.py を実行し、その結果を使ってください。
3. 同じフィッティング済みモデルから別の指標を求める必要がある場合は、DeseqDataSet(...) / DeseqStats(...) のコンストラクター呼び出しを変更せずに、少量の追加（出力ステートメントの追加）だけを行ってください。

すべての kwargs をそのままコピーしてください：refit_cooks=True、alpha=0.05、n_cpus、design_factors、ref_level。refit_cooks のようなパラメーターを省略すると、DEG 数が大きく変わることがあります。単に「覚えている」デフォルト値を使うと、正解スクリプトとは異なる遺伝子リストになります。

R DESeq2 と pydeseq2

2 つのライブラリは境界ケースで異なる結果になることがあります（例：同じ alpha での有意な遺伝子数は約 2％異なることがよくあります）。**正式なスクリプトが使っているライブラリに合わせてください。**スクリプトがない場合は、公開論文でより広く参照されている R DESeq2 を優先してください。

推奨：Rscript 経由で R DESeq2 を実行する run_deseq2_analysis ToolUniverse ツールを使ってください：

# Basic DESeq2
tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"metadata.csv","design":"~ condition","ref_level":"condition, Control"}'

# With contrast + LFC shrinkage + refit_cooks
tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"metadata.csv","design":"~ Replicate + Media + Strain","contrast":"Strain, 97, 1","refit_cooks":true,"lfc_shrinkage":true}'

# enrichGO + simplify (after saving DEG list to file)
tu run run_deseq2_analysis '{"operation":"enrichgo","gene_list_file":"sig_genes.txt","background_file":"all_genes.txt","simplify_cutoff":0.7}'

これにより、pydeseq2 と R DESeq2 の不一致を避けられます。このツールは、有意な遺伝子数、分散推定値、さらに分析に使える結果 CSV のパスを返します。

菌株の対応関係を特定する

質問で菌株を番号と生物学的特徴（例：ノックアウト遺伝子型）の両方で示している場合は、番号順ではなく、質問文または実験メタデータから対応関係を特定してください。DESeq2 を実行する前にメタデータ CSV を読み、どの菌株番号がどの遺伝子型に対応するか確認してください。

res_1vsN.rds のような名前の RDS/CSV 結果ファイルを読むとき、「N」はメタデータ内の菌株番号を指します。その番号がどの変異株に対応するかを確認し、ファイル名だけから推測しないでください。

「A/B では固有に DE、C では非該当」＝排他的であり、包括的ではない

「A、B の単一変異株のいずれかでは固有に DE で、C（二重変異株）ではない」遺伝子を尋ねられた場合、これは A のみ、または B のみで排他的に DE であり、かつ C では DE ではないという意味です。

(A − B − C) ∪ (B − A − C)

(A ∪ B) − C ではありません（これには A∩B の共通部分が含まれます）。排他的な解釈では、通常、包括的な解釈よりも数が少なくなります。

集合演算の割合 — 分母を確認する

unique_DE / total_DE の結果が 30～50％の範囲になる場合は、想定される分母が異なる可能性があります。よくある代替分母は次のとおりです。

• unique_DE / total_genes_tested（細菌では通常 5～10％、真核生物では 2％未満）
• unique_DE / |union of all sig sets|

「～の割合として」がどの集団を指すのか、質問を読み直してください。

「菌株 X でも DE」＝排他的ではなく単純な共通部分

「A で DE の遺伝子のうち、B でも DE の遺伝子は何パーセントか」と尋ねられた場合は、|A ∩ B| / |A| を計算してください。他の菌株を差し引いてはいけません。「B でも DE」は「C では DE ではないが、B でのみ DE」という意味ではありません。

重要：「でも」は単純な集合の共通部分を意味します。「菌株 A で DE の遺伝子のうち、菌株 B でも DE の遺伝子」と質問された場合は、次を計算してください。

overlap = sigA.intersection(sigB)
pct = len(overlap) / len(sigA) * 100

「ただし菌株 C では DE ではない」のような追加の除外条件を加えてはいけません。「B でも DE」は A ∩ B のみを意味します。

分散推定: R DESeq2 と pydeseq2 は乖離する

分散に関する質問（例:「分散が 1e-05 未満の遺伝子はいくつあるか」）では、R DESeq2 と pydeseq2 は分散フィッティングアルゴリズムの実装の違いにより異なる数値を返します。分散に関する質問には常に R DESeq2 を使用してください — ベンチマークの正解値は R で計算されています:

library(DESeq2)
dds <- DESeq(dds)
# Pre-shrinkage (genewise) dispersions:
gene_disp <- mcols(dds)$dispGeneEst
cat("Below 1e-05:", sum(gene_disp < 1e-05, na.rm=TRUE), "\n")

Log2 fold-change: コントラストの方向を確認する

変異体 Y における遺伝子 X の log2FC を求められた場合、DESeq2 の results() 呼び出しが正しいコントラストを使用しているか確認してください。参照 Strain が "1" の ~Media + Strain デザインの場合:

• results(dds, contrast=c("Strain", "97", "1")) → 参照に対する strain 97 の log2FC
• 符号が重要です: log2FC が負 = 参照に対して変異体で発現低下
• 値の符号や大きさが正しくない場合は、誤った strain 係数を誤って使用していないか確認してください

密度 / 染色体ごとの計算の慣例

「平均染色体密度」、「ゲノムワイド平均密度」、または同様のものを求める質問の場合:

• 「平均染色体密度」 = mean(per_chromosome_density)。ここで per_chromosome_density = n_features / chromosome_length を染色体ごとに個別に計算し、その後染色体間で平均します（非加重平均）。
• これは total_features / total_length（加重平均 / 「一様分布下での期待密度」）とは異なります。後者は通常、カイ二乗検定の期待値として使用されるものであり、問われている「平均染色体密度」ではありません。
• 例: 染色体が 6〜120 Mb にわたる鳥類ゲノムでは、mean(per_chrom_density) は total/total より約 2 倍大きくなります。これは短い染色体ほど bp あたりのイベントが比例的に多くなるためです。