RNA-seq数据分析概述：从测序数据到表达矩阵

做RNA-seq分析这么多年，我最大的体会就是：数据预处理的质量，直接决定了后续所有分析的成败。说白了，你后面用多牛的模型、画多漂亮的ROC曲线，如果表达矩阵本身有硬伤，那都是白搭。

这一章，咱们就把RNA-seq从原始测序数据到表达矩阵的全流程捋一遍。我会重点说清楚每个步骤到底是干什么的，以及我曾经踩过的坑。

1.1 RNA-seq分析的基本流程

RNA-seq的原始数据，说白了就是一堆几亿条的小短序列（reads）。核心任务就是搞清楚两件事：「这条read来自哪个基因？」 和 「那个基因到底有多少条read？」

我个人习惯把整个流程拆成四大块：

  1. 原始数据预处理——把脏数据洗洗干净（质量控制、去接头、过滤低质量数据）
  2. 比对或定量——把reads定位到基因组或转录组上
  3. 基因/转录本定量——统计每个基因或转录本上有多少reads
  4. 表达矩阵构建——把定量结果整理成行是基因、列是样本的表格

核心要点：不管是差异表达分析、聚类分析，还是咱们这门课要讲的ROC评估，输入数据都必须是高质量的、经过正确归一化的表达矩阵。这一步不要图快。

你想想看，如果原始数据里有残留的接头序列，或者比对参数设得不对，那下游模型评估时得到的「准确率」就是虚假的。我以前就吃过这个亏，后面会细说。

1.2 第一步：原始数据预处理

拿到下机数据后，千万别急着跑比对。先做质量评价。

1.2.1 质量评价（QC）

我通常用 FastQC 来做初步的质量评价。看一眼报告里的几个关键指标，基本就能判断这批数据能不能用。

fastqc sample_1.fastq.gz sample_2.fastq.gz -o ./qc_reports

重点关注这些地方：

  • Per base sequence quality：每个碱基位置的质量得分。如果尾巴上大面积低于Q20，建议直接截短或过滤。
  • Adapter Content：接头残留。如果有明显的接头峰，必须用工具去接头。
  • Overrepresented sequences：看看有没有可疑的序列污染。

注意： 我曾经遇到过一个项目，FastQC报告里显示样品间质量差异很大。当时赶进度没仔细看，直接全部合并处理，结果后续差异分析结果完全无法解释。后来排查才发现是其中一个样品的测序芯片有瑕疵。这件事之后，我养成了每个样品单独看FastQC报告的习惯。

1.2.2 去除接头与低质量序列

这一步我推荐用 Trimmomatic 或 cutadapt。这里给一个Trimmomatic的常见用法：

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_paired.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_paired.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 TRAILING:3 \
  SLIDINGWINDOW:4:15 MINLEN:36

参数简单解释一下：

  • ILLUMINACLIP：去接头的参数对。
  • SLIDINGWINDOW:4:15：4个碱基为一组，平均质量低于15就截断。嗯，这个值比较宽松，通常建议改成20。
  • MINLEN:36：过滤掉太短的reads。

小经验： 去接头这一步，不同测序平台用的接头序列不一样，一定要确认用对了适配器文件。我就犯过用错适配器文件导致大量数据被误删的错误。切记把平台确认好。

1.3 第二步：比对（Alignment）

这一步是把干净的reads映射到参考基因组或转录组上。比对工具的选型会影响后续很多细节。

1.3.1 常用比对工具

工具特性适用场景
STAR速度快，敏感度高，能处理剪切位点有参考基因组的标准RNA-seq
HISAT2基于图基因组，内存占用比STAR低大型基因组或复杂转录组
Salmon / kallisto比对+定量一体，不需要先比对再计数转录本级别的定量，轻量快速

我个人习惯用STAR，因为它稳定，而且输出格式比较规范。不过你要做好心理准备，STAR第一次建立基因组索引会比较慢，而且占内存。我曾经在服务器上建一个人类基因组索引，跑了将近2个小时。

1.3.2 STAR比对实战

先建立索引：

STAR --runMode genomeGenerate \
  --genomeDir ./star_index \
  --genomeFastaFiles GRCh38.fa \
  --sjdbGTFfile GRCh38.gtf \
  --sjdbOverhang 149

sjdbOverhang 设为读长减1，这个参数我建议别省。

然后比对：

STAR --genomeDir ./star_index \
  --readFilesIn sample_R1_paired.fastq.gz sample_R2_paired.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix sample_

输出会得到一个 sorted BAM 文件。这里有个细节——为了方便下游做基因计数，我通常会要求STAR输出直接按坐标排序的BAM，省去后面再排序的麻烦。

避坑指南： 对比对结果的QC。不要只看比对率大于80%就觉得万事大吉。我曾经遇到过一个样品，比对率高达95%，但是大部分reads都比对到线粒体基因组上，显然是样品污染了。所以记得要检查比对分布——看看reads是落在外显子、内含子还是基因间区。

1.4 第三步：基因定量（Gene Quantification）

拿到BAM文件后，需要统计每个基因上有多少条reads。这一步我强烈推荐用 featureCounts（来自Subread包），速度快且好用。

featureCounts -T 8 \
  -t exon \
  -g gene_id \
  -a GRCh38.gtf \
  -o counts.txt \
  sample1.bam sample2.bam sample3.bam

参数解析：

  • -T：线程数。
  • -t：指定GTF中的特征类型，通常用exon。
  • -g：分组依据，一般用gene_id。
  • -o：输出的计数矩阵文件。

你会得到一个类似这样的输出：

Geneid   Chr  Start  End  Strand  Length  sample1.bam  sample2.bam
ENSG000001  chr1  1000  2000  +  1000  45  32
ENSG000002  chr1  3000  4000  -  800  12  18

最后几列就是每个基因在每个样品中的原始read count。

注意： featureCounts默认是单端模式。如果你的数据是双端测序，记得加上 -p 参数，否则计数会翻倍，而且比对关系也会错。我刚开始做RNA-seq时就忘加这个参数，结果发现差异基因数量异常多，排查了好久才发现是计数翻倍了。

1.5 第四步：表达矩阵构建与归一化

featureCounts输出的counts.txt，其实就是一个原始计数矩阵。但直接用这个矩阵做下游分析是不行的。

为什么？因为不同样品间的测序深度不同、基因长度不同、GC含量不同，这些技术因素都会引入噪声。

常用的归一化方法：

方法适用场景说明
CPM（Counts Per Million）快速预览、文库大小标准化简单除以总read数，不能校正基因长度
RPKM/FPKM基因长度+文库大小双标准化适合基因内比较，不适合基因间比较
TPM基因间比较更公平先对基因长度归一化，再对文库大小归一化
DESeq2 中位数比值差异表达分析假设大部分基因不差异表达，自动校正

我个人习惯在构建表达矩阵时，同时保存原始count和TPM。原始count用于差异表达分析的统计检验，TPM用于后续机器学习的特征输入。

1.6 本章小结

好了，从测序数据到表达矩阵的完整流程就这些。我再用一张图帮你串起来：

从测序数据到表达矩阵流程图 从测序数据到表达矩阵：核心流程 原始测序数据 FASTQ文件 质量控制 (QC) FastQC + 去接头 比对/定量 STAR / HISAT2 Trimmomatic clean FASTQ 基因定量 featureCounts BAM文件 排序 + 索引 输出 表达矩阵 counts / TPM / RPKM 定量结果 归一化处理 DESeq2 / edgeR

嗯，上面这个流程图把关键环节和常用工具都标出来了。你跟着这个框架走，基本不会跑偏。从原始FASTQ到最终的表达矩阵，每一步踩实了，后面做模型评估才会心里有底。

这门课后续的交叉验证、ROC曲线绘制，全部是建立在「高质量的表达矩阵」之上的。所以这一章的内容，说白了就是课程的基石。你别嫌啰嗦，打好地基比什么都重要。

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