pydeseq2

PyDESeq2

概述

PyDESeq2 是 DESeq2 的 Python 实现，用于分析批量 RNA-seq 数据的差异表达。从数据加载到结果解读，设计并执行完整的工作流，包括单因素和多因素设计、带多重检验校正的 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.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. 拟合对数倍数变化
7. 计算 Cook 距离（异常值检测）
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()

对比指定方式：

• 格式：[变量, 测试水平, 参考水平]
• 示例：["condition", "treated", "control"] 用于检验处理组与对照组的差异
• 如果为 None，则使用设计中的最后一个系数

结果 DataFrame 列：

• baseMean：样本间的归一化计数均值
• log2FoldChange：条件之间的 log2 倍数变化
• lfcSE：LFC 的标准误
• stat：Wald 检验统计量
• pvalue：原始 p 值
• padj：校正后的 p 值（通过 Benjamini-Hochberg 方法进行 FDR 校正）

步骤 5：可选的 LFC 收缩

对倍数变化估计应用收缩，以减少噪声：

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)

常见分析模式

双组比较

标准的病例对照比较：

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)

常见问题排查

数据格式问题

问题：“计数与元数据的索引不匹配”

**解决方案：**确保样本名称完全一致

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]

问题：“所有基因的计数均为零”

**解决方案：**检查数据是否需要转置

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

设计矩阵问题

问题：“设计矩阵不满秩”

**原因：**变量混杂（例如，所有处理组样本都来自同一批次）

**解决方案：**移除混杂变量或添加交互项

# 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"))

可能原因：

• 效应量较小
• 生物变异较大
• 样本量不足
• 技术问题（批次效应、离群值）

参考文档

如需了解本工作流指南之外的详细信息，请参阅：

• 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. **基因筛选：**筛选低计数基因（例如，总 reads 数少于 10 的基因），以提高检验效能并缩短计算时间。

4. **设计公式顺序：**将调整变量放在目标变量之前（例如，使用 "~batch + condition"，而不是 "~condition + batch"）。

5. **LFC 收缩时机：**在统计检验后进行收缩，且仅用于可视化或排名。P 值仍基于未收缩的估计值。

6. **结果解读：**使用 padj < 0.05 判断显著性，而不是原始 P 值。Benjamini-Hochberg 程序可控制错误发现率。

7. **对比指定：**格式为 [variable, test_level, reference_level]，其中 test_level 与 reference_level 进行比较。

8. **保存中间对象：**使用 pickle 保存 DeseqDataSet 对象，以便稍后使用或执行其他分析，而无需重新运行耗时的拟合步骤。

安装和要求

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