序列比对：STAR / HISAT2 比对原理与实操归纳总结

好，进入第四章。前面我们拿到了经过质控和清洗的高质量 reads，下一步就是把它们“贴”回到参考基因组上。这个步骤叫序列比对，业内也叫 mapping。

其实说白了，RNA-seq 数据比对跟 DNA 比对不太一样。DNA 比对只要对齐就行，但 RNA 有内含子，reads 是跨外显子的剪接产物。你想想看，一条 read 前半段跟外显子 A 匹配，后半段突然跳到外显子 B，中间还隔着一个内含子——普通比对工具直接傻眼。

所以，专门给 RNA-seq 设计的比对器必须能处理剪接信号。今天我们就聊目前最主流的两个工具：STAR 和 HISAT2。我个人在十几个项目里都用过它俩，各有各的脾气。

为什么要用剪接感知比对器？

核心任务就一条：找到 reads 在基因组上的最佳比对位置，同时识别跨越内含子的剪接 junction。

这个过程简单拆解一下：

• 索引构建：把参考基因组建成可快速搜索的数据结构
• 种子搜索：用 read 的前半段在索引里找匹配位置
• 延伸与剪接发现：从种子向两端延伸，遇到不一致就尝试跨内含子跳跃
• 多比对处理：对于重复区或同源区，给每个比对打分

STAR 和 HISAT2 都做这三件事，但方法完全不同。

STAR 的核心原理

STAR，全称 Spliced Transcripts Alignment to a Reference。它最出名的就是快，而且比对率极高。怎么做到的？

STAR 先为参考基因组构建一个后缀数组 (Suffix Array) 或者近似索引。然后对每条 read，它找出最长的一段能精确匹配基因组的种子，再从种子两端向外延伸。延伸时一旦遇到 mismatch，它不直接放弃，而是试图跳过一段序列——这就是剪接检测。

我记得第一次跑 STAR 时，感觉它像一头猛兽。内存吃得很厉害，但出结果的速度让我直接改掉了以前用 TopHat 的习惯。后来我做过一个批量大型 RNA-seq 项目，一个样本平均 5000 万条 paired-end reads，STAR 单线程跑完只需要 4-6 小时，而且 junction 文件里的剪接位点准确率非常高。

HISAT2 的核心原理

HISAT2 是 HISAT 的升级版。它的设计初衷就是省内存。它用了层级 FM-index (Hierarchical Index)：一个全局索引覆盖整个基因组，再加上数万个局部索引覆盖外显子区域和剪接位点周围。

比对过程是这样的：先用全局 FM-index 快速定位 read 的大致位置，然后切换到局部索引完成精确比对。遇到跨内含子匹配时，HISAT2 利用图索引 (graph index) 来编码可能的剪接模式。

我一直觉得 HISAT2 在内存受限的环境下特别实用。有次我在一个只有 16 GB 内存的服务器上跑比对，STAR 直接报错说内存不够。换成 HISAT2 后，索引加载只用了一半内存，顺利跑完。那感觉，嗯，真香。

实操：构建参考基因组索引

不管你用 STAR 还是 HISAT2，第一步都是建索引。下面的命令会花几十分钟到几小时，看服务器性能。

STAR 索引构建

# 下载参考基因组和注释文件（假设已有）
# 示例：人类 hg38
STAR --runMode genomeGenerate \
     --genomeDir /path/to/star_index_hg38 \
     --genomeFastaFiles /path/to/hg38.fa \
     --sjdbGTFfile /path/to/gencode.v39.annotation.gtf \
     --sjdbOverhang 100 \
     --runThreadN 20

参数说明：

• --sjdbOverhang 100：通常设为 read 长度减 1，比如 150 bp 的 reads 就用 149
• --runThreadN：越大越好，但别超过物理核心数

HISAT2 索引构建

# HISAT2 需要额外构建剪接位点文件
extract_splice_sites.py /path/to/gencode.v39.annotation.gtf > splicesites.txt
extract_exons.py /path/to/gencode.v39.annotation.gtf > exons.txt

# 构建索引
hisat2-build -p 20 \
             --snp /path/to/common_snps.vcf \
             /path/to/hg38.fa \
             /path/to/hisat2_index_hg38 \
             --ss splicesites.txt \
             --exon exons.txt

实操：运行比对

建好索引后，就可以开始跑了。以下是典型命令。

STAR 比对命令

STAR --genomeDir /path/to/star_index_hg38 \
     --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
     --readFilesCommand zcat \
     --outFileNamePrefix sample1_ \
     --outSAMtype BAM SortedByCoordinate \
     --outBAMcompression 6 \
     --outFilterMultimapNmax 10 \
     --runThreadN 20

几个关键参数：

• --outFilterMultimapNmax 10：只保留比对到 10 个以下位置的 reads，过滤重复区干扰
• --outSAMtype BAM SortedByCoordinate：直接输出排序后的 BAM

HISAT2 比对命令

hisat2 -x /path/to/hisat2_index_hg38 \
       -1 sample_R1.fastq.gz \
       -2 sample_R2.fastq.gz \
       -S sample1.sam \
       --threads 20 \
       --known-splicesite-infile splicesites.txt

运行后别忘了转 BAM 和排序：

samtools view -bS sample1.sam > sample1.bam
samtools sort sample1.bam -o sample1_sorted.bam
samtools index sample1_sorted.bam

STAR vs HISAT2：对比清单

特性	STAR	HISAT2
索引结构	后缀数组 / SA-Index	层级 FM-index + 图索引
内存需求（人类基因组）	~30 GB	~8 GB
比对速度	极快（通常 2-4 小时/样本）	较快（通常 3-6 小时/样本）
剪接检测能力	灵敏度高，支持长内含子	灵敏度高，尤其依赖剪接位点文件
多比对处理	通过 --outFilterMultimapNmax 控制	通过 -k / --max-seeds 控制
推荐场景	服务器内存充足、追求速度	内存有限、多项目共享索引

比对后的文件解析

比对完后的 BAM 文件里，每条记录对应一条 read 的比对信息。重点关注以下字段：

• RNAME：比对到的染色体名称
• POS：比对起始位置（1-based）
• CIGAR：比对模式，比如 150M 表示完美匹配 150 bp，75M50000N75M 表示跨内含子的剪接
• NH / HI tag：记录该 read 的比对次数和排名

你想想看，如果 CIGAR 里出现大段的 N，比如 50000N，那就是一条跨内含子读段，N 的长度就是内含子长度。STAR 和 HISAT2 的 junction 文件其实就是把这些 N 信息汇总出来的。

常见问题与解决策略

• 比对率偏低（< 70%）？
  先查 RNA 完整性。我在一个项目中遇到过 RNA 降解严重的样本，比对只有 50%。重新建库后恢复正常。
• STAR 内存不够怎么办？
  降低 --genomeSAindexNbases 的值（比如从 14 降到 12），索引会更小，但搜索变慢。
• HISAT2 比对特别慢？
  检查剪接位点文件是否正确。如果没提供 --ss，HISAT2 会内部尝试推测剪接，耗时剧增。
• 多比对 reads (multimappers) 如何处理？
  我一直的做法是：保留唯一比对用于定量，多比对的 reads 按比例分配。这个后面课时再细说。