正文

测序数据比对流程详解 从原始fastq到变异检测全解析 新手也能看懂的BWA Bowtie2比对实战指南

/0 浏览量
0924

嘿，欢迎来到生物信息学的“翻译”现场。

如果你刚拿到测序仪吐出来的原始数据，看着满屏的 ACGT 感到头大，别慌。你现在的处境就像手里有一堆被打碎的拼图（原始测序片段），而你的任务是把这些碎片拼回原来那本书（参考基因组）里，然后找出这本书里哪里被画了叉、哪里漏了字。这一步，我们叫它序列比对（Sequence Alignment）。

很多新手一听就要跑 Linux 命令、调参数、看日志，直接劝退。但其实，比对的逻辑就像是你玩“找茬”游戏。今天，我们就把这套流程掰开了、揉碎了，用大白话带你走一遍从 fastq 到 VCF 的全过程。

一、 先搞懂：我们在比什么？

在敲任何代码之前，你得心里有数。测序仪出来的数据不是连续的长句子，而是一堆短的“词”。

想象一下，人类基因组有 30 亿个碱基对，但二代测序（NGS）一次只能读 150-300 个碱基。所以，你手里有几千亿条这种短序列。

比对的任务就是：

  1. 定位：这一小条序列，到底属于基因组的哪个位置？
  2. 对齐：它和参考基因组是完全匹配，还是有差异（突变）？

常用的“工具手”有两个：BWA 和 Bowtie2。它们都是将短序列快速“映射”到参考基因组上的神器。

小知识：BWA 适合长读长或基因组全测序（WGS/WES），而 Bowtie2 常用于 RNA-seq 或 ChIP-seq。但对于新手入门，两者原理相通，我们先以业界标准的 BWA-MEM 为主角，顺便聊聊 Bowtie2 的差别。

二、 战前准备：磨刀不误砍柴工

在开始之前，你需要三样东西：

  1. 参考基因组（Reference Genome）：比如人类的 GRCh38，通常是个 .fasta 文件。
  2. 测序数据：你的样本，通常是 .fastq 或 .fastq.gz 格式。
  3. 工具：BWA、Samtools、GATK 等。

1. 建立索引（Index）

这是最关键的一步，也是新手最容易漏的一步。参考基因组必须先生成索引，否则比对软件找不到北。

你可以把它想象成查字典前要先把字典编好目录。BWA 和 Bowtie2 的索引格式不一样，别搞混了。

使用 BWA 建立索引：

# 假设你的参考基因组是 hs38DH.fa (包含常用污染的灰度基因组)
# 1. 先转化为 bwa 可识别的索引
bwa index hs38DH.fa

# 运行后会生成这些文件，别慌，这是正常的：
# hs38DH.fa.amb
# hs38DH.fa.ann
# hs38DH.fa.bwt
# hs38DH.fa.pac
# hs38DH.fa.sa

使用 Bowtie2 建立索引（如果你决定用 Bowtie2）：

# Bowtie2 的索引命令不同，生成的是 .bt2 文件
bowtie2-build hs38DH.fa hs38DH_idx

专家提示：索引文件一旦生成，可以复用。下次再比对同一物种的样本时，直接调用索引即可，不需要重新构建。

三、 核心战场：开始比对

现在，我们要把样本的 FASTQ 文件“塞”进参考基因组里。这里我们以 BWA-MEM 算法为例，它是目前最主流的选择，精度高且支持长读长。

场景 A：单端测序（Single-end）

有些老数据或特定实验是单端的，读一条链。

bwa mem -t 16 hs38DH.fa sample_R1.fastq.gz > sample.sam
  • -t 16：使用 16 个线程加速（看你服务器配置，别把机器跑崩了）。
  • > sample.sam：输出结果为 SAM 格式（见下文解释）。

场景 B：双端测序（Paired-end）—— 最常见！

现在绝大多数 WGS/WES 都是双端测序，即从 DNA 片段的两端分别测序（R1 和 R2）。双端比对能大幅提高准确性，因为两端距离已知，就像你知道绳子的两头在哪。

bwa mem -t 16 -R '@RG\tID:Sample001\tSM:Sample001\tPL:ILLUMINA' \
  hs38DH.fa \
  sample_R1.fastq.gz \
  sample_R2.fastq.gz \
  > sample.sam
  • -R：这一串是读取标签（Read Group）。别嫌麻烦，这对后续的变异检测至关重要！它告诉后续流程：“这个样本叫 Sample001，平台是 Illumina”。如果漏掉这一步，GATK 可能会报错或忽略你的数据。
  • 注意 R1 和 R2 的顺序，必须先 R1 后 R2，否则结果会乱。

如果想用 Bowtie2 呢？

bowtie2 -p 16 -x hs38DH_idx -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
        --very-sensitive-local -S sample.sam
  • --very-sensitive-local：这是模式。Bowtie2 有 --end-to-end（强制两端都对齐）和 --local（允许软剪接，适合 RNA-seq 或存在大插入缺失的情况）。对于 DNA 变异检测，通常推荐 --very-sensitive 或 --local。

四、 垃圾清理：SAM 转 BAM，排序与去重

停！ 现在的 sample.sam 还只是原始的对齐记录，里面混杂着垃圾信息，而且顺序是乱的。直接拿去做变异检测会被老板（算法）骂死的。我们需要经过一道“精加工”流水线。

1. SAM 转 BAM 并排序

SAM 是文本格式，打不开也读不动；BAM 是二进制压缩格式，体积小、速度快。

# 转换并排序
samtools view -@ 16 -bS sample.sam | samtools sort -@ 16 -o sample.sorted.bam

# 建立 BAM 文件的索引，方便后续快速检索
samtools index sample.sorted.bam

2. 标记重复序列（Mark Dupicates）

这是一个极其重要但常被新手忽略的步骤。 测序过程中，同一个 DNA 片段可能被扩增多次（PCR 重复）。如果不剔除，这些重复片段的错误会误导变异识别，让你把测序错误当成真实突变。

# 使用 picard 工具标记重复
java -jar picard.jar MarkDuplicates \
  I=sample.sorted.bam \
  O=sample.dedup.bam \
  M=sample.metrics.txt \
  CREATE_INDEX=true

# 再次建立索引
samtools index sample.dedup.bam

注意：现在 GATK 更推荐用 GATK MarkDuplicates，但在很多老流程中 picard 依然通用。

五、 锦上添花：碱基质量重校准（BQSR）

这一步是 GATK 流程的精髓。测序仪给出的碱基质量值（Quality Score）往往不准，比如明明质量很低，却给了个 Q30。BQSR（Base Quality Score Recalibration）会通过统计已知位点（如 dbSNP 数据库中的常见变异）来修正这些系统误差。

# 1. 生成重校准表
gatk BaseRecalibrator \
  -I sample.dedup.bam \
  -R hs38DH.fa \
  --known-sites dbsnp.vcf \
  -O recal_data.table

# 2. 应用重校准
gatk ApplyBQSR \
  -I sample.dedup.bam \
  -R hs38DH.fa \
  --bqsr-recal-file recal_data.table \
  -O sample.recal.bam

经过这一步，你的 BAM 文件质量会显著提升，假阳性变异会大幅减少。

六、 终极目标：变异检测（Variant Calling）

准备好了吗？现在我们要把那些“不一样的地方”找出来。这里的“不一样”主要包括：

  • SNV（单核苷酸变异）：A 变成了 G。
  • Indel（插入缺失）：多了一个碱基或少了一个碱基。

使用 HaplotypeCaller（GATK 推荐）

GATK 的 HaplotypeCaller 不是简单地看每个位置，而是会在局部重新组装单倍型，这对处理 Indel 特别有效。

gatk HaplotypeCaller \
  -R hs38DH.fa \
  -I sample.recal.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF

关键点：这里我们生成的是 .g.vcf.gz（基因组 VCF），而不是直接的 .vcf。这是为了联合分析（Joint Calling）做准备。如果你只有这一个样本，可以直接把 GVCF 转成 VCF： gatk GenotypeGVCFs -R hs38DH.fa -V sample.g.vcf.gz -O sample.vcf.gz

如果是多样本分析呢？

这就是 GATK 强大的地方。你可以对 100 个、1000 个样本分别生成 gVCF，然后一起调用：

# 1. 合并所有样本的 gVCF
gatk CombineGVCFs -R hs38DH.fa -V sample1.g.vcf.gz -V sample2.g.vcf.gz -O combined.g.vcf.gz

# 2. 联合基因型调用
gatk GenotypeGVCFs -R hs38DH.fa -V combined.g.vcf.gz -O final_variants.vcf.gz

七、 过滤与注释：从海量结果中提炼金矿

final_variants.vcf.gz 里可能有几万个变异，但其中大部分是噪音。我们需要过滤。

1. 硬过滤（Hard Filtering）

对于非癌症体细胞突变，常用的规则如下：

gatk VariantFiltration \
  -R hs38DH.fa \
  -V final_variants.vcf.gz \
  -O filtered_variants.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 60.0" --filter-name "HighFS" \
  --filter-expression "MQ < 40.0" --filter-name "LowMQ" \
  --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRankSum" \
  --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowReadPosRankSum"

2. 功能注释

过滤完后，你得知道这些变异到底在基因哪里，有什么后果。

# 使用 SnpEff 进行注释（需要下载对应物种的数据库）
java -jar snpeff.jar GRCh38.99 filtered_variants.vcf.gz > annotated_variants.vcf

输出文件里，你会看到类似这样的信息：

  • Type：Missense（错义，氨基酸变了）、Synonymous（同义，没变）、StopGained（终止密码子，后果严重）。
  • Allele Frequency：变异频率。
  • Clinical Significance：如果连接了 ClinVar 数据库，还会告诉你这个变异是否致病。

八、 新手常见“坑”与避坑指南

  1. 参考基因组版本不一致：

    • 这是最致命的错误。确保你的 FASTA、BAM 索引、VCF 注释全部基于同一个版本（如 GRCh38 或 hg19）。混用会导致坐标错位，结果完全错误。

  2. 忽略 Read Group：

    • 再次强调，-R 参数必须加。没有 RG 信息，后续许多统计工具（如 GATK 的 BQSR 和 CNV 检测）会直接报错。

  3. 线程数设置不当：

    • 别为了快把所有核心都占满。留给操作系统和其他进程一点空间，通常使用物理核心的 70%-80% 比较稳妥。

  4. 文件路径含空格或特殊字符：

    • 编程界的一条铁律：路径中不要有空格。使用 my_analysis/sample_01.fastq 而不是 my analysis/sample 01.fastq。

  5. 盲目相信自动流程：

    • 虽然 nf-core 等流程很强大，但作为学习者，一定要亲手跑一遍基础命令，理解每一步在做什么。否则出了错，你连日志都看不懂。

结语：数据只是开始

从 FASTQ 到 VCF，你完成了一次从“数字”到“生物学意义”的跨越。但这只是第一步。接下来，你还需要进行家系分析、群体频率筛选、功能预测等，才能最终找到致病变异。

记住，生物信息学不仅是一套命令的堆砌，更是一门关于数据质量和生物学逻辑的艺术。当你在屏幕上看到那几个高置信度的致病突变时，那种成就感，值得你熬的所有夜。

保持好奇，保持严谨，祝你在比对的海洋里乘风破浪！

-- 展开阅读全文 --