NGS的变异检测高效工具Sentieon使用体验

2019-07-02 本文已影响124人 xuzhougeng

安装

以目前的20180806为例，下载并解压缩

wget https://s3.amazonaws.com/sentieon-release/software/sentieon-genomics-201808.06.tar.gz
tar xf sentieon-genomics-201808.06.tar.gz

之后需要将软件的安装位置加入到环境变量PATH中。

此外需要申请一个许可证文件，并且输出一个环境变量SENTIEON_LICENSE=license.lic

Sentieon读取文件的速度是 20-30Mb/每秒每核，因此为了获取极致性能，推荐使用SSD。

我的测试环境是：

Linux 3.10.0-862.14.4.el7.x86_64
Intel(R) Xeon(R) CPU E7-8890 v4 @ 2.20GHz
普通硬盘

bwa-mem测试

下面的测试基于200M基因组，213011760 条 read，约100X测序，80个线程

一些测试结果：

sentieon bwa建立的索引和原版的bwa索引可以共用
80个线程下，原版和sentieon版的运行用top看的时候都在8000%左右。

用time统计运行时间，结果如下：

# sentieon bwa比对
123772.31s user 3026.78s system 6601% cpu 32:00.84 total
# sentieon util排序
5775.55s user 659.32s system 319% cpu 33:35.59 total
# bwa比对
208744.12s user 2966.60s system 5942% cpu 59:22.51 total
# samtools sort排序
5457.31s user 377.73s system 84% cpu 1:55:41.84 total

最终输出的BAM，我用samtools idxstats比较结果时，两者输出一摸一样。

结论：是原版的bwa-mem的2倍以上速度。虽然现在有bwa-mem2了，但是据sentieon公司胡博士说，速度依旧比不上sentieon公司的优化版。同时，他们公司开发的配套排序和标记重复工具的速度也比目前开源软件快。

变异检测分析

单个样本

200M基因组，215050865条read，也是100X左右的深度，用了40个线程

标记重复: 49s + 139s

变异检测: 2402s

BSA项目

以我手头的几个BSA项目为例。

为参考序列(reference.fasta)建立索引文件

mkdir index && cd index
sentieon bwa index reference.fasta
samtools faidx reference.fasta
gatk CreateSequenceDictionary -R reference.fasta

然后创建输入文件

ls 0-raw-data | cut -d '_' -f 1 | uniq > samples.txt

编写批量运行的脚本

#!/bin/bash
# runing with bash align.sh
set -e
set -u
set -o pipefail

SAMPLES=samples.txt
REFERENCE=index/reference.fa
NUMBER_THREADS=80
SENTIEON_LICENSE=license.lic

export SENTIEON_LICENSE
export PATH=$PATH:/opt/biosoft/sentieon/sentieon-genomics-201808.06/bin

# read align
# \\t rather than \t
exec 0< $SAMPLES
while read sample;
do
(sentieon bwa mem -M -R "@RG\\tID:${sample}\\tSM:${sample}\\tPL:ILLUMINA" \
    -t $NUMBER_THREADS $REFERENCE 0-raw-data/${sample}_1.fq.gz 0-raw-data/${sample}_2.fq.gz  || echo -n 'error' ) \
    | sentieon util sort -r $REFERENCE -o 1-align/${sample}_sort.bam -t $NUMBER_THREADS --sam2bam -i - 

done

# Calculate data metrics
mkdir -p metrics
for BAM in 1-align/*_sort.bam;
do
    prefix=$(basename $BAM)
    # get metrics
    sentieon driver -t $NUMBER_THREADS -r $REFERENCE -i $BAM \
    --algo GCBias --summary metrics/${prefix}_GC_SUMMARY.txt metrics/${prefix}_GC_METRIC.txt \
    --algo MeanQualityByCycle metrics/${prefix}_MQ_METRIC.txt \
    --algo QualDistribution metrics/${prefix}_QD_METRIC.txt \
    --algo InsertSizeMetricAlgo metrics/${prefix}_IS_METRIC.txt \
    --algo AlignmentStat metrics/${prefix}_ALN_METRIC.txt
    # plotting
    sentieon plot QualDistribution -o metrics/${prefix}_QD_METRIC.pdf metrics/${prefix}_QD_METRIC.txt
    sentieon plot InsertSizeMetricAlgo -o metrics/${prefix}_IS_METRIC_PDF metrics/${prefix}_IS_METRIC.txt
done

# Remove duplicateds
for BAM in 1-align/*.bam;
do
    OUT=${BAM%_.*}
    sentieon driver -t $NUMBER_THREADS -i $BAM\
    --algo LocusCollector --fun score_info ${OUT}_SCORE.gz
    sentieon driver -t $NUMBER_THREADS -i $BAM \
    --algo Dedup --rmdup --score_info ${OUT}_SCORE.gz \
    --metrics ${OUT}_DEDUP_METRIC.txt ${OUT}_dup.bam
done

# Variant Calling
mkdir -p 3-variants
sentieon driver -t $NUMBER_THREADS -r $REFERENCE \
    $(for bam in 1-align/*_dup.bam; do echo "-i $bam"; done)\
    --algo Haplotyper 3-variants/final.vcf.gz

对于660M的两个亲本和两个子代重测序结果，换句话说就是要对4个样本进行变异建立。

时间记录：

比对+排序平均耗时: 40min（2390s，2463s，2000s，1808s）
标记重复平均耗时:6.4min (123 + 338, 63 + 193, 73 + 201, 74+ 272)
变异检测耗时: 2h (7252s)

此外我还用96个线程测试了基因组大小200M的5个样本的变异检测，运行时间是71分钟(4309)

200个群体GBS分析

我还分别测试了200个GBS样本和157个GBS样本从BAM文件到最终的VCF文件的时间，分别用的是40个线程和56个线程。

运行的脚本代码如下

#!/bin/bash
# runing with bash align.sh
set -e
set -u
set -o pipefail

REFERENCE=index/reference.fa
NUMBER_THREADS=40
SENTIEON_LICENSE=license.lic

# Variant Calling
mkdir -p 3-variants
sentieon driver -t $NUMBER_THREADS -r $REFERENCE \
    $(for bam in 1-align/*_dup.bam; do echo "-i $bam"; done)\
    --algo Haplotyper 3-variants/final.vcf.gz

对于56线程的157个GBS样本，时间约为26小时(94141 s)

对于40线程的200个GBS样本，时间约为19小时(71528 s)

差不多1天就解决了我原本需要一周时间运行的任务。

总结

GATK里用到的模型使用C语言实现而非Java，会提高最终的运行效率，其实我不怎么惊讶。但是同样是通过C语言编写的bwa，Sentieon的开发者居然还能让bwa速度提高了2倍以上，速度还优于目前bwa-mem2，这就体现出他们团队的专业了。

前段时间和Sentieon公司的胡晋南博士交流后还了解到几件事情：Sentieon公司其实比broad研究所更熟悉GATK，因为broad研究所维护GATK软件的人流动性大，有段时间所有人都到谷歌公司去上班了。这也是为啥我用GATK那么难受的原因，很多命令行参数随着版本更新不知不觉就变了；GATK4.0为了追求效率，做了很多取舍，因此反而不如3.8的准确性高；GATK速度慢并且准的原因是因为有一步local assembly，因此只对SNP感兴趣的话，我个人认为可以不用GATK。还有GATK其实还有一些BUG，Sentieon公司为了保证一致性，在其中一个版本中保留了所有的bug，此外还开发另一个优化版本。

当然Sentieon是一个收费软件，对于公司或者一个大的研究机构而言，如果追求效率那么可以考虑一下。但就针对我而言，速度可能并不是我的刚需，手头有那么多项目，我可以切换个项目。不过Sentieon还支持广州超算那边的按核时计费模式，如果以后能够按样本付费的话，对于一些小的课题组会更加好一点吧。