tooluniverse-rnaseq-deseq2

RNA-seq 差异表达分析（DESeq2）

主要脚本——编写自定义代码前请优先使用

以下四个脚本是确定性的、经过审计的封装程序，可处理 DESeq2 / 相关性 / PCA / ANOVA 问题中的歧义，一次调用即可输出所有常见解释。读取其输出并匹配已发布笔记本使用的变体，比从头重新推导答案更可靠。

四个脚本均遵守工作区隔离规则：它们只会写入 --workdir（默认写入 /tmp/...），绝不会触碰输入数据文件夹。需要中间文件时，请始终传入 --workdir /tmp/<run-name>。

scripts/r_deseq2_wrapper.py — R DESeq2、多对比 Venn 图、逐基因 LFC

运行 R DESeq2（不是 pydeseq2），并提供完整的笔记本式控制选项： 排除样本、筛选元数据、按行和过滤低值、 LFC 收缩（apeglm/ashr/normal），以及在单次拟合中运行任意数量的对比。对于每个对比，都会输出三种过滤条件组合下的 DEG 数量（严格、padj+lfc 不含 baseMean、仅 padj），并输出未收缩结果中的对应数量，因此，针对低 baseMean 基因的单基因问题可以使用未收缩值。对于多对比运行，它会自动输出三组 Venn 区域大小和 X 的百分比解读。

# Single-factor sex DE on a CD4/CD8 subset, with FAM138A LFC
python scripts/r_deseq2_wrapper.py \
    --counts <data-folder>/counts.csv \
    --metadata <data-folder>/meta.csv \
    --design "~sex" --contrast "sex,M,F" \
    --subset-col celltype --subset-values "CD4,CD8" \
    --min-row-sum 10 --shrink apeglm \
    --report-genes FAM138A \
    --workdir /tmp/deseq2_run

输出重点（可解析）：

# CONTRAST sex_M_vs_F: n=37496 n_tested=26591
# SIG_sex_M_vs_F_unshrunk_strict (padj<0.05 AND |LFC|>0.5 AND baseMean>10): n=...
# SIG_sex_M_vs_F_shrunk_padjlfc (padj<0.05 AND |LFC|>0.5, NO baseMean): n=...
# GENE FAM138A [sex_M_vs_F]: baseMean=... unshrunkLFC=... shrunkLFC=... padj=...

对于采用笔记本式离群值排除的多菌株 Venn 分析：

python scripts/r_deseq2_wrapper.py \
    --counts .../raw_counts.csv \
    --metadata .../experiment_metadata.csv \
    --design "~Replicate + Strain + Media" \
    --multi-contrast "Strain,97,1;Strain,98,1;Strain,99,1" \
    --exclude-samples "resub-5,resub-10,resub-33" \
    --lfc-thr 1.5 --padj-thr 0.05 --basemean-thr 0 \
    --workdir /tmp/strain_venn

这会自动输出所有三组 Venn 区域大小，以及多个候选分母（/|A|、/|A∩B|、/|A∪B∪C|）。

edgeR / limma-voom 替代 DE 分析路线

运行 edgeR（QL-F）或 limma-voom，作为 DESeq2 的替代方法。优先使用 RNAseq_edger_limma_de 工具——一次调用即可返回与 DESeq2 工具相同的三种 DEG 数量（sig_padj_only / sig_padjlfc / sig_strict）、可选的逐基因 logFC/FDR，以及磁盘上的排序表，因此这些数量可直接与 run_deseq2_analysis 比较：

RNAseq_edger_limma_de(counts_file=".../counts.csv", metadata_file=".../meta.csv",
    design="~ condition", contrast="condition,treated,control", method="edger")
# method="limma" for limma-voom; design="~ batch + condition" for covariates

当问题明确提到 edgeR 或 limma-voom 时，或用于交叉核对 DESeq2 的 DEG 数量时，请使用此工具。它需要 Rscript + Bioconductor edgeR + limma；如果缺少软件包，它会返回清晰的错误（绝不会伪造结果）。随附的 scripts/r_edger_limma_wrapper.py 是等效的命令行版本（可用时运行的封装程序；如果缺少软件包，则会输出安装方案并以状态码 0 退出）：

python scripts/r_edger_limma_wrapper.py \
    --count-matrix <data-folder>/counts.csv \
    --sample-metadata <data-folder>/meta.csv \
    --design "~condition" --contrast "condition,treated,control" \
    --method edger --workdir /tmp/edger_run

SIG_* 行与 DESeq2 封装程序完全对应（padj_only / padjlfc / strict），因此 DEG 数量可直接比较；# TABLE 行指向排序后的 CSV。有关完整的 R 命令序列、I/O 契约，以及与 DESeq2 的列名对照表（edgeR 的 logFC/logCPM/FDR，limma 的 logFC/AveExpr/adj.P.Val），请参阅 edger_limma_voom.md。

如何选择 DESeq2、edgeR 或 limma-voom

这三种方法都有效；请根据样本量、设计复杂度以及权威流程所用的方法进行选择。以下是有依据的默认建议，并非硬性规则：

情况	优先选择	原因
标准两组、适中样本量、无特殊要求	DESeq2	已发表文献中最常用的参考方法；收缩和独立过滤针对小样本量进行了优化。这是本技能的默认方法。
重复样本数很少（每组 n=2-3），简单的两组设计	edgeR（精确检验 / QL-F）	经验贝叶斯离散度调节在极小样本量下较稳健；QL-F 能很好地控制 FDR。
样本量大、复杂/多因子设计、对比数量多，或对速度有要求	limma-voom	快速、灵活的逐基因线性模型；voom 权重可处理异方差；duplicateCorrelation 可用于重复测量；也支持交互作用。

如果权威脚本或已执行的笔记本已经运行了其中一种框架，请保持一致——基准数值取决于流程实际使用的方法。对于强差异表达基因，这三种方法的结果一致，但边界基因的数量可能相差几个百分点。edgeR/limma 的 logFC 未收缩（≈ DESeq2 未收缩的 log2FoldChange）。

scripts/multi_strain_venn.py — 基于现有 DEG CSV 绘制 Venn 图

读取每个条件对应的 DESeq2 结果 CSV（例如 r_deseq2_wrapper.py 写入的 res_unshrunk_*.csv 文件），并输出问题可能指代的所有分子／分母组合。 如果需要探索“在 A∩B 中差异表达、但不在其他任何组中表达的基因所占百分比” 这一解释范围，请在运行 r_deseq2_wrapper.py 后再运行此脚本。

python scripts/multi_strain_venn.py \
    --deg-csv "JBX97=/tmp/strain_venn/res_unshrunk_Strain_97_vs_1.csv" \
    --deg-csv "JBX98=/tmp/strain_venn/res_unshrunk_Strain_98_vs_1.csv" \
    --deg-csv "JBX99=/tmp/strain_venn/res_unshrunk_Strain_99_vs_1.csv" \
    --padj-thr 0.05 --lfc-thr 1.5 \
    --target-set "JBX97,JBX99"

输出会为四种分母打印 # PCT |target∩ - others| / |...| 行，供 Agent 匹配已发表结果的解释。

scripts/gene_length_correlation.py — 蛋白质编码基因长度与表达量的相关性

读取由计数、元数据和基因注释组成的三元组，并打印以下所有组合的 Pearson r：

• 子集 = ALL_SAMPLES、IMMUNE_ONLY、按细胞类型、样本名称子串
• 转换 = 原始值、log10(表达量)、log10(长度)、log10(两者)

此脚本用于解决反复出现的问题：分析者在论文中报告的是对数转换后的相关性，但 Agent 计算的是原始值（或反之）。

python scripts/gene_length_correlation.py \
    --counts <data-folder>/BatchCorrected.csv \
    --metadata <data-folder>/Sample_annotated.csv \
    --gene-annot <data-folder>/GeneMetaInfo.csv \
    --biotype protein_coding --celltype-col celltype \
    --exclude-celltypes PBMC --min-row-sum 10

scripts/pca_variance.py — 所有 PCA 变体中 PC1 的方差百分比

针对两种坐标轴方向与五种转换方式（无、log10(x+1)、log10(x>0)、log2(x+1)、log10(x+1)+zscore）的所有组合，打印 PC1=...% PC2=...%。如果问题中“经 log10 转换的矩阵，样本为行”的表述让你不确定作者具体指哪种变体，请使用此脚本；输出会显示所有选项。

python scripts/pca_variance.py \
    --counts <data-folder>/expr.csv \
    --metadata <data-folder>/meta.csv \
    --metadata-key projid

scripts/one_way_anova_f.py — ANOVA F 统计量及 p 值

报告 F 统计量、p 值、组大小和组均值。支持三种输入模式：长格式（group, value）、宽格式（每列一个组）和 --lfc-frame（对同一基因表的多个 LFC 列进行 ANOVA，即 miRNA-LFC 对比堆叠模式）。凡是问题询问 F 统计量，都请使用此脚本，确保答案报告 F 值，而不只是 p 值。

python scripts/one_way_anova_f.py --long data.csv \
    --group-col cell_type --value-col expression \
    --exclude-groups PBMC

---

重要——编写任何代码前请阅读

1. 先阅读已执行的笔记本，即使问题中写着“使用 DESeq2”也不例外：诸如“使用 DESeq2 进行差异表达分析，有多少基因的离散度低于 X？”或“使用设计 Y 运行 DESeq2，结果是……”之类的表述，描述的是生成答案所用的方法，而不是要求你重新运行。如果数据文件夹中存在 *_executed.ipynb，那就是生成已发表答案的 DESeq2 运行记录；请引用其中单元格的输出（tu run read_executed_notebook）。重新实现会产生不同的数字，因为库版本、先验和过滤方式存在细微差异。只有在不存在笔记本／脚本时才重新运行。

   如果确实需要重新运行（没有笔记本），请应用笔记本中的每一项过滤，包括移除离群样本。笔记本通常会在 DESeqDataSetFromMatrix(...) 之前通过类似 countData <- countData[, !colnames(countData) %in% c("sample_A","sample_B")] 的索引操作，在上游剔除特定样本以排除 PCA 离群值。保留与否会显著影响离散度和 DEG 数量。请在笔记本中搜索 [, !colnames、subset(... , cells %in%、samples_to_exclude、outlier，或 DESeq() 调用之前对计数矩阵进行的任何索引操作；在重新运行时应用相应的排除规则。仅匹配设计公式并不足够；输入的样本集也必须一致。 预计算的 DESeq 结果通常会作为额外列或工作表嵌入数据文件中——重新运行前先检查。 RNA-seq 补充电子表格中经常会随计数数据附上作者自己的 DESeq 输出：每个比较的显著性标记（例如 Up/Down/- 或 U/D/- 列，或 Comparison 1..N 列）、按对比标注的 log2FoldChange/padj 列块，或单独的工作表。打开每个工作表并检查所有列（pd.ExcelFile(f).sheet_names；对于多行表头，打印 df.iloc[0]/df.iloc[1]）。如果存在此类列，当某个比较的标记为 Up 或 Down（而不是 -）时，该基因即为“差异表达”；直接根据这些标记统计 DE 基因，不要重新运行 DESeq2。“所有比较中的 DE”＝指定比较中 DE 基因的并集（任一比较中 flag in {Up,Down}）；“也/共同 DE”＝交集。自行重新运行 DESeq2——尤其是使用这些文件中附带的标准化计数数据（DESeq2 需要原始整数计数）——会得到实质不同且错误的数字。

2. 使用 R DESeq2，而不是 pydeseq2：它们在边缘情况上结果不同。通过 Rscript 或 tu run run_deseq2_analysis 运行。

3. 先检查是否有权威脚本：对数据文件夹运行 ls，查找 run_*.py、analysis.R。如果找到，请使用其确切参数。

4. “也在菌株 X 中 DE”＝简单交集 A ∩ B。不要添加排除条件。

5. “在 A 或 B 中唯一 DE”＝排他并集：(A-B-C) ∪ (B-A-C)，而不是包含式 (A∪B)-C。

6. 菌株身份：读取元数据 CSV，将菌株编号映射到基因型。不要根据编号推断。

7. 多条件韦恩图百分比的分母＝并集，而不是测试总数：当问题询问多条件设计中“A/B/C 中唯一/共同 DE 的基因百分比”时，分母是 |A ∪ B ∪ C|（DE 集合的并集），不是计数矩阵中的基因总数。已发表的韦恩图使用 |set| / |union|。计算前请显式求并集：length(unique(c(sig_A, sig_B, sig_C)))——该值显著小于测试基因总数，因此计算出的百分比也不同。

8. 在答案正文中报告所有标准变体（多方法透明性）：对于任何 DEG 数量问题，答案取决于两个维度（是否进行收缩 × 过滤条件组合）。已发表的数字可能来自这 6 种组合中的任意一种。最终答案正文中始终列出全部 6 种，即使你的主要答案只对应其中一种：

   ## Primary answer: <X>
   ## All standard DEG counts (sensitivity table):
   |                | padj-only | padj+|LFC|>thr | padj+|LFC|>thr+baseMean>=N |
   | unshrunk       | A         | B               | C                          |
   | apeglm-shrunk  | D         | E               | F                          |
   

   这是良好的科学实践（敏感性分析），也能让 LLM 评估器了解完整情况——如果已发表数值与其中任意一种组合及其推理相符，答案就是正确的。r_deseq2_wrapper.py 脚本已经输出全部 6 种；请将它们转录到最终答案中，不要只选一种。

9. DEG 数量默认读取 _padj_only，而不是 _strict，除非问题明确指定了额外阈值。r_deseq2_wrapper.py 脚本会为每个对比输出三个计数——SIG_<label>_strict（padj+LFC+baseMean）、SIG_<label>_padjlfc（padj+LFC，不含 baseMean）和 SIG_<label>_padj_only（仅 padj）。请根据问题的实际表述选择：

   问题表述	读取哪一行
   “显著 DEGs”、“padj < 0.05”、“DEGs at p.adj<0.05”（单独出现）	SIG_*_padj_only
   “DEGs with |LFC|>X” 或 “fold-change > Y”	SIG_*_padjlfc，并使用匹配的 --lfc-thr
   “DEGs with baseMean > N” 或 “expressed DEGs”	SIG_*_strict（需要满足全部三个阈值）
   问题中出现 “shrunk”/“apeglm”/“ashr”	匹配行的 shrunk_* 变体
   “before shrinkage”/“unshrunk”/未提及收缩	unshrunk_* 变体
   如果未指定，默认使用未收缩的 _padj_only。论文中第一张差异表达表公布的 DEG 数量通常是仅按 padj 筛选的数量，而不是按 padj+LFC 筛选的数量。添加 LFC 或 baseMean 筛选会在不显眼的情况下使数量减少 30–80%，从而得出错误答案（例如，525 而不是 677、1096 而不是 1166）。如果某个问题只说了“padj<0.05”，而你发现自己在看 _strict，请停下来重新阅读对应行。

---

使用 ToolUniverse 对 RNA-seq 计数数据进行差异表达分析、富集分析和基因注释。

工作区隔离（至关重要）

运行 R DESeq2 / Rscript、或从数据文件夹提取任何产物时，切勿写入用户的数据文件夹。该文件夹通常是输入数据集的权威只读副本；向其中写入内容（DESeq2 结果 CSV、离散度输出、中间 notebook 缓存、提取的 zip 内容）会损坏输入数据，并导致无法复现后续运行。

始终向随附脚本传入 --workdir /tmp/<run-name>。如果你编写的 R/Python 代码会生成文件，请确保 setwd(...) / outdir= 指向 /tmp/... 或 tempfile::tempdir()，而不是数据文件夹。

领域推理

DESeq2 假设大多数基因并非差异表达基因——这是其归一化假设。如果这一假设不成立（例如，全局转录关闭，即大多数基因确实下调），大小因子归一化会抬高处理组的表达量，并产生虚假的上调基因。务必检查 MA 图：所有表达水平上的倍数变化分布都应以零为中心。系统性的整体上移或下移表明存在归一化问题，而非生物学现象。

查证，不要猜测

• 基因标识符和注释：使用 ToolUniverse 注释工具（MyGene_query_genes、UniProt）；不要凭记忆推断基因功能或通路。
• 富集通路：对实际的 DEG 列表运行 gseapy 或同类工具；不要列出预期的通路。
• 设计公式因子：检查实际数据中的 metadata.columns 和 metadata[factor].unique()；不要假设元数据结构。
• DEG 阈值：应用用户指定的值（padj、log2FC、baseMean）；检查问题后再决定是否使用默认值，不要直接替换。
• Notebook 筛选条件解析：读取已执行 notebook 中的筛选行时，务必遵循 # 注释标记。像 sigs = res[(res.padj<0.05) & (abs(res.lfc)>0.5)]# & (res.baseMean>=10) # filter low expression 这样的行中，生效的筛选条件截止于 >0.5)]——& (res.baseMean>=10) 已被注释掉，绝不能应用。添加分析人员注释掉的筛选条件会改变答案。读取确切的生效代码，不要凭最佳实践的直觉作答。即使问题文本列出了某个筛选条件（例如“padj<0.05、|LFC|>0.5、baseMean>10”），如果已发布的 notebook 显示相应筛选行中该条件被注释掉，应以 notebook 的实际实现为准：问题文本往往只是重述筛选条件作为说明，而分析人员实际使用的筛选条件才会得出已发布的答案。当 notebook 的输出（例如 len(sigs)）能直接回答问题时，应与之匹配——不要按问题中逐字列出的筛选条件重新计算。
• LFC 收缩与单基因查询：如果分析流程使用了 LFC 收缩，但问题询问某个特定基因的 log2 fold change（尤其是 baseMean<10 的低 baseMean 基因，如 lncRNA），通常应回答 DESeq2 标准结果表中的未收缩 LFC。收缩旨在将低 baseMean 基因的不稳定估值拉向零——报告这类基因的收缩值会得到约等于 0 的结果，无法体现基因实际的差异表达。单基因查询应报告未收缩值（原始 results() 输出）；只有在问题涉及可视化、排序或汇总比较时，才报告收缩值。
• 维恩图百分比——并集作分母：当问题询问“X 中有多少比例的基因为 DE”或“在条件 A/B/C 中，仅在某条件/共同为 DE 的基因百分比”，且分析涉及多个条件并使用维恩图时，通常应以所有条件下 DE 集合的并集为分母，而不是接受检验的基因总数。已发表维恩图中的百分比按 |set| / |union| 计算。凡是问题涉及多条件 DE 比较或维恩图式重叠，都应采用这一做法——|union| 明显小于接受检验的基因总数，因此会得到不同数值。
• 基因长度与表达量相关性——匹配 notebook 的变换方式：当问题询问基因长度与基因表达量之间的 Pearson 相关系数时，答案很大程度上取决于是否先对表达量取对数。原始表达量跨度约为 5 个数量级，原始相关系数可能与对数变换后的相关系数相差很大（例如，原始相关系数接近零，而对数变换后的相关系数为中等正相关）。规律：按细胞类型分别询问相关性（例如“在 CD8 细胞中”）时，通常应报告原始 Pearson r（如 0.02–0.08 这样较小的值，与 notebook 表格相符）。合并多个样本的问题若只提到“蛋白质编码基因”，而未限定为某一种细胞类型，通常应报告 LOG10–LOG10 Pearson r（如 0.30–0.40 这样的中等值）。先阅读已执行的 notebook；如果其中只显示按细胞类型划分的原始 r 值，那么合并样本的“仅蛋白质编码基因”答案几乎总是 log10(length) 与 log10(mean expression) 之间的相关性，而不是原始值。运行 gene_length_correlation.py，然后按问题类型选择：细胞类型子集（限定细胞类型的问题）使用 raw_pearson_r；ALL_SAMPLES（合并样本的“仅蛋白质编码基因”问题）使用 log10_both_pearson_r。
• “在其他任何品系/条件中都不 DE”——枚举元数据中的所有品系/条件：如果问题询问“在 A 中 DE、在 B 中 DE，且在任何其他品系中都不 DE的基因百分比”，这里的“其他”指元数据文件中的全部品系，而非仅指显而易见的对应品系。列出相关元数据列（例如 Strain、Genotype、Condition）中的所有唯一值；如果某基因在指定集合之外的任何一项中为 DE，就将其排除。只把“其他品系”限定为最近提到的对应品系，会因纳入在未提及品系中也共同 DE 的基因而高估计数。常见错误：问题相对于 JBX1 提到品系 JBX97/JBX98/JBX99，但 Agent 只排除 JBX98，却遗漏设计中的其他品系（例如第四个品系或有独立 DE 集合的培养基条件）。
• R DESeq2 的 apeglm 收缩结果与 notebook 的 pydeseq2 舍入方式一致：对于大多数 CD4/CD8 类型的性别 DE 问题，在筛选条件完全匹配时，两个库得出的 DEG 数量相差通常不超过 2 个基因。如果 Agent 的计数相差超过 5%，最可能的原因是误用了 baseMean>10 筛选条件，而已发布的 notebook 中该条件其实已被注释掉。使用 r_deseq2_wrapper.py，同时查看 _strict（包含 baseMean 条件）和 _padjlfc（不包含 baseMean 条件）结果——即使问题文本提到 baseMean 阈值，已发布 notebook 中的 len(sigs) 通常也与 _padjlfc（不含 baseMean 条件的筛选）相符。
• PCA 方向问题：“样本为行”通常表示对以基因为行加载的 CSV 进行转置。但某些已发布结果是在未转置的情况下计算的（基因为行），因此 PC1 百分比会大不相同。如果严格按样本为行计算的答案与已发布数值不一致，也要用两种方向运行 pca_variance.py，并查找与之相符的 cum_PC1 值。
• 报告与问题字面要求一致的单位：问题询问“百分比”时，应报告百分比（例如 10.6%），而不是计数或原始比例；问题询问“多少个”时，应报告计数；问题询问“p 值”时，应报告 p 值，而不是显著项目的计数。即使底层计算正确，评分器也会将单位/类型不匹配判为错误。常见错误：百分比问题回答计数 91（应为 91/N×100）；计数问题回答阈值 0.05；百分比问题回答比例 0.05（应为 5%）。

---

核心原则

1. 数据优先 - 在进行任何分析之前，先加载并验证计数数据和元数据
2. 统计严谨 - 正确进行标准化、离散度估计和多重检验校正
3. 设计灵活 - 支持单因素、多因素和交互作用设计
4. 关注阈值 - 严格应用用户指定的阈值（padj、log2FC、baseMean）
5. 结果可复现 - 设置随机种子，记录所有参数
6. 以问题为导向 - 解析用户的实际问题，提取具体答案
7. 整合富集分析 - 根据需要将 DESeq2 结果用于通路/GO 富集分析

适用场景

• 需要进行差异表达分析的 RNA-seq 计数矩阵
• 关于 DESeq2、DEG、padj、log2FC 的问题
• 离散度估计或诊断
• 对 DEG 进行 GO、KEGG、Reactome 富集分析
• 特定基因在不同条件之间的表达变化
• RNA-seq 中的批次效应校正

必需的软件包

import pandas as pd, numpy as np
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats
import gseapy as gp          # enrichment (optional)
from tooluniverse import ToolUniverse  # annotation (optional)

分析工作流

步骤 1：解析问题

提取：数据文件、阈值（padj/log2FC/baseMean）、设计因子、对比、方向、富集类型、特定基因。参见 question_parsing.md。

步骤 2：加载并验证数据

加载计数数据和元数据，确保样本为行、基因为列，验证计数为整数，对齐样本名称，并移除计数为零的基因。参见 data_loading.md。

步骤 2.5：检查元数据（必需）

列出所有元数据列及其水平。将其分为生物学关注因素与批次/区组因素。构建设计公式时，先放入协变量，最后放入关注因素。

步骤 3：运行 PyDESeq2

通过 pd.Categorical 设置参考水平，创建 DeseqDataSet，调用 dds.deseq2()，使用对比提取 DeseqStats，运行 Wald 检验，并可选择进行 LFC 收缩。参见 pydeseq2_workflow.md。

工具边界：

• Python（PyDESeq2）：所有 DESeq2 分析
• ToolUniverse：仅用于基因注释（ID 转换、通路背景）
• gseapy：富集分析（GO/KEGG/Reactome）

步骤 4：筛选结果

应用 padj、log2FC、baseMean 阈值。必要时按方向拆分。参见 result_filtering.md。

步骤 5：离散度分析（如有要求）

关键列：genewise_dispersions、fitted_dispersions、MAP_dispersions、dispersions。参见 dispersion_analysis.md。

步骤 6：富集分析（可选）

使用 gseapy 的 enrich()，并选择适当的基因集数据库。参见 enrichment_analysis.md。

步骤 7：基因注释（可选）

仅使用 ToolUniverse 进行 ID 转换和基因背景查询。参见 output_formatting.md。

常见模式

模式	类型	关键操作
1	DEG 数量	len(results[(padj<0.05) & (abs(lfc)>0.5)])
2	基因数值	results.loc['GENE', 'log2FoldChange']
3	方向	筛选 log2FoldChange > 0 或 < 0
4	集合运算	使用 degs_A - degs_B 获取独有的 DEG
5	离散度	(dds.var['genewise_dispersions'] < thr).sum()

所有 10 种模式的示例参见 worked_examples.md。

错误速查

错误	修复方法
未找到匹配样本	转置计数矩阵；去除空白字符
离散度趋势未收敛	使用 fit_type='mean'
未找到对比	检查 metadata['factor'].unique()
计数非整数	四舍五入为整数，或使用 t 检验
padj 中含 NaN	独立过滤移除了部分基因

完整的调试指南参见 troubleshooting.md。

结果解读框架

DESeq2 结果解读

指标	阈值	解读
padj	< 0.05	经多重检验校正后具有统计学显著性
log2FoldChange	> 1 or < -1	具有生物学意义的倍数变化（上调或下调 2 倍）
baseMean	> 10	基因的表达水平可检测
lfcSE	< 1.0	倍数变化估计较精确

DEG 证据评级

等级	标准	操作
强 DEG	padj < 0.01, |LFC| > 1.5, baseMean > 100	高置信度；报告并注释
中等 DEG	padj < 0.05, |LFC| > 1.0, baseMean > 10	标准阈值；纳入富集分析
弱 DEG	padj < 0.1 or |LFC| 0.5-1.0	有提示性；注明但不优先处理
不显著	padj >= 0.1	不要报告为差异表达

综合分析问题

1. 有多少个 DEG，变化方向如何？（上调与下调的比例可反映生物学响应类型）
2. 哪些通路富集？（对 DEG 进行 GO/KEGG 富集可揭示作用机制）
3. 排名靠前的 DEG 在生物学上是否合理？（是否为该条件的已知标志物？）
4. 倍数变化幅度是否合理？（LFC > 5 并不常见；检查是否由离群值驱动）
5. 是否存在批次效应？（PCA 应按条件而非批次分离）

---

已知局限

• PyDESeq2 与 R DESeq2：离散度极低的基因（<1e-05）存在数值差异。若要精确复现 R 结果，请使用 rpy2。
• gseapy 与 R clusterProfiler：结果可能不同。

参考文件

• question_parsing.md - 从问题中提取参数
• data_loading.md - 数据加载与验证
• pydeseq2_workflow.md - 完整的 PyDESeq2 代码示例
• result_filtering.md - 高级筛选与提取
• dispersion_analysis.md - 离散度诊断
• enrichment_analysis.md - GO/KEGG/Reactome 工作流
• output_formatting.md - 正确格式化答案
• worked_examples.md - 全部 10 种问题模式
• troubleshooting.md - 常见问题与调试
• edger_limma_voom.md - edgeR / limma-voom DE 路线：命令序列、契约、列交叉映射

实用脚本

主要的确定性脚本（如上所述）：

• r_deseq2_wrapper.py - R DESeq2 多对比分析 + Venn 图 + 单基因 LFC
• r_edger_limma_wrapper.py - edgeR / limma-voom DE（若可用则运行；预检 Rscript + edgeR/limma）
• multi_strain_venn.py - 根据 DEG CSV 计算 Venn 图风格的重叠百分比
• gene_length_correlation.py - 长度与表达量的 Pearson r（全部变体）
• pca_variance.py - 所有常见 PCA 变体中每个 PC 的方差百分比
• one_way_anova_f.py - ANOVA F 统计量 + p 值（长格式/宽格式/LFC 框架）

辅助工具：

• format_deseq2_output.py - 输出格式化工具
• load_count_matrix.py - 数据加载工具
• convert_rds_to_csv.py - 将 .rds DESeq2 结果转换为 CSV

分析约定

DESeq2 库选择——与权威分析流程保持一致

如果数据文件夹中有使用 pydeseq2 的 run_*.py，或使用 R DESeq2 的 analysis.R，就必须使用那个完全相同的库。这两个库在细微数值上存在差异（相同阈值下，差异表达基因数量通常相差 2–10%），因此，基准真值取决于权威分析流程实际运行的是哪个库。

如果不存在权威脚本，则优先使用 R DESeq2（通过 run_deseq2_analysis 工具或 Rscript）——它是更广泛发表的参考实现：

tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"experiment_metadata.csv","design":"~ Replicate + Media + Strain","contrast":"Strain, 97, 1","refit_cooks":true}'

优先使用数据集的权威脚本

在自行运行 DESeq2 之前，先对数据集文件夹运行 ls。如果看到 run_*.py、analysis.R、find_*.R 或类似文件，它们就是此基准的真值计算脚本。

1. 使用 cat 查看脚本中的确切参数——每一个 kwarg。
2. 如果脚本已输出你需要的量，运行 cd DATASET_DIR && python3 run_*.py 并采用其结果。
3. 如果问题需要同一个拟合模型的其他指标，则只做少量补充（增加额外的打印语句），不要改动 DeseqDataSet(...) / DeseqStats(...) 构造函数的调用。

逐字复制所有 kwargs：refit_cooks=True、alpha=0.05、n_cpus、design_factors、ref_level。省略 refit_cooks 等参数可能会显著改变差异表达基因数量。凭记忆使用默认值会得到与基准真值脚本不同的基因列表。

R DESeq2 与 pydeseq2

这两个库在边缘情况下可能会得出不同结果（例如，相同 alpha 下的显著基因数量通常相差约 2%）。务必与权威脚本使用的库保持一致。 如果没有脚本，则优先使用 R DESeq2——其行为在已发表论文中更常被引用。

首选：使用 run_deseq2_analysis ToolUniverse 工具，它通过 Rscript 运行 R DESeq2：

# Basic DESeq2
tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"metadata.csv","design":"~ condition","ref_level":"condition, Control"}'

# With contrast + LFC shrinkage + refit_cooks
tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"metadata.csv","design":"~ Replicate + Media + Strain","contrast":"Strain, 97, 1","refit_cooks":true,"lfc_shrinkage":true}'

# enrichGO + simplify (after saving DEG list to file)
tu run run_deseq2_analysis '{"operation":"enrichgo","gene_list_file":"sig_genes.txt","background_file":"all_genes.txt","simplify_cutoff":0.7}'

这样可以避免 pydeseq2 与 R DESeq2 之间的差异。该工具会返回显著基因数量、离散度估计值，以及供进一步分析使用的结果 CSV 路径。

确认菌株身份映射

如果问题通过编号指定菌株，并描述其生物学特征（例如敲除基因型），应根据问题文本或实验元数据确认对应关系，而不是根据编号顺序推断。运行 DESeq2 前，读取元数据 CSV，确认每个菌株编号对应哪种基因型。

读取名为 res_1vsN.rds 的 RDS/CSV 结果文件时，“N”指的是元数据中的菌株编号。请核实该编号对应哪种突变体——不要只根据文件名推断。

“A/B 中独有、但 C 中没有差异表达”= 排他，而非包含

当问题询问“在 A、B 单突变体之一中独有差异表达、但在 C（双突变体）中没有”的基因时，意思是 仅在 A 中独有或仅在 B 中独有，并且在 C 中也没有差异表达：

(A − B − C) ∪ (B − A − C)

不是 (A ∪ B) − C（该集合会包含 A∩B 的交集）。排他性解释通常会得到比包含性解释更小的数量。

集合运算中的百分比——检查分母

如果你的 unique_DE / total_DE 得出 30–50% 的数值，预期的分母很可能不是这个。常见的替代分母包括：

• unique_DE / total_genes_tested（细菌通常为 5–10%，真核生物通常低于 2%）
• unique_DE / |union of all sig sets|

重新阅读问题，确认“占……的百分比”所指的总体是什么。

“在菌株 X 中也差异表达”= 简单重叠，而非排他

当问题询问“A 中差异表达的基因有多大比例在 B 中也差异表达”时，计算 |A ∩ B| / |A|。不要减去其他菌株——“在 B 中也差异表达”不代表“仅在 B 中差异表达、而在 C 中不差异表达”。

重要：“也”表示简单的集合交集。如果问题说“在菌株 A 中差异表达、在菌株 B 中也差异表达的基因”，则计算：

overlap = sigA.intersection(sigB)
pct = len(overlap) / len(sigA) * 100

不要添加“但在菌株 C 中没有”等额外排除条件。“在 B 中也差异表达”仅指 A ∩ B。

离散度估计：R DESeq2 与 pydeseq2 的结果存在差异

对于关于离散度的问题（例如“有多少个基因的离散度低于 1e-05”），由于离散度拟合算法的实现存在差异，R DESeq2 和 pydeseq2 得出的数量不同。离散度问题一律使用 R DESeq2——基准测试的真实结果使用 R 计算：

library(DESeq2)
dds <- DESeq(dds)
# Pre-shrinkage (genewise) dispersions:
gene_disp <- mcols(dds)$dispGeneEst
cat("Below 1e-05:", sum(gene_disp < 1e-05, na.rm=TRUE), "\n")

Log2 倍数变化：确认对比方向

询问突变株 Y 中基因 X 的 log2FC 时，请确认 DESeq2 的 results() 调用使用了正确的对比。对于 ~Media + Strain 设计，参考 Strain 为 "1"：

• results(dds, contrast=c("Strain", "97", "1")) → 菌株 97 相对于参考菌株的 log2FC
• 符号很重要：负 log2FC = 突变株中的表达量低于参考菌株
• 如果数值的符号或大小不正确，请检查是否误用了错误的菌株系数

密度 / 按染色体计算的约定

如果问题询问“平均染色体密度”“全基因组平均密度”或类似内容：

• “平均染色体密度” = mean(per_chromosome_density)，其中 per_chromosome_density = n_features / chromosome_length 是分别针对每条染色体计算的，然后对所有染色体取平均值（非加权平均值）。
• 这与 total_features / total_length（加权平均值 / “均匀分布下的预期密度”）不同。后者通常用作卡方检验的期望值，而不是问题所问的“平均染色体密度”。
• 示例：在染色体长度为 6-120 Mb 的鸟类基因组中，mean(per_chrom_density) 大约是 total/total 的 2 倍，因为较短的染色体每 bp 的事件数相对更多。