PyDESeq2 差分表达分析工作流完全指南从数据预处理到结果可视化与排错实战【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills本篇指南以当前仓库中 PyDESeq2 技能skill的完整工作流文档为核心系统讲解基于 pydeseq2 的 bulk RNA-seq 差分表达分析Differential Expression Analysis, DEA全流程。全文覆盖 12 步双阶段分析框架、数据加载与过滤、单因素与多因素实验设计、Wald 检验与多重检验校正、LFC 收缩、结果导出与火山图/MA 图绘制并结合仓库内的命令行脚本 run_deseq2_analysis.py 与测试用例 test_scripts.py 给出源码级佐证。读完本文你将能够独立完成从 counts.csv 到显著差异基因列表的完整分析闭环并掌握常见错误的定位与修复方法。一、分析框架总览12 步双阶段流程一个标准的 PyDESeq2 分析由两个阶段、12 个主要步骤构成理解这条主线有助于你把握每一步代码在整体流程中的位置。阶段一Read Counts 建模第 1–7 步归一化size factor 估计与离散度dispersion估计Log 折叠变化log fold-change, LFC拟合离群值检测基于 Cooks distance阶段二统计分析第 8–12 步Wald 检验多重检验校正Benjamini-Hochberg FDR 控制可选的 LFC 收缩apeGLM shrinkage仓库的核心步骤文档 core_workflow_steps.md 将deseq2()方法内部的拟合过程进一步细化为 8 个子步骤与上文框架一一对应计算 size factors归一化拟合逐基因离散度genewise dispersions拟合离散度趋势曲线dispersion trend curve计算离散度先验dispersion priors拟合 MAP 离散度即离散度收缩拟合 log fold changes计算 Cooks distance离群值检测若检测到离群值则可选地重新拟合refit_cooksTrue时而命令行脚本 run_deseq2_analysis.py 在执行dds.deseq2()前会打印这 7 个主要阶段的进度提示size factors → genewise dispersions → dispersion trend curve → dispersion priors → MAP dispersions → log fold changes → Cooks distances运行时可据此判断拟合进度。二、完整工作流代码可复制运行以下是一段覆盖数据加载、过滤、拟合、检验、收缩与结果访问的端到端代码是全篇文章的“总纲”后续各节会对其中的每个环节逐一深入import pandas as pd from pydeseq2.dds import DeseqDataSet from pydeseq2.default_inference import DefaultInference from pydeseq2.ds import DeseqStats # 加载数据 counts_df pd.read_csv(counts.csv, index_col0).T # 按需转置为 samples × genes metadata pd.read_csv(metadata.csv, index_col0) # 过滤低表达基因总计数 10 的基因被剔除 genes_to_keep counts_df.columns[counts_df.sum(axis0) 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] # 初始化 DeseqDataSet显式指定参考水平 metadata[condition] pd.Categorical( metadata[condition], categories[control, treated] ) inference DefaultInference(n_cpus4) dds DeseqDataSet( countscounts_df, metadatametadata, design~condition, refit_cooksTrue, inferenceinference, ) # 运行归一化与拟合 dds.deseq2() # 执行统计检验treated vs control ds DeseqStats( dds, contrast[condition, treated, control], alpha0.05, cooks_filterTrue, independent_filterTrue, inferenceinference, ) ds.summary() # 可选为可视化应用 LFC 收缩 ds.lfc_shrink(coeffcondition[T.treated]) # 访问结果 results ds.results_df print(results.head())几点关键说明metadata[condition] pd.Categorical(..., categories[control, treated])中categories列表的第一个元素control即参考水平reference level对比检验都以此为基准DefaultInference(n_cpus4)是可并行化的推断后端建议按机器可用核心数调整PyDESeq2 0.5.x 已不再支持“默认 contrast”必须显式传入contrast见 api_reference.md 中的版本兼容说明。三、数据加载与预处理3.1 从 CSV / TSV 加载计数数据常见的落盘格式是genes × samples基因行为而 PyDESeq2 要求samples × genes样本行为、基因列为列因此通常需要转置import pandas as pd # 加载计数矩阵genes × samples counts_df pd.read_csv(counts.csv, index_col0) # 转置为 samples × genes counts_df counts_df.T # 加载元数据本身就是 samples × variables 格式 metadata pd.read_csv(metadata.csv, index_col0)从 TSV 加载只需指定分隔符counts_df pd.read_csv(counts.tsv, sep\t, index_col0).T metadata pd.read_csv(metadata.tsv, sep\t, index_col0)3.2 从 AnnData / H5AD 加载import anndata as ad adata ad.read_h5ad(counts_and_metadata.h5ad) counts_df pd.DataFrame(adata.X, indexadata.obs_names, columnsadata.var_names) metadata adata.obs安全提醒不要从未受信任的来源加载 pickle 文件pickle 反序列化可能执行任意代码。跨工具、跨 Agent、跨协作者传递数据时优先使用 CSV/TSV 或.h5ad这类可移植格式。3.3 数据过滤过滤低计数基因总 reads 阈值即剔除阈值默认 10# 移除总 reads 少于 10 的基因 genes_to_keep counts_df.columns[counts_df.sum(axis0) 10] counts_df counts_df[genes_to_keep]过滤元数据缺失的样本# 移除 condition 列为 NA 的样本 samples_to_keep ~metadata.condition.isna() counts_df counts_df.loc[samples_to_keep] metadata metadata.loc[samples_to_keep]多条件组合过滤# 只保留同时满足所有条件的样本 mask ( ~metadata.condition.isna() (metadata.batch.isin([batch1, batch2])) (metadata.age 18) ) counts_df counts_df.loc[mask] metadata metadata.loc[mask]3.4 数据校验在进入模型拟合前务必确认数据形态与取值合法print(fCounts shape: {counts_df.shape}) # 应为 (samples, genes) print(fMetadata shape: {metadata.shape}) # 应为 (samples, variables) print(fIndices match: {all(counts_df.index metadata.index)}) # 检查负值负计数在负二项模型中没有意义 assert (counts_df 0).all().all(), Counts must be non-negative # 检查非整数值DESeq2 需要原始整数 read counts assert counts_df.applymap(lambda x: x int(x)).all().all(), Counts must be integers源码级印证命令行脚本的load_and_validate_data()函数见 run_deseq2_analysis.py把这一逻辑固化为自动化检查用counts_df.index.equals(metadata.index)做精确对齐判断该函数注释明确指出Index.equals同时覆盖长度不一致的情况而逐元素比较index_a index_b在长度不同时直接抛错会导致后续取交集的分支永远走不到检测到负值时抛出ValueError(Count matrix contains negative values)索引不完全匹配时自动取两个 DataFrame 的交集并对齐。对应测试 test_scripts.py 中的LoadAndValidateTests进一步验证了这些边界行为负计数被拒绝test_negative_counts_are_rejected、零计数合法test_zero_counts_are_accepted、元数据与计数索引乱序时会被重新对齐而不是错位匹配test_a_reordered_metadata_index_is_realigned_not_left_shuffled。四、单因素分析Single-Factor Analysis4.1 简单的两组比较最常见的情景处理组 vs 对照组。设计公式~condition表示“基因表达仅由 condition 解释”from pydeseq2.dds import DeseqDataSet from pydeseq2.default_inference import DefaultInference from pydeseq2.ds import DeseqStats # 设计把表达建模为 condition 的函数 inference DefaultInference(n_cpus4) dds DeseqDataSet( countscounts_df, metadatametadata, design~condition, inferenceinference, ) dds.deseq2() # 检验 treated vs control ds DeseqStats( dds, contrast[condition, treated, control], inferenceinference, ) ds.summary() # 结果 results ds.results_df significant results[results.padj 0.05] print(fFound {len(significant)} significant genes)contrast[condition, treated, control]的语义是“比较 condition 变量下 treated 相对 control 的变化”格式固定为[变量, 检验水平, 参考水平]。注意显著性判定应使用校正后的padj 0.05而非原始pvalueBenjamini-Hochberg 程序控制错误发现率参见 SKILL.md 的 Key Reminders。4.2 多个两两比较当存在多个处理组如 treated_A、treated_B、treated_C时可以对同一个dds反复构造DeseqStats分别与 control 对比# 每个处理组分别 vs control treatments [treated_A, treated_B, treated_C] all_results {} for treatment in treatments: ds DeseqStats( dds, contrast[condition, treatment, control] ) ds.summary() all_results[treatment] ds.results_df # 横向对比各组显著基因数 for name, results in all_results.items(): sig results[results.padj 0.05] print(f{name}: {len(sig)} significant genes)注意多次调用DeseqStats并各自summary()时每次检验都基于同一个拟合好的dds无需重复执行deseq2()。此模式在 analysis_patterns.md 的 “Multiple Comparisons” 一节中也有等价实现。五、多因素分析Multi-Factor Analysis5.1 双因素设计控制批次效应当样本来自不同批次batch时把batch作为调整变量放入设计公式再检验 condition 的效应# 设计公式同时包含 batch 与 condition dds DeseqDataSet( countscounts_df, metadatametadata, design~batch condition ) dds.deseq2() # 在控制 batch 的前提下检验 condition 效应 ds DeseqStats( dds, contrast[condition, treated, control] ) ds.summary()5.2 交互效应Interaction Effects当需要检验“处理效应是否在不同分组间有差异”时即交互项例如某药物只在特定基因型下有效使用冒号语法加入交互项# 设计包含交互项 group:condition dds DeseqDataSet( countscounts_df, metadatametadata, design~group condition group:condition ) dds.deseq2() # 使用与设计矩阵列数严格对应的 numpy 对比向量检验交互项 print(dds.obsm[design_matrix].columns) interaction_contrast_vector ... # 例如 np.array([...])长度等于设计矩阵列数 ds DeseqStats(dds, contrastinteraction_contrast_vector) ds.summary()交互项通常无法用[变量, 检验水平, 参考水平]三元组表达因此这里需要构造显式的数值对比向量。关键技巧先通过dds.obsm[design_matrix].columns查看 PyDESeq2 实际构建的设计矩阵列名例如[Intercept, group[T.g2], condition[T.treated], group[T.g2]:condition[T.treated]]再按列生成向量长度必须与设计矩阵列数完全一致。5.3 连续协变量Continuous Covariates年龄、剂量等连续变量可直接放入设计公式前提是元数据中该列确为数值类型# 确保 age 在元数据中是数值型 metadata[age] pd.to_numeric(metadata[age]) dds DeseqDataSet( countscounts_df, metadatametadata, design~age condition ) dds.deseq2()关于连续变量的自动识别core_workflow_steps.md 指出连续变量由设计公式自动检测分类变量可通过C(variable)语法或 pandas 的 category dtype 强制指定。当前仓库的 PyDESeq2 0.5.x 已弃用design_factors、continuous_factors、ref_level等旧构造参数新工作流一律使用 formulaic/Wilkinson 风格的设计字符串。六、结果导出与可视化6.1 保存结果导出为 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) # 按 padj 排序后保存 sorted_results ds.results_df.sort_values(padj) sorted_results.to_csv(sorted_results.csv)保存 DeseqDataSet 为可移植 AnnData# 保存为 AnnData/H5AD 供后续检查 dds.to_picklable_anndata().write_h5ad(dds_result.h5ad)读取已保存的结果# 读取结果 results pd.read_csv(deseq2_results.csv, index_col0) # 读取 AnnData import anndata as ad adata ad.read_h5ad(dds_result.h5ad)命令行脚本的save_results()函数run_deseq2_analysis.py把上述输出固化为四个标准产物deseq2_results.csv全量结果、significant_genes.csvpadj 0.05、results_sorted_by_padj.csv按 padj 升序、deseq_dataset.h5ad可移植数据集。对应测试SaveResultsTests见 test_scripts.py还验证了两个细节padj 恰好等于 0.05 的基因不计入显著显著性判定是严格的小于号以及 padj 为 NaN 的被过滤基因在排序时必须排到末尾而非最前。6.2 基础可视化火山图Volcano plot横轴为 log2 折叠变化纵轴为-log10(padj)直观展示显著性 vs 效应量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)) plt.scatter( results.log2FoldChange, results[-log10(padj)], alpha0.5, s10 ) plt.axhline(-np.log10(0.05), colorred, linestyle--, labelpadj0.05) plt.axvline(1, colorgray, linestyle--) plt.axvline(-1, colorgray, linestyle--) plt.xlabel(Log2 Fold Change) plt.ylabel(-Log10(Adjusted P-value)) plt.title(Volcano Plot) plt.legend() plt.savefig(volcano_plot.png, dpi300)MA 图展示折叠变化 vs 平均表达量可用于观察低表达基因的 LFC 波动plt.figure(figsize(10, 6)) plt.scatter( np.log10(results.baseMean 1), results.log2FoldChange, alpha0.5, s10, c(results.padj 0.05), cmapbwr ) plt.xlabel(Log10(Base Mean 1)) plt.ylabel(Log2 Fold Change) plt.title(MA Plot) plt.savefig(ma_plot.png, dpi300)实现细节仓库脚本的create_plots()函数run_deseq2_analysis.py在绘制前先用-np.log10(results.padj.fillna(1))处理被独立过滤independent filtering剔除的基因——它们的 padj 是 NaN若不填充则-log10(NaN)仍是 NaNmatplotlib 会静默丢弃这些点。测试 test_scripts.py 专门验证了“padj 全为 NaN 时火山图仍能正常生成”这一场景。七、常见模式与最佳实践7.1 数据预处理清单运行 PyDESeq2 之前逐项核对计数为非负整数矩阵方向为 samples × genescounts 与 metadata 的样本名完全一致缺失的元数据值已处理或移除已过滤低计数基因通常为总计数 10实验因素已正确编码分类变量用 category dtype 或C()语法7.2 设计公式最佳实践变量顺序很重要调整变量放在关注变量之前。# 正确先控制 batch再检验 condition design ~batch condition # 不够理想condition 排在前面 design ~condition batch离散变量使用分类类型metadata[condition] metadata[condition].astype(category) metadata[batch] metadata[batch].astype(category)使用 formulaic 设计语法PyDESeq2 0.5.x 推荐design ~batch condition # 避免使用已弃用的构造参数 # design_factors, continuous_factors, ref_level7.3 统计检验准则设置合适的 alpha显著性阈值# 标准显著性阈值 ds DeseqStats(dds, contrast[condition, treated, control], alpha0.05) # 探索性分析可用更严格的阈值 ds DeseqStats(dds, contrast[condition, treated, control], alpha0.01)启用独立过滤independent filtering过滤低检验功效的基因可以提升整体检验功效建议默认开启# 推荐过滤低功效检验 ds DeseqStats(dds, contrast[condition, treated, control], independent_filterTrue) # 仅在有特定理由时关闭 ds DeseqStats(dds, contrast[condition, treated, control], independent_filterFalse)DeseqStats的完整参数lfc_null、alt_hypothesis等可查阅 api_reference.md。其中lfc_null指定零假设下的 log2 折叠变化默认 0.0alt_hypothesis支持阈值化检验greaterAbs、lessAbs、greater、less。7.4 LFC 收缩Shrinkage的使用时机何时使用可视化火山图、热图按效应量排序基因为后续验证实验筛选候选基因何时不使用报告统计显著性应使用未收缩的 p 值基因集富集分析通常使用未收缩值# 同时保存收缩前后两个版本 ds.results_df.to_csv(results_unshrunken.csv) ds.lfc_shrink(coeffcondition[T.treated]) ds.results_df.to_csv(results_shrunken.csv)lfc_shrink(coeff, adaptTrue)使用 apeGLM 方法见 api_reference.mdcoeff必须与dds.obsm[design_matrix]中的实际列名一致如condition[T.treated]adapt控制是否从 MLE 估计自适应先验尺度。收缩只改变log2FoldChange列p 值与 padj 保持不变——因此收缩值仅用于可视化/排序显著性判定始终基于未收缩的统计检验结果。源码级印证脚本的infer_shrink_coeff()函数run_deseq2_analysis.py会根据 contrast 自动拼接f{变量}[T.{检验水平}]然后到真实的设计矩阵dds.obsm[design_matrix].columns中验证该系数是否存在若不存在例如把参考水平当成了检验水平则抛出带全部可用列名的 ValueError并提示用--shrink-coeff显式指定或--no-shrink跳过。测试 test_scripts.py 专门覆盖了这一行为。7.5 内存管理对于大规模数据集# 使用并行处理 inference DefaultInference(n_cpus4) dds DeseqDataSet( countscounts_df, metadatametadata, design~condition, inferenceinference, # 根据可用核心数调整 ) # 必要时分批处理 # 把基因拆成若干子集分别分析最后合并结果内存优化相关的更多参数还包括DeseqDataSet(low_memoryTrue)用后即弃中间结构以及fit_typeparametric/mean、size_factors_fit_typeratio/poscounts/iterative等建模参数完整说明见 api_reference.md。八、排错指南Troubleshooting8.1 报错counts 与 metadata 索引不匹配KeyError: Sample names in counts and metadata dont match排查与解决# 检查索引 print(Counts samples:, counts_df.index.tolist()) print(Metadata samples:, metadata.index.tolist()) # 需要时对齐 common_samples counts_df.index.intersection(metadata.index) counts_df counts_df.loc[common_samples] metadata metadata.loc[common_samples]8.2 报错所有基因的总计数为零ValueError: All genes have zero total counts通常是数据方向问题忘了转置。基因数大于样本数时基本可以判定需要转置# 检查数据方向 print(fCounts shape: {counts_df.shape}) # 如果基因 样本很可能需要转置 if counts_df.shape[1] counts_df.shape[0]: counts_df counts_df.T8.3 警告大量基因被过滤掉先观察基因计数的分布再决定是否调整阈值# 查看基因总计数分布 print(counts_df.sum(axis0).describe()) # 可视化 import matplotlib.pyplot as plt plt.hist(counts_df.sum(axis0), bins50, logTrue) plt.xlabel(Total counts per gene) plt.ylabel(Frequency) plt.show()必要时放宽过滤阈值genes_to_keep counts_df.columns[counts_df.sum(axis0) 5]注意测试 test_scripts.py 验证了阈值语义过滤判定是含边界即总计数恰好等于 10 的基因会被保留。8.4 报错设计矩阵不满秩not full rank典型原因是混淆设计confounded design例如所有处理组样本都来自同一个批次。排查# 检查设计混淆交叉表可直观看出 batch 与 condition 是否完全对应 print(pd.crosstab(metadata.condition, metadata.batch))解决要么删除混淆变量要么加入交互项design ~condition # 去掉 batch # 或 design ~condition batch condition:batch # 加入交互项8.5 没有找到显著基因可能原因效应量太小、生物学变异过高、样本量不足、技术问题批次效应、离群值。诊断步骤# 检查离散度估计 import matplotlib.pyplot as plt dispersions dds.var[dispersions] plt.hist(dispersions, bins50) plt.xlabel(Dispersion) plt.ylabel(Frequency) plt.show() # 检查 size factors应接近 1 print(Size factors:, dds.obs[size_factors]) # 即使不显著也看看按原始 p 值排名最靠前的基因 top_genes ds.results_df.nsmallest(20, pvalue) print(top_genes)额外诊断来自 SKILL.md 的 Result Interpretation 一节绘制 p 值分布直方图健康的结果应“大部分平坦、在 0 附近有尖峰”padj与|log2FoldChange|组合打分排序score -log10(padj) * abs(log2FoldChange)也有助于在弱信号中优先候选基因。8.6 大数据集上的内存错误解决方案按优先级# 1. 减少 CPU 核数看似矛盾但有时有效 inference DefaultInference(n_cpus1) dds DeseqDataSet(..., inferenceinference) # 2. 更激进地过滤 genes_to_keep counts_df.columns[counts_df.sum(axis0) 20] # 3. 分批处理 # 按基因子集拆分分析最后合并结果九、命令行脚本一键完成全流程除逐段编写 Python 代码外仓库还提供了开箱即用的 CLI 脚本 run_deseq2_analysis.py它整合了数据加载校验、过滤、拟合、检验、收缩、导出与可视化# 基本用法 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 \ --shrink-coeff condition[T.treated] \ --plots完整参数表与脚本argparse定义一致参数类型默认值说明--counts必填—计数矩阵 CSV 路径--metadata必填—元数据 CSV 路径--design必填—设计公式如~condition--contrast必填—对比规格变量 检验水平 参考水平三个值--output可选results输出目录--min-counts可选10基因过滤的最小总计数阈值--alpha可选0.05显著性阈值--no-transpose开关关闭若输入已是 samples × genes 则跳过转置--no-shrink开关关闭跳过 LFC 收缩--shrink-coeff可选自动推断要收缩的设计矩阵系数如condition[T.treated]--n-cpus可选1并行 CPU 核数--plots开关关闭生成火山图与 MA 图脚本会打印完整的分析摘要测试基因总数、显著基因数、上下调基因数、Top 10 显著基因方便直接判断结果质量。十、安装与运行环境本仓库的 PyDESeq2 技能面向 0.5.x 版本推荐使用 uv 安装uv pip install pydeseq20.5.4系统与依赖要求依据 SKILL.md 与 api_reference.mdPython 3.11PyDESeq2 0.5.3 起已放弃 Python 3.10 支持PyDESeq2 0.5.4pandas 2.2.0numpy 2.0.0scipy 1.12.0scikit-learn 1.4.0anndata 0.11.0formulaic 1.0.2 与 formulaic-contrasts 0.2.0可视化可选matplotlib、seaborn关于端到端正确性的验证测试 test_scripts.py 中的EndToEndFitTests构造了一个已知答案的合成数据集——将某个基因在处理组中人为放大 8 倍即真实 log2 折叠变化为 3随后验证(1) 推断出的收缩系数确实存在于 PyDESeq2 实际构建设计矩阵的列中(2) 拟合恢复出的 LFC 在 log2(8)3 附近允许 0.5 的收缩偏移(3) 该“种植”基因是全表唯一显著的基因(4) 每个被测基因在结果表中都有对应行。这一测试意味着本文中的工作流在“数据方向 → 对齐 → 阈值 → 收缩系数 → 显著性判定”每一个环节都有自动化回归保障可以直接作为你自己分析管线的正确性参照。【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考