4. 高通量测序数据处理流程

拿到下机数据的第一步，就是把那堆 FASTQ 文件处理干净。说白了，测序仪给我们的 raw data 就像一块刚挖出来的矿石——里面有好东西，但也混着不少石头渣子。你得先筛一遍，再磨一磨，才能真正拿去做差异表达分析。这套流程我跑了不下两百次，每次都会遇到点新花样，今天我把实战中踩过的坑和总结的经验掰开揉碎了讲给你。

4.1 原始数据：FASTQ 格式长什么样

FASTQ 文件每四条记录对应一条测序片段。记住这个节奏：

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

第一行是标识符（@开头），第二行是碱基序列，第三行是一个 + 号（有时还带点注释信息），第四行是质量分数，每个字符对应一个碱基的测序质量。

质量分数用的是 ASCII 编码，不同版本 Phred 分数偏移不一样。现在主流 Illumina 1.8+ 用 Phred+33，也就是 ASCII 码值减 33 得到质量值 Q。Q 和错误概率的关系简单粗暴：

Q 值	错误概率	碱基准确率
10	1/10	90%
20	1/100	99%
30	1/1000	99.9%
40	1/10000	99.99%

我个人习惯拿到数据先看一眼 Q20 和 Q30 的比例——低于 85% 的 Q30 就要警觉了。你想想看，如果一条 read 大部分碱基 Q 值不到 20，那它的比对和定量结果基本就是噪声。

4.2 质量控制：FastQC 是必备的第一道门

FastQC 是我每次跑 RNA-seq 必用的工具。别嫌它老，这东西稳得很。跑起来就一行命令：

fastqc sample_R1.fastq sample_R2.fastq -o fastqc_output/ -t 4

它会在输出目录生成一个 HTML 报告和一个 zip 压缩包。报告有 12 个模块，我最看重的是这几个：

• Per base sequence quality —— 看每个位置的质量箱线图。如果超过 Q20 的线开始往下掉，就得考虑修剪了。
• Per base sequence content —— 看看前几个碱基的 ATCG 比例是不是均匀。RNA-seq 里偶尔会出现接头残留，导致前几个碱基比例异常。
• Adapter Content —— 直接告诉你有没有接头污染。我遇到过接头上游能跑到 30% 的情况，那基本是建库时没有完全去除接头。
• Overrepresented sequences —— 出现太多相同序列要小心，可能是扩增偏好导致的重复。

我第一次做 single-end 数据时，没注意 FastQC 报告里有个模块叫 “Per tile sequence quality”。后来比对率死活上不去，找了好久才发现是某个特定 tile 区域出现了质量坍塌。所以，多花五分钟看一眼报告，能省一整天 debug 时间。

4.3 数据修剪：用 Trimmomatic 做精细活

FastQC 发现问题后，就该 Trimmomatic 上场了。它能去除接头、切除低质量碱基、过滤短 read。我这里讲一个双端数据的经典调用：

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

逐个参数拆开讲：

• ILLUMINACLIP——接头去除。冒号分隔的三个数字：允许的最大错配数（2）、比对得分阈值（30）、最小匹配长度（10）。数字越大，去接头越严格，但也可能误切。我一般用默认的 2:30:10，遇到 adapter 含量高的样本会调高到 2:32:12。
• LEADING:3——切除 read 开头质量低于 3 的碱基。3 这个值很低，基本只切掉那些完全极差的碱基。你可以提高到 5 或 10 试试。
• TRAILING:3——和 LEADING 对称，从尾部切。
• SLIDINGWINDOW:4:15——滑动窗口修剪。窗口大小 4 个碱基，窗口内平均质量低于 15 则切除该窗口及之后所有碱基。这是最核心的参数，直接影响保留 reads 的长度和数量。
• MINLEN:36——修剪后长度小于 36 bp 的 read 直接丢弃。建议设为 read 原始长度的 50% 左右。

我在项目中遇到过这样一个案例：某批次数据 Q30 只有 72%，但经过 SLIDINGWINDOW:4:20 修剪后，Q30 提升到了 94%，同时保留了 85% 的 reads。这个代价很划算，后期差异表达基因的数量多了两倍。

做完修剪后，一定要再跑一次 FastQC，看看接头是否彻底清除、质量是否有提升。这是闭环验证的思想——别自己骗自己，数据到底行不行，报告说了算。

上面这张图就是我每次做 RNA-seq 必走的流程，没有捷径。你可能会问：能不能跳过质控直接跑 Trimmomatic？我试过，结果发现比对率从 92% 掉到 78%，假阳性基因增加了三倍。说白了，质控和修剪花的时间，最终会以更高的数据质量和更少的假发现率回报你。