bulk-rnaseq 技能实战：从定量结果到 counts 矩阵的组装与 DE、富集分析交接-平芜编程栈

bulk-rnaseq 技能实战：从定量结果到 counts 矩阵的组装与 DE、富集分析交接

【免费下载链接】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

导读

本篇文章聚焦 scientific-agent-skills 仓库中 bulk-rnaseq 技能的核心衔接环节——把 Salmon / STAR / featureCounts 的定量输出，整理成 pydeseq2 技能所需的counts.csv（基因×样本整数矩阵）与metadata.csv，再把差异表达结果正确投喂给 pathway-enrichment 技能做 GSEA / ORA。读完你可以掌握scripts/build_counts_matrix.py三种输入模式的真实用法、整数计数与基因 ID 映射背后的统计陷阱，以及 "DE 结果 → 富集分析" 的完整交接配方。

参考主文档：counts-and-handoff.md，整体流程编排见 bulk-rnaseq SKILL.md。

为什么需要专门的 "counts 组装与交接" 环节

在 bulk-rnaseq 技能的管线（FastQC/trim → align/quant →counts→ DE → enrichment → figures）中，"定量 → counts 矩阵"是唯一没有上游或下游技能单独负责、由本技能自己兜底的关键缺口（见 SKILL.md）。整条链路的最终目标非常明确：产出pydeseq2 技能想要的两个文件，再对 DE 结果做阈值/排序，供pathway-enrichment 技能消费：

counts.csv——基因 × 样本的整数矩阵（原始计数或长度标定计数；绝不能是 TPM/FPKM）。
metadata.csv—— 每个样本一行（行索引 = 与 counts 列对应的样本 ID），列描述实验设计（condition、batch等）。

这两份文件分别对应 pydeseq2 技能在 SKILL.md 中的加载约定：pydeseq2 底层使用 pandas，先pd.read_csv("counts.csv", index_col=0).T把数据转置成"样本 × 基因"，再与 metadata 合并构造DeseqDataSet。换句话说，约定是传递出来的，不是可以随意更改的。

Orientation：矩阵方向与转置约定（先别在这里翻车）

PyDESeq2 最终需要样本 × 基因的矩阵。但本技能约定counts.csv写成基因 × 样本——这与 Salmon / STAR / featureCounts 及 nf-core 的输出方向一致，也方便你直接从上游拷贝矩阵列。真正把方向切换成 PyDESeq2 期望格式的动作，由 pydeseq2 技能的加载器用.T完成。

正确做法是：始终保持counts.csv为基因 × 样本，把转置交给 DE 步骤，不要自己转置两次。这也和 build_counts_matrix.py 源码里的注释一致："counts.csv stays genes x samples; the pydeseq2 loader transposes to samples x genes"（Salmon 分支内部用 AnnData 时是 obs=samples、X=samples×genes，因此用txi.to_df().T转回基因×样本再落盘）。

Salmon → 基因计数：pytximport

Salmon 的输出是转录本层面（transcript-level）的估计计数，需要先按基因汇总。做法是使用pytximport——tximport 的 Python 移植版，参数为counts_from_abundance="length_scaled_tpm"。

之所以选length_scaled_tpm，是因为这是基因层面差异表达的正确选择：它会校正跨样本的转录本长度差异与转录本使用差异，产出的计数可以直接喂给 DE 工具。核心调用方式：

from pytximport import tximport quant_files = ["quant/s1/quant.sf", "quant/s2/quant.sf", "quant/s3/quant.sf"] txi = tximport( quant_files, data_type="salmon", transcript_gene_map="tx2gene.tsv", # 两列: transcript_id, gene_id counts_from_abundance="length_scaled_tpm", output_type="xarray", ignore_transcript_version=True, # 去掉 Ensembl 版本号后缀 (.N) ) # txi 保存基因 x 样本的估计计数; 供给 PyDESeq2 前四舍五入为整数(见下)

仓库内封装的命令把这一逻辑做成了可复现的 CLI：

python scripts/build_counts_matrix.py --from salmon \ --quant-dir quant/ --tx2gene tx2gene.tsv --output-dir counts/

对照 build_counts_matrix.py 中build_from_salmon()的实现，脚本实际做四件事：

按quant_dir.glob("*/quant.sf")发现每个样本的quant.sf，以父目录名作为样本名（也兼容扁平布局*.sf，此时由_clean_sample_name()剥离.sf后缀命名）；
缺失--tx2gene时直接报错退出（--tx2gene is required for --from salmon），未安装 pytximport 时提示uv pip install pytximport；
以ignore_transcript_version=True调用 tximport（与参考文档示例保持一致的语义）；
counts.round().astype(int)四舍五入并转整数，counts.index.name = "gene_id"。

获取 tx2gene 映射表

tx2gene是一张两列表：transcript_id → gene_id。文档给出三种来源：

内置辅助函数：pytximport.utils.create_transcript_gene_map(species="human")，可直接取 human/mouse 等物种的映射。
从注释 GTF 自行生成（权威做法，与你的定量参考一致）：

awk -F'\t' '$3=="transcript"{ match($9,/transcript_id "([^"]+)"/,t); match($9,/gene_id "([^"]+)"/,g); print t[1]"\t"g[1] }' \ annotation.gtf | sort -u | sed '1i transcript_id\tgene_id' > tx2gene.tsv

复用 nf-core/rnaseq 实际使用的 tx2gene：Path A（nf-core 路线）在输出中写入了它真正使用的映射表，直接复用能保证与你定量所用参考版本严格一致（详见 upstream-nfcore.md）。

STAR → 基因计数：ReadsPerGene

STAR 在--quantMode GeneCounts下为每个样本产出*.ReadsPerGene.out.tab。每个文件共 4 列：gene_id、unstranded（非链特异）、forward-strand、reverse-strand。组装时必须：

跳过 STAR 文件开头的 4 行统计信息（N_unmapped、N_multimapping等汇总行）；
按链特异性选取对应列：列索引1/2/3 → unstranded/forward/reverse。

源码中把这一映射固化成了常量表：

STAR_STRAND_COL = {"unstranded": 1, "forward": 2, "reverse": 3}

build_from_star()用pd.read_csv(tab, sep="\t", header=None, skiprows=4)跳过头 4 行，按列索引取数并转为int64，最后用fillna(0)补齐跨样本缺失的基因。文档对应的封装命令：

python scripts/build_counts_matrix.py --from star \ --quant-dir star/ --strandedness reverse --output-dir counts/

注意--strandedness默认值就是reverse（这是 TruSeq 链特异性 mRNA 文库最常见的设置，但请务必与你的建库试剂盒确认）。STAR 产出的本身就是整数计数，无需取整。错误选择链方向会静默丢掉约一半读数——这是 bulk-rnaseq 技能列出的最典型错误之一，参考 SKILL.md 的 Common Pitfalls。

featureCounts → 基因计数

featureCounts一次写出所有样本的合并矩阵，头部有以#开头的命令行注释行，真正的列结构为：Geneid, Chr, Start, End, Strand, Length, <bam1>, <bam2>, …。组装时保留Geneid和每个 BAM 对应的计数列，并把 BAM 列重命名为样本 ID。

python scripts/build_counts_matrix.py --from featurecounts \ --counts-file counts/featurecounts.txt --output-dir counts/

源码实现build_from_featurecounts()用pd.read_csv(..., sep="\t", comment="#")吃掉注释行，用白名单去掉 6 个元数据列，剩下全部作为样本计数列并astype("int64")；样本名由_clean_sample_name()统一剥离.Aligned.sortedByCoord.out.bam/.bam等后缀。featureCounts 计数同样已是整数。

估计计数 / 整数的微妙之处

PyDESeq2 的负二项计数模型要求整数计数。这里要区分两件事：

STAR、featureCounts 直接给出整数；
Salmon、RSEM 给出的是估计计数（可能含小数）。

针对小数的情况，本文档定义的技能路线与 R 生态的 "标准路线" 不同：

本技能做法：使用length_scaled_tpm并四舍五入到最近整数。因为长度标定计数已经把所有文库大小与转录本长度信息折入数值本身，四舍五入后当作计数使用，是基因层面 DE 中成熟且站得住脚的近似。这条路线也正是 nf-core 面向下游用户暴露的方式（详见 upstream-nfcore.md：salmon.merged.gene_counts_length_scaled.tsv，文件本身就是非整数，官方说明"length-scaled counts → use for DESeq2"，拿到手 round 成整数即可）。
R 正统路线（tximport→DESeqDataSetFromTximport）：导入原始计数，同时传入每个基因的平均转录本长度偏移量（offset），让 DESeq2 在模型内部处理长度。PyDESeq2 目前不接受这种 offset，因此"长度标定 + 取整"是标准的 Python 等价方案。

两条路线有一条共同铁律：永远不要把 TPM / FPKM 喂给 DESeq2——它们是归一化后的值，会直接破坏基于计数的模型假设。bulk-rnaseq 技能在 SKILL.md 的 Common Pitfalls 中也明确把 "Feeding TPM/FPKM to DESeq2" 列为高危错误，而这个桥接脚本正是防线所在。

基因 ID 映射（进入富集分析之前必做）

DESeq2 输出通常以Ensembl 基因 ID为主键（如ENSG00000141510），经常还带版本后缀（如.17）。而 Enrichr / MSigDB / g:Profiler 的基因集库期望的是基因符号（人类基因全大写）。映射不匹配是 "富集结果为空" 的头号原因。

文档建议依次处理：

去掉版本后缀：ids.str.replace(r"\.\d+$", "", regex=True)。
Ensembl → symbol 映射：可用gget技能（gget info）、database-lookup、pybiomart或mygene。在 Path A 上，nf-core 输出的salmon.merged.gene_counts_length_scaled.tsv本身就带有gene_name列（往往还有gene_id列），把 symbol 列与 gene_id 一起保留即可——这也是 upstream-nfcore.md 中把gene_name去掉再set_index("gene_id")交给 PyDESeq2 的原因：DE 全程用 Ensembl ID 跑，symbol 只在最终基因列表上需要。
映射只服务于富集输入：DE 阶段可以全程保留 Ensembl ID，只在最终投喂富集分析之前映射符号。这样既避开 ID 漂移，也不牺牲 DE 阶段与参考注释的一致性。

DE → 富集交接配方

pydeseq2 技能产出deseq2_results.csv（关键列：log2FoldChange、pvalue、padj、stat）之后，两种富集方法对输入的要求完全不同：

方法	输入	做法	原因
GSEA（preranked）	完整排序的基因列表	用 Wald 统计量`stat`排序	符号表示方向、大小表示证据强度；比直接按`log2FoldChange`排序更稳定（低计数基因的 log2FC 噪声大）。先不要做阈值筛选
ORA	阈值化命中列表	`padj < 0.05`，可再加`\|log2FoldChange\| > 1`；建议上/下调基因分开跑	超几何/Fisher 检验只关心"列表中是否有基因"，阈值是它的前提

这条规则在 pathway-enrichment 技能侧得到了双重印证：其 SKILL.md 明示"Never threshold a list and then feed it to GSEA——那会丢弃 GSEA 赖以工作的排序信息"；其 run_enrichment.py 的_build_rank_from_deseq2()也正是按此约定实现的——优先取stat列排序，缺失时才回退到sign(log2FoldChange) * -log10(pvalue)，并对索引做符号清洗与去重。

pathway-enrichment 技能提供了直接读取 DESeq2 结果 CSV 的命令行工具，交接零代码：

# GSEA: 直接从 DE 表构建排序(自动用 stat 列) python ../pathway-enrichment/scripts/run_enrichment.py gsea \ --deseq2 deseq2_results.csv --organism human --outdir enrichment/ --seed 123 # ORA: 从符号命中列表出发 python ../pathway-enrichment/scripts/run_enrichment.py ora \ --genes sig_symbols.txt --organism human --outdir enrichment/

两处共用的可选参数值得留意（源码解析见 run_enrichment.py）：--libraries默认是MSigDB_Hallmark_2020 GO_Biological_Process_2023 KEGG_2021_Human Reactome_2022；--fdr默认 0.05；GSEA 还支持--min-size 15 --max-size 500 --permutations 1000 --threads 4，且 GSEA 必须显式给--seed保证 p 值可复现。脚本按库内 FDR 过滤后写出*_results.csv、*_significant.csv与 dotplot。

投喂前请确认deseq2_results.csv/sig_symbols.txt里的 ID 已经是基因符号（或先完成映射），之后再交给 scientific-visualization 技能出图。若走 nf-core Path A，直接用salmon.merged.gene_counts_length_scaled.tsv即可，桥接脚本不参与——两条路径在基因级 counts 矩阵处汇合后，下游交接完全一致。

一个完整的端到端示例

把三个上游来源、桥接脚本与下游交接串起来的典型流程（Path B 独立工具路线）：

# 1. 每样本定量（Salmon 示例，输出 quant/s1, quant/s2, ...） salmon quant -i salmon_index -l A -1 s1_R1.fq.gz -2 s1_R2.fq.gz \ --gcBias --seqBias -p 8 -o quant/s1 # 2. 组装基因 x 样本整数 counts 矩阵 + metadata 模板 python scripts/build_counts_matrix.py --from salmon \ --quant-dir quant/ --tx2gene tx2gene.tsv --output-dir counts/ # 3. 编辑 counts/metadata_template.csv，把 condition/batch 等列填成真实分组 # (模板中 condition 预置为 "CHANGE_ME", batch 为空串) # 4. pydeseq2 技能: counts.csv + metadata.csv -> DE 表 python ../pydeseq2/scripts/run_deseq2_analysis.py \ --counts counts/counts.csv --metadata counts/metadata_template.csv \ --design "~condition" --contrast condition treated control \ --output results/ # 5. pathway-enrichment 技能: stat 排序做 GSEA, 阈值列表做 ORA python ../pathway-enrichment/scripts/run_enrichment.py gsea \ --deseq2 results/deseq2_results.csv --organism human --outdir enrichment/ --seed 123 python ../pathway-enrichment/scripts/run_enrichment.py ora \ --genes sig_symbols.txt --organism human --outdir enrichment/

关于实验设计、重复数与批次模型的约束（≥3 生物学重复、~batch + condition设计公式），以及对齐后 PCA / 样本距离热图等质量门控，是这份交接配方能产出可信结果的前提，详见 design-and-qc.md；两种上游路线的完整命令与参数则分别见 upstream-nfcore.md（Path A）与 upstream-manual.md（Path B）。每一环节都建议固定工具与参考基因组/注释版本，以便在方法学部分完整复现。

交接配方中常见的五个坑（自检清单）

方向转置两次：counts.csv保持基因 × 样本，只让 pydeseq2 加载器.T一次。
把小数估计计数直接喂给 PyDESeq2：Salmon/RSEM 必须先length_scaled_tpm再 round 成整数（STAR/featureCounts 本来就是整数）。
用 TPM/FPKM 当 counts：归一化值会破坏 DESeq2 计数模型，宁可不做也不要这么投喂。
用 Ensembl ID 或带版本号 ID 直接跑富集：Enrichr/MSigDB 要人类全大写符号；先str.replace(r"\.\d+$", "")再去映射，否则大概率 "nothing is enriched"。
ORA/GSEA 输入用反：GSEA 喂完整stat排序列表、不设阈值；ORA 喂padj < 0.05（可叠加|log2FoldChange| > 1）的命中列表。

需要与本文联动阅读的仓库文件：桥接脚本 build_counts_matrix.py（三种来源的完整 CLI 参数与实现）、样本表校验工具 validate_samplesheet.py、下游的 pydeseq2 技能 SKILL.md（转置、过滤、contrast 约定）与 pathway-enrichment 技能的 SKILL.md（ORA vs GSEA 选型与库选择），以及本文的直接上一级编排文档 bulk-rnaseq SKILL.md。

创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考