pydeseq2

PyDESeq2

概要

PyDESeq2 は、バルク RNA-seq データの発現量差解析を行うための DESeq2 の Python 実装です。データの読み込みから結果の解釈までの完全なワークフローを設計・実行できます。これには、単一因子および多因子デザイン、多重検定補正を伴う Wald 検定、オプションの apeGLM シュリンケージ、pandas および AnnData との統合が含まれます。

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

このスキルは、次のような場合に使用します:

• バルク RNA-seq のカウントデータを発現量差解析のために解析する場合
• 実験条件間の遺伝子発現を比較する場合(例: 処理群 vs 対照群)
• バッチ効果や共変量を考慮した多因子デザインを実行する場合
• R ベースの DESeq2 ワークフローを Python に移行する場合
• Python ベースのパイプラインに発現量差解析を組み込む場合
• ユーザーが "DESeq2"、"differential expression"、"RNA-seq analysis"、"PyDESeq2" に言及する場合

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

標準的な発現量差解析を実行したいユーザー向けです:

import pandas as pd
from pydeseq2.dds import DeseqDataSet
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. Initialize and fit DESeq2
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True
)
dds.deseq2()

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

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

主要なワークフローの手順

ステップ 1: データの準備

入力要件:

• カウント行列: 非負整数のリード数を格納した、サンプル × 遺伝子の DataFrame
• メタデータ: 実験因子を格納した、サンプル × 変数の DataFrame

一般的なデータ読み込みパターン:

# From CSV (typical format: genes × samples, needs transpose)
counts_df = pd.read_csv("counts.csv", index_col=0).T
metadata = pd.read_csv("metadata.csv", index_col=0)

# From TSV
counts_df = pd.read_csv("counts.tsv", sep="\t", index_col=0).T

# From AnnData
import anndata as ad
adata = ad.read_h5ad("data.h5ad")
counts_df = pd.DataFrame(adata.X, index=adata.obs_names, columns=adata.var_names)
metadata = adata.obs

データのフィルタリング:

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

# Remove samples with missing metadata
samples_to_keep = ~metadata.condition.isna()
counts_df = counts_df.loc[samples_to_keep]
metadata = metadata.loc[samples_to_keep]

ステップ 2: デザインの指定

デザイン式は、遺伝子発現をどのようにモデル化するかを指定します。

単一因子デザイン:

design = "~condition"  # Simple two-group comparison

多因子デザイン:

design = "~batch + condition"  # Control for batch effects
design = "~age + condition"     # Include continuous covariate
design = "~group + condition + group:condition"  # Interaction effects

デザイン式のガイドライン:

• Wilkinson の式記法(R スタイル)を使用します
• 調整変数(例: バッチ)は、注目する主要な変数の前に置きます
• 変数がメタデータ DataFrame の列として存在することを確認します
• 適切なデータ型を使用します(離散変数にはカテゴリカル型)

ステップ 3: DESeq2 のフィッティング

DeseqDataSet を初期化し、完全なパイプラインを実行します:

from pydeseq2.dds import DeseqDataSet

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,  # Refit after removing outliers
    n_cpus=1           # Parallel processing (adjust as needed)
)

# Run the complete DESeq2 pipeline
dds.deseq2()

deseq2() が行うこと:

1. サイズファクターを計算する(正規化)
2. 遺伝子ごとの分散をフィッティングする
3. 分散のトレンド曲線をフィッティングする
4. 分散の事前分布を計算する
5. MAP 分散をフィッティングする(シュリンケージ)
6. log fold change をフィッティングする
7. Cook's 距離を計算する(外れ値の検出)
8. 外れ値が検出された場合は再フィッティングする(任意)

ステップ 4: 統計的検定

Wald 検定を実行し、発現量差のある遺伝子を特定します:

from pydeseq2.ds import DeseqStats

ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],  # Test treated vs control
    alpha=0.05,                # Significance threshold
    cooks_filter=True,         # Filter outliers
    independent_filter=True    # Filter low-power tests
)

ds.summary()

対比の指定:

• 形式: [variable, test_level, reference_level]
• 例: ["condition", "treated", "control"] は treated と control を比較します
• None の場合、デザイン内の最後の係数を使用します

結果 DataFrame の列:

• baseMean: サンプル間の正規化カウントの平均
• log2FoldChange: 条件間の log2 fold change
• lfcSE: LFC の標準誤差
• stat: Wald 検定統計量
• pvalue: 補正前の p 値
• padj: 調整済み p 値(Benjamini-Hochberg 法による FDR 補正)

ステップ 5: オプションの LFC シュリンケージ

fold change の推定値のノイズを減らすために、シュリンケージを適用します:

ds.lfc_shrink()  # Applies apeGLM shrinkage

LFC シュリンケージを使用する場合:

• 可視化のため(ボルケーノプロット、ヒートマップ)
• 効果量による遺伝子のランク付けのため
• 追試験の対象となる遺伝子に優先順位を付ける場合

重要: シュリンケージが影響するのは log2FoldChange の値のみで、統計検定の結果(p 値は変わりません)には影響しません。可視化にはシュリンク後の値を使用し、有意性の報告にはシュリンク前の p 値を使用してください。

ステップ 6: 結果のエクスポート

結果および中間オブジェクトを保存します:

import pickle

# Export results as CSV
ds.results_df.to_csv("deseq2_results.csv")

# Save significant genes only
significant = ds.results_df[ds.results_df.padj < 0.05]
significant.to_csv("significant_genes.csv")

# Save DeseqDataSet for later use
with open("dds_result.pkl", "wb") as f:
    pickle.dump(dds.to_picklable_anndata(), f)

一般的な解析パターン

2 群比較

標準的な症例対照比較:

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

results = ds.results_df
significant = results[results.padj < 0.05]

複数比較

複数の処理群を対照群と比較する場合:

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

treatments = ["treatment_A", "treatment_B", "treatment_C"]
all_results = {}

for treatment in treatments:
    ds = DeseqStats(dds, contrast=["condition", treatment, "control"])
    ds.summary()
    all_results[treatment] = ds.results_df

    sig_count = len(ds.results_df[ds.results_df.padj < 0.05])
    print(f"{treatment}: {sig_count} significant genes")

バッチ効果の考慮

技術的変動を制御します:

# Include batch in design
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~batch + condition")
dds.deseq2()

# Test condition while controlling for batch
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

連続共変量

年齢や投与量などの連続変数を含めます:

# Ensure continuous variable is numeric
metadata["age"] = pd.to_numeric(metadata["age"])

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~age + condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

分析スクリプトの使用

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

# 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 \
  --plots

スクリプトの機能:

• データの自動読み込みと検証
• 遺伝子およびサンプルのフィルタリング
• DESeq2 パイプラインの完全な実行
• カスタマイズ可能なパラメータによる統計検定
• 結果のエクスポート (CSV、pickle)
• オプションの可視化 (ボルケーノプロットおよび 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()
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.obsm["size_factors"])

# Examine dispersion estimates
import matplotlib.pyplot as plt
plt.hist(dds.varm["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)

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

データ形式の問題

問題: "Index mismatch between counts and metadata"

解決策: サンプル名が完全に一致していることを確認してください

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]

問題: "All genes have zero counts"

解決策: データの転置が必要かどうか確認してください

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

デザイン行列の問題

問題: "Design matrix is not full rank"

原因: 交絡変数 (例: 処理群のサンプルがすべて同一バッチに含まれる)

解決策: 交絡変数を削除するか、交互作用項を追加してください

# 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.varm["dispersions"], bins=50)
plt.show()

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

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

考えられる原因:

• 効果量が小さい
• 生物学的なばらつきが大きい
• サンプルサイズが不十分
• 技術的な問題 (バッチ効果、外れ値)

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

このワークフror指向のガイドを超える詳細については、以下を参照してください:

• 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 (トラブルシューティングのセクションを参照)

重要な注意事項

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

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

3. 遺伝子のフィルタリング: 低カウントの遺伝子 (例: 総リード数 < 10) をフィルタリングすると、検出力が向上し、計算時間を短縮できます。

4. デザイン式の順序: 調整変数は関心のある変数の前に置いてください (例: "~condition + batch" ではなく "~batch + condition")。

5. LFC シュリンケージの適用タイミング: シュリンケージは統計検定の後に、可視化やランキングの目的でのみ適用してください。P 値はシュリンクしていない推定値に基づいたままです。

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

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

8. 中間オブジェクトの保存: DeseqDataSet オブジェクトは pickle で保存しておくと、負荷の大きいフィッティングを再実行せずに、後で利用したり追加の分析を行ったりできます。

インストールと要件

uv pip install pydeseq2

システム要件:

• Python 3.10-3.11
• pandas 1.4.3+
• numpy 1.23.0+
• scipy 1.11.0+
• scikit-learn 1.1.1+
• anndata 0.8.0+

可視化に必要なオプション:

• matplotlib
• seaborn

Additional Resources

• Official Documentation: https://pydeseq2.readthedocs.io
• GitHub Repository: https://github.com/owkin/PyDESeq2
• Publication: Muzellec et al. (2023) Bioinformatics, DOI: 10.1093/bioinformatics/btad547
• Original DESeq2 (R): Love et al. (2014) Genome Biology, DOI: 10.1186/s13059-014-0550-8