← Blog

从 FASTQ 到差异表达：一个可复现的 RNA-seq 分析流程

RNA-seqSnakemake生物信息学DESeq2

在中大读博期间，RNA-seq 分析是我反复要做的工作。不同的课题、不同的物种、不同的实验设计，但分析流程总是那几步：下载数据、质控、 trimming、比对、定量、标准化、差异分析、功能富集。每次都要重新跑几十条命令，调参数，还得回忆上次到底怎么做的。

这种重复劳动正是我开发 RNA-seq-analysis-pipeline 的初衷。这是一套基于 Snakemake 的工作流，能从原始 FASTQ 文件一路跑到差异表达结果和功能富集分析。

流程包含哪些步骤

整个流程覆盖了 RNA-seq 分析的完整链路：

  1. 数据获取 — 自动从 SRA 下载原始测序数据。
  2. 质量控制与修剪 — 使用 fastp 去接头、过滤低质量 reads，并生成 QC 报告。
  3. 序列比对 — STAR 和 HISAT2 并行运行，提供两组独立的比对结果用于交叉验证。
  4. 表达定量 — 同时生成 featureCounts 和 STAR 的基因水平计数。
  5. 差异表达分析 — 使用 DESeq2 进行差异分析，配合 ashr 进行 log-fold-change 收缩，输出可直接用于发表的结果。
  6. 共表达网络分析 — 通过 WGCNA 识别基因模块及其与表型的关联。
  7. 功能富集分析 — 使用 genekitr 对差异基因集进行 GO 和 KEGG 富集分析。

每一步都会生成报告和可视化图表，方便你在任何阶段检查数据质量。

为什么选择 Snakemake

Snakemake 在灵活性和可复现性之间取得了很好的平衡。流程中的每条规则（rule）都有独立的 conda 环境，避免了不同步骤之间的依赖冲突。同时项目也提供了 Docker 容器支持，确保分析环境完全一致。

这个流程已经注册在 Snakemake Workflow Catalog 中，其他人可以直接发现和使用。

配置驱动的设计

我做了一个关键的设计决策：所有参数都通过配置文件控制。分析新数据集时，你不需要修改工作流本身。只需要准备两个文件：

# config.yaml
samples: "config/samples.tsv"
outdir: "results"
genome:
  source: "ensembl"
  organism: "homo_sapiens"
  release: 110
  fasta: "Homo_sapiens.GRCh38.dna.primary_assembly.fa"
  gtf: "Homo_sapiens.GRCh38.110.gtf"
params:
  fastp:
    qualified_quality_phred: 20
    length_required: 50
  star:
    outFilterMismatchNmax: 10
# samples.tsv
sample_id	condition	platform	layout
SRR1234567	treated	ILLUMINA	PAIRED
SRR1234568	control	ILLUMINA	PAIRED

然后一条命令运行：

snakemake --use-conda --cores 16

就这么简单。换一个课题，改一下配置文件就行了。流程同时支持双端和单端测序数据，参考基因组可以从 Ensembl 自由选择。

过程中的收获

搭建这套流程的过程，逼着我深入了很多以前只是浅尝辄止的领域。

RNA-seq 分析的理论基础 — 我必须理解每个步骤背后的原因：为什么 trimming 会影响比对率，为什么真核生物的 RNA 比对需要剪接感知的比对工具，不同的定量策略会引入什么偏差，标准化到底对计数数据做了什么。这些知识后来成了我做计算生物学的根基。

R 与 Bioconductor 生态 — 使用 DESeq2、WGCNA 和 genekitr 的过程让我深入了解了 R 社区处理基因组数据的方式。S4 对象系统、SummarizedExperiment 抽象、Bioconductor 包的设计哲学，这些都是在构建下游分析模块时逐渐理解的。

可复现研究的实践 — 每条规则独立 conda 环境、容器化、配置驱动，这些其实都是软件工程中的理念。把 DevOps 的思维方式引入科学分析是一个视角的转变。科研不仅要出结果，更要确保别人（或未来的自己）能复现同样的结果。

工作流管理 — Snakemake 基于规则的设计，通过输入输出声明自动构建有向无环图（DAG），是一种完全不同的数据处理思维方式。一开始确实需要时间适应，但现在让我回到手写 shell 脚本的方式，我已经做不到了。

写在最后

这套流程很大程度上是我博士经历的产物。它反映了我的日常工作内容和我想要解决的痛点。比如并行运行两组比对和定量工具的设计，就源于有一次我在排查不一致的结果，需要确认是不是比对工具引入的问题。

如果你也在频繁做 RNA-seq 分析，我建议花时间搭建（或者直接采用）一套标准化的分析流程。前期投入确实不小，但在一致性、可复现性和工作效率上的长期回报是值得的。

代码在 GitHub 上，欢迎反馈和贡献。