正文

测序数据比对流程详解从FASTQ到BAM格式转换的每一步操作指南

/0 浏览量
0909

嘿，我是 Agnes。今天咱们不聊那些枯燥的教科书定义，而是像剥洋葱一样，把高通量测序（NGS）数据处理中最核心、也是最“折腾人”的一环——从原始数据到可比对的 BAM 文件——彻底掰开揉碎讲清楚。

很多人看到服务器上的 FASTQ 文件，心里既兴奋又慌：兴奋的是数据质控过了，慌的是接下来该怎么让这些短序列“各归其位”。别担心，这篇指南就是为你准备的“保姆级”地图。我会结合实战场景，把每一步的原理、参数玄学和常见坑点都给你摆出来。

1. 起点：理解你的“原材料”——FASTQ 格式

在动手跑代码之前，你得先知道手里拿的是什么。FASTQ 是测序仪吐出的原始数据，它的结构看着简单，实则暗藏玄机。

FASTQ 的四个组成部分

每一段序列在 FASTQ 里占用 4 行：

  1. 序列头行：以 @ 开头，包含读段（Read）的唯一标识符和测序仪器信息。
  2. 序列行：实际的 DNA 碱基序列，由 A、T、C、G 组成（有时会有 N 表示未知）。
  3. 分隔行：以 + 开头，可以是空，也可以重复序列头，主要用于区分序列行和质量行。
  4. 质量行：长度与序列行完全一致，用 ASCII 码字符表示每个碱基的测序质量值。
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
AGCTTGGCGGTAGGCTC... (序列)
+
IIIHGIJIJJIH... (质量值，ASCII编码)

质量值的秘密：Phred 分数

这部分极其重要，直接决定你后面能不能过滤掉垃圾数据。质量值 \(Q\) 与错误概率 \(P\) 的关系是：

\[ Q = -10 \log_{10}(P) \]

常见的编码方式有两种：

  • Sanger/Illumina 1.8+：Phred+33。这是目前最主流的格式。例如，字符 I 对应的 ASCII 码是 73，减去 33 得到质量值 40，意味着错误概率为 \(10^{-4}\)，也就是 0.01% 的错误率。
  • Illumina 1.3-1.7 / Solexa：Phred+64。老版本 Illumina 测序仪使用，现在很少见了，但如果你的 FASTQ 质量值字符很密集（全是高 ASCII 字符），要小心是不是编码搞错了。

💡 专家提示：如果你拿到一批数据，发现质量值全是 ~ 或者特别低，别急着骂测序公司，先检查是不是 FASTQ 编码搞错了，或者接头污染严重。

2. 第一步：质控与清洗（QC & Trimming）

原始数据通常带着“杂质”：测序接头、低质量碱基、PCR 重复引物序列等。如果不处理干净，它们会干扰后续比对，导致大量读段被错误地丢弃或比对到错误位置。

常用工具：FastQC 和 Trimmomatic

FastQC 负责“体检”，看看数据整体情况：

  • Per base sequence quality：检查每个位置的碱基质量，通常在序列末尾质量会下降。
  • Per sequence quality scores：查看整个读段的平均质量分布。
  • Adapter Content：检测接头污染比例。

Trimmomatic（或 cutadapt、fastp）负责“整容”，去掉脏东西。

实战代码：使用 Trimmomatic 进行双端数据清洗

假设你有一对双端测序数据 sample_R1.fastq.gz 和 sample_R2.fastq.gz。

# 定义输入输出文件
INPUT_R1="sample_R1.fastq.gz"
INPUT_R2="sample_R2.fastq.gz"
OUTPUT_R1="sample_R1_clean.fastq.gz"
OUTPUT_R2="sample_R2_clean.fastq.gz"
UNPAIRED_R1="sample_R1_unpaired.fastq.gz"
UNPAIRED_R2="sample_R2_unpaired.fastq.gz"

# 定义适配器文件路径 (adapters.fa 需提前准备)
ADAPTERS="adapters.fa"

# 运行 Trimmomatic
java -jar /path/to/trimmomatic.jar \
    PE -phred33 \
    $INPUT_R1 $INPUT_R2 \
    $OUTPUT_R1 $UNPAIRED_R1 \
    $OUTPUT_R2 $UNPAIRED_R2 \
    ILLUMINACLIP:$ADAPTERS:2:30:10 \  # 去除接头，参数：种子匹配长度:最大不匹配数:最大不匹配比例
    LEADING:3 \                      # 去除头部质量低于 3 的碱基
    TRAILING:3 \                     # 去除尾部质量低于 3 的碱基
    SLIDINGWINDOW:4:20 \             # 滑动窗口，4bp 窗口平均质量低于 20 时截断
    MINLEN:36                        # 最短保留长度，低于 36bp 的读段丢弃

参数解读：

  • -phred33：明确指定编码方式，避免误判。
  • ILLUMINACLIP：这是最关键的一步。如果接头没切干净，比对软件可能会把接头序列强行比对到基因组上，产生大量假阳性。
  • SLIDINGWINDOW：这是最常用的质量控制策略。想象你在看一条生产线，每 4 个产品平均一下，一旦发现质量下滑，立刻裁掉剩下的。

💡 专家提示：对于 RNA-Seq 数据，如果你做的是转录组分析，建议保留 UNPAIRED 的文件，因为有些读段可能因为一端质量太差而丢失了配对关系，这些单端数据依然有价值。

3. 第二步：参考基因组准备

比对需要“地图”。这个地图就是参考基因组（Reference Genome）。你不能直接拿个 .fa 文件就去比对，必须先用工具把它索引化，建立数据结构（如 FM-index），否则比对速度会慢得像蜗牛。

以人类基因组 GRCh38 为例

3.1 下载基因组

通常从 UCSC、Ensembl 或 NCBI 下载。为了演示，我们假设你已经有了 Homo_sapiens.GRCh38.dna.primary_assembly.fa。

3.2 构建 Bowtie2 索引

Bowtie2 是常用的比对工具之一，它需要特定的索引格式。

# 使用 bowtie2-build 建立索引
bowtie2-build \
    Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    ./indexes/grch38

# 注意：
# 1. 输入是 fasta 文件
# 2. 输出前缀是 ./indexes/grch38
# 3. 构建完成后，你会看到目录下出现 .bt2 文件：
#    grch38.1.bt2, grch38.2.bt2 ... grch38.8.bt2

3.3 构建 BWA 索引（另一种选择）

如果你打算用 BWA-MEM（在长读段或变异检测中更常用），命令略有不同：

# 使用 bwa index 建立索引
bwa index \
    -p ./indexes/grch38_bwa \  # 输出前缀
    Homo_sapiens.GRCh38.dna.primary_assembly.fa

# 同时，如果你要做 SAMtools 的 faidx 索引（用于查看特定区域）
samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa

💡 专家提示：

  • 索引文件会占用很大磁盘空间（人类基因组索引约 10-20GB），但这是值得的，因为比对时内存访问速度极快。
  • 如果你比对的是线粒体 DNA 或非标准染色体（如嵌合体、病毒序列），确保这些序列也包含在参考基因组 FASTA 文件中，否则那些读段会被丢弃或比对错误。

4. 第三步：序列比对（Alignment）

这是整个流程的核心。我们要把清洗后的短读段“贴”回参考基因组上。这里我们重点介绍两种最主流的比对器：Bowtie2 和 BWA-MEM。

方案 A：使用 Bowtie2（适合小基因组、RNA-Seq 前处理、变异检测前的快速比对）

Bowtie2 速度快，内存占用低，但对长读段（>100bp）的支持不如 BWA。

# 定义变量
CLEAN_R1="sample_R1_clean.fastq.gz"
CLEAN_R2="sample_R2_clean.fastq.gz"
INDEX_PREFIX="./indexes/grch38"
OUTPUT_BAM="sample_aligned.sam"

# 运行 Bowtie2
bowtie2 -x $INDEX_PREFIX \
        -1 $CLEAN_R1 \
        -2 $CLEAN_R2 \
        -p 16 \                  # 使用 16 个线程加速
        --no-unal \              # 输出文件中不包含未比对的读段（可选，节省空间）
        --norc \                 # 不输出反向互补链的比对结果（可选）
        -S $OUTPUT_BAM \
        --rg-id sample1 \        # 读取组 ID
        --rg SM:sample1 \        # 样本名
        --rg PL:ILLUMINA \       # 测序平台
        --rg PU:unit1 \          # 平台单位/流动槽 ID

# 将 SAM 转换为 BAM 并排序
samtools view -@ 16 -bS $OUTPUT_BAM | samtools sort -@ 16 -o sample_aligned_sorted.bam
samtools index sample_aligned_sorted.bam

参数解读：

  • --rg-id, --rg SM 等：这些是“Read Group”信息。非常重要！ 许多下游分析工具（如 GATK 变异检测）强制要求 BAM 文件包含正确的 RG 信息，否则可能会报错或结果不准。
  • --no-unal：不输出未比对的读段。如果你需要计算比对率，建议去掉这个参数，或者单独统计未比对数。

方案 B：使用 BWA-MEM（适合 DNA 测序、大基因组、变异检测首选）

BWA-MEM 是目前体细胞突变检测和拷贝数变异分析的标准工具。它比 Bowtie2 更慢，但对 indel（插入缺失）和结构变异的处理能力更强。

# 使用 BWA-MEM
BWA_MEM="/path/to/bwa"
REF_PREFIX="./indexes/grch38_bwa"

# 1. 比对生成 SAM
$bwa_MEM mem \
    -t 16 \              # 线程数
    -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA\tPU:unit1' \  # Read Group 信息，写在参数里更简洁
    $REF_PREFIX \
    $CLEAN_R1 \
    $CLEAN_R2 > sample_raw.sam

# 2. 转换为 BAM
samtools view -@ 16 -bS sample_raw.sam > sample_raw.bam

# 3. 排序
samtools sort -@ 16 -o sample_sorted.bam sample_raw.bam

# 4. 建立索引
samtools index sample_sorted.bam

💡 专家提示：

  • BWA-MEM 默认会自动处理双端读段，并识别插入缺失。
  • 如果你发现比对率异常低（比如低于 70%），首先检查 FASTQ 质量，其次检查参考基因组物种是否匹配（别拿人类基因组比对细菌数据），最后考虑是否有大量重复序列区域。

5. 第四步：后处理与优化（Post-Processing）

比对完生成的 BAM 文件还不能直接用！原始 BAM 里有很多“噪音”和“错误”，必须经过严格的后处理，尤其是对于变异检测。

5.1 标记重复序列（Mark Duplicates）

PCR 扩增会导致同一个原始 DNA 片段被复制成多份，这些“重复读段”会给变异调用带来偏差（比如一个假阳性 SNP 因为来自多个 PCR 副本而被高估）。

工具：samtools markdup 或 Picard MarkDuplicates

# 使用 Picard（功能更全面，适合生信流程）
java -jar /path/to/picard.jar \
    MarkDuplicates \
    INPUT=sample_sorted.bam \
    OUTPUT=sample_dedup.bam \
    METRICS_FILE=sample_dedup_metrics.txt \
    CREATE_INDEX=true \
    VALIDATION_STRINGENCY=SILENT

# 或者使用 samtools（更轻量）
samtools markdup -r sample_sorted.bam sample_dedup.bam
samtools index sample_dedup.bam

METRICS_FILE 会告诉你有多少读段被标记为重复，这可以作为数据质量的一个指标。如果重复率超过 50%，可能说明测序深度过高或起始 DNA 量过低。

5.2 碱基质量重校准（BQSR）—— 可选但推荐

Illumina 测序仪给的质量值往往偏高，存在系统性偏差。GATK 的 BQSR（Base Quality Score Recalibration）可以校正这些偏差。

# 1. 生成 recalibration table
gatk BaseRecalibrator \
    -I sample_dedup.bam \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    --known-sites dbsnp.vcf \      # 已知的 SNP 位点，用于校准
    -O recal_data.table

# 2. 应用校准
gatk ApplyBQSR \
    -I sample_dedup.bam \
    --bqsr-recal-file recal_data.table \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -O sample_final.bam

samtools index sample_final.bam

💡 专家提示：

  • BQSR 需要已知的变异位点数据库（如 dbSNP）。如果没有这些信息，可以跳过此步，直接使用去重复后的 BAM。
  • 对于非人类物种（如小鼠、斑马鱼），如果缺乏高质量的已知位点数据库，BQSR 的效果可能有限，需谨慎使用。

6. 第五步：最终验证与质量控制

在将 BAM 文件交给下游分析（如 VarScan、GATK HaplotypeCaller）之前，务必做一次全面的质量检查。

使用 Samtools stats

samtools stats sample_final.bam > sample_stats.txt

查看关键指标：

  • Total sequences：总读段数
  • Mapped：比对上的读段数
  • Properly paired：正确比对的成对读段比例（理想情况下 >80%）
  • Insert size mean/std：插入片段大小的均值和标准差，用于评估文库构建质量。
  • Quality scores：质量值分布。

使用 QualiMap 或 IGV 进行可视化

# QualiMap 提供图形化报告和覆盖度图
qualimap bamqc -bam sample_final.bam -outdir qualimap_report

打开 IGV（Integrative Genomics Viewer），加载 sample_final.bam 和 .bai 索引，随机挑选几个已知基因位点（如 TP53、BRCA1），肉眼观察比对情况：

  • 读段是否整齐排列？
  • 是否有大量错配（绿色/红色框）？
  • 插入缺失（indel）区域是否正确处理？

7. 常见陷阱与解决方案

陷阱一：比对率极低（<50%）

  • 原因 1：参考基因组与样本物种不匹配。
    • 解决：确认 FASTA 文件的物种来源。
  • 原因 2：数据污染（如宿主 DNA 污染、接头污染）。
    • 解决：用 FastQC 检查序列组成，用 Kraken2 进行物种分类鉴定。
  • 原因 3：读段质量太差。
    • 解决：重新进行质控过滤，提高阈值。

陷阱二：比对率高，但变异检出率异常

  • 原因：未标记重复或未进行 BQSR。
    • 解决：严格执行后处理流程，特别是 MarkDuplicates。
  • 原因：比对参数过于宽松，导致大量错误比对。
    • 解决：收紧比对参数，如增加最小匹配长度、降低最大错配数。

陷阱三：BAM 文件损坏或索引丢失

  • 现象：IGV 无法加载，samtools 报错。
    • 解决：重新索引 samtools index，或使用 samtools fastq 转换回 FASTQ 重新比对。

8. 总结：从 FASTQ 到 BAM 的完整流程图

为了让你一目了然，我用一个简单的流程总结：

”` [原始 FASTQ] –> [FastQC 质控] –> [Trimmomatic 清洗] –> [干净 FASTQ]

                                                  |
                                                  v
                                           [参考基因组索引]
                                                  |
                                                  v
                                           [BWA-MEM/Bowtie2 比对]
                                                  |
                                                  v
                                           [SAM 文件] --> [SAMtools 转换 BAM]
                                                  |
                                                  v
                                           [BAM 排序] --> [Samtools Sort]
                                                  |
                                                  v
                                           [标记重复] --> [Picard MarkDuplicates]
                                                  |
                                                  v
                                           [BQSR 校准] (可选) --> [GATK ApplyBQSR]
                                                  |
                                                  v
                                           [最终 BAM + BAM.BAI 索引]
                                                  |
                                                  v
                                           [Samtools stats
-- 展开阅读全文 --