Pydeseq2 Differential Expression

何时使用

当你有批量 RNA-seq 的原始（未归一化）整数计数矩阵和对应样本元数据，想找出条件间差异表达基因时使用本条。典型场景：

• 两组比较（如 treated vs control），需要正经的统计检验而非简单 fold change。
• 多因子设计、校正批次/协变量（如 ~batch + condition）。
• 对一个模型跑多个 contrast（多个处理组 vs 共享对照）。
• 应用 apeGLM log2 fold change 收缩，用于排序与可视化。
• 产出火山图、MA 图、p 值分布/离散度 QC 图。
• 把 R 版 DESeq2 流程迁到纯 Python，或嵌入 scanpy/pandas 的 Python 生信管线。

不该用的边界：

• 上游比对/定量（FASTQ→counts）—— 用比对定量管线（如 nextflow-pipeline-builder）。
• 单细胞「每细胞」DE：细胞非独立观测会膨胀 p 值，先 pseudobulk 聚合（scanpy sc.get.aggregate）再喂本条。
• 通路/基因集富集 —— 属下游，交 gene-set-enrichment-analysis（用本条输出的 stat 列做预排序 GSEA）。
• 需要最广方法覆盖与社区验证的参考实现 —— 用 R/Bioconductor 的 DESeq2 或 edgeR。

步骤

PyDESeq2 期望 counts 为 samples × genes 的非负整数 DataFrame，metadata 为 samples × variables 且索引对齐。

1. 载入与校验：读 counts/metadata；CSV 常为 genes × samples，需转置；按索引交集对齐两表。
2. 基因过滤：去掉低表达基因（提统计功效、减多重检验负担）。
3. 建 DeseqDataSet 并拟合：指定 Wilkinson 设计公式，dds.deseq2() 一次跑完 size factor → 离散度 → 趋势 → MAP 收缩 → LFC 拟合。
4. Wald 检验：DeseqStats 指定 contrast=[变量, 检验水平, 参考水平]，ds.summary() 出结果。
5. （可选）LFC 收缩：ds.summary() 之后再 ds.lfc_shrink()；收缩值用于可视化/排序，不用于显著性判定。
6. 过滤与导出：按 padj + |log2FoldChange| 筛选，拆上调/下调，导 CSV。
7. 可视化 + QC：火山图、MA 图、p 值分布 + 离散度图。

指令

pip install pydeseq2 matplotlib seaborn
# 依赖：pydeseq2>=0.4, pandas>=1.4, numpy>=1.23, scipy>=1.11, scikit-learn>=1.1, anndata>=0.8；Python 3.10+

关键参数（建模时把协变量放在关注变量之前）：

示例

载入、对齐与过滤：

import pandas as pd
counts_raw = pd.read_csv("counts.csv", index_col=0)
metadata = pd.read_csv("metadata.csv", index_col=0)
counts_df = counts_raw.T if counts_raw.shape[0] > counts_raw.shape[1] else counts_raw  # → samples×genes
common = counts_df.index.intersection(metadata.index)        # 对齐索引，避免 Index mismatch
counts_df, metadata = counts_df.loc[common], metadata.loc[common]
counts_df = counts_df.loc[:, counts_df.sum(axis=0) >= 10]    # 过滤总计数<10 的基因

拟合 → 检验 →（可选）收缩 → 导出：

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=counts_df, metadata=metadata,
                   design="~condition", refit_cooks=True, n_cpus=4)
dds.deseq2()   # size factors → dispersions → trend → MAP → LFC

ds = DeseqStats(dds, contrast=["condition", "treated", "control"],
                alpha=0.05, cooks_filter=True, independent_filter=True)
ds.summary()
results = ds.results_df
ds.lfc_shrink()   # 必须在 summary() 之后；只改 log2FoldChange，供可视化/排序

sig = results[(results.padj < 0.05) & (results.log2FoldChange.abs() > 1.0)].copy()
sig[sig.log2FoldChange > 0].sort_values("padj").to_csv("deseq2_upregulated.csv")
sig[sig.log2FoldChange < 0].sort_values("padj").to_csv("deseq2_downregulated.csv")
results.to_csv("deseq2_all_results.csv")

火山图（MA 图把横轴换成 log10(baseMean+1)、纵轴 log2FoldChange）：

import matplotlib.pyplot as plt, numpy as np
res = results.dropna(subset=["padj"]).copy()
res["-log10_padj"] = -np.log10(res.padj)
is_sig = (res.padj < 0.05) & (res.log2FoldChange.abs() > 1.0)
fig, ax = plt.subplots(figsize=(10, 7))
ax.scatter(res.loc[~is_sig,"log2FoldChange"], res.loc[~is_sig,"-log10_padj"], c="grey", s=8, alpha=.3)
ax.scatter(res.loc[is_sig & (res.log2FoldChange>0),"log2FoldChange"],
           res.loc[is_sig & (res.log2FoldChange>0),"-log10_padj"], c="firebrick", s=12, alpha=.6)
ax.scatter(res.loc[is_sig & (res.log2FoldChange<0),"log2FoldChange"],
           res.loc[is_sig & (res.log2FoldChange<0),"-log10_padj"], c="steelblue", s=12, alpha=.6)
ax.axhline(-np.log10(0.05), ls="--", c="k", alpha=.4); ax.axvline(-1, ls="--", c="k", alpha=.4); ax.axvline(1, ls="--", c="k", alpha=.4)
plt.tight_layout(); plt.savefig("volcano_plot.png", dpi=300)

多 contrast / 批次校正：

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition"); dds.deseq2()
for name, c in {"A_vs_ctrl":["condition","treatment_A","control"],
                "B_vs_ctrl":["condition","treatment_B","control"]}.items():
    s = DeseqStats(dds, contrast=c, alpha=0.05); s.summary()
    s.results_df.to_csv(f"results_{name}.csv")

# 批次校正：协变量放公式最前；先查是否与条件混杂
print(pd.crosstab(metadata["batch"], metadata["condition"]))
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~batch + condition"); dds.deseq2()

注意事项

• counts 必须是原始未归一化整数；size factor >5 多半是误传了已归一化数据或文库异常。
• 样本名要对齐：Index mismatch 用 counts_df.index.intersection(metadata.index) 解决。
• 设计矩阵满秩：Design matrix is not full rank 说明变量混杂（如某批次全是 treated），用 pd.crosstab() 查，简化设计或剔除混杂变量。
• padj 全 NaN：基因零方差/全零计数 —— 加严过滤（提高 min_total_counts）。
• 没有显著基因：先看 p 值分布（健康分布应近似平坦 + 0 附近有峰）；可放宽 alpha 或 |LFC| 阈值。
• 收缩顺序：永远先 ds.summary() 再 ds.lfc_shrink()；收缩值仅供排序/作图，显著性看未收缩的 padj。
• 下游富集：用 results["stat"]（带方向与证据强度）构造预排序 GSEA，比按 log2FoldChange 更稳。
• 内存不足：先更激进地预过滤基因、减小 n_cpus。
• 断点续跑：pickle.dump(dds.to_picklable_anndata(), f) 存拟合结果，免重跑昂贵步骤。

互见

• combines_with：gene-set-enrichment-analysis —— 把 DE 基因 / stat 排序表交给富集分析做通路解读（标准下游）。
• combines_with：single-cell-rnaseq-analysis —— 单细胞跨条件严谨 DE 需先 pseudobulk 聚合再用本条。
• related：genomic-file-toolkit —— 处理上游基因组文件与计数表 I/O。
• related：nextflow-pipeline-builder —— 用管线产出 counts 矩阵作为本条输入。

本条采编自 jaechang-hits/SciAgent-Skills（CC-BY-4.0），适配重写而非逐字翻译。