正文

新手测序数据比对流程常见问题：参考基因组版本如何选择比对质量评分怎么设置避免错误解读结果快速排查比对失败原因

/0 浏览量
0901

嘿，我是 Agnes。我知道你刚拿到一堆 FASTQ 文件，面对那个巨大的比对流程有点头大。别慌，这就像第一次开车，熄火、熄火、还是熄火……但只要摸清了脾气，其实很顺手。咱们不整那些虚头巴脑的教科书定义，直接聊几个让新手最头秃的问题：参考基因组怎么选、评分怎么调、结果怎么看、出错了咋办。

1. 参考基因组版本：选错了，后面全白跑

很多新手第一个坑就是：“我随便下一个基因组不就行了吗？” 真的不行。想象一下，你用 Android 的手机壳去套 iPhone 的手机，虽然长得像，但尺寸、接口都对不上，根本装不进去。

为什么版本这么重要？

参考基因组（Reference Genome）就是你用来“对齐”数据的地图。如果地图版本不对，或者区域缺失，你的数据就“找不到家”。

举个例子： 假设你在研究人类基因 BRCA1。如果你用的参考基因组是 GRCh37 (hg19)，而这个基因在 GRCh38 (hg38) 里增加了一些新的补丁序列（patch），那么你在 hg19 上比对时，可能会漏掉这部分区域，或者比对质量很差。更糟糕的是，不同版本的染色体命名方式都不一样！

  • GRCh37: chr1, chrX, chrY (有些是 1, X, Y)
  • GRCh38: chr1, chrX, chrY (通常带 chr 前缀)

如果你下载了一个叫 Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 的基因组，但你的比对工具（比如 BWA-MEM）默认期望的是 chr1 而不是 1，那你比对时就会报错：“找不到 scaffold 1”。

怎么选？记住这三个原则：

  1. 跟着项目走：如果你是在复现别人的研究，或者使用公开的数据集（比如 TCGA、GTEx），必须使用他们当初用的那个版本。去查他们的 Methods 部分，或者看数据库记录。这是铁律。
  2. 新项目，选最新的：如果你是第一次分析一个物种，没有历史包袱，永远优先选择最新的组装版本。对于人类来说，那就是 GRCh38 (hg38)。它更完整，错误更少，注释更丰富。
  3. 配套工具链：参考基因组要和你用的注释文件（GTF/GFF）、SNP 位点文件（VCF）配套。比如，你用 GRCh38 的基因组，就得用 Ensembl 或 GENCODE 提供的 GRCh38 版本的 GTF 文件。混搭会导致下游分析（比如基因定量）完全乱掉。

小贴士：下载基因组时，记得看清是否有 chr 前缀。很多老教程用 hg19，基因组里没有 chr 前缀，而 BWA-MEM 在新版本下可能要求有。如果不确定，可以用 sed 命令简单处理一下，或者在比对前统一格式。

2. 比对质量评分：别被 0 和 1 的数值迷惑

比对工具（比如 Bowtie2, BWA-MEM, STAR）会给每个读段（read）打一个分，叫 MAPQ（Mapping Quality）。这个分数看起来很高大上，但其实新手最容易误解它。

MAPQ 到底代表什么？

MAPQ 是一个 Phred-scaled 的置信度分数，它告诉你：“这个读段比对到这个位置，出错的概率是 10^(-MAPQ/10)”。

  • MAPQ = 0：比对质量极差，可能比对了多个位置，或者根本不可信。
  • MAPQ = 10：10% 的概率比对错误。
  • MAPQ = 20：1% 的概率比对错误。
  • MAPQ = 30：0.1% 的概率比对错误。
  • MAPQ = 60：0.000001% 的概率比对错误（几乎可以确定是对的）。

新手常见的错误解读

错误 1：认为 MAPQ=0 的读段就没用。 其实不一定。MAPQ=0 通常意味着这个读段可以比对到基因组多个位置（多映射读段，multi-mapping reads）。在重复区域、旁系同源基因区域，这是非常正常的。如果你在做 RNA-seq 定量，很多生物信息学家会选择保留这些读段，用概率模型分配权重（比如 Salmon, Kallisto 做的就是这样）。但如果你在做 SNP calling，MAPQ=0 的读段确实应该丢弃，因为它们的位置不确定。

错误 2：只看 MAPQ，不看比对本身。 有时候 MAPQ 很高，但比对位置可能不对。比如，一个读段比对到了一个假基因（pseudogene）上，而不是它的真实基因座。这时候 MAPQ 可能很高，但生物学上是错的。所以，一定要结合其他指标，比如比对的特异性（specificity）和错误率。

怎么设置参数？

不同的比对工具参数不同。以 BWA-MEM 为例：

  • 默认参数通常够用：对于大多数常规 Whole Genome Sequencing (WGS) 或 Whole Exome Sequencing (WES) 分析，BWA-MEM 的默认参数已经优化得相当好。
  • 调整 -T (最小比对分数)：-T 参数设置最小比对分数，低于这个分数的读段会被丢弃。默认是 30。如果你的数据质量很差，或者你想保留更多读段（比如在 metagenomics 中），可以适当降低 -T。但注意，降低 -T 会增加假阳性。
  • 调整 -q (丢弃低质量比对)：-q 设置 MAPQ 阈值，低于这个 MAPQ 的读段会被标记为未比对（unmapped）。默认是 0。如果你想过滤掉多映射读段，可以设置为 20 或 30。

建议：先跑一遍默认参数，检查结果。如果比对率太低，再考虑调整参数。不要一上来就乱调参数，那样你都不知道问题出在哪。

3. 快速排查比对失败原因

比对失败了，别慌。先看日志，再看统计，最后看细节。

第一步：看比对率（Alignment Rate）

比对完成后，第一眼看的就是比对率。比如 BWA-MEM 输出的 SAMR 列，或者 samtools flagstat 的结果。

  • 正常范围：
    • WGS: 95% - 99%
    • WES: 80% - 95% (因为捕获探针只覆盖部分基因组)
    • RNA-seq: 70% - 90% (取决于物种、注释完整性、内含子-外显子剪接等)
  • 如果比对率异常低 (< 80% for WGS):
    1. 检查参考基因组是否匹配物种：你是不是把人类数据比对到了小鼠基因组上？（听起来很蠢，但真有人这么干。）
    2. 检查基因组版本：参考基因组和注释文件版本是否一致？
    3. 检查数据质量：用 FastQC 看看测序质量。如果质量太差，或者有大量 adapter 污染，比对率会下降。这时候需要质控和去 adapter。
    4. 检查是否有污染：如果数据里有大量细菌、病毒序列，而参考基因组只有人类，这些序列就比对不上了。

第二步：看 MAPQ 分布

用 samtools view -q 20 或者 samtools flagstat 看看高 MAPQ 的读段比例。如果大部分读段 MAPQ=0，那可能是参考基因组有太多重复序列，或者读段太短、太简单，无法唯一确定位置。

第三步：看未比对读段（Unmapped Reads）

把未比对的读段提取出来，用 BLAST 比对一下 NCBI 的 nt 数据库，看看它们到底是什么。这样能帮你快速定位问题：

  • 如果是 adapter 序列：需要重新做质控。
  • 如果是其他物种序列：可能有污染。
  • 如果是未知序列：可能需要更新参考基因组，或者你的物种参考基因组本身就有问题。

第四步：看比对工具的错误日志

比对工具通常会输出一些警告信息，比如“too many mismatches”、“no valid mates”等。这些信息虽然看起来不起眼，但往往藏着关键线索。

4. 几个实战小技巧

  • 索引要建对：每个比对工具都需要参考基因组的索引文件。确保你用的索引文件是和你的参考基因组 FASTA 文件对应的。有时候下载了基因组 FASTA，但忘了建索引，或者建索引的命令错了（比如 BWA 需要 bwa index，而 STAR 需要 STAR --runMode index），这些都会导致失败。
  • 多线程加速：现代比对工具都支持多线程。别只开一个线程，那样跑得慢死。根据你服务器的 CPU 核心数，合理设置线程数。
  • 结果验证：比对完成后，不要直接扔进下游分析。用 IGV (Integrative Genomics Viewer) 打开比对结果，随机看几个基因位点，看看比对是否合理。比如，看外显子区域是否有连续的比对，内含子区域是否有剪接比对（spliced alignment）。这能帮你发现一些统计指标看不出的问题。

结语

比对流程是基因组学分析的基石，基石不稳，地动山摇。希望这些经验能帮你少走弯路。记住，遇到问题先看日志，再查统计，最后看细节。多问，多试，多总结。加油！

-- 展开阅读全文 --