#PyDESeq2 ## 概述 PyDESeq2 是 DESeq2 的 Python 实现，用于使用批量 RNA-seq 数据进行差异表达分析。设计和执行从数据加载到结果解释的完整工作流程，包括单因素和多因素设计、具有多重测试校正的 Wald 测试、可选的 apeGLM 收缩以及与 pandas 和 AnnData 的集成。 ## 何时使用此技能该技能应该在以下情况下使用： * 分析大量 RNA-seq 计数数据以了解差异表达 * 比较实验条件之间的基因表达（例如，处理与对照） * 执行考虑批次效应或协变量的多因素设计 * 将基于 R 的 DESeq2 工作流程转换为 Python * 将差异表达分析集成到基于Python的管道中 * 用户提及“DESeq2”、“差异表达”、“RNA-seq 分析”或“PyDESeq2” ## 快速启动工作流程对于想要执行标准差异表达分析的用户：将 pandas 导入为 pd 从 pydeseq2.dds 导入 DeseqDataSet 从 pydeseq2.ds 导入 DeseqStats # 1.加载数据 counts_df = pd.read_csv("counts.csv", index_col=0).T # 转置为样本×基因元数据 = 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_df，元数据=元数据，设计=“〜条件”， refit_cooks=真） dds.deseq2() # 4. 执行统计测试 ds = DeseqStats(dds,对比度=[“条件”,“治疗”,“对照”]) ds.summary() # 5. 访问结果结果= ds.results_df 显着=结果[结果.padj < 0.05] print(f"找到{len(significant)}显着基因") ## 核心工作流程步骤 ### 第 1 步：数据准备 **输入要求：** * **计数矩阵：** 具有非负整数读取计数的样本 × 基因 DataFrame * **元数据：** 样本 × 变量 DataFrame 与实验因素 **常见的数据加载模式：** # 来自CSV（典型格式：基因×样本，需要转置） counts_df = pd.read_csv("counts.csv",index_col=0).T 元数据 = pd.read_csv("metadata.csv",index_col=0) # 来自 TSV counts_df = pd.read_csv("counts.tsv", sep="\t", index_col=0).T # 来自安数据将 anndata 导入为广告 adata = ad.read_h5ad("data.h5ad") counts_df = pd.DataFrame(adata.X, index=adata.obs_names, columns=adata.var_names) 元数据 = 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[要保留的样本] 元数据=metadata.loc[要保留的样本] ### 第 2 步：设计规范设计公式指定了如何对基因表达进行建模。 **单因素设计：** design = "~condition" # 简单的两组比较 **多因素设计：** design = "~batch + condition" # 控制批次效果 design = "~age + condition" # 包含连续协变量 design = "~group + condition + group:condition" # 交互效果 **设计公式指南：** * 使用威尔金森公式符号（R 式） * 将调整变量（例如批次）放在感兴趣的主要变量之前 * 确保变量作为元数据 DataFrame 中的列存在 * 使用适当的数据类型（离散变量的分类） ### 步骤 3：DESeq2 拟合初始化 DeseqDataSet 并运行完整的管道：从 pydeseq2.dds 导入 DeseqDataSet dds = DeseqDataSet( 计数=counts_df，元数据=元数据，设计=“〜条件”， refit_cooks=True, # 去除异常值后重新拟合 n_cpus=1 # 并行处理（根据需要调整）） # 运行完整的 DESeq2 管道 dds.deseq2() **`deseq2()` 的作用：** 1. 计算尺寸因子（标准化） 2. 适合 Genewise 分散体 3.拟合离散趋势曲线 4. 计算色散先验 5. 适合MAP分散体（收缩） 6. 适应日志折叠变化 7. 计算库克距离（异常值检测） 8. 如果检测到异常值则进行重新调整（可选） ### 步骤 4：统计测试执行 Wald 检验来识别差异表达基因：从 pydeseq2.ds 导入 DeseqStats ds = DeseqStats( dds, 对比=[“条件”，“治疗”，“对照”]，#测试治疗与对照 alpha=0.05, # 显着性阈值 cooks_filter=True, # 过滤异常值 Independent_filter=True # 过滤低功耗测试） ds.summary() **对比规格：** * 格式：`[变量，测试级别，参考级别]` * 示例：`["condition", "treatment", "control"]` 测试处理与对照 * 如果“无”，则使用设计中的最后一个系数 **结果数据框列：** * `baseMean`：样本的平均归一化计数 * `log2FoldChange`: Log2 条件之间的折叠变化 * `lfcSE`: LFC的标准误 * `stat`: Wald 检验统计量 * `pvalue`：原始 p 值 * `padj`：调整后的 p 值（通过 Benjamini-Hochberg 进行 FDR 校正） ### 步骤 5：可选的 LFC 收缩应用收缩来减少倍数变化估计中的噪声： ds.lfc_shrink() # 应用 apeGLM 收缩 **何时使用 LFC 收缩：** * 用于可视化（火山图、热图） * 按效应大小对基因进行排序 * 当对后续实验的基因进行优先排序时 **重要提示：** 收缩仅影响 log2FoldChange 值，而不影响统计测试结果（p 值保持不变）。使用缩小的值进行可视化，但报告未缩小的 p 值以显示显着性。 ### 第 6 步：结果导出保存结果和中间对象：进口泡菜 # 将结果导出为 CSV ds.results_df.to_csv("deseq2_results.csv") # 只保存重要基因显着 = ds.results_df[ds.results_df.padj < 0.05] important.to_csv("significant_genes.csv") # 保存 DeseqDataSet 供以后使用将 open("dds_result.pkl", "wb") 作为 f： pickle.dump(dds.to_picklable_anndata(), f) ## 常见分析模式 ### 两组比较标准病例对照比较： dds = DeseqDataSet(counts=counts_df, 元数据=元数据, design="~condition") dds.deseq2() ds = DeseqStats(dds,对比度=[“条件”,“治疗”,“对照”]) ds.summary() 结果= ds.results_df 显着=结果[结果.padj < 0.05] ### 多重比较对照对照测试多个治疗组： dds = DeseqDataSet(counts=counts_df, 元数据=元数据, design="~condition") dds.deseq2() 治疗 = [“治疗_A”，“治疗_B”，“治疗_C”] 所有结果 = {} 用于治疗中的治疗： ds = DeseqStats(dds,对比度=[“条件”,治疗,“对照”]) ds.summary() all_results[治疗] = 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, 元数据=元数据, design="~batch + 条件") dds.deseq2() # 控制批次时测试条件 ds = DeseqStats(dds,对比度=[“条件”,“治疗”,“对照”]) ds.summary() ### 连续协变量包括连续变量，如年龄或剂量： # 确保连续变量是数字元数据[“年龄”] = pd.to_numeric(元数据[“年龄”]) dds = DeseqDataSet(counts=counts_df, 元数据=元数据, design="~年龄 + 条件") dds.deseq2() ds = DeseqStats(dds,对比度=[“条件”,“治疗”,“对照”]) ds.summary() ## 使用分析脚本该技能包括用于标准分析的完整命令行脚本： # 基本用法 python 脚本/run_deseq2_analysis.py \ --counts counts.csv \ --元数据元数据.csv \ --设计“~条件”\ --对比条件处理对照\ --输出结果/ # 带有附加选项 python 脚本/run_deseq2_analysis.py \ --counts counts.csv \ --元数据元数据.csv \ --design "~batch + 条件" \ --对比条件处理对照\ --输出结果/ \ --最小计数 10 \ --阿尔法 0.05 \ --n-CPU 4 \ --绘图 **脚本特点：** * 自动数据加载和验证 * 基因和样本过滤 * 完整的DESeq2管道执行 * 可自定义参数的统计测试 * 结果导出（CSV、pickle） * 可选的可视化（火山图和 MA 图）当用户需要独立的分析工具或想要批量处理多个数据集时，请参考“scripts/run_deseq2_analysis.py”。 ## 结果解释 ### 识别重要基因 # 按调整后的 p 值过滤显着 = 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) ] # 上调和下调分开上调=显着[significant.log2FoldChange > 0] 下调=显着[significant.log2FoldChange < 0] print(f"上调：{len(上调)}") print(f"下调：{len(下调)}") ### 排名和排序 # 按调整后的 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", 升序=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("分数", 升序=False).head(20) ### 质量指标 # 检查标准化（大小因子应接近 1） print("尺寸系数：", dds.obsm["size_factors"]) # 检查离散度估计将 matplotlib.pyplot 导入为 plt plt.hist(dds.varm["分散度"], 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() ## 可视化指南 ### 火山图可视化显着性与效应大小：将 matplotlib.pyplot 导入为 plt 将 numpy 导入为 np 结果= ds.results_df.copy() 结果["-log10(padj)"] = -np.log10(results.padj) plt.figure(figsize=(10, 6)) 显着 = 结果.padj < 0.05 plt.分散（ results.loc[~显着，“log2FoldChange”]， results.loc[~显着, "-log10(padj)"], alpha=0.3, s=10, c='灰色', label='不显着' ） plt.分散（ results.loc[显着, "log2FoldChange"], results.loc[显着, "-log10(padj)"], alpha=0.6，s=10，c='红色'，标签='padj < 0.05' ） plt.axhline(-np.log10(0.05)，颜色='蓝色'，线型='--'，alpha=0.5) plt.xlabel("Log2 折叠变化") plt.ylabel("-Log10(调整后的 P 值)") plt.title("火山图") plt.图例() plt.savefig("volcano_plot.png", dpi=300) ### MA 情节显示倍数变化与平均表达： plt.figure(figsize=(10, 6)) plt.分散（ np.log10(results.loc[~显着, "baseMean"] + 1), results.loc[~显着，“log2FoldChange”]， alpha=0.3, s=10, c='灰色' ） plt.分散（ np.log10(results.loc[显着性, "baseMean"] + 1), results.loc[显着, "log2FoldChange"], alpha=0.6，s=10，c='红色' ） plt.axhline(0, 颜色='蓝色', 线型='--', alpha=0.5) plt.xlabel("Log10(基平均值 + 1)") plt.ylabel("Log2 折叠变化") plt.title("MA 图") plt.savefig("ma_plot.png", dpi=300) ## 常见问题故障排除 ### 数据格式问题 **问题：**“计数和元数据之间的索引不匹配” **解决方案：** 确保样本名称完全匹配 print("样本计数：", counts_df.index.tolist()) print("元数据样本：",metadata.index.tolist()) # 如果需要的话取交集常见 = counts_df.index.intersection(metadata.index) counts_df = counts_df.loc[通用] 元数据=metadata.loc[通用] **问题：**“所有基因的计数为零” **解决办法：** 检查数据是否需要转置 print(f"计算形状：{counts_df.shape}") # 如果基因 > 样本，则需要转置如果 counts_df.shape[1] < counts_df.shape[0]： counts_df = counts_df.T ### 设计矩阵问题 **问题：**“设计矩阵不是满秩的” **原因：** 混杂变量（例如，一批中所有处理过的样品） **解决方案：** 删除混杂变量或添加交互项 # 检查混杂情况打印（pd.crosstab（元数据.条件，元数据.批）） # 要么简化设计，要么添加交互 design = "~condition" # 删除批次 # 或 design = "~condition + batch + condition:batch" # 模型交互 ### 没有显着基因 **诊断：** # 检查色散分布 plt.hist(dds.varm["分散度"], bins=50) plt.show() # 检查大小因素打印（dds.obsm[“尺寸因子”]） # 通过原始 p 值查看顶级基因打印（ds.results_df.nsmallest（20，“pvalue”）） **可能的原因：** * 效应量小 * 生物变异性高 * 样本量不足 * 技术问题（批次效应、异常值） ## 参考文档有关此面向工作流程的指南之外的全面详细信息： * **API 参考** (`references/api_reference.md`)：PyDESeq2 类、方法和数据结构的完整文档。当需要详细参数信息或了解对象属性时使用。 * **工作流程指南** (`references/workflow_guide.md`)：涵盖完整分析工作流程、数据加载模式、多因素设计、故障排除和最佳实践的深入指南。在处理复杂的实验设计或遇到问题时使用。当用户需要时将这些引用加载到上下文中： * 详细的API文档：`阅读references/api_reference.md` * 全面的工作流程示例：`阅读参考文献/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 值。本杰明尼-霍赫伯格程序控制错误发现率。 7. **对比规范：** 格式为“[变量，测试级别，参考级别]”，其中测试级别与参考级别进行比较。 8. **保存中间对象：** 使用pickle保存DeseqDataSet对象以供以后使用或其他分析，而无需重新运行昂贵的拟合步骤。 ## 安装和要求 uv pip 安装 pydeseq2 **系统要求：** * Python 3.10-3.11 * 熊猫 1.4.3+ * numpy 1.23.0+ * scipy 1.11.0+ * scikit学习1.1.1+ * 安达 0.8.0+ **可视化可选：** * matplotlib * 海上生 ## 其他资源 * **官方文档：** * **GitHub 存储库：** * **出版物：** Muzellec 等人。（2023）生物信息学，DOI：10.1093/生物信息学/btad547 * **原始 DESeq2 (R)：** Love 等人。（2014）基因组生物学，DOI：10.1186/s13059-014-0550-8

pydeseq2

包含技能