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及之后的releaseSTAR版本选用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与DedupRNA-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第二步根据分数表标记并输出去重后的BAMsentieon 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拿VCF4.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内存。整个流程各环节耗时如下表所示环节工具线程数运行时间内存峰值FASTQ质控FastQC4约15分钟约4GBSTAR索引构建STAR32约50分钟约100GBSTAR比对STAR32约60分钟约64GB排序加RGSentieon util sort32约25分钟约40GB去重Sentieon Dedup32约20分钟约50GB变异检测Sentieon Haplotyper40约3小时约80GBVCF过滤与注释bcftools snpEff4约30分钟约12GB可以看到变异检测是绝对的耗时大头。如果用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变异检测常见问题速查表我把实际踩过的坑整理成下面这个速查表每次做新项目时我都会复查一遍。现象可能原因处理方式STAR比对率很低参考基因组版本或物种错误检查fasta来源和GTF一致性BAM文件没有RG标签排序时漏加--RGID用sentieon util sort重新添加VCF里没有样本名--RGSM未正确指定确保排序加RG时使用规范样本ID高表达基因假阳性多未去重或重复率过高检查Dedup步骤的metrics文件剪接位点附近变异密集比对软裁剪的影响过滤ReadPosRankSum异常位点大量A到G变异可能为RNA编辑使用RNA编辑数据库辅助判断内存不足导致driver崩溃STAR或Sentieon线程过多降低线程数并限制单样本内存上限VCF后续工具报错BAM和VCF的样本ID不一致统一RGSM、VCF样本名实际项目中最常被忽略的就是样本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这些数据在后续做融合基因分析或调试异常比对时非常有用我当时临时想回头补一个融合检测分析靠的就是这套保留文件省去了重跑比对的麻烦。