2、RNA-seq 分析概述：流程步骤、数据格式、质控、比对、定量、差异表达

RNA-seq 分析说白了就是“看基因在说什么”。我们拿到的是细胞里所有 RNA 的碎片，通过测序仪读成几千万条短序列，然后一步步把它还原成基因表达量，再找出不同条件下哪些基因被“调高”或“调低”。

我刚开始做这个分析时，总觉得每一步都挺简单，结果一跑起来各种报错。后来发现，搞懂每步的核心逻辑比背命令更重要。下面这张图就是整个 RNA-seq 标准流程的主干——你跟着它走，基本不会跑偏。

RNA-seq 分析标准流程 原始数据 质量控制 比对 定量 差异表达分析 说明： 原始数据 (FASTQ)→ 质量控制 (FastQC, Trimmomatic)→ 比对 (STAR, Hisat2)→ 定量 (featureCounts, Salmon)→ 差异表达 (DESeq2, edgeR) 每一步都会产出中间文件，但真正用于生物学解释的最后一步的结果。

原始数据格式：.fastq 文件里到底有什么？

RNA-seq 的测序结果通常是 FASTQ 格式。说白了，每条记录包含四行：序列标识符、碱基序列、分隔符（一个加号）、质量字符串。质量字符串每个字符对应一个碱基的测序质量，用 ASCII 码编码。

举个例子，一条 read 长这样：

@HWUSI-EAS100R:6:73:941:1973#0/1
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

那个第二行的字母就是碱基，第四行的符号代表质量得分。符号越“靠后”，测序质量越高。为什么要关心这个？因为低质量的碱基会影响后续比对和定量。我做过一个项目，一批数据里大约 15% 的碱基质量低于 Q20，结果差异表达基因的数量明显偏多——后来发现是假阳性，吃了大亏。

个人习惯：拿到新数据后第一件事不是直接跑比对，而是先用 less 瞄几眼 FASTQ 文件的前 1000 行。看看质量字符分布是否均匀，序列长度是否一致。这比任何统计报告都直观。

质量控制：把脏数据洗干净

这一步我通常会拆成两个子任务：质量评估和数据清洗。

质量评估——FastQC

FastQC 是黄金标准工具。它会生成一个 HTML 报告，从碱基质量、GC 含量、N 含量、接头污染等多个维度打分。

你直接在终端里跑：

fastqc sample_1.fastq sample_2.fastq -o qc_reports/

然后打开 HTML 文件，看每个模块的“红绿灯”。如果某个模块出现红色叉号，就说明有问题。

注意：FastQC 的红色标记不一定要全部清除。比如“Per base sequence content”在很多 RNA-seq 数据中可能是红色，因为文库制备的随机引物会引入偏差。我曾经被这个红色吓到，花了三天做二次清洗，结果发现它不影响下游差异分析。要学会区分“需要处理的警示”和“可以忽略的警示”。

数据清洗——Trimmomatic / Cutadapt

如果原始数据里有接头残留或低质量尾部，就得切掉。我常用 Trimmomatic，因为它支持双端数据还能自动检测接头。

一条典型的命令：

java -jar trimmomatic-0.39.jar PE \
  sample_1.fastq sample_2.fastq \
  sample_1_trimmed.fastq sample_1_unpaired.fastq \
  sample_2_trimmed.fastq sample_2_unpaired.fastq \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

参数含义其实很简单：去掉接头（ILLUMINACLIP），去掉头尾低质量碱基（LEADING/TRAILING），用 4 碱基窗口扫描，平均质量低于 15 就切掉（SLIDINGWINDOW），最后长度小于 36 的扔掉。

我个人建议把清洗后的数据再跑一次 FastQC 做对比，确保质量有改善。

比对：把 reads 放到参考基因组上

质控之后，我们会得到干净的 reads。下一步就是把它们“贴”回参考基因组上。RNA-seq 的比对和 DNA 比对不一样，它必须处理剪接位点——因为内含子被剪掉了，一条 read 可能跨越两个外显子。

所以不能用 BWA 或 Bowtie，得用剪接感知比对器。目前主流是 STAR 和 Hisat2。

STAR 的使用要点

STAR 的速度很快，但内存消耗大。它第一步需要建索引：

STAR --runMode genomeGenerate \
  --genomeDir star_index/ \
  --genomeFastaFiles hg38.fa \
  --sjdbGTFfile hg38.gtf \
  --sjdbOverhang 100

然后比对：

STAR --genomeDir star_index/ \
  --readFilesIn sample_1_trimmed.fastq sample_2_trimmed.fastq \
  --outFileNamePrefix sample_ \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode TranscriptomeSAM GeneCounts

注意那个 --quantMode GeneCounts，它可以直接输出每个基因的原始计数，省去后面再跑 featureCounts 的步骤。不过我个人更习惯单独定量，方便灵活调整参数。

踩坑经验：我刚开始用 STAR 时遇到过 --genomeDir 路径写错，结果跑了一晚上才发现索引不对。后来我在 Nextflow 里把索引路径写进配置文件，再也不用担心手误。建议你也把所有输入路径参数化，不要硬编码。

定量：从比对结果到基因表达计数

比对后我们会得到 SAM/BAM 文件。定量就是统计每个基因被多少个 reads 覆盖。有两种常见思路：

  • 基于比对的结果：用 featureCounts 或 HTSeq-count 直接从 BAM 文件中计数。
  • 基于转录本的定量：用 Salmon/Sailfish/kallisto，不需要完整比对，直接对 reads 做“伪比对”，速度快得多。
工具 方法 优点 缺点
featureCounts 基于 BAM 比对计数 准确，内存低 需要先比对，产生大文件
HTSeq-count 基于 BAM 比对计数 可自定义计数模式 速度慢，单线程
Salmon 转录组伪比对 极快，无需参考基因组索引 对剪接复杂样本准确率略低

用 featureCounts 的典型命令：

featureCounts -T 8 -t exon -g gene_id \
  -a hg38.gtf \
  -o counts.txt \
  sample_Aligned.sortedByCoord.out.bam

输出是一个矩阵：每一列是一个样品，每一行是一个基因，单元格里的数字就是原始计数。

差异表达分析：找出哪些基因变了

得到计数矩阵后，我们终于可以回答生物学问题了：两组样本（比如处理组 vs 对照组）之间，哪些基因的表达量有统计学差异？

差异表达分析的核心假设是：计数分布符合负二项分布（不是正态分布，因为 RNA-seq 计数有过度离散）。常用的工具有 DESeq2（R 包）和 edgeR（R 包）。它们在统计模型上略有差异，但基本流程相似：

  1. 输入计数矩阵 + 实验设计（哪几个样本属于哪个组）。
  2. 内部进行归一化（比如 DESeq2 的中位数比值归一化），消除文库大小和组成带来的偏差。
  3. 拟合负二项模型，计算每个基因的 fold change 和 p 值。
  4. 多重假设检验校正（通常用 BH 方法），得到调整后的 p 值。

输出结果是一张表，包含 log2FoldChange、pvalue、padj 等列。常用的筛选标准是 |log2FoldChange| > 1 且 padj < 0.05。

别只看 p 值：我记得有一次发现一个基因 padj 小到 10⁻¹²，但 log2FoldChange 只有 0.3。这种基因虽然统计显著但变化幅度很小，生物学意义可能不大。我一般会结合 fold change 和 padj 两个维度做火山图来选候选基因。

小汇

RNA-seq 分析从原始数据到差异表达，每步都不复杂，但环环相扣。质控做不好，后续全是垃圾；比对参数错了，定量结果直接偏掉；差异表达如果不理解统计假设，很容易筛选出假阳性。Nextflow 的好处就是能把这几步串成一条自动化的流水线，而且每一步都能做容器化封装，避免环境冲突。

← 上一章 📖 返回目录 下一章 →