RNA-seq变异检测实战：Sentieon全流程详解

前段时间接了一个RNA-seq项目，需求不是常规的找差异基因，而是直接做变异检测，说白了就是要在转录组数据里把样本的SNV和InDel从BAM里挖出来，然后输出一份带注释的突变列表。做之前我查了一圈方案，最后落到了Sentieon上。之所以选它，是因为这套流程跑RNA-seq变异检测时，在去重和变异识别阶段的速度优势非常明显，而且输出的VCF可以直接套用GATK体系下的注释和过滤生态，不需要额外折腾格式转换。这篇文章就把整个流程从头到尾拆一遍，适合手里有转录组FASTQ、还没想清楚怎么挖变异的人，也适合被GATK HaplotypeCaller跑得怀疑人生的同行。

1. 为什么RNA-seq变异检测要单独讲一套流程

1.1 RNA-seq变异检测与DNA-seq的典型差异

不少朋友刚接触这个需求时，第一反应是“我平时DNA测序的变异检测流程很熟，把reads比上去然后call variant不就行了”。实际上RNA-seq变异检测还真不能照搬DNA-seq那套，核心差异集中在两个地方：比对环节和变异过滤环节。

DNA-seq的reads基本都来自基因组连续区段，比对工具只需要考虑线性匹配。RNA-seq的reads则大量跨越外显子连接区域，一段150bp的reads可能前80bp落在某个外显子，后70bp落在下一个外显子，中间还要跨过内含子。这时候如果用普通的BWA或Bowtie2去做比对，跨剪接位点的reads会直接丢失或比对到错误位置，变异检测也就无从谈起。所以第一步就必须换成STAR这类支持剪接联配的比对器，后续所有结果都建立在这个基础上。

过滤环节差异更明显。RNA-seq的覆盖度由基因表达水平决定，高表达基因可能有几百上千X深度，低表达基因可能只有几个X，某些外显子区域甚至完全没覆盖。直接用DNA-seq那套基于深度的过滤阈值，会把大量真实变异当成“低深度”误删。同时RNA-seq还会引入RNA编辑事件，最典型的是ADAR导致的A到I编辑，在测序数据里看起来和A到G的SNV几乎一样，这一点在过滤时要特别留意，不然报告里会多出一批“假突变”。

1.2 为什么用Sentieon来做这条流程

Sentieon本质上是一套基于GATK工具语义重构的生物信息分析软件，它的命令和GATK高度兼容，但在底层算法和工程实现上做了大量优化。对RNA-seq项目来说，最直观的收益就是速度。同样一份转录组数据，用GATK HaplotypeCaller跑可能要几十小时，换成Sentieon Haplotyper往往几个小时就能完成，而且输出VCF的SNV和InDel位点与GATK一致性非常高，可以无缝使用GATK后续的过滤、注释工具链。

另一个实际优势是内存占用更可控。RNA-seq的BAM文件通常很大，200M条reads的BAM动辄六七十GB，用GATK处理这类文件时对服务器内存很敏感。Sentieon的driver模式支持流式读取和分区域处理，在实际项目中我常用40线程跑，内存峰值大约在60GB到100GB之间，比起GATK那种动不动就要200GB以上内存的场景，友好很多。

1.3 一个典型的RNA-seq变异检测项目设定

为了方便后面展开命令和参数，我先明确一下这篇文章对应的典型项目假设：人类转录组样本，Illumina双端150bp测序，FASTQ总量约60M对reads，目标是检测SNV和小于50bp的InDel。参考基因组使用GRCh38，加上对应的GTF注释文件。软件环境为Linux服务器，Sentieon版本建议用202112及之后的release，STAR版本选用2.7.x。

之所以选这个设定，是因为它覆盖了绝大多数转录组变异检测场景：既不是超大规模队列，也不是单基因局部测序。在这个量级下，Sentieon的处理耗时大约在3到6小时之间，STAR比对约 1小时，资源分配合理，非常适合作为一篇文章的演示基准。如果你的数据是肿瘤样本做DNA/RNA联合分析，整体流程会稍复杂一些，需要额外考虑配对的正常样本，但RNA-seq这一侧的核心逻辑和本文一致。

2. 开工前的准备：参考基因组、索引与输入数据

2.1 参考基因组与GTF注释文件的选型

RNA-seq变异检测对参考基因组的要求比DNA-seq更“挑剔”，因为STAR在构建索引时不仅需要基因组序列，还需要知道基因结构信息，也就是外显子和剪接位点位置。这个信息来自GTF/GFF3注释文件。常见选型有两种方案：Ensembl的GTF，或者GENCODE的GTF。我建议优先用GENCODE，因为它整合了Ensembl和RefSeq的注释，且对Havana手工注释基因覆盖得更好，对变异检测时过滤基因间区假阳性很有帮助。

版本上需要注意GRCh38的GTF要和参考基因组fasta对应，不能GRCh37的基因组配GRCh38的注释文件，否则STAR索引构建时会出现大量转录本比对失败或基因模型错位。下载完参考序列之后，建议用samtools faidx建一下fai索引，同时用picard CreateSequenceDictionary生成dict文件，后面Sentieon处理BAM时如果缺少这些伴随文件，个别环节会直接报错。

2.2 FASTQ数据质控与预处理

RNA-seq变异检测的FASTQ预处理，基本原则是“能不做就尽量少做”。很多人习惯上来就做一轮Trimmomatic，把低质量碱基和接头都切掉。但要注意，STAR本身对低质量末端的容忍度不错，并且变异检测需要尽可能保留真实变异位点的支持碱基，过度修剪反而会降低低深度区域的检出率。我的建议是先用FastQC看一眼数据质量，如果平均Q30比例低于80%或者明显存在接头污染，再做针对性的修剪；如果数据质量没问题，直接进STAR比对就行。

这里有一个容易被忽略的细节：如果项目建库时用了UMI分子标签，预处理阶段就不能只做简单修剪，必须先根据建库说明提取UMI，并且在后续比对排序时把UMI信息写入BAM的RX标签。Sentieon在去重环节能读取UMI信息做单分子一致性去重，这对RNA-seq中PCR重复比例偏高的情况非常有效。没有UMI的普通转录组数据，直接跳过这一步，不要画蛇添足。

2.3 用STAR构建基因组索引及参数选择

STAR索引构建是整个流程中最耗时的前置步骤之一，GRCh38加上GENCODE注释大约需要40分钟到一个小时，内存建议给足80GB。命令我一般这样写：

STAR \ --runMode genomeGenerate \ --genomeDir /ref/star_grch38 \ --genomeFastaFiles /ref/GRCh38.fa \ --sjdbGTFfile /ref/gencode.v42.annotation.gtf \ --sjdbOverhang 149 \ --runThreadN 32

这里--sjdbOverhang是一个需要认真讲的参数。它表示构建剪接位点数据库时，考虑的内含子两侧外显子延伸长度，官方推荐值是读长减1，也就是双端150bp测序就填149。这个值决定了STAR在比对上能否有效支持跨越剪接位点的reads：如果填得太小，长reads跨剪接位点时比对信息不足，后续变异检测在剪接点附近容易出问题。如果填得明显大于读长，又会在运行时多耗费内存且索引构建时间变长。所以“读长减1”是唯一的正确选择。

索引构建完会出现一个Genome目录，里面有SA、SAindex等文件。注意这些文件不能和别的参考基因组版本混用，换参考版本后必须重新构建索引，这是STAR使用里最基础也最容易踩的环境问题。

3. 从FASTQ到BAM：比对与数据清洗实践

3.1 STAR比对完整命令与参数解析

索引准备好之后，进入正式比对阶段。STAR比对命令我通常写成这样：

STAR \ --runMode alignReads \ --genomeDir /ref/star_grch38 \ --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM Unsorted \ --outSAMattributes NH HI AS NM MD \ --outSAMunmapped Within \ --outFileNamePrefix sample_ \ --runThreadN 32

几个参数我逐个说下选择理由。--outSAMtype BAM Unsorted表示直接输出未排序的BAM，排序和加RG放到后面统一用Sentieon处理，这样能减少一次中间文件的写入。--outSAMattributes里我保留了NH、HI、AS、NM、MD等标签，这些在变异检测阶段会被用于判断reads的多重比对情况、比对质量、错配数量和碱基质量校准，缺了NM或MD，后续变异检测的过滤会少掉几个重要维度。

--outSAMunmapped Within是把未比对的reads保留在BAM文件中，而不是单独输出一个Unmapped文件。有些场景下，比如后续想用这些unmapped reads做融合基因或病原体筛查，保留下来会非常方便。代价是BAM体积会略微增加，我通常保留这个选项。

比对完成后，STAR目录下会出现sample_Aligned.out.bam和sample_Log.final.out。后者记录了总体比对率、唯一比对reads比例、多位点比对reads比例等关键指标。RNA-seq项目里唯一比对率如果低于75%，就要警惕数据污染、参考基因组不匹配或者建库问题，建议先排查再继续后续流程。

3.2 Read Group添加与排序转换的细节

STAR输出的BAM默认没有Read Group信息，而Sentieon和GATK在处理多样本合并、重复标记和变异检测时，都严重依赖RG信息来区分样本来源。所以拿到比对BAM后，第一个操作就是加RG、排序并转成坐标排序BAM。这一步可以用Sentieon一条命令完成：

sentieon util sort \ -i sample_Aligned.out.bam \ -o sample.sorted.bam \ -t 32 \ --sam2bam \ --RGID sample \ --RGSM sample \ --RGLB lib1 \ --RGPL ILLUMINA

注意--RGSM的值会作为后续VCF中每个样本的标识，建议直接使用样本编号，不要包含特殊符号。--RGID一般和样本名保持一致即可，如果同一份FASTQ被拆分过多次，可以给不同批次不同的RGID，但RGSM保持一致，这样Sentieon去重时能正确处理。

排序这一步会生成sample.sorted.bam和对应的sample.sorted.bam.bai，后续所有操作都基于这个文件。我在实际项目中通常会在这一步后跑一次samtools flagstat做快检，确认比对率、配对率都正常，再继续向下。

3.3 用Sentieon做去重：LocusCollector与Dedup

RNA-seq去重这件事，实际执行时和DNA-seq略有区别。DNA-seq里我们认定同一位点、同一条链上的reads大概率来自PCR扩增，直接标记重复；RNA-seq里同一基因的高表达区域确实会产生大量天然相同的转录本片段，这类reads严格意义上并不能算“技术重复”。但在变异检测实践中，如果不做任何去重，高表达基因区域的覆盖深度会被严重高估，导致变异检测在这些区域对假阳性的判断失效。Sentieon的Dedup模块做的就是标记并移除这种冗余reads。

Sentieon的去重分两步。第一步是收集位点重复信息：

sentieon driver \ -t 32 \ -i sample.sorted.bam \ --algo LocusCollector \ --fun score_info \ sample.score

第二步根据分数表标记并输出去重后的BAM：

sentieon driver \ -t 32 \ -i sample.sorted.bam \ --algo Dedup \ --score_info sample.score \ --dedup sample.dedup.bam \ --metrics sample.dedup_metrics.txt

这里--metrics输出的txt文件里记录了去重前后的reads条数、重复率和估计的文库复杂度。RNA-seq样本的重复率通常在20%到50%之间，如果你看到重复率超过70%，很可能是建库PCR循环数过多，这种情况下即使去重，后续变异检测的可靠性也会受影响，最好回到建库源头去排查。

去重后的BAM文件，可以顺手再用samtools index建一次索引，然后就可以进入变异检测阶段了。

4. 变异检测：用Sentieon Haplotyper拿VCF

4.1 Sentieon Haplotyper调用方式与关键参数

Sentieon的Haplotyper模块在算法上对标GATK的HaplotypeCaller，但执行速度和资源消耗都优化了不少。这一步是整条流程里最吃CPU的阶段，但也是Sentieon优势最明显的地方。我的常见命令写法如下：

sentieon driver \ -t 40 \ -i sample.dedup.bam \ -r /ref/GRCh38.fa \ --algo Haplotyper \ --annotation ClippingRankSumTest \ --annotation DepthPerAlleleBySample \ --annotation DepthPerSample \ --annotation FisherStrand \ --annotation MappingQualityRankSumTest \ --annotation MappingQualityZero \ --annotation QualByDepth \ --annotation RMSMappingQuality \ --annotation StrandOddsRatio \ sample.vcf.gz

这些--annotation参数对应变异位点的多个质控字段，后续过滤都要用到，不能省。比如QualByDepth就是QD，衡量变异位点质量值除以深度的比值；FisherStrand表示链偏倚；MappingQualityRankSumTest是比对质量秩和检验。这些字段在肿瘤样本的变异筛选，或者转录组高假阳性场景下特别重要。

--emit_mode参数这里我刻意没有加。如果要做多样本联合分析，建议输出gVCF模式，也就是加上--emit_mode gvcf，后续再通过GenomicsDBImport或Sentieon的GVCF联合工具做群体call。如果只是单样本出结果，直接输出VCF即可，避免生成中间文件占用空间。

4.2 RNA-seq变异结果的过滤指标选择

RNA-seq相比DNA-seq最大的过滤难点，是覆盖度不均匀带来的一系列连锁反应。DNA-seq里常见的过滤阈值比如“DP < 10过滤”、“GQ < 20过滤”，在RNA-seq中不能机械套用。低表达基因的真实位点可能DP只有4到6，直接把DP过滤阈值拉高就会把好东西丢掉。

我自己一般会在拿到VCF后，先用bcftools view做一轮基础筛选：

bcftools view \ -i 'QUAL > 30 && FMT/DP >= 5 && FMT/AF >= 0.2 && FMT/AD[1] >= 2' \ sample.vcf.gz \ -o sample.filtered.vcf.gz

这个命令的意思是：位点QUAL大于30，样本深度至少5，变异等位基因频率不低于0.2，支持变异的reads至少2条。这几个参数比DNA-seq的阈值都放宽了不少，原因很简单：转录组低表达基因的覆盖度就是达不到DNA测序的标准，硬卡深度会漏掉真实突变。

还要补充一个针对剪接位点侧翼区域的决策。跨剪接位点的reads比对时，外显子边界附近往往会被软裁剪，这会导致变异检测在剪接位点上游1到3bp处产生大量低质量候选位点。对这类位点，我会额外检查ReadPosRankSum和ClippingRankSumTest这两个指标，如果偏差过大，一般建议过滤掉，除非有很强的生物学证据支持该位点为真实变异。

4.3 变异注释：snpEff/VEP/ANNOVAR的实践思路

变异注释这一步的目的是把VCF中的基因组坐标信息转化为基因和转录本层面的功能后果。常用工具包括ANNOVAR、snpEff、VEP。这里我以snpEff为例做一个快速实践。

snpEff使用前需要先下载对应参考版本的数据库：

java -Xmx16g -jar snpEff.jar download GRCh38.105

然后对过滤后的VCF执行注释：

java -Xmx16g -jar snpEff.jar \ -v GRCh38.105 \ sample.filtered.vcf.gz \ > sample.ann.vcf.gz

注释后的VCF里可以看到每个变异位点在基因上的位置、影响类型（如missense、nonsense、synonymous）、涉及的转录本ID、氨基酸变化等信息。如果项目规模不大，或者你不想额外安装很多软件，也可以直接用VEP在线或离线版本，输出HTML和VCF两种结果。

RNA-seq变异检测有一点和DNA-seq不一样：注释时经常要小心“伪基因”和“同源基因高相似区”。转录组reads往往来自基因家族中高度相似的成员，这些区域比对时天然存在多映射，变异检测结果容易在多个基因的相同位置同时call出相同变异。注释出来后，如果发现同一变异位点映射到多个基因，基本可以判定是多映射reads导致的假阳性，需要在后续报告中单独标记，不能当成多个基因的真实突变分别解读。

5. 运行实录与常用问题排查

5.1 一份典型的运行时间和资源占用记录

我用一份真实项目数据做个参考：人类转录组，60M对reads，双端150bp，服务器配置为Intel Xeon 48核、512GB内存。整个流程各环节耗时如下表所示：

可以看到，变异检测是绝对的耗时大头。如果用GATK HaplotypeCaller跑同样的数据，这个步骤通常要8到12小时甚至更久，Sentieon在这个环节能省下将近三分之二的时间，而且内存峰值也没有失控，这是它在转录组这种大BAM场景下最实用的优势。

5.2 比对率低和变异假阳性的排查路径

比对率低是RNA-seq变异检测中最先暴露的问题。如果STAR的Log.final.out里uniquely mapped占比低于75%，我一般按以下顺序排查。第一步检查参考基因组是否与样本物种一致，尤其是人的样本却误用了小鼠参考；第二步查看数据是否混入接头或rRNA序列，这在质量报告中往往有明显表现；第三步检查GTF注释版本和参考基因组是否配套，不配套时STAR在剪接位点数据库构建阶段就会产生偏差，直接影响跨外显子reads的比对率。

变异假阳性偏高则通常发生在以下几种场景：比对质量差的区域、序列同源性高的基因家族区域、剪接位点侧翼区、以及RNA编辑事件密集区域。前面提到过，ADAR介导的A到I编辑在VCF中看起来就是A到G的SNV，这类位点往往质量值不低，但生物学上并不是来源于基因组突变。如果项目只需关注真正的胚系或体细胞突变，建议在最终报告阶段对这类候选变异单独增加一个“RNA编辑可能性”备注。另一个实用技巧是，如果同一项目里有DNA-seq配对数据，直接用DNA-seq的变异结果去过滤RNA-seq结果，只保留两个平台都检出的位点，能被这个策略筛掉的RNA-seq候选位点数量通常非常可观。

5.3 RNA-seq变异检测常见问题速查表

我把实际踩过的坑整理成下面这个速查表，每次做新项目时我都会复查一遍。

实际项目中，最常被忽略的就是样本ID一致性。很多人比完对、跑完变异检测，最后在注释或者临床解读阶段才发现，VCF中的样本名和原始编号对不上，只能重新跑一遍排序和变异检测，浪费时间不说还容易引入新的错误。所以所有流程一开始，我就建议把样本命名规则定死，比如用“项目编号-样本编号”的格式统一写入--RGSM，后续所有步骤都沿用，不要随便换。

5.4 提升RNA-seq变异检出的一个补充思路

如果你的项目重点不是单纯拿到一份VCF，而是想提高真正有意义变异的检出率，可以再考虑两种策略。第一种是多样性本联合call，也就是多个RNA-seq样本合并分析，用统一的gVCF模式生成各样本gVCF后再联合变异检测。转录组覆盖度不均匀的弱点，在多样本联合分析中会被部分稀释，同一个位点只要有几个样本支持，可信度就会明显提升。第二种是在变异检测前做BAM的拆分再合并，比如把STAR输出的唯一比对reads单独抽出来，去掉多比对reads再进Sentieon Haplotyper，这一步在很多公开教程里提得不多，但实际操作里能有效降低同源基因区域的假阳性。

当然，这些策略需要根据项目本身的数据规模和研究目的来权衡。如果是单样本的医学报告，不希望漏掉罕见突变，单纯追求高召回率，那就不必过度过滤，更多依赖注释和人工复核来解决假阳性问题；如果是多队列筛选候选易感基因，过滤策略则可以更严格一些，确保最终候选列表的精确度。

我在实际使用Sentieon跑RNA-seq变异检测时，最深的体会是：这款工具确实把变异识别和去重的效率拉到了一个新高度，但这并不代表前面的比对和后续的过滤可以被轻视。STAR的比对参数、RG信息的规范、去重策略的选择，任何一个环节偷懒，都会在最终VCF里以成百上千条噪声位点的方式还回来。最后再分享一个小技巧：STAR比对结束后保留sample_Chimeric.out.junction文件和未比对reads，这些数据在后续做融合基因分析或调试异常比对时非常有用，我当时临时想回头补一个融合检测分析，靠的就是这套保留文件，省去了重跑比对的麻烦。