2、RNA-seq 分析流程回顾：从原始数据到差异表达的关键步骤

聊审计追踪之前，咱们先把 RNA-seq 分析流程再过一遍。说实话，我见过不少团队把审计追踪做得特别复杂，结果发现流程本身就有漏洞——那就本末倒置了。你今天要是能把这套标准流程刻在脑子里，后面做审计点设计会顺手得多。

RNA-seq 分析，说白了就是一趟「从原始读段到生物学结论」的流水线。每个环节都有它的目的，也都有它容易出问题的地方。嗯，咱们今天就按实际跑项目的顺序，把这几个关键步骤捋一遍。

2.1 数据预处理阶段

所有分析的第一步，都从 FASTQ 文件开始。这步要是出了问题，后面再怎么调参数都没用。我个人习惯是：拿到数据的第一件事，先做质量检查。

2.1.1 原始数据质量评估

推荐用 FastQC 看几个核心指标：

  • 碱基质量分数（Q-score）：一般看 Q30 比例，低于 80% 的样本我建议直接标记警告
  • GC 含量分布：正常应该在 40-60%，偏离太离谱的可能存在污染
  • 重复序列比例：文库扩增过度会导致重复率升高

执行质量评估的命令很简单：

# 单样本 FastQC 检查
fastqc sample_R1.fastq.gz -o qc_report/

# 批量检查并生成汇总
multiqc qc_report/ -o multiqc_report/

我的经验：一定要保留 FastQC 生成的 HTML 报告。我在项目中遇到过两次「实验室说数据没问题，结果我一跑 FastQC 发现 adapter 残留高达 15%」的情况。审计时要能拿出来说「我检查过质量了，这是当时的报告」。

2.1.2 接头去除与低质量过滤

质量不达标的数据必须处理掉。我最常用的工具是 cutadapt 和 Trimmomatic。具体步骤见下表：

处理步骤常用参数目的
去除接头-a AGATCGGAAGAGC移除测序接头序列
切除低质量碱基SLIDINGWINDOW:4:154bp 窗口平均质量低于 Q15 则切除
长度过滤MINLEN:36去掉太短的序列，避免比对出错

一个标准的 Trimmomatic 命令长这样：

java -jar trimmomatic.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_trimmed.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_trimmed.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

注意：剪切的参数不能随意改。我曾经因为某个项目想「多保留一些数据」而放松了质量阈值，结果下游差异表达分析出来的显著基因有 60% 都是假阳性。审计追踪里必须记录每次运行的参数版本。

2.2 比对与定量阶段

比对就是把清洗后的 reads 定位到参考基因组上。这步决定了基因定量的准确性，是整个分析的核心节点。

2.2.1 比对工具选择

现在主流的比对工具有 STAR、HISAT2、Salmon 这几个。我个人偏好 STAR，因为它的比对速度快，而且对剪接点识别做得比较好。当然，你也可以根据项目特点来选：

  • STAR：适用性最广，需要做融合基因检测的推荐用它
  • HISAT2：内存占用低，适合服务器配置一般的场景
  • Salmon：直接做定量，不需要提前比对，适合只看表达量的项目

一个典型的 STAR 比对命令：

STAR --genomeDir /path/to/genome_index \
     --readFilesIn sample_R1_trimmed.fastq.gz \
                  sample_R2_trimmed.fastq.gz \
     --readFilesCommand zcat \
     --outFileNamePrefix sample_ \
     --outSAMtype BAM SortedByCoordinate \
     --outSAMattributes NH HI AS nM

2.2.2 基因表达定量

比对之后拿到 BAM 文件，接下来就是定量。这里有两个策略：

  1. 基于基因注释的定量：使用 featureCounts 或 htseq-count，直接统计每个基因覆盖的 reads 数
  2. 转录本级别的定量：使用 Salmon 或 RSEM，先预估转录本丰度再汇总到基因

featureCounts 的用法很直接：

featureCounts -T 8 -t exon -g gene_id \
  -a annotation.gtf \
  -o counts.txt \
  *.sorted.bam

关键审计点：定量结果要保留完整 log。featureCounts 会在输出文件最后一行记录「Assigned」和「Unassigned」的 reads 数量。我一般要求项目日志里必须记录这两项——如果未分配的 reads 超过 30%，比对参数一定有问题。

2.3 差异表达分析阶段

终于到了大家最关心的环节。前面的步骤都是在准备数据，现在要回答「哪些基因在两个条件之间有显著变化？」。

2.3.1 表达矩阵与分组设置

把 featureCounts 输出的 counts 矩阵整理到一个表格里，格式大概是：

GeneID  Control_1  Control_2  Treatment_1  Treatment_2
ENSG001  120        145        320          298
ENSG002  55         48         12           9
ENSG003  2000       1980       2100         2050

然后定义分组信息。我习惯用 DESeq2，它有一个很清晰的 colData 格式：

# R 代码示例
colData <- data.frame(
  row.names = c("Control_1","Control_2","Treatment_1","Treatment_2"),
  condition = c("Control","Control","Treatment","Treatment")
)

2.3.2 差异基因筛选标准

差异表达分析不是随便跑个脚本就完事的。筛选标准直接决定了结果的可信度：

常用标准阈值说明
padj（校正后 p 值）< 0.05多重假设检验校正，控制 FDR
log2FoldChange|log2FC| > 1相当于表达量变化 2 倍以上
基础表达量baseMean > 10排除低表达基因，减少噪音

完整的 DESeq2 分析流程：

library(DESeq2)
dds <- DESeqDataSetFromMatrix(
  countData = counts_matrix,
  colData = colData,
  design = ~ condition
)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition","Treatment","Control"))
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
write.csv(as.data.frame(res_sig), "differential_genes.csv")

避坑指南：千万记住，DESeq2 的 design 公式必须包含批次信息（如果有）。我曾经接手过一个项目，分析的人忘记把批次配成协变量，结果差异基因列表完全被批次效应主导了——那重新跑了三天三夜才把结果纠正过来。

2.4 流程总结：一张图讲清楚

下面这张图把 RNA-seq 分析的核心步骤串起来了。你看一眼就能明白，每个环节之间有明确的数据流向，而每个步骤都留下了可供审计的「痕迹」——比如 FastQC 报告、比对日志、参数配置文件。这些就是咱们后续做审计追踪的原材料。

RNA-seq 分析核心流程 ① 原始 FASTQ ② 质量控制 Clean FASTQ ③ 比对 ④ 基因定量 ⑤ 差异表达分析 ⑥ 差异基因列表 每个步骤留下的审计线索 • FastQC HTML 报告 • Trimmomatic 参数日志 • STAR 比对统计 审计追踪核心记录 参数版本 · 环境信息 · 时间戳

从 FASTQ 到差异基因列表，每一步都环环相扣。你可能注意到了，我特意在流程图两侧加了「审计线索」的标注——FastQC 报告、比对日志、参数文件这些，每一个环节都需要记录下来。这也是为什么我总跟团队说：好的分析习惯，就是好的审计基础。

2.5 关键参数汇总表格

最后，我把整个流程中出现的关键参数整理成一个对照表。这个表格在做审计追踪时特别有用——你看一眼就知道哪些参数版本必须记录：

分析步骤工具/参数是否必须审计记录内容示例
质量评估FastQC 版本建议记录v0.11.9
数据清洗Trimmomatic 参数必须SLIDINGWINDOW:4:15 MINLEN:36
比对STAR 版本 + --outSAMattributes必须v2.7.10a, NH HI AS nM
定量featureCounts -t -g必须-t exon -g gene_id
差异表达DESeq2 版本 + design formula必须v1.34.0, ~ condition
筛选标准padj, log2FC 阈值建议记录padj<0.05, |log2FC|>1

小技巧：我写了一个自动化脚本，每次运行完一步分析，就把当前命令行、工具版本和输出文件的 MD5 值一起写进项目日志里。这样审计的时候，我能 10 分钟之内把所有记录翻出来——而不是花两天去回忆「当时用的哪个版本」。

好了，RNA-seq 分析流程就回顾到这里。你应该已经清楚了从原始数据到差异表达需要经过哪些关卡。下一块内容咱们就要聊正题了——审计追踪到底是什么，它为什么在 RNA-seq 项目里那么重要。

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