pydeseq2

PyDESeq2

概要

PyDESeq2 は、バルク RNA-seq データを用いた差次的発現解析のための DESeq2 の Python 実装です。データの読み込みから結果の解釈まで、単一因子および多因子の formulaic デザイン、複数検定補正を伴う Wald 検定、オプションの apeGLM による縮小、pandas および AnnData との統合を含む、完全なワークフローを設計して実行できます。

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

次の場合にこのスキルを使用してください：

• バルク RNA-seq のカウントデータを解析して差次的発現を調べる
• 実験条件間（例：処理群と対照群）の遺伝子発現を比較する
• バッチ効果や共変量を考慮した多因子デザインを実施する
• R ベースの DESeq2 ワークフローを Python に変換する
• 差次的発現解析を Python ベースのパイプラインに統合する
• ユーザーが「DESeq2」、「差次的発現」、「RNA-seq 解析」、「PyDESeq2」に言及している

クイックスタートのワークフロー

標準的な差次的発現解析を実施する場合：

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats

# 1. Load data
counts_df = pd.read_csv("counts.csv", index_col=0).T  # Transpose to samples × genes
metadata = pd.read_csv("metadata.csv", index_col=0)

# 2. Filter low-count genes
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# 3. Make the reference level explicit and fit DESeq2
metadata["condition"] = pd.Categorical(
    metadata["condition"], categories=["control", "treated"]
)
inference = DefaultInference(n_cpus=4)
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,
    inference=inference,
)
dds.deseq2()

# 4. Perform statistical testing
ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],
    inference=inference,
)
ds.summary()

# 5. Access results
results = ds.results_df
significant = results[results.padj < 0.05]
print(f"Found {len(significant)} significant genes")

ワークフローの主な手順

コード付きの 6 つの手順については、こちらを参照してください。 references/core_workflow_steps.md：

1. データの準備 — 遺伝子を列、サンプルを行とする生の整数カウントと、対応するメタデータを用意します。正規化または変換済みの値を DESeq2 に渡してはいけません。
2. デザインの指定 — デザイン因子と、それぞれの参照レベルを指定します。
3. DESeq2 のフィッティング — サイズファクター、分散、および GLM のフィッティングを行います。
4. 統計検定 — 名前付きコントラストに対して Wald 検定を行います。
5. オプションの LFC 縮小 — ランキングと可視化に使用します。
6. 結果のエクスポート — 調整済み p 値を含む結果テーブルを出力します。

多因子デザイン、コントラスト、交互作用項については、こちらを参照してください。 references/analysis_patterns.md。

解析スクリプトの使用

このスキルには、標準的な解析を行うための完全なコマンドラインスクリプトが含まれています：

# Basic usage
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~condition" \
  --contrast condition treated control \
  --output results/

# With additional options
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~batch + condition" \
  --contrast condition treated control \
  --output results/ \
  --min-counts 10 \
  --alpha 0.05 \
  --n-cpus 4 \
  --shrink-coeff "condition[T.treated]" \
  --plots

スクリプトの機能：

• データの自動読み込みと検証
• 遺伝子とサンプルのフィルタリング
• DESeq2 パイプライン全体の実行
• カスタマイズ可能なパラメーターによる統計検定
• 結果のエクスポート（CSV および移植性の高い AnnData/H5AD）
• PyDESeq2 0.5.x 向けの LFC 縮小係数の明示的な指定
• オプションの可視化（ボルケーノプロットと MA プロット）

単体の解析ツールが必要な場合や、複数のデータセットを一括処理する場合は、scripts/run_deseq2_analysis.py を案内してください。

結果の解釈

有意な遺伝子の特定

# Filter by adjusted p-value
significant = ds.results_df[ds.results_df.padj < 0.05]

# Filter by both significance and effect size
sig_and_large = ds.results_df[
    (ds.results_df.padj < 0.05) &
    (abs(ds.results_df.log2FoldChange) > 1)
]

# Separate up- and down-regulated
upregulated = significant[significant.log2FoldChange > 0]
downregulated = significant[significant.log2FoldChange < 0]

print(f"Upregulated: {len(upregulated)}")
print(f"Downregulated: {len(downregulated)}")

ランキングと並べ替え

# Sort by adjusted p-value
top_by_padj = ds.results_df.sort_values("padj").head(20)

# Sort by absolute fold change (use shrunk values)
ds.lfc_shrink(coeff="condition[T.treated]")
ds.results_df["abs_lfc"] = abs(ds.results_df.log2FoldChange)
top_by_lfc = ds.results_df.sort_values("abs_lfc", ascending=False).head(20)

# Sort by a combined metric
ds.results_df["score"] = -np.log10(ds.results_df.padj) * abs(ds.results_df.log2FoldChange)
top_combined = ds.results_df.sort_values("score", ascending=False).head(20)

品質指標

# Check normalization (size factors should be close to 1)
print("Size factors:", dds.obs["size_factors"])

# Examine dispersion estimates
import matplotlib.pyplot as plt
plt.hist(dds.var["dispersions"], bins=50)
plt.xlabel("Dispersion")
plt.ylabel("Frequency")
plt.title("Dispersion Distribution")
plt.show()

# Check p-value distribution (should be mostly flat with peak near 0)
plt.hist(ds.results_df.pvalue.dropna(), bins=50)
plt.xlabel("P-value")
plt.ylabel("Frequency")
plt.title("P-value Distribution")
plt.show()

可視化のガイドライン

ボルケーノプロット

有意性と効果量を可視化します：

import matplotlib.pyplot as plt
import numpy as np

results = ds.results_df.copy()
results["-log10(padj)"] = -np.log10(results.padj)

plt.figure(figsize=(10, 6))
significant = results.padj < 0.05

plt.scatter(
    results.loc[~significant, "log2FoldChange"],
    results.loc[~significant, "-log10(padj)"],
    alpha=0.3, s=10, c='gray', label='Not significant'
)
plt.scatter(
    results.loc[significant, "log2FoldChange"],
    results.loc[significant, "-log10(padj)"],
    alpha=0.6, s=10, c='red', label='padj < 0.05'
)

plt.axhline(-np.log10(0.05), color='blue', linestyle='--', alpha=0.5)
plt.xlabel("Log2 Fold Change")
plt.ylabel("-Log10(Adjusted P-value)")
plt.title("Volcano Plot")
plt.legend()
plt.savefig("volcano_plot.png", dpi=300)

MA プロット

フォールドチェンジと平均発現量を示します：

plt.figure(figsize=(10, 6))

plt.scatter(
    np.log10(results.loc[~significant, "baseMean"] + 1),
    results.loc[~significant, "log2FoldChange"],
    alpha=0.3, s=10, c='gray'
)
plt.scatter(
    np.log10(results.loc[significant, "baseMean"] + 1),
    results.loc[significant, "log2FoldChange"],
    alpha=0.6, s=10, c='red'
)

plt.axhline(0, color='blue', linestyle='--', alpha=0.5)
plt.xlabel("Log10(Base Mean + 1)")
plt.ylabel("Log2 Fold Change")
plt.title("MA Plot")
plt.savefig("ma_plot.png", dpi=300)

よくある問題のトラブルシューティング

データ形式の問題

問題：「カウントとメタデータのインデックスが一致しません」

**解決策：**サンプル名が正確に一致していることを確認します

print("Counts samples:", counts_df.index.tolist())
print("Metadata samples:", metadata.index.tolist())

# Take intersection if needed
common = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common]
metadata = metadata.loc[common]

問題：「すべての遺伝子のカウントが 0 です」

**解決策：**データの転置が必要かどうかを確認します

print(f"Counts shape: {counts_df.shape}")
# If genes > samples, transpose is needed
if counts_df.shape[1] < counts_df.shape[0]:
    counts_df = counts_df.T

デザイン行列の問題

問題：「デザイン行列のランクが不足しています」

**原因：**交絡した変数（例：処理群のサンプルがすべて 1 つのバッチに含まれている）

**解決策：**交絡変数を削除するか、交互作用項を追加します

# Check confounding
print(pd.crosstab(metadata.condition, metadata.batch))

# Either simplify design or add interaction
design = "~condition"  # Remove batch
# OR
design = "~condition + batch + condition:batch"  # Model interaction

有意な遺伝子がない場合

診断：

# Check dispersion distribution
plt.hist(dds.var["dispersions"], bins=50)
plt.show()

# Check size factors
print(dds.obs["size_factors"])

# Look at top genes by raw p-value
print(ds.results_df.nsmallest(20, "pvalue"))

考えられる原因：

• 効果量が小さい
• 生物学的変動が大きい
• サンプルサイズが不十分
• 技術的な問題（バッチ効果、外れ値）

リファレンスドキュメント

このワークフロー中心のガイドで扱う内容を超える詳細については、以下を参照してください。

• API リファレンス（references/api_reference.md）：PyDESeq2 のクラス、メソッド、データ構造に関する完全なドキュメントです。パラメーターの詳細情報が必要な場合や、オブジェクトの属性を理解する場合に使用してください。

• ワークフローガイド（references/workflow_guide.md）：分析ワークフロー全体、データ読み込みのパターン、多因子デザイン、トラブルシューティング、ベストプラクティスを扱う詳細なガイドです。複雑な実験デザインを扱う場合や、問題が発生した場合に使用してください。

以下が必要な場合は、これらのリファレンスをコンテキストに読み込んでください。

• 詳細な API ドキュメント：Read references/api_reference.md
• 包括的なワークフロー例：Read references/workflow_guide.md
• トラブルシューティングのガイダンス：Read references/workflow_guide.md（Troubleshooting セクションを参照）

重要な注意事項

1. データの向きが重要です： カウント行列は通常、遺伝子 × サンプルとして読み込まれますが、サンプル × 遺伝子にする必要があります。必要に応じて、必ず .T で転置してください。

2. サンプルのフィルタリング： エラーを避けるため、分析前にメタデータが欠損しているサンプルを削除してください。

3. 遺伝子のフィルタリング： 検出力を高め、計算時間を短縮するため、カウント数の少ない遺伝子（例：総リード数が 10 未満）をフィルタリングしてください。

4. デザイン式の順序： 調整変数は、注目する変数より前に配置してください（例："~batch + condition"。"~condition + batch" ではありません）。

5. LFC の縮小を行うタイミング： 縮小は統計検定後に、可視化や順位付けの目的に限って適用してください。P 値は縮小前の推定値に基づいたままです。

6. 結果の解釈： 有意性の判定には、生の P 値ではなく padj < 0.05 を使用してください。Benjamini-Hochberg 法により、偽発見率が制御されます。

7. コントラストの指定： 形式は [variable, test_level, reference_level] です。test_level を reference_level と比較します。

8. 中間オブジェクトの保存： 移植可能な出力には、dds.to_picklable_anndata().write_h5ad("dds_result.h5ad") の使用を推奨します。自分で作成し、信頼できる pickle ファイルのみを読み込んでください。

インストールと要件

uv pip install pydeseq2==0.5.4

システム要件：

• Python 3.11+
• PyDESeq2 0.5.4
• pandas 2.2.0+
• numpy 2.0.0+
• scipy 1.12.0+
• scikit-learn 1.4.0+
• anndata 0.11.0+
• formulaic 1.0.2+ および formulaic-contrasts 0.2.0+

可視化に必要なオプションのライブラリ：

• matplotlib
• seaborn

その他のリソース

• 公式ドキュメント： https://pydeseq2.readthedocs.io
• GitHub リポジトリ： https://github.com/scverse/PyDESeq2
• 論文： Muzellec et al. (2023) Bioinformatics, DOI: 10.1093/bioinformatics/btad547
• オリジナルの DESeq2（R）： Love et al. (2014) Genome Biology, DOI: 10.1186/s13059-014-0550-8

エージェントの作業手順

1. 環境を確認します。 ツールのバージョン、参照ゲノム／アノテーションのビルド、入力形式（FASTQ、BAM、VCF、h5ad）を確認します。
2. 入力を明確にします。 データまたはユーザーから、形式、識別子、パラメーターを確認します。結果に影響する値は推測せず、質問してください。
3. まず小規模なバージョンを実行します。 最初に、パイプラインを小さなサブセット（1サンプル、1染色体、数千個の細胞）で実行します。
4. 上記の手順と参照情報を使って、タスク全体を実行します。
5. 結果を検証します。 結果を解釈する前に、QC指標、サンプルの同一性、ゲノムビルドの整合性、バッチ効果を確認します。
6. 報告します。 実行した内容（バージョン、コマンド、パラメーター）、確認した内容、未確定の事項を明記します。

次のような場合	対応方法
入力間でゲノムビルドや識別子が一致しない	いったん停止し、続行する前に調整します（リフトオーバー、ID マッピング）。
インストール済みバージョンに、これらの手順に記載された関数、フラグ、エンドポイントがない	インストール済みバージョンのドキュメント（help()、--help、公式ドキュメント）を確認して適宜対応し、ユーザーに伝えます。API を決して捏造しないでください。
必須の入力、識別子、パラメーターが曖昧	ユーザーに質問するか、実行前に仮定を明示します。

完全性に関するルール

• 結果、パラメーター、識別子、引用、統計値を決して捏造しないでください。実行または検証できない場合は、その旨を明確に伝えます。
• 統計的な裏付けを超えて生物学的な意義を解釈しないでください。多重検定補正について報告します。
• ここに記載されたバージョン固有の情報は古くなっている可能性があるため、インストール済みバージョンの公式ドキュメントで確認してください。
• 費用が発生する操作、共有 GPU やクラウドクォータを消費する操作、個人データや患者データに触れる操作、取り消せない操作を行う前に、許可を求めてください。

関連スキル

• bulk-rnaseq：バルク RNA-seq のエンドツーエンドのオーケストレーター — 生の FASTQ リードから QC とトリミング（FastQC、fastp/Trim Galore）、アライメント、定量までを実行します…
• scanpy：標準的なシングルセル RNA-seq 解析パイプライン。
• arboreto：スケーラブルなアルゴリズム（GRNBoost2、GENIE3）を使って、遺伝子発現データから遺伝子制御ネットワーク（GRN）を推定します。