pydeseq2-differential-expression

PyDESeq2による差次的発現解析

概要

PyDESeq2は、バルクRNA-seqのカウントデータを用いた差次的遺伝子発現解析のために、RのDESeq2パッケージをPythonで再実装したものです。遺伝子ごとに負の二項一般化線形モデルを当てはめ、経験ベイズ縮小によって分散を推定し、Benjamini-Hochberg法によるFDR補正を適用したWald検定を行います。このスキルでは、未加工のカウントから論文掲載に適した結果テーブルや可視化まで、パイプライン全体を扱います。

使用する場面

• バルクRNA-seqを用いて、2つ以上の実験条件間で差次的に発現する遺伝子を特定する場合
• 適切な統計検定を用いて、2群間の比較（例：処理群と対照群）を行う場合
• バッチ効果や共変量（例：~batch + condition）を考慮した多因子デザインを実行する場合
• ランキングや可視化のために、log2 fold changeの縮小（apeGLM）を適用する場合
• DE後の差次的発現結果について論文掲載品質のプロットを作成するには、omics-plotting SKILLを使用してください
• RベースのDESeq2ワークフローを純粋なPython環境に移行する場合
• DE解析をより大規模なPythonバイオインフォマティクスパイプライン（例：scanpy、pandas）に統合する場合
• 最も幅広い手法への対応とコミュニティによる検証を備えたRの参照実装を使用するには、代わりにDESeq2（R/Bioconductor）またはedgeRを使用してください

前提条件

• Pythonパッケージ：pydeseq2>=0.4、pandas>=1.4、numpy>=1.23、scipy>=1.11、scikit-learn>=1.1、anndata>=0.8
• データ要件：未加工（正規化されていない）の整数カウント行列（サンプル x 遺伝子）＋サンプルメタデータのDataFrame
• 環境：Python 3.10以降。可視化にはmatplotlib、seabornを任意で使用できます

pip install pydeseq2 matplotlib seaborn

事前確認

解析コードを書く前に、ユーザーと以下の点を確認してください。

decisions:
  - id: D1
    param: design
    kind: required
    source: data
    ask: "Which column of the sample sheet separates the groups you want to compare?"
    default: null

  - id: D2
    param: contrast
    kind: required
    source: user
    depends_on: [D1]
    ask: "Within that column, which group is the baseline the others are measured against?"
    default: null

  - id: D3
    param: design
    kind: required
    source: data
    depends_on: [D1]
    ask: "Is there a nuisance variable to adjust for (batch, donor, sex, sequencing run)?"
    default: "none - compare groups without covariates"
    skip_if: "sample sheet carries no column besides the grouping variable"

  - id: D4
    param: alpha
    kind: optional
    source: user
    ask: "How strong must the evidence be before a gene counts as changed (false-discovery rate)?"
    default: 0.05

  - id: D5
    param: lfc_cutoff
    kind: optional
    source: user
    ask: "How large a fold change is worth reporting, on top of the statistical call?"
    default: "|log2FC| > 0.58 (1.5x change)"

  - id: D6
    param: min_total_counts
    kind: optional_conditional
    source: data
    ask: "Genes with very few reads cannot be tested - where should the floor sit?"
    default: 10

  - id: D7
    param: lfc_shrink
    kind: optional
    source: user
    ask: "Shrink fold-change estimates so low-count genes stop dominating the rankings and plots?"
    default: "applied for ranking and visualization, not for significance calls"

  - id: D8
    param: refit_cooks, cooks_filter, independent_filter
    kind: optional
    source: user
    ask: "Should DESeq2 retain its default outlier replacement and adaptive filtering behavior?"
    default: true

  - id: D9
    param: n_cpus
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the result"
    default: 4

designが2回登場するのは、D1とD3の両方から数式が組み立てられるためです。 次にD2で、当てはめたモデルのどの係数を検定するかを選択します。

ワークフロー

ステップ1：データの読み込みと検証

カウント行列とメタデータを読み込みます。PyDESeq2では、非負の整数を含むサンプル x 遺伝子のDataFrameとしてカウントを、インデックスが一致するサンプル x 変数のDataFrameとしてメタデータを指定します。

import pandas as pd

# Load data — typical CSV has genes as rows, samples as columns
counts_raw = pd.read_csv("counts.csv", index_col=0)
metadata = pd.read_csv("metadata.csv", index_col=0)

# Transpose if needed: PyDESeq2 requires samples x genes
if counts_raw.shape[0] > counts_raw.shape[1]:
    counts_df = counts_raw.T  # genes x samples → samples x genes
else:
    counts_df = counts_raw

# Validate alignment
common_samples = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common_samples]
metadata = metadata.loc[common_samples]

print(f"Samples: {counts_df.shape[0]}, Genes: {counts_df.shape[1]}")
print(f"Metadata columns: {list(metadata.columns)}")
print(f"Condition counts:\n{metadata['condition'].value_counts()}")

ステップ2：遺伝子のフィルタリング

発現量の低い遺伝子を除外し、統計的検出力を高め、多重検定の負担を軽減します。

# Filter genes with total counts below threshold
min_total_counts = 10
gene_counts = counts_df.sum(axis=0)
genes_to_keep = gene_counts[gene_counts >= min_total_counts].index
counts_df = counts_df[genes_to_keep]

# Optional: require minimum counts in a minimum number of samples
min_count_per_sample = 5
min_samples = 3
genes_expressed = (counts_df >= min_count_per_sample).sum(axis=0) >= min_samples
counts_df = counts_df.loc[:, genes_expressed]

print(f"Genes after filtering: {counts_df.shape[1]}")

ステップ3：DeseqDataSetの初期化とフィッティング

DESeqデータセットオブジェクトを作成し、デザイン数式を指定して、パイプライン全体（サイズファクターの推定、分散の推定、モデルのフィッティング）を実行します。

from pydeseq2.dds import DeseqDataSet

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",   # Wilkinson-style formula
    refit_cooks=True,      # Refit after Cook's outlier removal
    n_cpus=4               # Parallel threads
)

# Run: size factors → dispersions → trend → MAP shrinkage → LFC fitting
dds.deseq2()

# Inspect normalization
print(f"Size factors (first 5): {dds.obsm['size_factors'][:5]}")
print(f"Size factor range: {dds.obsm['size_factors'].min():.2f} - {dds.obsm['size_factors'].max():.2f}")

ステップ4：統計検定（Wald検定）

差次的に発現する遺伝子を特定するため、Wald検定を行います。コントラストは[variable, test_level, reference_level]の形式で指定します。

from pydeseq2.ds import DeseqStats

ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],
    alpha=0.05,              # FDR threshold
    cooks_filter=True,       # Filter Cook's outliers
    independent_filter=True  # Independent filtering for power
)

ds.summary()

# Access full results
results = ds.results_df
print(f"Total genes tested: {len(results)}")
print(f"Significant (padj < 0.05): {(results.padj < 0.05).sum()}")

ステップ5：LFCの縮小（任意）

log2 fold changeの推定値のノイズを減らすために、apeGLMによる縮小を適用します。縮小後の値は可視化やランキングに使用し、有意性の判定には使用しないでください。

# Apply shrinkage — modifies results_df.log2FoldChange in place
ds.lfc_shrink()

# Compare pre/post shrinkage effect
print(f"Max |LFC| after shrinkage: {results.log2FoldChange.abs().max():.2f}")
print(f"Genes with |LFC| > 2: {(results.log2FoldChange.abs() > 2).sum()}")

ステップ6：結果のフィルタリングとエクスポート

有意な遺伝子をフィルタリングし、下流解析用に結果をエクスポートします。

import numpy as np

# Significance + effect size filter
significant = results[
    (results.padj < 0.05) &
    (results.log2FoldChange.abs() > 1.0)
].copy()

# Separate up/down-regulated
up = significant[significant.log2FoldChange > 0].sort_values("padj")
down = significant[significant.log2FoldChange < 0].sort_values("padj")

print(f"Upregulated: {len(up)}, Downregulated: {len(down)}")

# Export
results.to_csv("deseq2_all_results.csv")
significant.to_csv("deseq2_significant.csv")
up.to_csv("deseq2_upregulated.csv")
down.to_csv("deseq2_downregulated.csv")
print("Results exported to CSV files")

ステップ7：可視化 — 火山プロット

DEGテーブルを用意し、エクスポートしたCSVに対して**skills/data-visualization/omics-plotting/SKILL.mdを読み、「Volcano」のレシピに従ってください**（→ figures/volcano_plot.png）。

# Volcano input = log2FoldChange vs padj; saved in Step 6 as deseq2_all_results.csv
volcano_input = results[["log2FoldChange", "padj"]].dropna()
print(f"Volcano input: {len(volcano_input)} genes -> deseq2_all_results.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Volcano" recipe.

ステップ8：可視化 — MAプロット

DEGテーブルを用意し、エクスポートしたCSVに対して**skills/data-visualization/omics-plotting/SKILL.mdを読み、「MA plot」のレシピに従ってください**（→ figures/ma_plot.png）。

# MA input = log2FoldChange vs mean normalized count (baseMean); from deseq2_all_results.csv
ma_input = results[["baseMean", "log2FoldChange", "padj"]].dropna(subset=["padj"])
print(f"MA input: {len(ma_input)} genes -> deseq2_all_results.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "MA plot" recipe.

主なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
design	(必須)	Wilkinson の式	モデル式。共変量を関心対象の変数より前に指定します
contrast	None	[var, test, ref]	検定する比較対象です。None の場合は最後の係数を使用します
alpha	0.05	0.01–0.10	有意性判定の FDR 閾値です
refit_cooks	True	True/False	Cook の距離の外れ値を除外した後、モデルを再適合します
cooks_filter	True	True/False	検定中に Cook の距離によるフィルタリングを適用します
independent_filter	True	True/False	検出力を最適化するための独立フィルタリングです
n_cpus	1	1–N	分散の適合に使用する並列スレッド数です
min_total_counts (ユーザー指定)	10	5–50	遺伝子フィルタリング：全サンプルにおけるリード数の最小合計値です
lfc_shrink()	オフ	summary() の後に呼び出す	apeGLM による縮小。ノイズの多い LFC 推定値を抑えます

よく使うレシピ

レシピ：1 つのモデルから複数のコントラストを取得する

複数の処理群を共通の対照群と比較する場合に使用します。

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

contrasts = {
    "A_vs_ctrl": ["condition", "treatment_A", "control"],
    "B_vs_ctrl": ["condition", "treatment_B", "control"],
    "C_vs_ctrl": ["condition", "treatment_C", "control"],
}

for name, contrast in contrasts.items():
    ds = DeseqStats(dds, contrast=contrast, alpha=0.05)
    ds.summary()
    n_sig = (ds.results_df.padj < 0.05).sum()
    ds.results_df.to_csv(f"results_{name}.csv")
    print(f"{name}: {n_sig} significant genes")

レシピ：バッチ補正済みの解析

複数のバッチまたはシーケンシング実行から得たサンプルを扱う場合に使用します。

# Verify batch is not confounded with condition
print(pd.crosstab(metadata["batch"], metadata["condition"]))

dds = DeseqDataSet(
    counts=counts_df, metadata=metadata,
    design="~batch + condition"  # batch first, then variable of interest
)
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()
print(f"Significant genes (batch-corrected): {(ds.results_df.padj < 0.05).sum()}")

レシピ：P 値分布の QC

診断チェックです。健全な解析では、0 付近にスパイクがあり、その他は平坦なヒストグラムになります。

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# P-value histogram
axes[0].hist(ds.results_df.pvalue.dropna(), bins=50, edgecolor="black")
axes[0].set_xlabel("P-value")
axes[0].set_ylabel("Frequency")
axes[0].set_title("P-value Distribution")

# Dispersion plot
axes[1].scatter(
    np.log10(dds.varm["_normed_means"] + 1),
    np.log10(dds.varm["dispersions"]),
    alpha=0.3, s=5
)
axes[1].set_xlabel("Log10(Mean Expression + 1)")
axes[1].set_ylabel("Log10(Dispersion)")
axes[1].set_title("Dispersion vs Mean")

plt.tight_layout()
plt.savefig("qc_diagnostics.png", dpi=300)
print("Saved qc_diagnostics.png")

レシピ：DESeq データセットの保存と再読み込み

計算負荷の高い適合処理を再実行せずに解析を再開できます。

import pickle

# Save
with open("dds_fitted.pkl", "wb") as f:
    pickle.dump(dds.to_picklable_anndata(), f)
print("Saved fitted DESeqDataSet")

# Reload
with open("dds_fitted.pkl", "rb") as f:
    adata = pickle.load(f)
# Note: reload requires re-constructing DeseqDataSet from the AnnData

出力されるファイル

• deseq2_all_results.csv — 検定したすべての遺伝子の完全な結果テーブル（baseMean、log2FoldChange、lfcSE、stat、pvalue、padj）
• deseq2_significant.csv — フィルタリング済みの結果（padj < 0.05 かつ |LFC| > 1）
• deseq2_upregulated.csv — padj 順に並べた、有意に上方制御された遺伝子
• deseq2_downregulated.csv — padj 順に並べた、有意に下方制御された遺伝子
• figures/volcano_plot.png — 有意性とフォールド変化の閾値を示すボルケーノプロット
• figures/ma_plot.png — フォールド変化と平均発現量を示す MA プロット
• figures/qc_diagnostics.png — P 値分布と分散プロット

トラブルシューティング

問題	原因	解決策
ValueError: Index mismatch	カウントとメタデータでサンプル名が異なります	counts_df.index.intersection(metadata.index) を使用して揃えます
すべての遺伝子で padj = NaN になる	遺伝子の分散がゼロ、またはカウントがすべてゼロです	より厳しい遺伝子フィルタリングを適用します（min_total_counts を増やします）
Design matrix is not full rank	交絡した変数があります（例：すべての処理群が 1 つのバッチに含まれる）	pd.crosstab() で確認し、デザインを単純化するか、交絡した変数を削除します
有意な遺伝子が見つからない	効果量が小さい、変動が大きい、またはサンプル数が少ない	P 値分布を確認し、alpha または
適合中に MemoryError が発生する	遺伝子数が多すぎるか、データセットが非常に大きい	事前フィルタリングを強め、n_cpus を減らすか、RAM の多いマシンを使用します
サイズファクターが非常に大きい（>5）	ライブラリーサイズに極端な差があります	生カウントが正規化されていないことを確認し、汚染やライブラリー調製の失敗がないか確認します
縮小によって予期しない LFC が生じる	summary() より前に lfc_shrink() を呼び出しています	必ず先に ds.summary() を呼び出し、その後 ds.lfc_shrink() を呼び出します

参考資料

• PyDESeq2 のドキュメント — 公式 API リファレンスとチュートリアル
• Muzellec et al. (2023) — PyDESeq2：DESeq2 によるバルク差次的発現解析のための Python パッケージ。Bioinformatics。
• Love et al. (2014) — DESeq2 を用いた RNA-seq データのフォールド変化と分散のモデレート推定。Genome Biology。
• DESeq2 ビネット（Bioconductor） — デザイン式の例を含む、R の DESeq2 総合ガイド