⚡

PyDESeq2差异基因表达分析Skill pydeseq2

这是一个用于批量RNA测序数据差异基因表达分析的Python技能。基于DESeq2算法，提供完整的分析工作流程，包括数据预处理、统计测试、结果可视化和导出。适用于生物信息学研究和RNA测序数据分析。关键词：差异基因表达，RNA测序，DESeq2，Python，统计分析。

数据分析 0 次安装 0 次浏览 更新于 3/16/2026

name: pydeseq2 description: “差异基因表达分析（Python DESeq2）。从批量RNA测序计数中识别差异表达基因，进行Wald检验、FDR校正，生成火山图和MA图，用于RNA测序分析。”

PyDESeq2

概述

PyDESeq2 是用于批量RNA测序数据差异表达分析的DESeq2的Python实现。设计和执行完整的工作流程，从数据加载到结果解释，包括单因素和多因素设计、带多重测试校正的Wald检验、可选的apeGLM收缩，以及与pandas和AnnData的集成。

何时使用此技能

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

  • 分析批量RNA测序计数数据以进行差异表达
  • 比较实验条件之间的基因表达（例如，处理组 vs 对照组）
  • 执行多因素设计，考虑批次效应或协变量
  • 将基于R的DESeq2工作流程转换为Python
  • 将差异表达分析集成到基于Python的管道中
  • 用户提到“DESeq2”、“差异表达”、“RNA测序分析”或“PyDESeq2”

快速入门工作流程

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

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# 1. 加载数据
counts_df = pd.read_csv("counts.csv", index_col=0).T  # 转置为样本 × 基因
metadata = pd.read_csv("metadata.csv", index_col=0)

# 2. 过滤低计数基因
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# 3. 初始化和拟合DESeq2
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True
)
dds.deseq2()

# 4. 执行统计测试
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

# 5. 访问结果
results = ds.results_df
significant = results[results.padj < 0.05]
print(f"找到 {len(significant)} 个显著基因")

核心工作流程步骤

步骤 1: 数据准备

输入要求：

  • 计数矩阵： 样本 × 基因的DataFrame，包含非负整数读取计数
  • 元数据： 样本 × 变量的DataFrame，包含实验因素

常见数据加载模式：

# 从CSV（典型格式：基因 × 样本，需要转置）
counts_df = pd.read_csv("counts.csv", index_col=0).T
metadata = pd.read_csv("metadata.csv", index_col=0)

# 从TSV
counts_df = pd.read_csv("counts.tsv", sep="\t", index_col=0).T

# 从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

数据过滤：

# 移除低计数基因
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# 移除缺失元数据的样本
samples_to_keep = ~metadata.condition.isna()
counts_df = counts_df.loc[samples_to_keep]
metadata = metadata.loc[samples_to_keep]

步骤 2: 设计规范

设计公式指定基因表达如何建模。

单因素设计：

design = "~condition"  # 简单的两组比较

多因素设计：

design = "~batch + condition"  # 控制批次效应
design = "~age + condition"     # 包括连续协变量
design = "~group + condition + group:condition"  # 交互效应

设计公式指南：

  • 使用Wilkinson公式表示法（R风格）
  • 将调整变量（例如，批次）放在主要关注变量之前
  • 确保变量存在于元数据DataFrame的列中
  • 使用适当的数据类型（离散变量用分类）

步骤 3: DESeq2拟合

初始化DeseqDataSet并运行完整管道：

from pydeseq2.dds import DeseqDataSet

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,  # 移除异常值后重新拟合
    n_cpus=1           # 并行处理（根据需要调整）
)

# 运行完整的DESeq2管道
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"],  # 测试处理组 vs 对照组
    alpha=0.05,                # 显著性阈值
    cooks_filter=True,         # 过滤异常值
    independent_filter=True    # 过滤低功率测试
)

ds.summary()

对比规范：

  • 格式：[变量, 测试级别, 参考级别]
  • 示例：["condition", "treated", "control"] 测试处理组 vs 对照组
  • 如果为 None，则使用设计中的最后一个系数

结果DataFrame列：

  • baseMean: 样本间平均归一化计数
  • log2FoldChange: 条件间对数2倍数变化
  • lfcSE: LFC的标准误差
  • stat: Wald检验统计量
  • pvalue: 原始p值
  • padj: 调整后的p值（通过Benjamini-Hochberg进行FDR校正）

步骤 5: 可选LFC收缩

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

ds.lfc_shrink()  # 应用apeGLM收缩

何时使用LFC收缩：

  • 用于可视化（火山图、热图）
  • 用于按效应大小对基因排序
  • 当优先考虑后续实验的基因时

重要： 收缩仅影响log2FoldChange值，不影响统计测试结果（p值保持不变）。使用收缩值进行可视化，但报告非收缩p值以表示显著性。

步骤 6: 结果导出

保存结果和中间对象：

import pickle

# 将结果导出为CSV
ds.results_df.to_csv("deseq2_results.csv")

# 仅保存显著基因
significant = ds.results_df[ds.results_df.padj < 0.05]
significant.to_csv("significant_genes.csv")

# 保存DeseqDataSet供以后使用
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} 个显著基因")

考虑批次效应

控制技术变异：

# 在设计中包括批次
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~batch + condition")
dds.deseq2()

# 在控制批次的情况下测试条件
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

连续协变量

包括连续变量如年龄或剂量：

# 确保连续变量是数字
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()

使用分析脚本

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

# 基本用法
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~condition" \
  --contrast condition treated control \
  --output results/

# 带额外选项
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。

结果解释

识别显著基因

# 按调整后p值过滤
significant = ds.results_df[ds.results_df.padj < 0.05]

# 按显著性和效应大小过滤
sig_and_large = ds.results_df[
    (ds.results_df.padj < 0.05) &
    (abs(ds.results_df.log2FoldChange) > 1)
]

# 分离上调和下调基因
upregulated = significant[significant.log2FoldChange > 0]
downregulated = significant[significant.log2FoldChange < 0]

print(f"上调: {len(upregulated)}")
print(f"下调: {len(downregulated)}")

排序和排序

# 按调整后p值排序
top_by_padj = ds.results_df.sort_values("padj").head(20)

# 按绝对倍数变化排序（使用收缩值）
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)

# 按组合指标排序
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)

质量指标

# 检查归一化（大小因子应接近1）
print("大小因子:", dds.obsm["size_factors"])

# 检查离散度估计
import matplotlib.pyplot as plt
plt.hist(dds.varm["dispersions"], bins=50)
plt.xlabel("离散度")
plt.ylabel("频率")
plt.title("离散度分布")
plt.show()

# 检查p值分布（应主要平坦，峰值接近0）
plt.hist(ds.results_df.pvalue.dropna(), bins=50)
plt.xlabel("P值")
plt.ylabel("频率")
plt.title("P值分布")
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='不显著'
)
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("对数2倍数变化")
plt.ylabel("-Log10(调整后P值)")
plt.title("火山图")
plt.legend()
plt.savefig("volcano_plot.png", dpi=300)

MA图

显示倍数变化 vs 平均表达：

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(基均值 + 1)")
plt.ylabel("对数2倍数变化")
plt.title("MA图")
plt.savefig("ma_plot.png", dpi=300)

故障排除常见问题

数据格式问题

问题： “计数和元数据之间的索引不匹配”

解决方案： 确保样本名称完全匹配

print("计数样本:", counts_df.index.tolist())
print("元数据样本:", metadata.index.tolist())

# 如果需要，取交集
common = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common]
metadata = metadata.loc[common]

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

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

print(f"计数形状: {counts_df.shape}")
# 如果基因 > 样本，则需要转置
if counts_df.shape[1] < counts_df.shape[0]:
    counts_df = counts_df.T

设计矩阵问题

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

原因： 混杂变量（例如，所有处理样本在一个批次中）

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

# 检查混杂
print(pd.crosstab(metadata.condition, metadata.batch))

# 简化设计或添加交互项
design = "~condition"  # 移除批次
# 或
design = "~condition + batch + condition:batch"  # 模型交互

无显著基因

诊断：

# 检查离散度分布
plt.hist(dds.varm["dispersions"], bins=50)
plt.show()

# 检查大小因子
print(dds.obsm["size_factors"])

# 查看原始p值最高的基因
print(ds.results_df.nsmallest(20, "pvalue"))

可能原因：

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

参考文档

如需超出此工作流程导向指南的综合详细信息：

  • API参考 (references/api_reference.md): PyDESeq2类、方法和数据结构的完整文档。需要详细参数信息或理解对象属性时使用。

  • 工作流程指南 (references/workflow_guide.md): 深入指南，涵盖完整分析工作流程、数据加载模式、多因素设计、故障排除和最佳实践。处理复杂实验设计或遇到问题时使用。

当用户需要时，将这些参考加载到上下文中：

  • 详细API文档：阅读 references/api_reference.md
  • 综合工作流程示例：阅读 references/workflow_guide.md
  • 故障排除指导：阅读 references/workflow_guide.md（见故障排除部分）

关键提醒

  1. 数据方向重要： 计数矩阵通常加载为基因 × 样本，但需要为样本 × 基因。如果需要，总是用 .T 转置。

  2. 样本过滤： 在分析前移除缺失元数据的样本以避免错误。

  3. 基因过滤： 过滤低计数基因（例如，< 10 总读取）以提高功率并减少计算时间。

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

  5. LFC收缩时机： 在统计测试后应用收缩，仅用于可视化/排序目的。P值基于非收缩估计。

  6. 结果解释： 使用 padj < 0.05 表示显著性，而不是原始p值。Benjamini-Hochberg程序控制错误发现率。

  7. 对比规范： 格式是 [变量, 测试级别, 参考级别]，其中测试级别与参考级别比较。

  8. 保存中间对象： 使用pickle保存DeseqDataSet对象供以后使用或额外分析，无需重新运行昂贵的拟合步骤。

安装和需求

PyDESeq2可以通过pip或conda安装：

# 通过pip
pip install pydeseq2

# 通过conda
conda install -c bioconda 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

额外资源