专栏名称: 生信媛
生信媛,从1人分享,到8人同行。坚持分享生信入门方法与课程,持续记录生信相关的分析pipeline, python和R在生物信息学中的利用。内容涵盖服务器使用、基因组转录组分析以及群体遗传。
目录
相关文章推荐
51好读  ›  专栏  ›  生信媛

Sentieon:NGS下游变异检测高效工具使用体验

生信媛  · 公众号  · 生物  · 2019-07-04 12:14

正文

请到「今天看啥」查看全文


安装

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

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

  2. 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个线程

一些测试结果:

  1. sentieon bwa建立的索引和原版的bwa索引可以共用

  2. 80个线程下,原版和sentieon版的运行用top看的时候都在8000%左右。

用time统计运行时间,结果如下:

  1. # sentieon bwa比对

  2. 123772.31s user 3026.78s system 6601% cpu 32:00.84 total

  3. # sentieon util排序

  4. 5775.55s user 659.32s system 319% cpu 33:35.59 total

  5. # bwa比对

  6. 208744.12s user 2966.60s system 5942% cpu 59:22.51 total

  7. # samtools sort排序

  8. 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)建立索引文件

  1. mkdir index && cd index

  2. sentieon bwa index reference.fasta

  3. samtools faidx reference.fasta

  4. gatk CreateSequenceDictionary -R reference.fasta

然后创建输入文件

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

编写批量运行的脚本

  1. #!/bin/bash

  2. # runing with bash align.sh

  3. set -e

  4. set -u

  5. set -o pipefail


  6. SAMPLES=samples.txt

  7. REFERENCE=index/reference.fa

  8. NUMBER_THREADS=80

  9. SENTIEON_LICENSE=license.lic


  10. export SENTIEON_LICENSE

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


  12. # read align

  13. # \\t rather than \t

  14. exec 0< $SAMPLES

  15. while read sample;

  16. do

  17. (sentieon bwa mem -M -R "@RG\\tID:${sample}\\tSM:${sample}\\tPL:ILLUMINA" \

  18. -t $NUMBER_THREADS $REFERENCE 0-raw-data/${sample}_1.fq.gz 0-raw-data/${sample}_2.fq.gz || echo -n 'error' ) \

  19. | sentieon util sort -r $REFERENCE -o 1-align/${sample}_sort.bam -t $NUMBER_THREADS --sam2bam -i -


  20. done


  21. # Calculate data metrics

  22. mkdir -p metrics

  23. for BAM in 1-align/*_sort.bam;

  24. do

  25. prefix=$(basename $BAM)

  26. # get metrics

  27. sentieon driver -t $NUMBER_THREADS -r $REFERENCE -i $BAM \

  28. --algo GCBias --summary metrics/${prefix}_GC_SUMMARY.txt metrics/${prefix}_GC_METRIC.txt \

  29. --algo MeanQualityByCycle metrics/${prefix}_MQ_METRIC.txt \

  30. --algo QualDistribution metrics/${prefix}_QD_METRIC.txt \

  31. --algo InsertSizeMetricAlgo metrics/${prefix}_IS_METRIC.txt \

  32. --algo AlignmentStat metrics/${prefix}_ALN_METRIC.txt

  33. # plotting

  34. sentieon plot QualDistribution -o metrics/${prefix}_QD_METRIC.pdf metrics/${prefix}_QD_METRIC.txt

  35. sentieon plot InsertSizeMetricAlgo -o metrics/${prefix}_IS_METRIC_PDF metrics/${prefix}_IS_METRIC.txt

  36. done


  37. # Remove duplicateds

  38. for BAM in 1-align/*.bam;

  39. do

  40. OUT=$ {BAM%_.*}

  41. sentieon driver -t $NUMBER_THREADS -i $BAM\

  42. --algo LocusCollector --fun score_info ${OUT}_SCORE.gz

  43. sentieon driver -t $NUMBER_THREADS -i $BAM \

  44. --algo Dedup --rmdup --score_info ${OUT}_SCORE.gz \

  45. --metrics ${OUT}_DEDUP_METRIC.txt ${OUT}_dup.bam

  46. done


  47. # Variant Calling

  48. mkdir -p 3-variants

  49. sentieon driver -t $NUMBER_THREADS -r $REFERENCE \

  50. $(for bam in 1-align /*_dup.bam; do echo "-i $bam"; done)\

  51. --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个线程。

运行的脚本代码如下

  1. #!/bin/bash

  2. # runing with bash align.sh

  3. set -e

  4. set -u

  5. set -o pipefail


  6. REFERENCE=index/reference.fa

  7. NUMBER_THREADS=40

  8. SENTIEON_LICENSE=license.lic


  9. # Variant Calling

  10. mkdir -p 3-variants

  11. sentieon driver -t $NUMBER_THREADS -r $REFERENCE \

  12. $(for bam in 1-align/*_dup.bam; do echo "-i $bam"; done)\

  13. --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。

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








请到「今天看啥」查看全文