pydeseq2-differential-expression

PyDESeq2 差异表达分析

概述

PyDESeq2 是 R DESeq2 软件包的 Python 重实现，用于分析 bulk RNA-seq 计数数据中的差异基因表达。它为每个基因拟合负二项广义线性模型，使用经验贝叶斯收缩估计离散度，并通过 Benjamini-Hochberg FDR 校正执行 Wald 检验。本技能涵盖从原始计数到可用于发表的结果表和可视化的完整流程。

适用场景

• 根据 bulk RNA-seq 数据，识别两个或更多实验条件之间的差异表达基因
• 使用适当的统计检验进行两组比较（例如，处理组与对照组）
• 运行多因素设计，纳入批次效应或协变量（例如，~batch + condition）
• 应用 log2 倍数变化收缩（apeGLM），用于排序和可视化
• DE 后使用 omics-plotting SKILL 为差异表达结果绘制适用于发表的图表
• 将基于 R 的 DESeq2 工作流转换为纯 Python 环境
• 将差异表达分析整合到更大的 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 会出现两次，因为公式由 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 检验以识别差异表达基因。将 contrast 指定为 [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 收缩（可选）

应用 apeGLM 收缩，以降低 log2 倍数变化估计中的噪声。收缩后的数值用于可视化和排序，不用于判定显著性。

# 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 表格，然后阅读 skills/data-visualization/omics-plotting/SKILL.md 并按照其中的“Volcano”配方，对导出的 CSV 作图（→ 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 表格，然后阅读 skills/data-visualization/omics-plotting/SKILL.md 并按照其中的“MA plot”配方，对导出的 CSV 作图（→ 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	基因筛选：所有样本中的最小总 reads 数
lfc_shrink()	关闭	在 summary() 后调用	apeGLM 收缩；降低 LFC 估计中的噪声

常见用法

用法：从一个模型进行多重对比

比较多个处理组与同一个对照组时使用。

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 值分布质控

诊断检查——健康的分析会显示平坦的直方图，并在 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 与 metadata 中的样本名称不同	使用 counts_df.index.intersection(metadata.index) 进行对齐
所有基因的 padj = NaN	基因方差为零或计数全为零	应用更严格的基因筛选（提高 min_total_counts）
Design matrix is not full rank	变量混杂（例如，所有处理组样本都在同一个批次）	使用 pd.crosstab() 检查；简化设计或移除混杂变量
未发现显著基因	效应量小、变异性高或样本量少	检查 P 值分布；放宽 alpha 或
拟合时出现 MemoryError	基因过多或数据集过大	更严格地进行预筛选；减少 n_cpus；使用 RAM 更多的计算机
size factor 非常大（>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 Vignette（Bioconductor） — 全面的 R DESeq2 指南，附设计公式示例