4 高通量测序数据下载与质量控制

4.1 数据下载——从公共库里捞数据

嗯，拿到测序数据的套路是什么？说白了，就是从 NCBI 的 SRA 数据库或者 ENA 里，把别人或者自己上传的原始数据拽下来。这个东西，坑不少，我先说说我自己的习惯。

我最常用的是 SRA Toolkit 里的 fastq-dump 或者 prefetch。prefetch 先下载 .sra 文件，然后 fastq-dump 把它转成 fastq。但这里有个坑——fastq-dump 默认会有内存问题，尤其是大文件。我建议你用 --split-files 参数，特别是双端测序的时候，它会把 R1 和 R2 分开。

# 下载单个 SRA 文件
prefetch SRR12345678

# 转为 fastq，双端拆分
fastq-dump --split-files SRR12345678

# 或者直接一条命令
fastq-dump --split-files --gzip SRR12345678

你想想看，为什么要加 --gzip？因为测序数据量大，很多项目都有几十个 G。压缩后传输快，而且 FastQC 也能直接读压缩文件。我在项目中遇到过一个案例——没加压缩，结果磁盘直接撑爆了，复盘时才发现白白浪费了 30% 的 IO 时间。

4.2 为什么要做质量控制？

测序仪也是精密仪器，总有状态不好的时候。我见过一批数据，读长看起来没问题，跑完 FastQC 才发现接头根本没被切干净。结果下游差异分析一塌糊涂，白白跑了三天流程。所以——数据下载完，第一件事不是做 mapping，而是做质控。

质量控制的核心点有三个：

• 碱基质量分数 — Q30 以下的比例太高，你就等着看到处都是 unreliable 的比对吧。
• GC 含量分布 — 如果 GC 出现双峰，很可能有外源污染（比如细菌混进了你的细胞系样本）。
• Duplication levels — 时间序列 RNA-seq 经常出现技术重复偏差，高 duplication 会误导差异分析的统计模型。
• Adapter content — 接头的残留也是老问题。我建议直接看“Overrepresented sequences”那栏。

我曾经在坑里躺过——有些 SRA 记录用的还是 Illumina 的老版本接头，你用默认的 Trimmomatic 参数可能切不干净。最好先跑 FastQC 看看过表达序列，再决定接头类型。

4.3 FastQC——你的测序体检报告

FastQC 就是给每条 fastq 文件做个“体检报告”。我个人习惯一开始就跑所有样本，然后拿 MultiQC 汇总，一目了然。

# 单样本
fastqc SRR12345678_1.fastq.gz SRR12345678_2.fastq.gz

# 批量跑
for f in *.fastq.gz; do fastqc "$f" & done
wait

它会生成一个 .html 和 .zip 的报告。我建议你看下面这些关键指标：

指标	正常绿勾	黄/红警告信号	我的经验
Per base sequence quality	Q30 覆盖整条 read 中位线 > 30	尾巴突然跌到 Q20 以下	直接裁掉后面低质量部分，不用心疼
GC content	接近理论分布（单峰）	双峰，或偏离理论值超过 15%	双峰大概率是污染，我曾经追查过是建库时样本交叉了
Adapter content	前 20bp 无过量	某类 adapter 富集超过 5%	直接用 Trimmomatic 或 Cutadapt 暴力裁剪
Duplication levels	约 10-20% 重复	> 60% 重复	如果文库复杂度本来就低（如一些单细胞建库），需要留意影响定量

4.4 自适应裁剪与过滤——我说的是“够用就好”

Trimmomatic 和 Cutadapt 是两把常用的剪刀。我的原则是：够用就好，不要过度裁剪。你想想看，RNA-seq 的 read 本来就短，多裁 5bp 可能就丢了一个外显子交界的信息。

下面是我在项目里常用的 Trimmomatic 双端参数：

# Trimmomatic 双端
trimmomatic PE \
  SRR12345678_1.fastq.gz SRR12345678_2.fastq.gz \
  SRR12345678_1_clean.fastq.gz SRR12345678_1_unpaired.fastq.gz \
  SRR12345678_2_clean.fastq.gz SRR12345678_2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

参数解释一下——MINLEN:36 是底线，低于 36bp 的直接丢掉；SLIDINGWINDOW:4:15 是滑窗：每 4bp 算一个窗口，平均质量低于 Q15 就裁掉。我建议对于高覆盖的时间序列数据，可以适当调低到 Q10，保住更多 read。

4.5 MultiQC——把一堆报告变成一张大表

说实话，单个 FastQC 报告看 20 个样本还行，如果手上有 60 个样本（比如一个 3 个时间点 × 5 个处理 × 4 replicate 的设计），翻报告会翻到崩溃。

MultiQC 就是用来汇总的。一行命令搞定：

# 在当前目录下找所有 fastqc_data 文件夹
multiqc .

它会生成一个 multiqc_report.html，里面有热图、箱线图、表格，一眼就能看出哪些样本有问题。我建议你看看“Sequence Quality Histograms”这个模块——如果某个时间点的样本整体质量偏低，很可能说明那批试剂的保存条件出了问题。

4.6 整合一下：本章的核心流程

下面这张 SVG 图总结了我个人在数据下载与质控环节的流向逻辑：

4.7 最后的检查清单

做完以上步骤，我习惯列一个清单核对：

1. 所有样本的 FastQC 报告有没有红灯？如果没有，说明你已经规避了最常见的坑。
2. MultiQC 汇总表里每个时间点的数据量是否一致？我曾遇到过某一时间点所有样本数据量骤减 30%，后来发现是那批样本在测序时某个通道堵了。
3. 裁剪后的 read 长度分布是否稳定？如果某个样本的 read 长度中位数跟其他样本差太多，考虑要不要剔除它。
4. 有没有保留未配对的 read？时间序列分析时，我建议至少保留一份未配对的作为备份，有些拼接软件用得着。