DeepVariant — step3 DeepVariant、

找了一组数据，是100个基因的小panel的bam文件，用DeepVariant、GATK、Sentieon三个软件分别call SNP/INDEL，最后做一下比较。
bed文件：100.bed
samplename_N.recal.bam
samplename_N.recal.bam.bai

一、GATK

GATK call 胚系突变的流程参考：GATK Best Practices — step2 胚系突变 SNP/INDEL（Germline SNPs + Indels）

命令：

time   gatk --java-options "-Xmx6G -XX:GCTimeLimit=50 -XX:GCHeapFreeLimit=10"   \
HaplotypeCaller -R   ucsc.hg19.fasta  \
-I   samplename_N.recal.bam -L 100.bed  \
-O  GATK/ samplename_N.gatk.vcf -contamination 0  \
--dbsnp   dbsnp_138.hg19.vcf

得到：

samplename_N.gatk.vcf
samplename_N.gatk.vcf.idx

所用时间：

real    14m7.563s
user    31m37.614s
sys     0m18.182s

二、Sentieon

命令：

time  sentieon driver -t 32 -r   ucsc.hg19.fasta  \
--interval 100.bed -i  bam/samplename_N.recal.bam  \
--algo Haplotyper --bam_output   Sentieon/samplename_N.Sentieon.bam  \
-d dbsnp_138.hg19.vcf --min_base_qual 20  \
--prune_factor 2 --emit_conf 30 --call_conf 30  \
--emit_mode variant --pcr_indel_model conservative  \
Sentieon/ samplename_N.Sentieon.vcf

得到：

samplename_N.Sentieon.vcf
samplename_N.Sentieon.vcf.idx
samplename_N.Sentieon.bam
samplename_N.Sentieon.bam.bai

所用时间（简直太快了）：

real    0m43.262s
user    8m32.122s
sys     0m6.830s

三、Deepvariant

命令：

#进入镜像
docker run -it -v /tmp/:/tmp/ -v /deepvariant/:/deepvariant/ docker.io/dajunluo/deepvariant:latest  bash

/opt/deepvariant/bin/run_deepvariant --model_type=WGS \
--ref=  ucsc.hg19.fasta --reads=samplename_N.recal.bam \
--regions 100.bed --output_vcf=samplename_N.deepvariant.vcf  \
--output_gvcf=samplename_N.deepvariant.gvcf --num_shards=4

得到：

samplename_N.deepvariant.vcf
samplename_N.deepvariant.gvcf

运行时间：

#三步分别时间：
#make_examples
real  1m38.409s
user  5m17.508s
sys  0m7.040s
#call_variants
real  0m23.316s
user  6m52.188s
sys  1m22.548s
#postprocess_variants
real  0m26.652s
user  0m21.988s
sys  0m3.616s

四、三者的比较

1.运行时间

不过后来发现，在运行这三个软件的时候，所给定的线程并不一样，所以其实在这里的可比性不能完全依赖这个时间... ...

GATK 14min
Sentieon 40s （名不虚传，但是收费）
Deepvariant 2min30s （时间上可以接受，比GATK要快，而且开源免费）

2.得到的突变个数

deepvariant得到的突变个数最多，gatk和Sentieon的突变个数十分接近。

$ grep -v "#"  samplename_N.gatk.vcf | wc
   2334   23340  517745

$ grep -v "#"  samplename_N.Sentieon.vcf | wc
   2325   23250  559982

$ grep -v "#"  samplename_N.deepvariant.vcf | wc
   3216   32160  272788

用这三个vcf的1，2，4，5列作为输入画一个韦恩图，更直观：

Deepvariant、GATK和Sentieon得到的突变结果韦恩图

3.准确性

随机找几个位点，从IGV里看一下：

1. chr1 110883569 A->T
   这个位点是在三个软件里都找到的位点，IGV里显示也认为这个位点可靠。

chr1 110883569 A->T

但是，这个位点的深度在三个软件不同，其中deepvariant的和IGV中一致，GATK和Sentieon一致，GATK中的深度和真实深度不同的原因可能在于，它会过滤掉一些reads再计算深度：

$ grep -w 110883569  GATK/ samplename_N.gatk.vcf
chr1    110883569       rs3738752       A       T       4219.77 .       AC=1;AF=0.500;AN=2;BaseQRankSum=7.551;DB;DP=279;ExcessHet=3.0103;FS=6.425;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=15.18;ReadPosRankSum=-0.318;SOR=0.432     GT:AD:DP:GQ:PL  0/1:142,136:278:99:4248,0,4204


$ grep -w 110883569  Deepvariant/ samplename_N.deepvariant.vcf
chr1    110883569       .       A       T       53.5    PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:53:280:142,138:0.492857:53,0,64


$ grep -w 110883569  Sentieon/ samplename_N.Sentieon.vcf
chr1    110883569       rs3738752       A       T       4209.77 .       AC=1;AF=0.500;AN=2;BaseQRankSum=7.520;ClippingRankSum=0.000;DB;DP=279;ExcessHet=3.0103;FS=6.425;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=15.14;ReadPosRankSum=-1.010;SOR=0.432       GT:AD:DP:GQ:PL  0/1:142,136:278:99:4238,0,4214

2. chr15 28518114 TC->T
   这个位点是在三个软件里都找到的位点且IGV验证成功

chr15 28518114 TC->T

3. chr1 120539687 T->C
   这个位点是在GATK和Sentieon中call到 但是没有在deepvariant中得到。
   这个位点的突变频率较低，IGV也不能确定这是一个可靠的突变。
   这个位点理应被展示出来。

chr1 120539687 T->C

4. chr19 8999535 ACT->A
   这个位点是在GATK和Sentieon中call到 但是没有在deepvariant中得到。
   从IGV可以看出，这应该是一个假突变。

chr19 8999535 ACT->A

5. chr1,116916118 A->C

这个位点是deepvariant中得到 但是没有在GATK和Sentieon中call到
深度很低，不能被认为是一个突变位点。

chr1,116916118 A->C

但是Deepvariant给出的判定也不是PASS，而是RefCall（Genotyping model thinks this site is reference）

$ grep 116916118  Deepvariant/samplename_N.deepvariant.vcf
chr1    116916118       .       A       C       0       RefCall .       GT:GQ:DP:AD:VAF:PL      0/0:38:8:6,2:0.25:0,37,58

6. chr11 3742058 G->T
   这个位点是deepvariant中得到 但是没有在GATK和Sentieon中call到
   这个位点的突变频率较低，IGV也不能确定这是一个可靠的突变。

chr11 3742058 G->T

Deepvariant给出的判定也不是PASS，而是RefCall（Genotyping model thinks this site is reference）

$ grep 3742058  Deepvariant/samplename_N.deepvariant.vcf
chr11   3742058 .       G       T       0       RefCall .       GT:GQ:DP:AD:VAF:PL      0/0:37:239:203,36:0.150628:0,38,41

位点太多，没办法一一去核对，要是有个矫正好的数据就好了~

五、总结

综上，可以看到，Deepvariant能够call到GATK和SENTIEON得到的大多数突变，并且call到很多其他突变（可能是假阳性或者深度过低无法确定的情况）

速度上可以了~看过一篇博客说Deepvariant速度比GATK要慢很多，不过我测试的结果反而更快呢。。。

自从2017年底后，Deepvariant就火热了一小阵，2018-2019年很少听到这个软件的进展了~看这个测试结果，Deepvariant得到的vcf再经过矫正，还是能够再精确的，而且速度也可以，是应该被推广的。