1. 内容整体设计与思路拆解1.1 RNA-seq 变异检测为什么绕不开 Sentieon做生信的人应该都有这种体会RNA-seq 数据里挖变异比 DNA 数据麻烦得多。DNA 变异检测的流程早就成熟了GATK 的 Best Practices 一套组合拳打下来大家照着做就行。可 RNA-seq 不一样splice junction 的存在让比对环节就变得很棘手reads 跨过外显子接头时如果不做特殊处理比对软件很容易把它们当成“垃圾”序列丢掉或者硬生生比对到错误的位置上。比对错了后面的变异检测就是空中楼阁。Sentieon 这一套 RNA-seq 变异检测流程解决的正是这个痛点。它提供的 RNA-seq 短序列比对方案在比对阶段就引入了对 splice junction 的感知能力再配合它那套出了名的、与 GATK 数学原理高度一致的变异检测算法能把 RNA-seq 数据的变异检测精度拉到一个非常可用的水平。很多人一提到 Sentieon 就只想到速度快其实它的 RNA-seq 流程在设计上更值得琢磨速度只是表象真正厉害的是它对 RNA-seq 数据特有噪声的处理方式。这条流程适用的场景很广既有肿瘤 RNA-seq 样本的体细胞变异筛查也有遗传病研究中基于转录组的罕见变异发现还包括一些无法获取 DNA 样本、只能做转录组测序的回顾性研究。只要你手里的数据是链特异性或者非链特异性的 RNA-seq 短读长数据想从中挖出 SNV 和 InDel这套流程就是一个相当稳妥的参考基准。可能有人会问RNA-seq 本身就不是做变异检测的首选技术为什么还要费劲去跑这套流程RNA-seq 的优势在于它能同时反映基因表达和变异信息对于某些特殊类型样本比如 FFPE 保存多年的组织DNA 可能已经降解得不成样子但 RNA 反而还能提取出可用的信息。另外RNA-seq 变异检测可以帮助验证 DNA 水平发现的剪接区域变异是否真的影响了转录本这种“转录组验证”的视角是纯 DNA 测序给不了的。1.2 项目选型时我为什么没用 GATK 而选了 Sentieon先说实话GATK 的 RNA-seq 变异检测流程本身没毛病它能用但用起来确实有点憋屈。RNA-seq 数据量通常比 WES 大不少动辄上百 G 的 fastqGATK 那套流程跑下来耗时经常以“天”为单位计算。如果手头有几十个样本光等结果就能等到怀疑人生。而且 GATK 官方对 RNA-seq 变异检测的支持力度一直不如 DNA 流程很多参数需要自己反复调文档里语焉不详的地方也不少。Sentieon 打动我的有两点。第一点自然是速度这没得洗。Sentieon 把 HaplotypeCaller 的核心算法用底层优化重新实现了一遍在相同硬件条件下加速比通常在 5 到 15 倍之间。我实测过一批 100X 左右的 RNA-seq 肿瘤样本GATK 跑一个样本的 HaplotypeCaller 大概要 8 小时Sentieon 同等配置下 40 分钟出头就跑完了。这还只是单样本的差距大规模队列跑下来省出来的时间够做很多事了。第二点是它在 RNA-seq 流程上做的针对性优化。Sentieon 不需要像 GATK 那样强制你用 STAR 比对后再做一系列繁琐的预处理它对 Input 比对结果的要求更灵活而且它对 splice junction 信息的使用方式也更直接。它有一套推荐参数是针对 RNA-seq 的比如在比对时推荐使用特定工具生成带有NMtag 的 BAM变异检测时能更准确地利用转录组特征。当然选型还要考虑成本。Sentieon 是商业软件授权费不便宜GATK 免费且开源。如果是个人学习、发小文章GATK 完全够用但如果项目有工期压力、样本量又大把时间成本折算进去Sentieon 其实并不算贵。我自己在项目里倾向于混合使用小样本探索用 GATK批量正式跑用 Sentieon两边结果做个交叉验证心里更有底。1.3 这套流程能解决的核心问题与适合人群这套 RNA-seq 变异检测全流程能解决的最核心问题是在 RNA-seq 数据中可靠地检测出 SNV 和短 InDel尤其是在涉及剪接位点附近的变异时减少比对错误带来的假阳性。RNA-seq 数据的变异检测假阳性率天然高于 DNA 数据这是转录组本身的结构决定的可变剪接、RNA 编辑、序列比对的多重映射都会在突变位点上制造噪声。Sentieon 的流程通过比对阶段的 splice-aware 处理和变异检测阶段的局部重组装把假阳性控制在一个可接受的范围内。适合参考这套流程的人我总结下来有三类。第一类是肿瘤方向的生信工程师需要从 RNA-seq 数据中寻找潜在的驱动突变特别是那些在 DNA 水平上被漏掉的、只在转录组水平表达的融合或点突变。第二类是遗传病诊断方向的研究人员当 DNA 测序没有找到明确致病位点时补充分析 RNA-seq 数据看是否有异常剪接或等位基因特异性表达的问题。第三类是刚刚接触 RNA-seq 分析、想建立一套标准流程的入门者这套流程的逻辑清晰、参数明确照着做一遍就能理解 RNA-seq 变异检测的各个环节。我特别想提醒的一点是这套流程虽然是目前 RNA-seq 变异检测的标杆但它并不能替代 DNA 测序的变异检测。RNA-seq 只能检测到在转录组中有表达的基因表达量低的基因或者被 NMD 途径降解的突变转录本在 RNA 水平上根本看不到。所以 RNA-seq 变异检测的结果更准确地说是“转录组层面的表达变异”阴性结果并不能排除 DNA 水平存在变异。2. 核心细节解析与实操要点2.1 从 fastq 到 BAM比对环节的“知易行难”RNA-seq 流程的第一步也是最关键的一步就是把原始测序数据比对到参考基因组上。这一步做得不好后面再怎么调参都救不回来。DNA 测序比对时只需要考虑 reads 在基因组线性坐标上找位置就行但 RNA-seq 的 reads 经常跨过内含子一条 read 的 5 端比对外显子 A3 端比对到外显子 B中间隔着几千甚至几万 bp 的内含子。如果比对软件不认识剪接位点这种 read 就会被打上“无法比对”的标签丢掉或者错误地比对到其它有相似序列的区域。所以我强烈建议RNA-seq 的比对环节一定要选择支持 splice junction 的比对软件。Sentieon 官方推荐的比对策略中STAR 是首选。STAR 的比对速度极快它的 seed-search 和 stitching 策略非常契合 RNA-seq 数据的特点而且它输出的 BAM 里包含NH、HI等 tag这些 tag 对后续 Sentieon 变异检测的去重环节非常有用。比对时有一个参数非常重要--outSAMattributes。如果用 STAR 比对建议在参数里加上NH HI AS NM MD这样输出的 BAM 文件里就有了比对质量、错配信息等关键属性。Sentieon 的LocusCollector去重步骤会依赖这些 tag如果缺失可能会影响去重的准确性。我自己踩过这个坑一开始没注意 STAR 的输出属性配置结果下游去重时总是报一些莫名其妙的信息加上这几个属性之后立刻清爽了。还要提醒的是RNA-seq 比对时普遍存在 multi-mapping 的问题。转录组里高度同源的基因家族、假基因、重复序列都会导致一条 read 能同时比对到基因组的多个位置。对于变异检测来说multi-mapping 的 reads 是很大的噪声来源。我推荐在比对后过滤掉某些比对质量过低或标记为 multi-mapping 的 reads但过滤尺度要把握好如果过滤太狠一些真正来自转录本高表达区域的 reads 也会被误伤。2.2 数据预处理里容易被忽略的细节比对完成后接下来的数据预处理环节看似简单但细节决定了最终结果的质量。这个环节主要包括排序、去重复、碱基质量校正三步。传统流程里这三步要分别调用不同的工具Sentieon 的优势在于它用一条命令把多个工具串联起来了省去了频繁读写中间文件的 I/O 开销这也是它速度快的一个重要原因。排序这一步没什么好说的按坐标排序是必须的关键是SortSam之后要建立 index否则很多下游工具会拒绝工作。Sentieon 的做法比较聪明它通过--sort_order coordinate参数直接在流程中完成排序和建索引不用中间环节额外处理。去重复这一步要特别讲一下它可能是整个预处理里对变异检测结果影响最大的一步。RNA-seq 测序时由于反转录和 PCR 扩增的偏好性某些高表达基因的片段会被大量重复测到这些重复 reads 如果不去掉会在变异检测时被当作独立的证据堆高假阳性。但 RNA-seq 的去重和 DNA 不太一样同一个基因的不同转录本异构体可能共享相同的外显子区域导致它们的 reads 在基因组坐标上看起来完全一样但实际上来自不同的 RNA 分子。如果按照 DNA 测序的标准去重这些 reads 会被错误地当成 PCR 重复去掉造成表达量和变异频率的偏差。Sentieon 的去重算法考虑到了这一点通过LocusCollector收集 reads 的坐标和标签信息再利用Dedup进行去重。在实际使用中我会把--rmdup相关参数逐个检查一遍确保去重策略符合自己项目的需求。如果是做表达量分析有人会选择不去重但做变异检测我还是建议去重否则重复序列区域的变异检测结果基本没法看。碱基质量校正BQSR是另一个容易被忽视但很重要的环节。测序仪给出的碱基质量值往往跟真实的错误率有偏差比如某些 motif 后面容易出现特定类型的错误。BQSR 就是通过对已知变异位点dbSNP 等的实际观察建立错误率模型重新校准碱基质量值。Sentieon 的BaseQualityScoreRecalibration和ApplyBQSR对应 GATK 的 BQSR 流程但速度更快。RNA-seq 数据里有些位点的碱基组成有偏好性BQSR 的纠偏效果在 RNA-seq 上尤其明显。2.3 变异检测阶段的参数选择与调优思路变异检测是整个流程的核心环节。Sentieon 的Haplotyper算法和 GATK 的 HaplotypeCaller 思路一致在目标区域发现变异信号后提取相关的 reads对这些 reads 做局部 de Bruijn 图的组装构造出可能的单倍型再把 reads 比对回单倍型上最后用隐马尔可夫模型计算每个位点的基因型概率。这种局部重装配的策略在处理 RNA-seq 数据时非常关键因为它能处理 splice junction 区域的复杂变异模式。RNA-seq 变异检测时有几个参数需要特别注意。第一个是--min_base_qual默认值是 10但 RNA-seq 数据在剪接位点附近经常会由比对软件引入一些质量偏低的碱基如果这个阈值太高可能会漏掉真实变异。我一般会把这个参数设为 15 左右然后结合具体数据表现再调整。第二个是--max_reads_per_alignment_start控制每个位点最多参与计算的 reads 数量RNA-seq 数据在高表达区域会出现大量 reads 堆叠这个参数默认值往往不够用导致高覆盖区域的变异检测效果变差。我会根据数据的平均深度适当调高这个值。还有一个非常实际的问题RNA-seq 数据里覆盖度在基因之间的差异非常大。高表达基因的覆盖深度可能达到几千甚至上万 X低表达基因可能只有个位数 X。这种极端的覆盖度不均一性对变异检测算法的碱基质量模型是一个巨大的考验。Sentieon 在默认参数下对高覆盖区域做了一定的降采样处理避免计算资源被单一区域耗尽但如果你的研究重点恰好是高表达基因的变异可能需要手动调整采样参数。在变异检测结果的过滤上我倾向于不过度依赖单一的过滤标准。RNA-seq 的变异检测结果需要用多种信息交叉验证比如比对质量、覆盖深度、链特异性、是否位于外显子-内含子边界等。Sentieon 的介绍里也强调了这一点它生成的 VCF 文件里包含了丰富的注释字段充分利用这些字段做过滤比死记硬背一套“万能过滤参数”要靠谱得多。3. 实操过程与核心环节实现3.1 环境准备与测试数据获取先把环境跑通这是所有后续工作的基础。Sentieon 提供了多种安装方式最常见的是下载官方编译好的安装包。安装本身不难但要注意版本匹配问题Sentieon 版本更新很快不同版本对 GATK 兼容模式的支持略有差异建议直接用最新的稳定版本。安装好之后测试数据可以从 Sentieon 官网的教程目录下载里面有配套的 RNA-seq 示例数据和参考基因组的对应区域。如果没有现成的测试数据也可以用公开数据库里的 RNA-seq 数据比如 ENCODE 项目或者 GTEx 项目里的人类 RNA-seq 数据。第一次跑流程强烈建议先用小规模的测试数据把全流程跑通确认每一步的输出都符合预期再上真实的大规模数据。这个习惯帮我避免过很多低级错误比如路径写错、参考基因组版本不对之类的问题。参考基因组的选择也是一个需要提前确认的点。RNA-seq 变异检测推荐使用完整的参考基因组包括所有染色体 contig这样比对软件才能正确处理来自线粒体等区域的 reads。另外参考基因组的版本必须和比对软件及 Sentieon 的配置一致混合使用不同版本会导致染色体命名对不上、变异位点坐标错位等一系列问题。3.2 用 STAR 完成 splice-aware 比对STAR 比对 RNA-seq 数据之前需要先建立基因组索引。这一步不是特别耗时但需要注意构建索引时加上--sjdbGTFfile参数指定基因注释文件GTF/GFF3。这样 STAR 在比对时就能参考已知的剪接位点显著提高 splice junction reads 的比对准确性。如果省略这一步STAR 就只能依赖从头发现的剪接位点对已知剪接位点的支持会弱很多。STAR 建索引的命令大致是这样STAR --runMode genomeGenerate \ --genomeDir /path/to/star_index \ --genomeFastaFiles /path/to/reference.fa \ --sjdbGTFfile /path/to/annotation.gtf \ --sjdbOverhang 149 \ --runThreadN 20sjdbOverhang这个参数需要根据测序读长来设置一般推荐设为读长-1比如 150bp 读长就设为 149。这个值是给 STAR 在已知剪接位点两侧延伸的序列长度设得太小可能覆盖不了完整的剪接位点特征太大则会浪费内存。建索引时还需要注意内存要够用人类基因组加 GTF 注释的索引构建建议至少准备 30G 内存。比对命令的写法如下STAR --genomeDir /path/to/star_index \ --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM Unsorted \ --outSAMattributes NH HI AS NM MD \ --outFilterMultimapNmax 20 \ --outFilterMismatchNmax 10 \ --runThreadN 20这里面有几个参数值得细说。--outFilterMultimapNmax 20表示允许 read 最多比对到 20 个位置超过就会被过滤掉。转录组里 repeats 区域产生的 multi-mapping reads 非常多这个阈值设得太高会让大量多映射 reads 进入下游分析反而增加噪声设得太低又可能丢掉来自同源基因家族的有用 reads我常用的经验值是 10 到 20 之间根据项目需求调整。--outSAMattributes列出的这些属性很重要特别是NM错配数和MD错配字符串后续 Sentieon 的质量校正环节会用到。STAR 默认输出 Unsorted BAM这个后面要交给 Sentieon 做排序。这一步输出的 BAM 文件非常大建议在比对时开启--outBAMcompression 6之类的压缩选项可以节省不少磁盘空间。3.3 Sentieon 流程的完整命令与参数解读比对完成后的 Sentieon 流程可以用一条 slog 串联起来。这里给出一个比较完整的 RNA-seq 变异检测命令模板同时写了注释方便对照sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo LocusCollector \ --fun score_info \ sample.sorted.bam \ sample.score sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo Dedup \ --score_info sample.score \ sample.sorted.bam \ sample.dedup.bam sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo BaseQualityScoreRecalibration \ sample.dedup.bam \ sample.recal.table sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo ApplyBQSR \ -i sample.recal.table \ sample.dedup.bam \ sample.recaled.bam sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo Haplotyper \ --genotype_model multinomial \ sample.recaled.bam \ sample.vcf拆开来看每一步的作用。LocusCollector是在统计每个基因组位点上的 read 起始位置信息为去重做准备。它的输出是一个文本格式的 score 文件不是 BAM所以用重定向而不是-o。这一步的核心是标记 PCR 重复 reads 的候选位置它不需要真实输出 BAM所以速度很快。Dedup读取上一步的 score 文件真正执行去重操作。输出文件是去除重复后的 BAM。这里可以加上--rmdup参数表示在去除重复的同时也物理性地删除标记为重复的 reads否则只是给 reads 加上duplicateflag实际 reads 还保留在 BAM 里。从节省空间的角度我一般会开启--rmdup不过如果是后续还要做表达量分析保留重复 reads 也有道理具体看需求。BaseQualityScoreRecalibration生成碱基质量校正表它不会直接修改 BAM而是把校正参数输出到recal.table里。这一步需要提供已知变异位点数据库可以通过--known_sites参数传入 dbSNP 的 VCF 文件。如果没有已知位点数据库也可以不传校正效果会打折但流程仍然能跑通。ApplyBQSR是真正应用校正表的步骤把重校准后的质量值写回 BAM 文件。从这一步开始BAM 就已经是可用的“干净数据”了。Haplotyper是变异检测的核心命令。这里有一个小细节--genotype_model multinomial是适用于 RNA-seq 数据的模型选择它会用多项分布模型来评估基因型概率这个模型在杂合位点的检测上表现更好。如果这里不加这个参数默认模型也能跑但对 RNA-seq 数据的适配性差一些建议显式指定。3.4 从 VCF 到可用结果过滤、注释与解读Haplotyper输出的原始 VCF 文件直接拿去用是不行的。原始 VCF 里包含了很多低质量的候选位点需要一个过滤步骤把这些噪声去掉。Sentieon 官方推荐使用VariantFilter配合一组 RNA-seq 特异的过滤条件但实际使用中我更喜欢把 VCF 导出来自己在 Python 或者 R 里按照项目的具体需求做针对性过滤灵活性更高。常用的过滤思路包括过滤掉深度小于 10 的位点、过滤掉 QUAL 值小于 30 的位点、过滤掉链偏差显著的位点。RNA-seq 数据的链偏差问题尤其严重因为转录组测序本身就有链特异性如果两条链的覆盖度差异过大变异检测结果很容易出现假阳性。过滤完的 VCF 还需要做功能注释最常用的是SnpEff或VEP把变异位点映射到基因和转录本上标出外显子区、剪接位点区、UTR 区等位置信息。注释这一步特别重要因为 RNA-seq 检测到的变异最终要跟基因功能联系起来才有意义。我用过一个比较省力的组合Sentieon Haplotyper出原始 VCF接着用SnpEff做注释再把注释结果丢进ANNOVAR做进一步的数据库交叉比对。这个流程跑熟之后从 raw fastq 到最终可解释的变异列表一个样本大约需要 2 到 3 小时其中大头时间是 STAR 比对和 Haplotyper 这两步。4. 常见问题与排查技巧实录4.1 STAR 比对阶段的高频报错与解决方案STAR 比对报错是 RNA-seq 流程里最常碰到的问题。我总结了几个高频问题大家遇到可以直接对照排查。第一个是内存不足。构建人类基因组索引时默认参数下大概需要 30G 内存如果服务器内存不够会出现FATAL ERROR: not enough memory之类的提示。解决办法是降低--sjdbOverhang的值或者用--genomeSAindexNbases参数减少索引的规模。另外把--runThreadN调小一度内存压力也会降一点虽然速度会慢一些。第二个是zcat命令找不到。在读入压缩 fastq 时--readFilesCommand zcat依赖系统里有 gzip 解压工具有些精简版系统里没有装。解决办法是把zcat换成gzip -dc或者干脆先解压 fastq 再输入前提是磁盘空间足够。第三个是 GTF 文件和参考基因组版本不匹配。这个问题比较隐蔽表现在比对结果上就是剪接位点识别率异常低、比对率也偏低。解决办法是建立索引前严格确认 GTF 和 fasta 来自同一个版本比如都来自 Ensembl release 110或者都来自 GENCODE v44。还有一个容易被忽视的问题是STAR 比对率正常但后续 Sentieon 流程跑出来的变异数少得离谱。这种情况多半是 BAM 文件的NMtag 缺失导致的。如果比对时没有加上--outSAMattributes NM MDSentieon 在质量校正阶段就可能出错。解决办法是回到 STAR 比对那一步把参数补全重新比对。如果已经生成了大量 BAM没法重新比对也可以用samtools fillmd之类的工具把NM和MDtag 补上但效果还是不如一开始就加参数来得干净。4.2 Sentieon 运行时的常见错误与处理Sentieon 运行时的报错多数是输入文件或参数的问题。最常见的是参考基因组字典文件缺失。在跑Haplotyper之前实际上不强制要求提供 dict 文件但许多辅助功能比如按染色体并行、区间切分会依赖它。如果报错提示找不到.dict文件用picard CreateSequenceDictionary或者samtools dict生成一个放在参考基因组同目录下即可。另一种常见问题是 BAM 文件的 header 不完整。当 STAR 输出的 BAM 用samtools sort处理过而没有保留原始 header 时Sentieon 可能因为读不到正确的SQ信息而报错。解决办法是不要手工对 BAM 做多余处理直接把 STAR 的 Unsorted BAM 交给 Sentieon由它内部完成排序和索引这样最稳妥。如果是跑队列任务还容易遇到内存分配不足的问题。Sentieon 的driver -t 20表示用 20 个线程但线程数和内存是两回事如果每个线程分配的内存太少会出现Cannot allocate memory之类的报错。解决办法是在提交任务时给足内存常规人类 RNA-seq 样本建议至少 16G高深度样本建议 32G 以上。4.3 RNA-seq 变异结果的假阳性排查思路跑完流程拿到 VCF才是真正考验功力的开始。RNA-seq 变异检测的假阳性问题比 DNA 数据严重得多所以一定要有一套系统的排查思路。第一个要排查的是链偏差。打开 BAM 文件检查变异位点的 supporting reads 是否集中在某一条链上。如果是大概率是链特异性建库导致的技术偏差不是真实的体细胞突变。判断标准是看 the reads supporting the variant 的XStagSTAR 输出的链方向 tag如果所有突变 reads 都来自同一条链就要高度警惕。第二个要排查的是剪接位点附近的变异。位于外显子-内含子边界附近的变异很容易因为比对错误造成假阳性。特别是内含子侧翼 1~2bp 的位置那里是剪接体的识别核心区域如果检测到变异必须回到 IGV 里人工检查比对情况。我见过很多次看似“重要”的剪接位点突变实际上是 reads 比对到假剪接位点造成的假象。第三个要排查的是 RNA 编辑事件。RNA-seq 数据里常见的 A-to-I 编辑会把基因组上的 A 检测成 G。如果发现变异集中在已知的 RNA 编辑位点附近特别是 Alu 元件区域要警惕这是 RNA 编辑信号而不是真正的基因组变异。目前有一些 RNA 编辑位点数据库可以做过滤比如 REDIportal但要注意它对不同组织类型的覆盖度并不均衡。最后如果条件允许强烈建议对 RNA-seq 变异结果做 Sanger 测序验证或者用 DNA 测序数据做个交叉验证。RNA-seq 变异检测的结论尤其是那些影响临床决策的重要位点一定要有独立的技术验证才能下结论。这不是对流程的不信任而是分子生物学实验的基本严谨性。4.4 处理多个样本时的内存与运行时间管理单个样本跑通之后真实项目通常面临的是几十甚至上百个样本的批量运行。这时候需要考虑的是如何高效地管理资源而不是盲目地堆线程。先看时间。RNA-seq 样本的比对和变异检测单样本大约占用 2~4 小时的计算时间具体取决于数据量和服务器配置。如果有 50 个样本按照单样本 3 小时算串行跑需要 150 小时这显然是不可接受的。解决办法是并行化同一时刻跑多个样本的流程同时留出足够的 CPU 和内存资源。再看内存。Sentieon 的Haplotyper在高覆盖区域会消耗较多内存如果多个样本同时运行内存峰值会有叠加。我的经验是每同时运行 4 个样本至少准备 64G 内存。如果内存紧张可以通过--interval参数把基因组拆分成多个区间并行处理跑完再合并 VCF这样单样本的峰值内存可以降下来同时整体吞吐量也有提升。还有一个小技巧是监控中间文件大小。STAR 输出的 BAM 动辄几十 G如果一组样本同时跑磁盘空间很容易被打爆。建议在处理前先规划好目录结构给每个样本独立子目录并定时清理中间文件。流程跑完后保留最终的 VCF 和重校准后的 BAM 就够了原始比对产生的超大 BAM 可以先用samtools view -b -F 1024过滤掉重复 reads 再压缩存档能省下不少空间。5. 流程验证与结果质控5.1 用已知变异样本验证流程的准确性跑流程之前最好先用一个已知变异位点的样本做验证。Sentieon 官网和公开数据库都提供了一些带有真实变异注释的 RNA-seq 测试数据比如 NA12878 的 RNA-seq 数据它对应的高置信度变异集可以从基因组学权威资源中获取。验证的方法是把流程跑完得到 VCF 文件然后和已知变异集做比较。比较的指标包括 Recall召回率、Precision精确率和 F1 score。如果某个指标明显偏低说明流程中某些环节有问题需要回头检查。我见过比较多的坑是向量子集比较时参考基因组版本不一致导致同一变异位点的坐标对不上这会造成召回率虚低排查时往往浪费很多时间。实际操作里我是这么做的先用bcftools isec或者RTG Tools vcfeval把检测到的变异与已知变异集比较重点看常见变异位点是否都能检测到以及检测到的变异里有多少是已知集里没有的。vcfeval在处理复杂位点时比bcftools isec更准确因为它会做等位基因级别的比对而不是简单的坐标比较。5.2 覆盖度、比对率和样本一致性检查RNA-seq 样本之间差异很大同一个批次的数据也可能有不同的质量表现。所以在正式做变异检测之前先检查数据的几个核心质控指标比对率STAR 输出日志里会报告唯一比对 reads 的比例一般要求大于 80%如果低于这个值说明样本可能存在降解或者污染。基因体覆盖均匀度检查 reads 在转录本上的覆盖是否均匀如果 3 端严重偏向提示 RNA 质量可能有问题。重复率去重后剩余 reads 的比例高重复率可能说明起始 RNA 量不足或者 PCR 循环过多。链特异性如果建库是链特异性需要确认链方向是否和预期一致。这些指标可以用Picard CollectRnaSeqMetrics和MultiQC汇总检查。Per-sample 都通过质控后再做变异检测就能减少由于样本质量差异引入的假信号。有个情况特别提醒一下肿瘤样本可能会有较大的拷贝数变异或者杂合性缺失区域导致某些染色体区域的 reads 覆盖度显著偏离平均水平。这种区域在变异检测时经常会出现大片的连续变异信号看起来像是突变热点其实就是拷贝数变化导致的技术信号。分析这类样本时最好把拷贝数信息和变异检测结果联合看别被假象带偏了。5.3 变异结果的生物学合理性评估最后的质控关卡是评估变异结果的生物学合理性。即使所有技术指标都正常也不能保证变异结果是真实可用的。先从突变频谱入手。人类基因组中CT 转换通常是最常见的突变类型。如果检测到某个样本的突变频谱明显异常比如 CA 或 AT 颠换占主导有可能是测序错误或者样本处理过程中引入了氧化损伤。这种频谱异常在 FFPE 样本中特别常见是人工假象的典型标志。再看等位基因频率分布。胚系变异的等位基因频率一般集中在 0.5 和 1.0 附近体细胞变异则呈现低频分布。如果大量变异的等位基因频率集中在一个奇怪的档位比如 0.4 左右可能是样本污染也可能是拷贝数异常。这个判断需要结合样本类型做综合考量。最后把变异映射到通路和基因功能上。RNA-seq 变异检测的结果应该跟样本的表型或者疾病背景有逻辑关联。如果发现一个正常组织样本里出现了大量高频的已知驱动基因突变那先别高兴多半是样本标签搞错了或者比对时发生了样本间的交叉污染。这种“做出来结果太完美反而可疑”的场景在实际项目中并不少见。6. 项目总结与实操心得6.1 我跑完这套流程后的几点体会前前后后用 Sentieon 跑了不少 RNA-seq 样本积累了一些体会。第一个体会是RNA-seq 变异检测的快和准本质上是“比对阶段打底、变异检测阶段兜底”。STAR 的 splice-aware 比对决定了基础比对质量而 Sentieon 的 Haplotyper 通过局部重装配把变异检测精度拉高了一个档次。两者缺一不可。如果比对阶段就丢掉了大量剪接位点 reads后面的检测能力再强也无济于事。第二个体会是RNA-seq 变异检测的结果一定要做技术验证。Sentieon 流程给出的 VCF 文件在统计意义上是可信的但它毕竟是从转录组数据反推基因组变异受到表达水平、降解程度、建库偏好等多种因素影响。批量分析时可以用算法过滤掉大部分假阳性但真正的“金标准”位点比如临床报告里要写的位点绝对要用 Sanger 或者靶向重测序验证。这不是流程的缺陷而是 RNA-seq 数据本身的边界。第三个体会是流程的工程化能力。Sentieon 之所以适合大规模队列是因为它可以把比对之外的步骤全部用一条流水线串起来避免了中间文件的反复读写。在一个 100 样本的肿瘤 RNA-seq 队列中我用这套流程跑完发现耗时主要花在 STAR 比对和 Haplotyper 上去重和质量校正的速度几乎感觉不到存在。这种工程效率的提升是实际项目里能实打实感受到的。6.2 几个值得提前规划的点如果是从零开始搭建这套流程有几点值得提前规划好。第一磁盘空间一定要提前摸排。STAR 输出的大 BAM、Sentieon 去重后的BAM、重校准后的 BAM每个都是原始 fastq 的 2~3 倍大小三个版本叠加起来磁盘占用非常夸张。我建议流程跑完后立刻清理不需要的中间 BAM只保留最下游的重校准 BAM 和 VCF 文件。如果项目周期长中间文件可以压缩存档到冷存储。第二参考基因组和注释文件的版本管理一定要严格。我在多个项目里都吃过版本不一致的亏最典型的例子是参考基因组用了 GRCh37而 dbSNP 用了 GRCh38 的版本导致已知位点数据库在比对时大范围错位质量校正形同虚设。建议在项目启动时就把参考基因组、GTF 注释、已知变异数据库的版本固定下来用统一的配置记录在项目文档里。第三流程的可重复性要靠脚本化保证。不要手工一句一句地敲命令把整个流程写成 shell script 或者用流程管理工具Snakemake、Nextflow封装起来参数、路径、版本全部集中在一个配置文件里。这样不仅自己能复用团队协作时别人也能轻松复现你的结果。6.3 后续功能扩展的思路这套流程跑通之后可以沿着几个方向做扩展。最直接的扩展是融合基因检测。RNA-seq 数据除了点突变还能用来找融合基因。可以在 Sentieon 变异检测的基础上额外跑一遍STAR-Fusion或者Arriba把融合基因的检测结果和点突变结果合并分析为肿瘤研究提供更全面的变异图景。这两个工具的输入正好就是 STAR 比对后的 BAM不需要额外重新比对增加的时间成本很小。另一个扩展方向是等位基因特异性表达分析。RNA-seq 数据里本身就包含了表达量的信息如果某个基因的体细胞突变导致该等位基因的表达发生了明显偏移这种 ASE 信号是很有价值的生物学线索。可以在 VCF 变异位点的基础上用ASEReadCounter或类似工具统计每个变异位点的等位基因特异常 reads 数分析变异对表达的影响。第三个方向是把变异检测和剪接事件分析联动起来。RNA-seq 变异检测中发现的剪接位点突变可以通过rMATS或者MAJIQ这类工具做差异剪接分析直接看到该变异是否导致了异常的剪接产物。这个“突变-剪接-表达”的联动分析是 RNA-seq 数据独有的优势DNA 测序完全无法替代。我在实际工作中最常用的组合是Sentieon 做变异检测STAR-Fusion 做融合基因筛选再用rMATS做剪接分析三套结果拼在一起基本能覆盖转录组层面的主要变异类型。每套结果单独看都有局限但放在一起互相印证就能拼出一幅比较完整的肿瘤转录组变异图景。如果你也在用这套流程不妨往这个方向试试应该能挖出不少有意思的东西。
RNA-seq变异检测实战：基于Sentieon的全流程解析与调优
1. 内容整体设计与思路拆解1.1 RNA-seq 变异检测为什么绕不开 Sentieon做生信的人应该都有这种体会RNA-seq 数据里挖变异比 DNA 数据麻烦得多。DNA 变异检测的流程早就成熟了GATK 的 Best Practices 一套组合拳打下来大家照着做就行。可 RNA-seq 不一样splice junction 的存在让比对环节就变得很棘手reads 跨过外显子接头时如果不做特殊处理比对软件很容易把它们当成“垃圾”序列丢掉或者硬生生比对到错误的位置上。比对错了后面的变异检测就是空中楼阁。Sentieon 这一套 RNA-seq 变异检测流程解决的正是这个痛点。它提供的 RNA-seq 短序列比对方案在比对阶段就引入了对 splice junction 的感知能力再配合它那套出了名的、与 GATK 数学原理高度一致的变异检测算法能把 RNA-seq 数据的变异检测精度拉到一个非常可用的水平。很多人一提到 Sentieon 就只想到速度快其实它的 RNA-seq 流程在设计上更值得琢磨速度只是表象真正厉害的是它对 RNA-seq 数据特有噪声的处理方式。这条流程适用的场景很广既有肿瘤 RNA-seq 样本的体细胞变异筛查也有遗传病研究中基于转录组的罕见变异发现还包括一些无法获取 DNA 样本、只能做转录组测序的回顾性研究。只要你手里的数据是链特异性或者非链特异性的 RNA-seq 短读长数据想从中挖出 SNV 和 InDel这套流程就是一个相当稳妥的参考基准。可能有人会问RNA-seq 本身就不是做变异检测的首选技术为什么还要费劲去跑这套流程RNA-seq 的优势在于它能同时反映基因表达和变异信息对于某些特殊类型样本比如 FFPE 保存多年的组织DNA 可能已经降解得不成样子但 RNA 反而还能提取出可用的信息。另外RNA-seq 变异检测可以帮助验证 DNA 水平发现的剪接区域变异是否真的影响了转录本这种“转录组验证”的视角是纯 DNA 测序给不了的。1.2 项目选型时我为什么没用 GATK 而选了 Sentieon先说实话GATK 的 RNA-seq 变异检测流程本身没毛病它能用但用起来确实有点憋屈。RNA-seq 数据量通常比 WES 大不少动辄上百 G 的 fastqGATK 那套流程跑下来耗时经常以“天”为单位计算。如果手头有几十个样本光等结果就能等到怀疑人生。而且 GATK 官方对 RNA-seq 变异检测的支持力度一直不如 DNA 流程很多参数需要自己反复调文档里语焉不详的地方也不少。Sentieon 打动我的有两点。第一点自然是速度这没得洗。Sentieon 把 HaplotypeCaller 的核心算法用底层优化重新实现了一遍在相同硬件条件下加速比通常在 5 到 15 倍之间。我实测过一批 100X 左右的 RNA-seq 肿瘤样本GATK 跑一个样本的 HaplotypeCaller 大概要 8 小时Sentieon 同等配置下 40 分钟出头就跑完了。这还只是单样本的差距大规模队列跑下来省出来的时间够做很多事了。第二点是它在 RNA-seq 流程上做的针对性优化。Sentieon 不需要像 GATK 那样强制你用 STAR 比对后再做一系列繁琐的预处理它对 Input 比对结果的要求更灵活而且它对 splice junction 信息的使用方式也更直接。它有一套推荐参数是针对 RNA-seq 的比如在比对时推荐使用特定工具生成带有NMtag 的 BAM变异检测时能更准确地利用转录组特征。当然选型还要考虑成本。Sentieon 是商业软件授权费不便宜GATK 免费且开源。如果是个人学习、发小文章GATK 完全够用但如果项目有工期压力、样本量又大把时间成本折算进去Sentieon 其实并不算贵。我自己在项目里倾向于混合使用小样本探索用 GATK批量正式跑用 Sentieon两边结果做个交叉验证心里更有底。1.3 这套流程能解决的核心问题与适合人群这套 RNA-seq 变异检测全流程能解决的最核心问题是在 RNA-seq 数据中可靠地检测出 SNV 和短 InDel尤其是在涉及剪接位点附近的变异时减少比对错误带来的假阳性。RNA-seq 数据的变异检测假阳性率天然高于 DNA 数据这是转录组本身的结构决定的可变剪接、RNA 编辑、序列比对的多重映射都会在突变位点上制造噪声。Sentieon 的流程通过比对阶段的 splice-aware 处理和变异检测阶段的局部重组装把假阳性控制在一个可接受的范围内。适合参考这套流程的人我总结下来有三类。第一类是肿瘤方向的生信工程师需要从 RNA-seq 数据中寻找潜在的驱动突变特别是那些在 DNA 水平上被漏掉的、只在转录组水平表达的融合或点突变。第二类是遗传病诊断方向的研究人员当 DNA 测序没有找到明确致病位点时补充分析 RNA-seq 数据看是否有异常剪接或等位基因特异性表达的问题。第三类是刚刚接触 RNA-seq 分析、想建立一套标准流程的入门者这套流程的逻辑清晰、参数明确照着做一遍就能理解 RNA-seq 变异检测的各个环节。我特别想提醒的一点是这套流程虽然是目前 RNA-seq 变异检测的标杆但它并不能替代 DNA 测序的变异检测。RNA-seq 只能检测到在转录组中有表达的基因表达量低的基因或者被 NMD 途径降解的突变转录本在 RNA 水平上根本看不到。所以 RNA-seq 变异检测的结果更准确地说是“转录组层面的表达变异”阴性结果并不能排除 DNA 水平存在变异。2. 核心细节解析与实操要点2.1 从 fastq 到 BAM比对环节的“知易行难”RNA-seq 流程的第一步也是最关键的一步就是把原始测序数据比对到参考基因组上。这一步做得不好后面再怎么调参都救不回来。DNA 测序比对时只需要考虑 reads 在基因组线性坐标上找位置就行但 RNA-seq 的 reads 经常跨过内含子一条 read 的 5 端比对外显子 A3 端比对到外显子 B中间隔着几千甚至几万 bp 的内含子。如果比对软件不认识剪接位点这种 read 就会被打上“无法比对”的标签丢掉或者错误地比对到其它有相似序列的区域。所以我强烈建议RNA-seq 的比对环节一定要选择支持 splice junction 的比对软件。Sentieon 官方推荐的比对策略中STAR 是首选。STAR 的比对速度极快它的 seed-search 和 stitching 策略非常契合 RNA-seq 数据的特点而且它输出的 BAM 里包含NH、HI等 tag这些 tag 对后续 Sentieon 变异检测的去重环节非常有用。比对时有一个参数非常重要--outSAMattributes。如果用 STAR 比对建议在参数里加上NH HI AS NM MD这样输出的 BAM 文件里就有了比对质量、错配信息等关键属性。Sentieon 的LocusCollector去重步骤会依赖这些 tag如果缺失可能会影响去重的准确性。我自己踩过这个坑一开始没注意 STAR 的输出属性配置结果下游去重时总是报一些莫名其妙的信息加上这几个属性之后立刻清爽了。还要提醒的是RNA-seq 比对时普遍存在 multi-mapping 的问题。转录组里高度同源的基因家族、假基因、重复序列都会导致一条 read 能同时比对到基因组的多个位置。对于变异检测来说multi-mapping 的 reads 是很大的噪声来源。我推荐在比对后过滤掉某些比对质量过低或标记为 multi-mapping 的 reads但过滤尺度要把握好如果过滤太狠一些真正来自转录本高表达区域的 reads 也会被误伤。2.2 数据预处理里容易被忽略的细节比对完成后接下来的数据预处理环节看似简单但细节决定了最终结果的质量。这个环节主要包括排序、去重复、碱基质量校正三步。传统流程里这三步要分别调用不同的工具Sentieon 的优势在于它用一条命令把多个工具串联起来了省去了频繁读写中间文件的 I/O 开销这也是它速度快的一个重要原因。排序这一步没什么好说的按坐标排序是必须的关键是SortSam之后要建立 index否则很多下游工具会拒绝工作。Sentieon 的做法比较聪明它通过--sort_order coordinate参数直接在流程中完成排序和建索引不用中间环节额外处理。去重复这一步要特别讲一下它可能是整个预处理里对变异检测结果影响最大的一步。RNA-seq 测序时由于反转录和 PCR 扩增的偏好性某些高表达基因的片段会被大量重复测到这些重复 reads 如果不去掉会在变异检测时被当作独立的证据堆高假阳性。但 RNA-seq 的去重和 DNA 不太一样同一个基因的不同转录本异构体可能共享相同的外显子区域导致它们的 reads 在基因组坐标上看起来完全一样但实际上来自不同的 RNA 分子。如果按照 DNA 测序的标准去重这些 reads 会被错误地当成 PCR 重复去掉造成表达量和变异频率的偏差。Sentieon 的去重算法考虑到了这一点通过LocusCollector收集 reads 的坐标和标签信息再利用Dedup进行去重。在实际使用中我会把--rmdup相关参数逐个检查一遍确保去重策略符合自己项目的需求。如果是做表达量分析有人会选择不去重但做变异检测我还是建议去重否则重复序列区域的变异检测结果基本没法看。碱基质量校正BQSR是另一个容易被忽视但很重要的环节。测序仪给出的碱基质量值往往跟真实的错误率有偏差比如某些 motif 后面容易出现特定类型的错误。BQSR 就是通过对已知变异位点dbSNP 等的实际观察建立错误率模型重新校准碱基质量值。Sentieon 的BaseQualityScoreRecalibration和ApplyBQSR对应 GATK 的 BQSR 流程但速度更快。RNA-seq 数据里有些位点的碱基组成有偏好性BQSR 的纠偏效果在 RNA-seq 上尤其明显。2.3 变异检测阶段的参数选择与调优思路变异检测是整个流程的核心环节。Sentieon 的Haplotyper算法和 GATK 的 HaplotypeCaller 思路一致在目标区域发现变异信号后提取相关的 reads对这些 reads 做局部 de Bruijn 图的组装构造出可能的单倍型再把 reads 比对回单倍型上最后用隐马尔可夫模型计算每个位点的基因型概率。这种局部重装配的策略在处理 RNA-seq 数据时非常关键因为它能处理 splice junction 区域的复杂变异模式。RNA-seq 变异检测时有几个参数需要特别注意。第一个是--min_base_qual默认值是 10但 RNA-seq 数据在剪接位点附近经常会由比对软件引入一些质量偏低的碱基如果这个阈值太高可能会漏掉真实变异。我一般会把这个参数设为 15 左右然后结合具体数据表现再调整。第二个是--max_reads_per_alignment_start控制每个位点最多参与计算的 reads 数量RNA-seq 数据在高表达区域会出现大量 reads 堆叠这个参数默认值往往不够用导致高覆盖区域的变异检测效果变差。我会根据数据的平均深度适当调高这个值。还有一个非常实际的问题RNA-seq 数据里覆盖度在基因之间的差异非常大。高表达基因的覆盖深度可能达到几千甚至上万 X低表达基因可能只有个位数 X。这种极端的覆盖度不均一性对变异检测算法的碱基质量模型是一个巨大的考验。Sentieon 在默认参数下对高覆盖区域做了一定的降采样处理避免计算资源被单一区域耗尽但如果你的研究重点恰好是高表达基因的变异可能需要手动调整采样参数。在变异检测结果的过滤上我倾向于不过度依赖单一的过滤标准。RNA-seq 的变异检测结果需要用多种信息交叉验证比如比对质量、覆盖深度、链特异性、是否位于外显子-内含子边界等。Sentieon 的介绍里也强调了这一点它生成的 VCF 文件里包含了丰富的注释字段充分利用这些字段做过滤比死记硬背一套“万能过滤参数”要靠谱得多。3. 实操过程与核心环节实现3.1 环境准备与测试数据获取先把环境跑通这是所有后续工作的基础。Sentieon 提供了多种安装方式最常见的是下载官方编译好的安装包。安装本身不难但要注意版本匹配问题Sentieon 版本更新很快不同版本对 GATK 兼容模式的支持略有差异建议直接用最新的稳定版本。安装好之后测试数据可以从 Sentieon 官网的教程目录下载里面有配套的 RNA-seq 示例数据和参考基因组的对应区域。如果没有现成的测试数据也可以用公开数据库里的 RNA-seq 数据比如 ENCODE 项目或者 GTEx 项目里的人类 RNA-seq 数据。第一次跑流程强烈建议先用小规模的测试数据把全流程跑通确认每一步的输出都符合预期再上真实的大规模数据。这个习惯帮我避免过很多低级错误比如路径写错、参考基因组版本不对之类的问题。参考基因组的选择也是一个需要提前确认的点。RNA-seq 变异检测推荐使用完整的参考基因组包括所有染色体 contig这样比对软件才能正确处理来自线粒体等区域的 reads。另外参考基因组的版本必须和比对软件及 Sentieon 的配置一致混合使用不同版本会导致染色体命名对不上、变异位点坐标错位等一系列问题。3.2 用 STAR 完成 splice-aware 比对STAR 比对 RNA-seq 数据之前需要先建立基因组索引。这一步不是特别耗时但需要注意构建索引时加上--sjdbGTFfile参数指定基因注释文件GTF/GFF3。这样 STAR 在比对时就能参考已知的剪接位点显著提高 splice junction reads 的比对准确性。如果省略这一步STAR 就只能依赖从头发现的剪接位点对已知剪接位点的支持会弱很多。STAR 建索引的命令大致是这样STAR --runMode genomeGenerate \ --genomeDir /path/to/star_index \ --genomeFastaFiles /path/to/reference.fa \ --sjdbGTFfile /path/to/annotation.gtf \ --sjdbOverhang 149 \ --runThreadN 20sjdbOverhang这个参数需要根据测序读长来设置一般推荐设为读长-1比如 150bp 读长就设为 149。这个值是给 STAR 在已知剪接位点两侧延伸的序列长度设得太小可能覆盖不了完整的剪接位点特征太大则会浪费内存。建索引时还需要注意内存要够用人类基因组加 GTF 注释的索引构建建议至少准备 30G 内存。比对命令的写法如下STAR --genomeDir /path/to/star_index \ --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM Unsorted \ --outSAMattributes NH HI AS NM MD \ --outFilterMultimapNmax 20 \ --outFilterMismatchNmax 10 \ --runThreadN 20这里面有几个参数值得细说。--outFilterMultimapNmax 20表示允许 read 最多比对到 20 个位置超过就会被过滤掉。转录组里 repeats 区域产生的 multi-mapping reads 非常多这个阈值设得太高会让大量多映射 reads 进入下游分析反而增加噪声设得太低又可能丢掉来自同源基因家族的有用 reads我常用的经验值是 10 到 20 之间根据项目需求调整。--outSAMattributes列出的这些属性很重要特别是NM错配数和MD错配字符串后续 Sentieon 的质量校正环节会用到。STAR 默认输出 Unsorted BAM这个后面要交给 Sentieon 做排序。这一步输出的 BAM 文件非常大建议在比对时开启--outBAMcompression 6之类的压缩选项可以节省不少磁盘空间。3.3 Sentieon 流程的完整命令与参数解读比对完成后的 Sentieon 流程可以用一条 slog 串联起来。这里给出一个比较完整的 RNA-seq 变异检测命令模板同时写了注释方便对照sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo LocusCollector \ --fun score_info \ sample.sorted.bam \ sample.score sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo Dedup \ --score_info sample.score \ sample.sorted.bam \ sample.dedup.bam sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo BaseQualityScoreRecalibration \ sample.dedup.bam \ sample.recal.table sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo ApplyBQSR \ -i sample.recal.table \ sample.dedup.bam \ sample.recaled.bam sentieon driver -t 20 \ -r /path/to/reference.fa \ --algo Haplotyper \ --genotype_model multinomial \ sample.recaled.bam \ sample.vcf拆开来看每一步的作用。LocusCollector是在统计每个基因组位点上的 read 起始位置信息为去重做准备。它的输出是一个文本格式的 score 文件不是 BAM所以用重定向而不是-o。这一步的核心是标记 PCR 重复 reads 的候选位置它不需要真实输出 BAM所以速度很快。Dedup读取上一步的 score 文件真正执行去重操作。输出文件是去除重复后的 BAM。这里可以加上--rmdup参数表示在去除重复的同时也物理性地删除标记为重复的 reads否则只是给 reads 加上duplicateflag实际 reads 还保留在 BAM 里。从节省空间的角度我一般会开启--rmdup不过如果是后续还要做表达量分析保留重复 reads 也有道理具体看需求。BaseQualityScoreRecalibration生成碱基质量校正表它不会直接修改 BAM而是把校正参数输出到recal.table里。这一步需要提供已知变异位点数据库可以通过--known_sites参数传入 dbSNP 的 VCF 文件。如果没有已知位点数据库也可以不传校正效果会打折但流程仍然能跑通。ApplyBQSR是真正应用校正表的步骤把重校准后的质量值写回 BAM 文件。从这一步开始BAM 就已经是可用的“干净数据”了。Haplotyper是变异检测的核心命令。这里有一个小细节--genotype_model multinomial是适用于 RNA-seq 数据的模型选择它会用多项分布模型来评估基因型概率这个模型在杂合位点的检测上表现更好。如果这里不加这个参数默认模型也能跑但对 RNA-seq 数据的适配性差一些建议显式指定。3.4 从 VCF 到可用结果过滤、注释与解读Haplotyper输出的原始 VCF 文件直接拿去用是不行的。原始 VCF 里包含了很多低质量的候选位点需要一个过滤步骤把这些噪声去掉。Sentieon 官方推荐使用VariantFilter配合一组 RNA-seq 特异的过滤条件但实际使用中我更喜欢把 VCF 导出来自己在 Python 或者 R 里按照项目的具体需求做针对性过滤灵活性更高。常用的过滤思路包括过滤掉深度小于 10 的位点、过滤掉 QUAL 值小于 30 的位点、过滤掉链偏差显著的位点。RNA-seq 数据的链偏差问题尤其严重因为转录组测序本身就有链特异性如果两条链的覆盖度差异过大变异检测结果很容易出现假阳性。过滤完的 VCF 还需要做功能注释最常用的是SnpEff或VEP把变异位点映射到基因和转录本上标出外显子区、剪接位点区、UTR 区等位置信息。注释这一步特别重要因为 RNA-seq 检测到的变异最终要跟基因功能联系起来才有意义。我用过一个比较省力的组合Sentieon Haplotyper出原始 VCF接着用SnpEff做注释再把注释结果丢进ANNOVAR做进一步的数据库交叉比对。这个流程跑熟之后从 raw fastq 到最终可解释的变异列表一个样本大约需要 2 到 3 小时其中大头时间是 STAR 比对和 Haplotyper 这两步。4. 常见问题与排查技巧实录4.1 STAR 比对阶段的高频报错与解决方案STAR 比对报错是 RNA-seq 流程里最常碰到的问题。我总结了几个高频问题大家遇到可以直接对照排查。第一个是内存不足。构建人类基因组索引时默认参数下大概需要 30G 内存如果服务器内存不够会出现FATAL ERROR: not enough memory之类的提示。解决办法是降低--sjdbOverhang的值或者用--genomeSAindexNbases参数减少索引的规模。另外把--runThreadN调小一度内存压力也会降一点虽然速度会慢一些。第二个是zcat命令找不到。在读入压缩 fastq 时--readFilesCommand zcat依赖系统里有 gzip 解压工具有些精简版系统里没有装。解决办法是把zcat换成gzip -dc或者干脆先解压 fastq 再输入前提是磁盘空间足够。第三个是 GTF 文件和参考基因组版本不匹配。这个问题比较隐蔽表现在比对结果上就是剪接位点识别率异常低、比对率也偏低。解决办法是建立索引前严格确认 GTF 和 fasta 来自同一个版本比如都来自 Ensembl release 110或者都来自 GENCODE v44。还有一个容易被忽视的问题是STAR 比对率正常但后续 Sentieon 流程跑出来的变异数少得离谱。这种情况多半是 BAM 文件的NMtag 缺失导致的。如果比对时没有加上--outSAMattributes NM MDSentieon 在质量校正阶段就可能出错。解决办法是回到 STAR 比对那一步把参数补全重新比对。如果已经生成了大量 BAM没法重新比对也可以用samtools fillmd之类的工具把NM和MDtag 补上但效果还是不如一开始就加参数来得干净。4.2 Sentieon 运行时的常见错误与处理Sentieon 运行时的报错多数是输入文件或参数的问题。最常见的是参考基因组字典文件缺失。在跑Haplotyper之前实际上不强制要求提供 dict 文件但许多辅助功能比如按染色体并行、区间切分会依赖它。如果报错提示找不到.dict文件用picard CreateSequenceDictionary或者samtools dict生成一个放在参考基因组同目录下即可。另一种常见问题是 BAM 文件的 header 不完整。当 STAR 输出的 BAM 用samtools sort处理过而没有保留原始 header 时Sentieon 可能因为读不到正确的SQ信息而报错。解决办法是不要手工对 BAM 做多余处理直接把 STAR 的 Unsorted BAM 交给 Sentieon由它内部完成排序和索引这样最稳妥。如果是跑队列任务还容易遇到内存分配不足的问题。Sentieon 的driver -t 20表示用 20 个线程但线程数和内存是两回事如果每个线程分配的内存太少会出现Cannot allocate memory之类的报错。解决办法是在提交任务时给足内存常规人类 RNA-seq 样本建议至少 16G高深度样本建议 32G 以上。4.3 RNA-seq 变异结果的假阳性排查思路跑完流程拿到 VCF才是真正考验功力的开始。RNA-seq 变异检测的假阳性问题比 DNA 数据严重得多所以一定要有一套系统的排查思路。第一个要排查的是链偏差。打开 BAM 文件检查变异位点的 supporting reads 是否集中在某一条链上。如果是大概率是链特异性建库导致的技术偏差不是真实的体细胞突变。判断标准是看 the reads supporting the variant 的XStagSTAR 输出的链方向 tag如果所有突变 reads 都来自同一条链就要高度警惕。第二个要排查的是剪接位点附近的变异。位于外显子-内含子边界附近的变异很容易因为比对错误造成假阳性。特别是内含子侧翼 1~2bp 的位置那里是剪接体的识别核心区域如果检测到变异必须回到 IGV 里人工检查比对情况。我见过很多次看似“重要”的剪接位点突变实际上是 reads 比对到假剪接位点造成的假象。第三个要排查的是 RNA 编辑事件。RNA-seq 数据里常见的 A-to-I 编辑会把基因组上的 A 检测成 G。如果发现变异集中在已知的 RNA 编辑位点附近特别是 Alu 元件区域要警惕这是 RNA 编辑信号而不是真正的基因组变异。目前有一些 RNA 编辑位点数据库可以做过滤比如 REDIportal但要注意它对不同组织类型的覆盖度并不均衡。最后如果条件允许强烈建议对 RNA-seq 变异结果做 Sanger 测序验证或者用 DNA 测序数据做个交叉验证。RNA-seq 变异检测的结论尤其是那些影响临床决策的重要位点一定要有独立的技术验证才能下结论。这不是对流程的不信任而是分子生物学实验的基本严谨性。4.4 处理多个样本时的内存与运行时间管理单个样本跑通之后真实项目通常面临的是几十甚至上百个样本的批量运行。这时候需要考虑的是如何高效地管理资源而不是盲目地堆线程。先看时间。RNA-seq 样本的比对和变异检测单样本大约占用 2~4 小时的计算时间具体取决于数据量和服务器配置。如果有 50 个样本按照单样本 3 小时算串行跑需要 150 小时这显然是不可接受的。解决办法是并行化同一时刻跑多个样本的流程同时留出足够的 CPU 和内存资源。再看内存。Sentieon 的Haplotyper在高覆盖区域会消耗较多内存如果多个样本同时运行内存峰值会有叠加。我的经验是每同时运行 4 个样本至少准备 64G 内存。如果内存紧张可以通过--interval参数把基因组拆分成多个区间并行处理跑完再合并 VCF这样单样本的峰值内存可以降下来同时整体吞吐量也有提升。还有一个小技巧是监控中间文件大小。STAR 输出的 BAM 动辄几十 G如果一组样本同时跑磁盘空间很容易被打爆。建议在处理前先规划好目录结构给每个样本独立子目录并定时清理中间文件。流程跑完后保留最终的 VCF 和重校准后的 BAM 就够了原始比对产生的超大 BAM 可以先用samtools view -b -F 1024过滤掉重复 reads 再压缩存档能省下不少空间。5. 流程验证与结果质控5.1 用已知变异样本验证流程的准确性跑流程之前最好先用一个已知变异位点的样本做验证。Sentieon 官网和公开数据库都提供了一些带有真实变异注释的 RNA-seq 测试数据比如 NA12878 的 RNA-seq 数据它对应的高置信度变异集可以从基因组学权威资源中获取。验证的方法是把流程跑完得到 VCF 文件然后和已知变异集做比较。比较的指标包括 Recall召回率、Precision精确率和 F1 score。如果某个指标明显偏低说明流程中某些环节有问题需要回头检查。我见过比较多的坑是向量子集比较时参考基因组版本不一致导致同一变异位点的坐标对不上这会造成召回率虚低排查时往往浪费很多时间。实际操作里我是这么做的先用bcftools isec或者RTG Tools vcfeval把检测到的变异与已知变异集比较重点看常见变异位点是否都能检测到以及检测到的变异里有多少是已知集里没有的。vcfeval在处理复杂位点时比bcftools isec更准确因为它会做等位基因级别的比对而不是简单的坐标比较。5.2 覆盖度、比对率和样本一致性检查RNA-seq 样本之间差异很大同一个批次的数据也可能有不同的质量表现。所以在正式做变异检测之前先检查数据的几个核心质控指标比对率STAR 输出日志里会报告唯一比对 reads 的比例一般要求大于 80%如果低于这个值说明样本可能存在降解或者污染。基因体覆盖均匀度检查 reads 在转录本上的覆盖是否均匀如果 3 端严重偏向提示 RNA 质量可能有问题。重复率去重后剩余 reads 的比例高重复率可能说明起始 RNA 量不足或者 PCR 循环过多。链特异性如果建库是链特异性需要确认链方向是否和预期一致。这些指标可以用Picard CollectRnaSeqMetrics和MultiQC汇总检查。Per-sample 都通过质控后再做变异检测就能减少由于样本质量差异引入的假信号。有个情况特别提醒一下肿瘤样本可能会有较大的拷贝数变异或者杂合性缺失区域导致某些染色体区域的 reads 覆盖度显著偏离平均水平。这种区域在变异检测时经常会出现大片的连续变异信号看起来像是突变热点其实就是拷贝数变化导致的技术信号。分析这类样本时最好把拷贝数信息和变异检测结果联合看别被假象带偏了。5.3 变异结果的生物学合理性评估最后的质控关卡是评估变异结果的生物学合理性。即使所有技术指标都正常也不能保证变异结果是真实可用的。先从突变频谱入手。人类基因组中CT 转换通常是最常见的突变类型。如果检测到某个样本的突变频谱明显异常比如 CA 或 AT 颠换占主导有可能是测序错误或者样本处理过程中引入了氧化损伤。这种频谱异常在 FFPE 样本中特别常见是人工假象的典型标志。再看等位基因频率分布。胚系变异的等位基因频率一般集中在 0.5 和 1.0 附近体细胞变异则呈现低频分布。如果大量变异的等位基因频率集中在一个奇怪的档位比如 0.4 左右可能是样本污染也可能是拷贝数异常。这个判断需要结合样本类型做综合考量。最后把变异映射到通路和基因功能上。RNA-seq 变异检测的结果应该跟样本的表型或者疾病背景有逻辑关联。如果发现一个正常组织样本里出现了大量高频的已知驱动基因突变那先别高兴多半是样本标签搞错了或者比对时发生了样本间的交叉污染。这种“做出来结果太完美反而可疑”的场景在实际项目中并不少见。6. 项目总结与实操心得6.1 我跑完这套流程后的几点体会前前后后用 Sentieon 跑了不少 RNA-seq 样本积累了一些体会。第一个体会是RNA-seq 变异检测的快和准本质上是“比对阶段打底、变异检测阶段兜底”。STAR 的 splice-aware 比对决定了基础比对质量而 Sentieon 的 Haplotyper 通过局部重装配把变异检测精度拉高了一个档次。两者缺一不可。如果比对阶段就丢掉了大量剪接位点 reads后面的检测能力再强也无济于事。第二个体会是RNA-seq 变异检测的结果一定要做技术验证。Sentieon 流程给出的 VCF 文件在统计意义上是可信的但它毕竟是从转录组数据反推基因组变异受到表达水平、降解程度、建库偏好等多种因素影响。批量分析时可以用算法过滤掉大部分假阳性但真正的“金标准”位点比如临床报告里要写的位点绝对要用 Sanger 或者靶向重测序验证。这不是流程的缺陷而是 RNA-seq 数据本身的边界。第三个体会是流程的工程化能力。Sentieon 之所以适合大规模队列是因为它可以把比对之外的步骤全部用一条流水线串起来避免了中间文件的反复读写。在一个 100 样本的肿瘤 RNA-seq 队列中我用这套流程跑完发现耗时主要花在 STAR 比对和 Haplotyper 上去重和质量校正的速度几乎感觉不到存在。这种工程效率的提升是实际项目里能实打实感受到的。6.2 几个值得提前规划的点如果是从零开始搭建这套流程有几点值得提前规划好。第一磁盘空间一定要提前摸排。STAR 输出的大 BAM、Sentieon 去重后的BAM、重校准后的 BAM每个都是原始 fastq 的 2~3 倍大小三个版本叠加起来磁盘占用非常夸张。我建议流程跑完后立刻清理不需要的中间 BAM只保留最下游的重校准 BAM 和 VCF 文件。如果项目周期长中间文件可以压缩存档到冷存储。第二参考基因组和注释文件的版本管理一定要严格。我在多个项目里都吃过版本不一致的亏最典型的例子是参考基因组用了 GRCh37而 dbSNP 用了 GRCh38 的版本导致已知位点数据库在比对时大范围错位质量校正形同虚设。建议在项目启动时就把参考基因组、GTF 注释、已知变异数据库的版本固定下来用统一的配置记录在项目文档里。第三流程的可重复性要靠脚本化保证。不要手工一句一句地敲命令把整个流程写成 shell script 或者用流程管理工具Snakemake、Nextflow封装起来参数、路径、版本全部集中在一个配置文件里。这样不仅自己能复用团队协作时别人也能轻松复现你的结果。6.3 后续功能扩展的思路这套流程跑通之后可以沿着几个方向做扩展。最直接的扩展是融合基因检测。RNA-seq 数据除了点突变还能用来找融合基因。可以在 Sentieon 变异检测的基础上额外跑一遍STAR-Fusion或者Arriba把融合基因的检测结果和点突变结果合并分析为肿瘤研究提供更全面的变异图景。这两个工具的输入正好就是 STAR 比对后的 BAM不需要额外重新比对增加的时间成本很小。另一个扩展方向是等位基因特异性表达分析。RNA-seq 数据里本身就包含了表达量的信息如果某个基因的体细胞突变导致该等位基因的表达发生了明显偏移这种 ASE 信号是很有价值的生物学线索。可以在 VCF 变异位点的基础上用ASEReadCounter或类似工具统计每个变异位点的等位基因特异常 reads 数分析变异对表达的影响。第三个方向是把变异检测和剪接事件分析联动起来。RNA-seq 变异检测中发现的剪接位点突变可以通过rMATS或者MAJIQ这类工具做差异剪接分析直接看到该变异是否导致了异常的剪接产物。这个“突变-剪接-表达”的联动分析是 RNA-seq 数据独有的优势DNA 测序完全无法替代。我在实际工作中最常用的组合是Sentieon 做变异检测STAR-Fusion 做融合基因筛选再用rMATS做剪接分析三套结果拼在一起基本能覆盖转录组层面的主要变异类型。每套结果单独看都有局限但放在一起互相印证就能拼出一幅比较完整的肿瘤转录组变异图景。如果你也在用这套流程不妨往这个方向试试应该能挖出不少有意思的东西。