FASTQ格式及质量控制

RNA-seq数据分析的第一步，就是拿到测序仪吐出来的原始数据。

嗯，最常见的格式就是FASTQ。说白了，它就是一个文本文件，每四条记录代表一条测序片段。我当年刚接触时，总觉得这格式很简单，不就是序列加上质量嘛。但实际处理中，质量值解析、编码方式、过滤策略……任何一个地方出问题，后续分析都会崩掉。

这一章我们就来聊透它：从FASTQ格式本身，到质量评估，再到经典的去污工具Trimmomatic。我个人习惯，拿到任何RNA-seq数据，第一件事就是跑FastQC——先看数据质量，再决定怎么清。

1. FASTQ格式解析

每条read在FASTQ里长这样：

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

第一行以@开头，是序列ID，通常包含仪器、flowcell、lane等信息。第二行是碱基序列。第三行是一个+，可选项，后面可能跟着同样的ID。第四行是质量字符串，每个字符对应一个碱基的测序质量。

核心：质量值编码

质量值用ASCII码表示。最常见的两种：

  • Phred+33（Sanger格式）：ASCII码 = Phred质量分 + 33。字符范围 ! 到 ~，对应分值0~93。
  • Phred+64（旧Illumina 1.3+）：ASCII码 = Phred质量分 + 64。字符范围 @ 到 h，对应分值0~62。

现代数据几乎全是Phred+33。但如果你拿到的是老数据（比如2008-2010年的），要小心。我见过有人把Phred+64当+33解析，结果质量值全偏了64分，后续过滤完全失准。

怎么快速判断？看质量字符串的第一个字符。如果小于@（ASCII 64），那就肯定是Phred+33。如果全是@以上，有可能是Phred+64，但更可能是高质量reads的Phred+33（因为分值≥31时字符就是?以上）。稳妥的做法：让FastQC自动检测，它有内置的编码猜测逻辑。

2. 质量值解读

Phred质量分Q = -10 × log₁₀(P)，其中P是碱基错误的概率。我常用这个对照表快速判断：

Q值 错误概率 碱基准确率 字符（Phred+33）
10 1/10 90% +
20 1/100 99% 5
30 1/1000 99.9% ?
40 1/10000 99.99% I

Q30是公认的“高质量”门槛。RNA-seq中，我一般要求平均质量≥30，低于20的碱基就要警惕。

实战技巧：写一个小脚本，打印每条read的平均质量，快速扫一眼整体水平。我曾经在项目中遇到过一批数据，平均质量只有Q25——以为是批次问题，结果发现是测序仪上样量偏高。提前发现能省下很多后续排查时间。

3. FastQC质量评估

FastQC是我工具箱里最常用的质检工具。装起来很简单：

# 下载并运行
wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip
unzip fastqc_v0.12.1.zip
./FastQC/fastqc sample_R1.fastq.gz -o ./fastqc_output/

它会生成一个HTML报告，里面包含11个模块。我每次必看前三个：

  • Per base sequence quality：每个位置的碱基质量分布。理想情况是箱线图都落在绿色区（Q28以上）。如果末尾明显下降，说明reads长度过长或测序质量衰减。
  • Per sequence quality scores：每条read的平均质量分布。如果有个明显的双峰，说明混入了低质量reads。
  • Adapter Content：接头污染。RNA-seq测序片段一般短于read长度时，就会读到接头。这个指标如果超过5%，就要用Trimmomatic去接头。

常见陷阱：FastQC报告里出现“Failed”并不一定代表数据不可用。比如Per base sequence quality在末尾变黄甚至红，但你的后续分析有质量过滤步骤，完全没问题。重点是看有没有系统性偏差——比如前10个碱基质量异常低，可能暗示文库构建问题。

另外，FastQC会报告Overrepresented sequences。这个指标在RNA-seq中经常有“假阳性”——比如高表达的核糖体RNA（rRNA）片段会被标记为过代表。你想想看，如果你的数据是rRNA depletion文库，这个模块可以忽略。

4. Trimmomatic过滤

看完FastQC报告后，就该动手清数据了。Trimmomatic是我首选的过滤工具，支持多线程、能处理paired-end和single-end数据。下面是我在实际项目里用的典型命令：

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_clean_R1.fastq.gz sample_unpaired_R1.fastq.gz \
  sample_clean_R2.fastq.gz sample_unpaired_R2.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

参数拆开讲：

  • ILLUMINACLIP：去除接头。需要提供接头序列文件（如TruSeq3-PE.fa），三个数字分别代表：允许的最大错配数（2）、比对分数阈值（30）、最小匹配长度（10）。
  • LEADING:3：从read开头切掉质量低于3的碱基。
  • TRAILING:3：从read末尾切掉质量低于3的碱基。
  • SLIDINGWINDOW:4:15：滑动窗口4个碱基，计算平均质量，如果低于15就切掉。这是最关键的过滤策略。我一般设窗口大小4，阈值15~20，具体取决于数据质量。
  • MINLEN:36：过滤后长度小于36的read直接丢弃。

重要细节：对于paired-end数据，Trimmomatic会同时处理两条reads，保证配对。如果一个read被丢弃，对应的另一端也会被丢掉（或者输出到_unpaired文件）。我建议保留_unpaired文件，有时它们还能用上，比如在组装中。

另外，Trimmomatic还有一个TOPHRED33参数，用来强制转换质量编码。如果你确定数据是Phred+64，加上TOPHRED33会让后续处理统一。

我自己的经验：第一次跑Trimmomatic时，可以先设SLIDINGWINDOW阈值15跑一轮，然后再次对过滤后的数据跑FastQC。如果发现Per base sequence quality在read末尾依然很高（比如Q30以上），可以尝试把阈值调到20，进一步剔除质量稍差的碱基。但别调太高，不然数据量损失太大。平衡样本数和质量，是RNA-seq分析中最常遇到的权衡。

5. 质量控制流程全景

这一套下来，核心流程可以用下面的图记住：

RNA-seq 质量控制流程 原始 FASTQ FastQC 质量评估报告 Trimmomatic 去接头 + 质量裁剪 过滤后 FASTQ 再次 FastQC 验证清理效果 提示：清理后务必再次运行FastQC 对比前后报告，确保问题已解决 质量合格？

图上显示得很清楚：原始数据先送FastQC，生成报告后根据问题用Trimmomatic处理，然后一定要再跑一次FastQC验证。我如果看到清理后报告里的Per base sequence quality绿油油一片，Adapter Content降到0%，心里就踏实了。

你想想看，如果跳过第二次FastQC，你永远不知道Trimmomatic有没有把问题解决掉。我曾经遇到过数据质量极差，第一次跑了Trimmomatic后没检查，直接往下走比对——结果比对率只有40%。后来补了二次FastQC才发现，原来是接头序列文件选错了，根本没去掉接头。白白浪费了两天时间。

好了，这一章的核心就是：读得懂FASTQ，看得懂FastQC报告，用Trimmomatic精准过滤。把这套流程练熟，后面建模的数据质量就有保障了。

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