序列比对与表达定量

好，数据质控做完了，拿到了一堆clean reads。下一步该干嘛？说白了，就是把这个几千万条的短序列，放到参考基因组上，看看它们分别来自哪个基因，每个基因上有多少条。这个过程，我习惯分成两个大步骤：序列比对和表达定量。

你可能觉得，比对不就是BLAST一下吗？其实没那么简单。RNA-seq的比对，要比基因组重测序复杂得多——因为有剪接事件。一个read可能跨越两个外显子，中间还有内含子。你想想看，如果比对软件不支持剪接，那这条read就直接对不上了。

核心要点：RNA-seq比对的本质，就是把片段化的RNA序列，重新映射到基因组上，同时识别出剪接位点。表达定量，就是统计每个基因/转录本上覆盖了多少read。

从read到表达矩阵，流程长这样

我个人喜欢把整个流程画出来，带着团队一起走一遍。你看下面的图，基本就是标准流程了。

嗯，这里要注意。图里画了两条路。一条是先比对到基因组，生成BAM文件，然后再跑定量。另一条是直接用Salmon这种工具，跳过BAM生成环节。哪种好？我后面会聊。

比对工具那么多，怎么选？

市面上的比对工具，少说也有十几个。但最常用的，我排个序：STAR > HISAT2 > TopHat2（已过时）。我曾经拿同一个数据集，对比过STAR和HISAT2的时间跟内存。

工具	速度	内存需求	剪接支持	我的推荐场景
STAR	极快（约2-4小时/样本）	高（约30-50GB）	优秀（基于MAXIC）	人、鼠等大型基因组，批量跑
HISAT2	较快（约3-5小时/样本）	中等（约5-10GB）	优秀（基于FM-index）	笔记本或服务器内存不足时
Salmon	极快（约30分钟/样本）	低（约4-8GB）	基于转录组（需要注释）	想跳过比对，直接定量

避坑指南：STAR跑之前，一定要预先生成基因组索引（genomeGenerate）。我第一次做的时候，忘了这一步，直接跑到报错，白白浪费了一个小时去查问题。另外，索引目录记得留够磁盘空间，人类基因组索引大约需要80GB左右。

我个人习惯，大项目（超过50个样本）优先用STAR。小项目或者只想看看几个基因的表达量，就用Salmon，省时省力。但要注意，Salmon是直接定量到转录本水平，跟基于比对的count方法，结果可能有一点点偏差。我见过案例，同一个基因，FPKM和TPM算出来的差异倍数方向相反，最后发现是定量策略不一致导致的。

定量这一步，到底在做什么？

比对完成后，你会得到一个BAM文件。这个文件记录了每个read在基因组上的位置。定量，就是统计每个基因的外显子区域上，覆盖了多少条read。

最简单的方法，就是用HTSeq-count或featureCounts。它们做的事情很直白：给定一个BAM文件和GTF注释文件，数一数每个基因被多少条read命中。代码长这样：

# 用featureCounts进行基因水平定量
featureCounts -T 8 -p -t exon -g gene_id \
  -a annotation.gtf \
  -o counts_output.txt \
  sample1.bam sample2.bam sample3.bam

这个命令跑完，会生成一个counts矩阵。第一列是基因ID，后面每一列是一个样本的count数。

Geneid	sample1.bam	sample2.bam	sample3.bam
ENSG00000141510	1256	345	789
ENSG00000123456	2345	678	123
ENSG00000012345	789	234	567

但这里有个坑：count数是绝对数值，受测序深度影响很大。同一个基因，样本A测了5000万条read，样本B测了500万条，A的count数肯定高。所以，后续必须要做标准化，比如FPKM、TPM、DESeq2的中位数比值法等。

经验之谈：我建议一开始就决定好你的分析终点。如果只是比较组间差异，用DESeq2或edgeR的count数据就够了，它们内部会自己处理文库大小。但如果你要跨样本比较表达水平，或者做热图展示，建议用TPM。FPKM虽然经典，但它的计算方式有点过时了——FPKM先对基因长度归一化再对总read数归一化，TPM是反过来，所以TPM在不同样本之间更有可比性。

质量控制与常见陷阱

比对和定量做完，不能直接就上差异表达分析。你得先看看比对质量，确认数据是OK的。

比对率：正常情况下，人类样本应该有80%-95%的read能比对到基因组。低于70%，就要怀疑样品污染或建库问题了。
唯一比对率：也就是只比对到一个位置的read比例。multi-mapping reads（多比对read）在重复序列区域常见，定量时通常只保留unique reads。
rRNA污染：如果比对到rRNA基因的read比例超过5%，说明建库时rRNA去除不彻底。我遇到过一次，rRNA污染高达20%，最后只能重新做。
3' bias：检查read覆盖是否集中在转录本3'端。如果是，可能是RNA降解导致的。这种样本，差异分析时要谨慎。

怎么检查这些指标？STAR输出日志里就有。路径是：

# STAR运行日志，关键内容摘录
cat sample1_Log.final.out

# 你会看到类似这样的关键行：
# Uniquely mapped reads % | 91.23%
# % of reads mapped to multiple loci | 4.56%
# % of reads mapped to too many loci | 0.34%
# % of reads unmapped: too many mismatches | 2.11%

嗯，这个文件一定要看。我见过有人跑完STAR直接删掉日志，结果发现一个样本的比对率只有30%，还不知道哪儿出了问题。

实战：我的一次定量对比经历

以前有个项目，我同时用了三种定量方法：STAR+featureCounts、Salmon、kallisto。我想看看哪种方法在差异基因检测上最稳定。结果很有意思：

高表达基因：三种方法结果高度一致，差异倍数几乎一致。
低表达基因：分歧就大了。基于比对的count方法检测到很多低表达基因，但Salmon/kaillsto因为直接处理read的ambiguity，对低表达基因更保守一些。
转录本水平定量：Salmon的优势明显，它直接统计每个转录本的表达量。如果你关注可变剪接，Salmon比基因水平的featureCounts更合适。

核心结论：没有绝对的“最好”工具。你的选择取决于研究目标和可用的计算资源。我个人比较倾向基于比对+count的策略，因为它更经典，后续差异分析工具（DESeq2, edgeR）都是基于count设计的。如果你想节省时间，或者关注转录本表达，可以直接跳过比对，用Salmon。

最后，记住一点：定量结果的质量，很大程度上取决于参考基因组和注释文件的版本。不同版本的GTF，基因ID、染色体命名可能都不一样。我建议，同一个项目内，所有样本都共用同一个版本，不要混用。我以前遇到过这种情况：重复样本用了不同版本的GTF，结果定量到同一个基因的ID不同，排查了半天才发现是注释版本不一致。