从GATK流程转为Sentieon流程需要注意的细节
上两期，我分别测试了Sentieon的两个分析模块，分别是DNAseq流程和TNseq流程，两个模块都是基于等同GATK最佳流程的数学模型，分析速度让人很满意。
我在上两期只是写了“飞起来”的速度,没有去展开比较：GATK、MuTect/MuTect2与Sentieon命令行有哪些不一样？GATK、MuTect/MuTect2的默认参数和Sentieon是一致的吗？为了得到有效的比较，GATK需要调节哪些特殊的参数？
带上这些问题，我找到了Sentieon的技术专家2017年三月份在biorxiv发表的文章（见阅读原文）。这一篇文章比较了GATK、MuTect/MuTect2和Sentieon的运行效率和结果的比对，同时我在文章的附录找到了它们各自相对应的命令行。总结下来，也想分享给需要的朋友。
如果你也用GATK、MuTect/MuTect2，想把自己的GATK等的命令，转换成Sentieon 命令，这篇文章的附录绝对不要错过。它可以帮你进行命令行的调整。
（原文附录GATK和Sentieon运行命令比对的详细列表请参见：
https://www.biorxiv.org/highwire/filestream/33653/field_highwire_adjunct_files/0/115717-1.pdf）
注：上图为附录的截图
我们先来回忆一下我上一期为大家整理的Sentieon的基本命令行格式：
Sentieon driver是用于pipeline的命令。单次调用driver程序可以运行多个算法比如Metrics就调用了多个algo，--algo就是告诉软件采用什么算法。
总结下来，大的区别有以下四点：
多线程
Sentieon在每一个计算环节都可以用多线程，默认值为服务器所有核数。所以我在编辑命令行的时候，就设定了默认所有的分析环境要利用服务器所有的线程去运算。
而GATK默认是1个线程，有些地方可以用到多个线程，但在很多计算步骤里面不能用到多线程。如果有的填了大于1线程数（比如HC环节，nct不能填大于1）就会报错。详细请参考（https://software.broadinstitute.org/gatk/documentation/article?id=1975）GATK4的Beta版本我也试过，很多地方比如MuTect2摒弃了以前的单线程，加速不少。但是毕竟是测试版本，没有正式发布，一旦正式发布，我会比一比。
无降采样
什么是降采样呢？降采样即在特定位置或区域内降低reads深度的过程 （具体设置见：https://gatkforums.broadinstitute.org/gatk/discussion/1323/downsampling）。
GATK流程默认是有降采样的环节。降采样的优点是可以节省计算成本，让reads深度达到平衡。但缺点是损失信息，有可能会出现批次做出来的结果不一致。如果希望得到GATK最精确的结果， -dcov设置为一个很大的值，比如1000000.
在这里Sentieon没有降采样过程，因此也没有损失信息，在这个环节就不用去做设定。
Sentieon可以实现同时对100,000个文件的联合基因分型
(见joint genotyping环节)，避免产生中间文件合并的冗长过程。
命令行界面不同：
对于GATK来说，如果要告诉软件将要采用什么算法，需要用-T告诉GATK。而使用Sentieon的时候，需要记住--algo就等同于-T。Sentieon会直接出log文件，不需要在命令行里面另外加-log。Sentieon每个算法后直接输出文件，不需要在命令行外加-o。
用惯了GATK的朋友，改用Sentieon的时候，每一步需要注意的细节是什么？这里针对DNAseq流程，和TNseq肿瘤正常配对样本流程进行说明。
DNAseq – Single Sample（单个样本举例）
分析流程：
我们一个个步骤开始看。
GATK对原始数据进行mapping，生成sam格式的比对文件，对sam文件进行重排序，排好以后再对sam文件转换成bam文件，再对bam文件进行sort排序，创建索引。这个过程中，GATK流程需要额外安装picard，BWA， SAMtools等。
Sentieon优化并简化了操作步骤，不需要再安装其他工具，因为都封装到Sentieon流程中去了。与GATK的BWA alignment过程是相同的。不同的是，Sentieon把结果从BWA-MEM直接传送到sorting阶段，从而直接产生bam文件。而GATK是把结果从BWA-MEM输出到sam文件，再经过转换，产生bam文件。
除此以外，Sentieon将自动为排序的bam文件创建一个索引文件，而GATK不会自动创建索引，需要额外的命令产生索引文件。
更重要的是Sentieon对BWA-MEM算法和内存管理进行了优化，使其运行速度不同的数据测试有1.0 – 3.9倍的加速（参考biobiorxiv文章）。
具体命令行我这里就不列了，大家可以去链接中下载和比较（https://www.biorxiv.org/highwire/filestream/33653/field_highwire_adjunct_files/0/115717-1.pdf）
Sentieon在该步骤的计算过程和picard格式一致。在这里，它使用了五个 algo: MeanQualityByCycle, QualDistribution, GCBias, AlignmentStat, InsertSizeMetricAlgo 对每个样本进行统计。
基于得到的bam文件，Sentieon有两步计算：--algo LocusCollector和--algo Dedup。第一个是收集Read 比对的信息，第二个用来执行去重复。会自动产生deduped bam的index文件，不需要再敲输出index的命令。
GATK需要-T RealignerTargetCreator和-T IndelRealigner两个命令步骤来执行Realign Indels，Sentieon需要--algo Realigner一个命令步骤即可。
这个过程分为3步，分别是计算、执行碱基质量得分的校正（calculate recalibration，apply recalibration），以及为碱基质量得分的校正出具图和报告（plot recalibration）。
GATK和Sentieon的计算方法和步骤是一样的。区别在：
GATK在执行这些步骤的时候，使用了这些命令：
-T BaseRecalibrator：Detect systematic errors in base quality scores（第一，二步）
-T PrintReads：Write out sequence read data (for filtering, merging, subsettingetc)（第二步）
-T AnalyzeCovariates：Create plots to visualize base recalibration results（第三步）
而Sentieon在执行这三步（计算、执行、报告）时，都用了-- algo QualCal 。-- algo QualCal: for the base quality recalibration stage. （第一~三步）。
此外，Sentieon如果要对碱基质量得分校正以后的BAM文件的合并这个步骤，使用--algo ReadWriter来操作。可以认为它等同于GATK的-T PrintReads。
对于UnifiedGenotyper ，GATK调用-T UnifiedGenotyper来执行，而Sentieon则为 --algo Genotyper。对于HaplotypeCaller，GATK调用-T HaplotyperCaller 来执行，而Sentieon则为 --algo Haplotyper。
针对联合基因分型（joint genotyping）阶段，GATK和Sentieon之间的执行命令也是很相似的。其中，GATK是通过-T GenotypeGVCFs来执行，而Sentieon是通过--algo GVCFtyper来执行。
此外，为了本次的分析，GATK需要重新设定-maxAltAlleles。GATK软件对此默认值为6, 在这里，需要改为比较大的值，比如20。
maxAltAlleles的意思为Maximum numberof alternate alleles to genotype。也就是说，软件只接受这么多的alternate alleles到genotyping。默认值之所以设为6，是因为GATK在这个步骤其实需要吃CPU和内存，而且是指数级别（见GATK --max_alternate_alleles的说明）。这个值设定得越大，GATK跑起来越慢。设得低，会让GATK在这个步骤不至于太慢。但是当测序的深度和测序的样本数量上去了以后(尤其是biorxiv的这篇文章这么大的数据量)，如果还是默认值，会损失很有用的信息，而且还会造成每一次运行的结果不一致。而Sentieon软件没有这些顾虑， 更不需要对这个参数有任何设置，因此Sentieon没有这个参数。
此外，Sentieon也可以进行体细胞变异的分析。
Tumor-Normal Pipeline（肿瘤-对照配对样本）
分析流程：
在分析操作的时候，从“read alignment”到“Base quality score recalibration”和DNAseq都是一样的。当分别对肿瘤和正常样本做过Base quality score recalibration以后，会获得tumor_recaled.bam和normal_recaled.bam。这两个文件可作为Indel joint realignment的输入文件。
肿瘤和正常样本的组合数据进行indel附近的局部调整。和各自进行indel realignment一样，GATK在这里也要调用-T RealignerTargetCreator和-T IndelRealigner，Sentieon仍然需要--algo Realigner来执行。
这样得到的结果，GATK还是两个文件（参数设置为-nWayOut _corealigned.bam，其中_corealigned.bam是输出文件的后缀，最后生成tumor_recaled_corealigned.bam 和normal_recaled_corealigned.bam，这两个bam作为MuTect/Mutect2的输入），Sentieon处理以后会形成一个文件（tn_corealigned.bam作为后续Variant calling的唯一输入文件）。
MuTect调用了-T MuTect来执行，Sentieon则为--algo TNsnv。MuTect2调用了-T MuTect2来执行，Sentieon则为--algo TNhaplotyper。
Sentieon一开始我设置默认调用所有的核数去处理，到这一步也是如此，它能运用到所有的核数。MuTect只支持单线程，GATK3.x MuTect2可以设为多线程，但是运行起来容易跑死，而如果设为单线程，又会特别慢。GATK4的beta版本我也测试过，核数我设置调用所有的核数，但是只能带起来个位数。原因暂时未知。
对于GATK和Sentieon出来的结果，biorxiv的那篇文章最后利用了RTG的vcfeval工具，对各自得到的VCF结果进行比较。结果在Fig4中体现了。
在biorxiv发表的文章的Fig4中，使用了F1-Score这个值作为Sentieon与GATK的结果比较参考。这个值越大，差异越小。这个值如果为1，则二者结果完全相同了。
这个值的计算意义是什么？由于篇幅关系，我下次再给大家展开说明。
关于bioRxiv那个用于比较的文章，大家先看原文（可以通过点击下方原文链接来跳转），再把附录里面的命令行（开头给出的链接）下载下来多琢磨一下。另外，原文还附有用于分析的几百例原始数据的链接，对于我这种常常找不到合适数据练手的人，帮助很大。
猜你喜欢
菜鸟入门
数据分析
ChIP-seq（上）| ChIP-seq（下）| RNA-seq | miRNA
WGS,WES,RNA-seq组与ChIP-seq之间的异同
编程实践
直播基因组分析
生信技能树编辑：吃瓜群众