4. 序列比对原理：参考基因组选择、索引构建、比对算法简介

拿到测序仪吐出来的 reads 之后，第一步是什么？

我个人的习惯，是先安静地想清楚一个问题：我要把这些短片段贴到哪张地图上？

这个“地图”就是参考基因组。没有它，你手里的几千万条 reads 就是一盘散沙。几年前我第一次做 RNA-seq 时，随便下了一个参考基因组版本就跑流程，结果发现基因注释老对不上。后来才意识到，参考基因组的选择直接影响后面所有分析的质量。

好了，咱们一步步拆解。序列比对的核心流程可以概括为三步：选基因组 → 建索引 → 做比对。

序列比对核心流程 序列比对核心流程 参考基因组选择 索引构建 比对算法执行 SAM/BAM 输出 核心：索引构建是比对效率的关键，相当于给参考基因组建立一本“电话簿”

4.1 参考基因组选择——选对地图，事半功倍

参考基因组说白了就是一个物种的“标准模板”。通常以 FASTA 格式存放，里面是一堆染色体序列和 scaffolds。选择时我主要看三个指标：组装质量、版本号、文件格式。

考量维度说明我的建议
组装质量基因组是否完整，gap 多不多优先选染色体级别 (chromosome-level) 的组装
版本号GRCh38 vs GRCh37 / hg38 vs hg19新项目一律用最新主要版本
文件格式FASTA + GTF/GFF 注释是否配套检查注释文件的基因名格式是否一致
我曾经踩过的坑： 下载了参考基因组 FASTA，但注释文件是从另一渠道拿的。比对时发现染色体命名方式不一致——有些叫“chr1”，有些叫“1”。结果大量 reads 比对不上，白白浪费了两天时间。所以记住：基因组和注释文件要来自同一来源。

常见的参考基因组来源有几个：Ensembl、NCBI RefSeq、UCSC Genome Browser。我个人更习惯用 Ensembl，因为它的注释比较完整，且版本命名清晰。下载时一般这样操作：

# 以人类 GRCh38 为例，从 Ensembl 下载
wget ftp://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 同时下载配套注释
wget ftp://ftp.ensembl.org/pub/release-110/gtf/homo_sapiens/Homo_sapiens.GRCh38.110.gtf.gz
如果你做的是非模式生物，比如某种细菌或者植物，建议选 RefSeq 版本。因为 RefSeq 对非模式生物的注释质量做得更均匀。

4.2 索引构建——给基因组建立“电话簿”

为什么要有索引？你想想看，几千万条 reads 逐条去跟全长 3Gbp 的基因组比对，那太慢了。索引的作用就是把基因组序列转化成一个可以快速检索的数据结构。

最流行的索引方案是 Burrows-Wheeler Transform (BWT) 配合 FM-index。说白了，就是先把基因组序列变换一下，再压缩存储，让搜索过程更快、更省内存。这个领域的代表作就是 BWA 和 Bowtie2 这两个工具。

不同比对工具的索引构建命令略有不同：

# BWA 索引
bwa index Homo_sapiens.GRCh38.dna.primary_assembly.fa

# 运行后你会得到 5 个文件：
# .amb, .ann, .bwt, .pac, .sa

# Bowtie2 索引
bowtie2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa GRCh38

# 运行后生成 .bt2 系列文件

需要注意，索引一旦建好就可以重复使用。我通常会在服务器上建立一个共享目录，比如 /ref/genomes/，把常用物种的索引都放进去，这样不同项目之间就不用重复建索引了。建一次索引的时间取决于基因组大小，人类基因组用 BWA 大概需要 2-3 小时。

索引文件体积通常比原始 FASTA 大 1.5-2 倍。人类 GRCh38 的原始 FASTA 大约 3GB，建完 BWA 索引后总大小在 5-6GB 左右。磁盘空间不够的话，可以考虑用压缩格式或选择更轻量的索引方式。

4.3 比对算法——三种思路，各有千秋

市面上主流的比对工具虽然多，但核心算法思路基本可以归为三类。我根据实际使用经验整理了一个对比表：

算法类型代表工具适用场景核心原理速度和准确率
BWT + FM-indexBWA-MEM, Bowtie2常规 WGS、WES、RNA-seq基于 Burrows-Wheeler 变换的精确匹配，再允许少量错配和 gap速度快，内存低，准确率不错
哈希表（Hash-based）BWA-backtrack, MAQ短 reads（≤100bp）、早期数据分析将 reads 的 k-mer 或基因组位置映射到哈希表中短 reads 时速度较快，长 reads 时性能下降明显
Seed-and-extendSTAR, HISAT2RNA-seq 剪接比对、长 reads先找小段精确匹配（seed），再向两端延伸并合并能处理剪接位点和结构变异，但内存消耗大

我最常用的是 BWA-MEM。它能处理 70bp 到 1Mbp 的 reads，准确率也很稳。RNA-seq 场景下我会换用 STAR，因为它能高效处理剪接位点的比对。操作起来很简单：

# BWA-MEM 比对（单端）
bwa mem -t 8 Homo_sapiens.GRCh38.dna.primary_assembly.fa sample_R1.fastq.gz > aligned.sam

# BWA-MEM 比对（双端）
bwa mem -t 8 Homo_sapiens.GRCh38.dna.primary_assembly.fa sample_R1.fastq.gz sample_R2.fastq.gz > aligned.sam

# 参数说明
# -t 8：使用 8 个线程
# -M：短 split reads 标记为 secondary（Picard 流程常用）

比对完成后你会得到一个 SAM 文件。这里面每一行记录了一条 read 的比对信息：比对位置、CIGAR 字符串、MAPQ（比对质量分数）等等。我一般会先看一下 MAPQ 的分布，大致评估比对质量。

关于 CIGAR 字符串：比如 "3M1D5M" 意思就是前 3 个碱基匹配，然后 1 个碱基缺失，再 5 个碱基匹配。它记录了 read 和参考基因组之间的对齐细节。理解 CIGAR 能帮你快速诊断比对中的问题——比如大量 insert 出现，可能说明 reads 中有污染序列。

4.4 实操中的几个小经验

  • 检查重复率：比对跑完后，用 samtools flagstat 看一下有多少 reads 比对不上。如果比对率低于 70%，先检查参考基因组版本是否正确。
  • 内存监控：BWA-MEM 用 8 线程时内存大概 5-6GB，STAR 起步就在 30GB 以上。跑之前用 top 或 htop 看一眼。
  • 输出格式转换：SAM 文件很大，通常转成 BAM（二进制）来节省空间。这一步用 samtools view 就能搞定。
# 查看比对统计
samtools flagstat aligned.sam

# 转换 SAM 到 BAM，并排序
samtools view -bS aligned.sam | samtools sort -o aligned.sorted.bam
samtools index aligned.sorted.bam

我记得有一次做微生物组分析，参考基因组用的是从 NCBI 下的一个不完整组装版本。跑完比对后一看，70% 的 reads 都没比对上去。换了一个更完整的参考版本后，比对率直接升到 95%。这给我的教训就是：花 10 分钟选择合适的参考基因组，能省掉 10 小时的排错时间。

好，这章的内容就到这儿。核心点你们记一下：选对参考基因组版本 → 建好索引 → 用对的比对工具。这三步走稳了，后面的特征矩阵才能建得扎实。

← 上一章 📖 返回目录 下一章 →