📚 课程目录第3章数据预处理概览

3、数据预处理概览：从原始FASTQ到清洁读段

数据预处理这步，我见过太多人栽跟头了。说白了，你后面做拼接、分箱、基因预测，结果好不好，八成取决于你前期的数据清洗。原始测序数据里啥都有——接头、低质量碱基、宿主DNA污染，甚至测序仪本身的错误。不把这些清掉，你后面花再多GPU算力也是白搭。

3.1 原始测序数据格式：FASTQ

FASTQ是目前最通用的原始测序数据格式。每一条read由四行组成，缺一不可。

我记得第一次接触宏基因组数据时，拿到一个几百G的FASTQ文件，差点没把我电脑撑爆。后来才明白，格式得先搞透。

第一行：以@开头，后面跟着序列标识符和可选描述
第二行：碱基序列（A、T、C、G、N）
第三行：以+开头，通常只保留@行内容或为空
第四行：质量得分，每个碱基对应一个ASCII字符

举个例子，一条典型的FASTQ记录长这样：

@SRR12345678.1 1 length=150
ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

质量得分那行，用的是ASCII码。早期有个坑——不同版本的质量编码方式不一样。早些年Sanger格式用ASCII 33-126，而Illumina 1.3+用ASCII 64-126。现在绝大多数数据都是Phred+33标准，但万一你拿到老数据，质量解读就全错了。

注意：质量得分Q = -10 * log10(P)，P是碱基错误概率。Q30意味着错误率千分之一，也就是每1000个碱基可能有1个错误。宏基因组分析中，我建议至少保留Q20以上的碱基。

3.2 质量控制：FastQC与Trimmomatic

3.2.1 FastQC — 一眼看穿数据质量

FastQC是个经典的质量评估工具。我个人的习惯是，不管数据来源是哪里，先跑一遍FastQC看看整体面貌。

安装很简单：

wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip
unzip fastqc_v0.12.1.zip
cd FastQC
chmod +x fastqc
./fastqc sample_R1.fastq.gz -o fastqc_output/

运行结束后，每个样本会生成一个HTML报告。常见的问题我总结在下面：

模块	红灯含义	我的处理建议
Per base sequence quality	某个位置碱基质量明显下降	裁剪低质量末尾
Per sequence GC content	GC分布偏离正态	检查是否混入非目标物种
Adapter Content	接头残留比例高	必须去除接头
Overrepresented sequences	某些序列出现频率异常	可能是引物二聚体或污染

小技巧：如果你处理的是双端测序数据，建议把R1和R2分别跑FastQC。我习惯先看R2的质量，因为大部分情况下R2质量更差。

3.2.2 Trimmomatic — 数据清洗的主力

Trimmomatic是目前最主流的清洗工具之一。它支持多线程，处理大规模数据时效率很高。

我常用的Trimmomatic命令长这样：

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

参数解释一下：

PE：双端模式。单端用SE。
ILLUMINACLIP：去掉接头序列。后面的2:30:10分别表示允许的最少匹配碱基数、匹配得分阈值、以及单端匹配得分阈值。
LEADING:3：去掉开头质量低于3的碱基。
TRAILING:3：去掉末尾质量低于3的碱基。
SLIDINGWINDOW:4:15：滑动窗口裁剪。窗口大小4个碱基，平均质量低于15就裁剪。这个参数我一般现场调。
MINLEN:36：保留的最短读长。

我曾经踩过的坑：有一次把SLIDINGWINDOW的窗口设成了10:20，结果大部分数据都被裁成了碎片。后来才发现，宏基因组数据的读长本身就参差不齐，窗口太大反而误伤。我现在一般用4:15，对大多数项目都适用。

3.3 去宿主污染

宏基因组样本里最烦人的就是宿主DNA污染。你明明想研究肠道菌群，结果测出来的数据大半是人基因组。这些宿主reads会显著降低拼接效率，还会让分箱结果一团糟。

去宿主污染的常见策略有两种：

基于比对的方法：把reads比对到宿主参考基因组上，保留比对不上的reads。
基于k-mer的方法：构建宿主基因组的k-mer索引，快速过滤匹配的reads。

我个人倾向于前者。Bowtie2或BWA都是好选择，尤其Bowtie2在宏基因组数据上表现不错。

一个去宿主的典型流程：

# 先构建宿主基因组索引
bowtie2-build human_genome.fa human_index

# 比对宿主基因组
bowtie2 -x human_index \
  -1 sample_trimmed_R1.fastq.gz \
  -2 sample_trimmed_R2.fastq.gz \
  -S sample.sam \
  --very-sensitive \
  -p 32

# 提取未比对上的reads
samtools view -bS sample.sam > sample.bam
samtools view -b -f 12 sample.bam > sample_unmapped.bam
samtools sort -n sample_unmapped.bam -o sample_unmapped_sorted.bam
bedtools bamtofastq -i sample_unmapped_sorted.bam \
  -fq sample_clean_R1.fastq \
  -fq2 sample_clean_R2.fastq

-f 12这个参数意思很明确：保留R1和R2都没比对上的pair。我在实际项目中遇到过，有时候只检查R1，结果R2配对但R1没比对，数据对不齐就炸了。

个人经验：如果你的数据量特别大，比如超过100G的FASTQ，可以考虑用Kraken2这类基于k-mer的工具来去宿主。它在内存利用率和速度上比Bowtie2有优势，尤其在CPU资源有限时。

3.4 整个预处理流程的样子

下面这张图把今天讲的内容串起来了。从原始FASTQ开始，经过质量评估、质量裁剪、去宿主污染，最终得到清洁读段。

  核心要点回顾：
  FASTQ是原始数据格式，注意质量编码方式（Phred+33 vs Phred+64）。
FastQC做全面质量评估，Trimmomatic做精细裁剪。
去宿主污染是宏基因组数据分析的第一步硬门槛，跳过这一步后面全是噪声。
预处理后的清洁读段质量直接影响后续GPU加速分析的效率。

嗯，说到这里你应该明白了。数据预处理不是可有可无的步骤，而是整个分析流程的地基。地基不稳，楼再高也是危楼。我这些年做过的项目，但凡预处理阶段做得糙的，后期调参调到头秃。反过来，预处理做到位了，后面的拼接、分箱几乎是一马平川。