📚 课程目录第3章 RNA-seq数据预处理流程概览

RNA-seq数据预处理流程概览

说实话，很多同学一上来就急着跑差异表达分析，结果折腾半天发现结果根本不可信。为什么？因为原始数据里藏着太多坑。我个人的习惯是，不管数据来自哪家测序公司，拿到手必须先过一遍预处理流程。这步做扎实了，后续的分析才有意义。

RNA-seq数据预处理，说白了就是把测序仪吐出来的原始信号，转换成我们可以用的基因表达定量值。整个过程大致可以分为六个核心步骤。下面这张图可以帮你快速建立整体印象。

原始数据长什么样？

测序公司交付给你的，通常是一个或多个 FASTQ 文件。这个文件格式没什么神秘的，就是四行一组的文本记录。

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
GTGCTTGCCTGAAAATGGAACGAGTGCGATCCACGGGAGGTTCTGG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

第一行是序列标识符，包含测序仪器、flowcell 位置等信息。第二行是碱基序列。第三行是个加号，有时候也带信息。第四行是质量分数，每个字符对应一个碱基的质量，用的是 Phred 打分系统。

我在项目中遇到过，有些同学忽略了质量分数，直接拿序列去比对。结果下游差异分析时发现一堆假阳性。说白了，质量分数是你判断数据好坏的「照妖镜」。

质量评估：先看数据懂不懂事

拿到数据第一件事是什么？我个人习惯先跑 FastQC。这工具虽然老，但确实管用。它会输出一个 HTML 报告，告诉你数据的基本情况。

小提示：批量处理时建议用 MultiQC 汇总多个样本的 FastQC 报告，一目了然。我一般这么写命令：

# 对每个样本跑 FastQC
fastqc sample_1.fastq sample_2.fastq -o qc_reports/

# 汇总所有报告
multiqc qc_reports/ -o multiqc_report/

快去看看报告里的这些指标：

每个碱基位置的测序质量：看箱线图是否大部分在绿色区域（Q30 以上）。如果 3' 端质量掉得厉害，后面就得考虑截断。
碱基含量分布：理论上 A/T/G/C 四条线应该平行。如果有明显偏移，可能是有接头污染或建库偏向。
GC 含量：正常分布是单峰。出现双峰可能提示有污染或数据来源有问题。
接头污染：报告里会直接标记是否找到过多的 adapter 序列。
重复序列比例：有些重复是正常的，但过高可能提示 PCR 扩增偏差。

注意：如果 Q30 比例低于 80%，我个人会建议重测序。我曾经为了省这点钱硬着头皮往下做，结果验证实验全部对不上，白白浪费了两周时间。所以该砍就砍，别侥幸。

数据清洗：把垃圾倒掉

质量评估完了，就该动手清洗了。你想想看，如果原始数据里有 5% 的垃圾，你的表达量就会偏 5%。对于低表达基因来说，这点误差可能直接让结果翻转。

清洗主要包括三件事：

去除接头：用 Cutadapt 或 Trimmomatic 把 adapter 序列切掉。
质量修剪：把 3' 端质量差的碱基切掉。
长度过滤：扔掉那些太短的序列，比如低于 36 bp 的读段。

我常用的 Trimmomatic 命令长这样：

java -jar trimmomatic-0.39.jar SE \
  sample.fastq \
  sample_trimmed.fastq \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

参数解释一下：

ILLUMINACLIP：接头文件，允许 2 个错配，比对分值阈值 30，palindrome 模式分值 10。
LEADING:3：切掉 5' 端质量低于 3 的碱基。
TRAILING:3：切掉 3' 端质量低于 3 的碱基。
SLIDINGWINDOW:4:15：4 个碱基一个窗口，平均质量低于 15 就切掉。
MINLEN:36：保留长度不低于 36 的读段。

清洗完之后，记得再用 FastQC 跑一次，看看前后对比。嗯，这一步叫「清洗前后对比」，是很多审稿人重点关注的内容。

序列比对：把 reads 放回基因组上

数据洗干净了，接下来就得知道这些 reads 到底来自基因组的哪个位置。这步叫比对。比对工具有很多，但有参考基因组和没参考基因组的策略完全不一样。

我先说有参考基因组的情况。我个人最常用的是 STAR，它跑得快，准确率也不错。

STAR 的基本用法：

# 第一步：建立基因组索引
STAR --runThreadN 10 \
     --runMode genomeGenerate \
     --genomeDir star_index \
     --genomeFastaFiles genome.fa \
     --sjdbGTFfile annotation.gtf

# 第二步：比对
STAR --genomeDir star_index \
     --readFilesIn sample_trimmed.fastq \
     --runThreadN 10 \
     --outFileNamePrefix sample_ \
     --outSAMtype BAM SortedByCoordinate

比对完你会得到一个 BAM 文件。这文件是二进制格式，不能直接看，但你可以用 samtools 查看它的统计信息。

samtools flagstat sample_Aligned.sortedByCoord.out.bam

输出里重点看两个数字：

Mapping rate：就是比对率。正常样本应该在 70%-90% 之间。低于 60% 就值得警惕了。
Unique mapping rate：比对到唯一位置的 reads 比例。用于定量的一般只取这部分。

我曾经踩过一个坑：有一个样本的比对率只有 40%，我当时以为是比对参数问题，调了半天还是不行。后来一查，才发现是样本标签搞混了，那根本就是个细菌样本。所以说，比对率异常时先别怀疑工具，先确认样本来源。

对于没有参考基因组的物种，就得用 Trinity 或 StringTie 做 de novo 组装。这步计算量会大很多，后面课程再细讲。

从片段到计数：基因表达定量

比对完了，怎么知道每个基因有多少表达？这就用到了基因定量。最经典的工具是 featureCounts，它能把比对到基因区域内的 reads 数出来。

featureCounts -T 10 \
  -p \
  -t exon \
  -g gene_id \
  -a annotation.gtf \
  -o counts.txt \
  sample_Aligned.sortedByCoord.out.bam

这里解释一下参数：

-p：如果测的是 paired-end，加上这个参数。
-t exon：在 exon 水平上进行计数。
-g gene_id：按基因 ID 汇总。
-a annotation.gtf：提供基因注释文件。

输出的 counts.txt 长这样：

Geneid  sample_Aligned.sortedByCoord.out.bam
ENSG00000223972  12
ENSG00000227232  45
ENSG00000278267  0
ENSG00000243485  89

每一行是一个基因，第二列是对应样本中比对到这个基因的 reads 数。这就是原始的 count 数据，也是后续归一化的输入。

小建议：如果样本很多，可以把所有样本的 BAM 文件一起输入 featureCounts，它会自动生成一个矩阵形式的 counts 文件，便于后续合并分析。

表达矩阵的最终形态

定量完成后，我们就得到了最原始的 count 矩阵。它长这样：

GeneID	Sample1	Sample2	Sample3	Sample4
GeneA	120	94	137	105
GeneB	23	19	28	21
GeneC	0	2	5	1
GeneD	400	435	390	412

为什么要强调「原始」？因为这些 count 数据还没做任何归一化。样本间的测序深度不同、基因长度不同，直接比较会出大问题。比如 Sample3 的 GeneA 是 137 条，Sample1 是 120 条，是不是 Sample3 表达更高？不一定，也许只因为 Sample3 测得更深。

这也就是为什么我们在进行差异分析之前，必须做归一化处理。后面的章节会详细讲 TPM、FPKM、DESeq2 的归一化策略，以及蛋白质组学里那些「套路相似但细节不同」的方法。

回顾一下，整个预处理流程的核心目标就是：从 FASTQ 到 count 矩阵。每一步都要留下中间文件和质量报告，方便回溯。我之前吃过一次亏，跑完了所有样本才发现有一个样本的接头没去掉，所有定量结果都得重来。从那以后，我每做完一步都会检查并保存中间结果。

好，预处理流程就讲到这。下面该进入归一化环节了。