pydeseq2

PyDESeq2

概述

PyDESeq2 是 DESeq2 的 Python 实现，用于分析 bulk RNA-seq 数据的差异表达。它支持设计并执行完整的工作流，从数据加载到结果解读，包括公式化的单因素和多因素设计、带多重检验校正的 Wald 检验、可选的 apeGLM 收缩，以及与 pandas 和 AnnData 的集成。

何时使用此技能

在以下情况下应使用此技能：

• 分析 bulk RNA-seq 计数数据以检测差异表达
• 比较不同实验条件下的基因表达（例如，处理组与对照组）
• 执行可考虑批次效应或协变量的多因素设计
• 将基于 R 的 DESeq2 工作流转换为 Python
• 将差异表达分析集成到基于 Python 的流水线中
• 用户提到“DESeq2”“差异表达”“RNA-seq 分析”或“PyDESeq2”

快速入门工作流

对于希望执行标准差异表达分析的用户：

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
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. Make the reference level explicit and fit DESeq2
metadata["condition"] = pd.Categorical(
    metadata["condition"], categories=["control", "treated"]
)
inference = DefaultInference(n_cpus=4)
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,
    inference=inference,
)
dds.deseq2()

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

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

核心工作流步骤

六个步骤及其代码见 references/core_workflow_steps.md：

1. 数据准备 — 原始整数计数，基因作为列、样本作为行， 并提供匹配的元数据。切勿向 DESeq2 输入归一化或转换后的值。
2. 设计指定 — 设计因子以及各因子的参考水平。
3. DESeq2 拟合 — 大小因子、离散度和 GLM 拟合。
4. 统计检验 — 对指定对比执行 Wald 检验。
5. 可选的 LFC 收缩 — 用于排序和可视化。
6. 结果导出 — 包含校正后 p 值的结果表。

多因素设计、对比和交互项见 references/analysis_patterns.md。

使用分析脚本

此技能包含用于标准分析的完整命令行脚本：

# 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 \
  --shrink-coeff "condition[T.treated]" \
  --plots

脚本功能：

• 自动加载和验证数据
• 过滤基因和样本
• 执行完整的 DESeq2 流水线
• 使用可自定义参数进行统计检验
• 导出结果（CSV 和可移植的 AnnData/H5AD）
• 明确支持 PyDESeq2 0.5.x 的 LFC 收缩系数
• 可选的可视化（火山图和 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(coeff="condition[T.treated]")
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.obs["size_factors"])

# Examine dispersion estimates
import matplotlib.pyplot as plt
plt.hist(dds.var["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.var["dispersions"], bins=50)
plt.show()

# Check size factors
print(dds.obs["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. **基因过滤：**过滤低计数基因（例如，总读取数 < 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. **保存中间对象：**建议使用 dds.to_picklable_anndata().write_h5ad("dds_result.h5ad") 保存可移植的输出。仅加载由你自己创建且信任的 pickle 文件。

安装和要求

uv pip install pydeseq2==0.5.4

系统要求：

• Python 3.11+
• PyDESeq2 0.5.4
• pandas 2.2.0+
• numpy 2.0.0+
• scipy 1.12.0+
• scikit-learn 1.4.0+
• anndata 0.11.0+
• formulaic 1.0.2+ 和 formulaic-contrasts 0.2.0+

可选的可视化依赖：

• matplotlib
• seaborn

更多资源

• 官方文档： https://pydeseq2.readthedocs.io
• GitHub 仓库： https://github.com/scverse/PyDESeq2
• 论文： Muzellec et al. (2023) Bioinformatics，DOI: 10.1093/bioinformatics/btad547
• 原始 DESeq2（R）： Love et al. (2014) Genome Biology，DOI: 10.1186/s13059-014-0550-8

Agent 操作流程

1. 检查环境。 确认工具版本、参考基因组／注释版本和输入格式（FASTQ、BAM、VCF、h5ad）。
2. 明确输入。 根据数据或用户提供的信息确认格式、标识符和参数。任何会改变结果的值都应询问确认，不要猜测。
3. 先运行小规模版本。 先在小型子集上运行流水线（一个样本、一个染色体、几千个细胞）。
4. 按照上述说明和参考资料执行完整任务。
5. 验证结果。 在解读结果前，检查 QC 指标、样本身份、基因组版本一致性和批次效应。
6. 汇报。 说明运行了什么（版本、命令、参数）、检查了什么，以及还有哪些不确定之处。

如果发生这种情况	请执行此操作
输入之间的基因组版本或标识符不匹配	停止并进行协调（坐标转换、ID 映射）后再继续。
已安装版本中缺少这些说明里的函数、标志或端点	查阅已安装版本自身的文档（help()、--help、官方文档），进行调整并告知用户。切勿虚构 API。
必需的输入、标识符或参数不明确	询问用户，或在运行前明确说明所作的假设。

完整性规则

• 切勿捏造结果、参数、标识符、引用或统计数据。如果无法运行或验证某项内容，请直言说明。
• 不要对超出统计结果支持范围的生物学意义作出解读；报告多重检验校正结果。
• 此处与版本相关的细节可能已经过时：请对照已安装版本的官方文档进行确认。
• 对于会产生费用、占用共享 GPU 或云配额、涉及个人或患者数据，或无法撤销的操作，请先征求同意。

相关技能

• bulk-rnaseq：端到端 bulk RNA-seq 编排器——从原始 FASTQ reads 开始，完成 QC 和剪切（FastQC、fastp/Trim Galore）、比对和定量…
• scanpy：标准的单细胞 RNA-seq 分析流水线。
• arboreto：使用可扩展算法（GRNBoost2、GENIE3）从基因表达数据推断基因调控网络（GRN）。