bulk-rnaseq

Overview

本技能把 bulk RNA-seq 的原始 FASTQ 一路处理到可发表的差异表达结论：先做质控与接头去除，再用 Salmon 做转录本定量并汇总成基因水平整数 counts，随后运行 DESeq2 风格负二项 Wald 检验得到差异基因列表，接着做 GO/KEGG 超几何富集，最后输出火山图、MA 图、PCA 图、热图与文字报告。

差异表达统计提供 双后端：--backend python 使用内置的纯 NumPy/SciPy 实现（size factor 归一化、Cox-Reid 离散度估计、趋势拟合、MAP 收缩、Wald 检验、BH 校正），--backend r 调用 Rscript 上的 Bioconductor DESeq2。两者输出同一套列名的结果表，可直接互换对比。

When to Use This Skill

• 用户要做 bulk 转录组差异表达分析（关键词：RNA-seq、bulk RNA-seq、转录组、差异表达、DEG、DESeq2、Salmon、counts matrix）
• 用户已有 Salmon/STAR 的定量结果，需要构建 counts 矩阵并做差异表达
• 用户需要 GO / KEGG 富集分析（关键词：GO 富集、KEGG 富集、pathway enrichment、功能富集）
• 用户需要 RNA-seq 的标准发表级图表（火山图、MA 图、PCA、热图）
• 用户想在没有 R / Bioconductor 的环境里得到可靠的 DESeq2 风格结果

不适用：单细胞数据（用 single-cell-rnaseq）、仅做原始 reads 可视化、蛋白组与代谢组分析。

上游衔接：如果还没有 FASTQ 文件，可先用 ena-download 技能按 SRA/ENA accession 下载原始测序数据（支持断点续传与 MD5 校验），再交给本技能做定量与差异分析。

Quick Start

# 完整流程（需要 fastqc / fastp / salmon 在 PATH 中）
export PYTHONPATH="$PWD/packages"
python skills/bulk-rnaseq/scripts/run_pipeline.py \
    --input  /data/fastq \
    --output /results/rnaseq \
    --transcriptome /ref/transcriptome.fa \
    --tx2gene /ref/tx2gene.tsv \
    --metadata /data/metadata.csv \
    --go-annotation /ref/go_annotation.csv \
    --group-a treated --group-b control

# 跳过外部工具，直接从已有 Salmon 结果开始
python skills/bulk-rnaseq/scripts/run_pipeline.py \
    --input /data/fastq --output /results/rnaseq \
    --skip-qc --skip-trimming --skip-quant \
    --quant-dir /results/rnaseq/03_salmon \
    --tx2gene /ref/tx2gene.tsv \
    --metadata /data/metadata.csv \
    --group-a treated --group-b control --backend python

# 单独构建 counts 矩阵
python skills/bulk-rnaseq/scripts/build_counts_matrix.py \
    --quant-dir /results/rnaseq/03_salmon \
    --tx2gene /ref/tx2gene.tsv \
    --output /results/counts

Workflow

1. FastQC 质控 — 对每个 .fastq.gz 生成 HTML/ZIP 报告，可选 MultiQC 汇总。关注 per-base quality、adapter content、GC 分布。
2. fastp 过滤 — 质量过滤（默认 Phred Q20）、PE 自动检测接头、最短长度 36 bp；同时输出 JSON/HTML 报告。
3. Salmon 定量 — salmon index 构建转录组索引，salmon quant 带 --gcBias --seqBias --validateMappings 输出转录本级 TPM / NumReads / Length。
4. counts 矩阵构建 — 用 tx2gene.tsv 把转录本汇总到基因；counts 取整（NumReads 四舍五入），TPM/FPKM 同步输出；生成元数据模板。
5. DESeq2 差异表达 — size factor 归一化 → 基因离散度估计 → 趋势拟合 → MAP 离散度 → Wald 检验 → BH FDR。按 padj 与 |log2FC| 双阈值标注 Up / Down / NS。
6. GO/KEGG 富集 — 单侧超几何检验 + BH FDR，按 BP / CC / MF 分别输出。
7. 可视化 — 火山图、MA 图、富集条形图、样本 PCA、Top DEG 热图（PDF + PNG，300 DPI）。
8. 报告 — 汇总实验设计、DEG 列表、富集结果、方法说明与产出清单。

Parameters

Outputs

Validation

• counts 矩阵：检查 counts.csv 各样本列和是否接近（差异过大说明测序深度差异大，需检查 size factor）；确认 gene ID 与注释文件 ID 空间一致。
• size factor：应集中在 1 附近（0.5–2）。出现极端值通常意味着某样本建库或测序异常。
• 火山图：显著点应呈对称的双翼分布；MA 图应大致平坦。大量点集中在 log2FC = 0 附近是正常现象。
• p 值分布：期望绝大多数基因 padj > 0.05；若几乎全部显著，多半是分组标签写反或存在批次混杂。
• 离散度：dispersions.csv 中落在下限（1e-8 附近）说明该基因数据量不足，其 p 值不可用。
• 富集：检查 Overlap_Genes 与 Term_Genes 的比例与背景基因数；min_gene_set 太小会得到大量无意义的极小 term。
• 交叉验证：同一数据用 --backend python 与 --backend r 各跑一次，比较 log2FC 排序与显著基因集合是否基本一致。

References

• references/methodology.md — 各步骤方法学原理、统计推导与参数选择依据
• references/troubleshooting.md — 常见报错、结果异常与排错手册