GT calls where DP=0

[PfeiferLab]

Hi @dancooke,

We have been calling variants in the trio mode using three different settings - each with different results - and are wondering what to make out of GT calls where the DP=0.

Below is an example for chromosome 2 position 168,349,309.

Setting 1 using octopus v0.7.4 (installed via conda):

octopus 
-R reference.fa 
-I sire.bam dam.bam offspring.bam 
-M dam 
-F sire 
--forest germline.v0.7.4.forest 
-o triovcf.gz

In this setting, the position is called as fixed alternative in the trio - however, please note that individual 1 (GT 1/1) has a DP of 0. chr2 168349309 . T C 5508.5 PASS AC=6;AN=6;DP=42;MP=81.61;MQ=60;NS=2;RFGQ_ALL=12.6 GT:GQ:DP:MQ:PS:PQ:RFGQ:FT 1|1:10:0:0:168348962:80:8.56:PASS 1|1:10:15:60:168348962:94:16.99:PASS 1|1:96:27:60:168348962:96:100:PASS

Setting 2 using octopus v0.7.4 (installed via conda) with "--assemble-all" and "--denovos-only":

octopus 
-R reference.fa 
-I sire.bam dam.bam offspring.bam 
-M dam 
-F sire 
--forest germline.v0.7.4.forest 
-o triovcf.gz
--assemble-all yes
--denovos-only yes

In this setting, the position is called as DENOVO with individual 1 now having GT 0|1 - however, still with a DP of 0. chr2 168349309 . T C 37.17 PASS AC=5;AN=6;DENOVO;DP=42;MP=86.56;MQ=60;NS=2;PP=37.17;REVERSION;RFGQ_ALL=5.99 GT:GQ:DP:MQ:PS:PQ:RFGQ:FT 0|1:37:0:0:168348962:80:3.18:PASS 1|1:55:15:60:168348962:100:8.24:PASS 1|1:135:27:60:168348962:100:13.01:PASS

Setting 3 using octopus v0.7.4 with patch (https://github.com/luntergroup/octopus/commit/09ebd28945026556e77d73f1dcd8f0212265183c) and "--assemble-all" and "--denovos-only"

octopus 
-R reference.fa 
-I sire.bam dam.bam offspring.bam 
-M dam 
-F sire 
--forest germline.v0.7.4.forest 
-o triovcf.gz
--assemble-all yes
--denovos-only yes

Similar to setting 2, the position is still called as DENOVO despite individual 1 having DP=0 (which the patch was meant to fix): chr2 168349309 . T C 70.35 PASS AC=5;AN=6;DENOVO;DP=42;MP=121.34;MQ=60;NS=2;PP=52.13;REVERSION;RFAQ_ALL=6.36 GT:GQ:DP:MQ:PS:PQ:RFAQ:RFGQ:FT 0|1:68:0:0:168348962:100:3.19:.:PASS 1|1:52:15:60:168348962:100:10.46:.:PASS 1|1:139:27:60:168348962:100:13.98:.:PASS

Here is an IGV screenshot:

[dancooke]

Hi @PfeiferLab, thanks for the report. That does look a little strange - also the QUAL decreasing considerably. Could you provide BAMs for the trio around this region?

[PfeiferLab]

@dancooke: Sure, for setting 3?

[dancooke]

I assumed it was the same set of BAMs for all three settings?

[PfeiferLab]

Yes, the input BAMs are the same though we also have the re-assembled/re-aligned BAMs (--bamout) for each setting in case they might be helpful.

[dancooke]

Ah. Just the input BAMs (subsetted) please.

[PfeiferLab]

Thank you for looking into this! I emailed you the subsetted input bams for the example above.

[dancooke]

Hey @PfeiferLab, I no longer have access to my .well.ox.ac.uk email. Please could you re-send to daniel.cooke@me.com?

[PfeiferLab]

Sure, re-sent.

[dancooke]

Thanks for sending the data and sorry it's taken so long to get back to you. The problem is that there's a 1521bp deletion spanning this region that isn't getting called because it doesn't get proposed as a candidate with the settings you've tried - both genotype calls are incorrect. We can rescue the deletion by increasing the assembly window size:

$ octopus \
  -R papAnu4.fa \
  -I 10173_chr2_168349309_1kb.bam \
     1X1765_chr2_168349309_1kb.bam \
     1X2124_chr2_168349309_1kb.bam \
  -T chr2:168,347,848-168,351,135 \
  -M 1X1765 -F 1X2124 \
  --assemble-all \
  --max-assembly-region-size 2000  \
  --disable-call-filter \
  -o octopus.vcf \
  --bamout realigned

Sure enough, the deletion gets called and the correct genotype is called at chr2:168349309:

chr2    168348545   .   AAAGTTTTATATATATATATAAAAAATATATATATGTATGACTGAATCATTGCAATCTCACAATAAAACTTGAATAGATGAGGAGTTGCTTCTTCTGGATGAGCAAAGAAATAGGTATATATATATATATATATATATATATATATATATATAGTATTTATATATAGTGAATTAGGCCACCTCTCATCACCTCCACTGAAGTTATAGTTTTGAAAAGTGTTTACCAGTTCTTTAACATTCCTCCCTTCAGGAAGTGAATCTTAATTATCCTTCTCTTATGTGTGGATGAGACTTGGCAACTCATGGCTAATGAACTGAATTATAGCGGAACTGATGAGTTGCTACTTCCAAGGTTAGGTTATAAAGAGACTGCGCCTTCCAATTGTTTTGCTCTCTCTTGCTCACTGACTCTCTCTCAGAGCCCAGCCACCATGTAGAATAGTAAAGGTTATTCTGCTGAAGAGGTTACATAGAGAAGAGCAGCCCTCAAGGATAAAAGGCCAAATGTAAGAGAAATGAGGTGCTCACTCCAGCCAAAAGCCAAAACTCACTGCCAGCTATGGTAGTGAGTCTTCTTGGAAGCTGACCCTCCAGCCCCAGTTCAGGCTTCAGATGATGCATCCCTGGCTGACAGTTTGACTACAGCCTCATGAGAAACCCTGAGTCAGAATCACCCAGATAAACTTCTCCCGGACTCCTGACCCTCAAAAACTATAAGAAAATAAATGTTTGTTATTTCAACCCTCCAGATTTCTGGGGAATTTTTACCCAGCAATGGATAATGAATACAACTGCTATCATCCTAGTCCATGCTGTTGTGACCTTTCACTGAGGCTACTGAGCTATTCCTCCTGAACAGACTCCCTCCTTCCATCTCTCCTCATGAAATATAATTCCCCACTTCAGATCCCAAATGATCTTTTAAATATCTAAAGTGGCCTTCACGATTGTCCATTTACTCATTTGTCACTACTTACAAGGTTAAAGTTCTACAAAGACTTTTGATAAATACTAAGAATAAAAGCCAAACTACTACCAGGTCTATACGACCGTACATGATCTGGGCATTCAGAATGACTTTTCTAGTACATACTATTTTAAAAAATGCCCTCTCCTTCAATATTCTTCTATTATATTTATTTAATTCTTAACATATTATCTGAAATTACTTTATGTATTATTTAGTATTTTTTCAATCATTTAGGGTCTGTCATTTCTAACTAGTATGAAAGATCCATGAGAGGACTTTGTATTGCTTATCCCTGTAGACTGTTGGATACATATCTTTTAATACATGTATATGTACATATATATGCATAACTATATGTATATCCATCCACCTCTCTGTCTATCTATCCATCCATGTATCTATCTGTCTATCTATCTGTCTGTCTAACCTCTCCAACTAGTAGTACTACTCACATATAGGACATATTAGAAGACTAAATCATATAAGGTGGTGTTAACAACTATTCACAAACTAGATGACTTGCTATCTTTCCTTTCCCCTTTTAGGTACCGG  A,* 29707.2 .   AC=4,2;AN=6;DP=123;MQ=59;NS=3   GT:GQ:DP:MQ:PS:PQ:MP    1|1:289:46:60:168348388:100:68.3    2|1:569:33:59:168348388:100:0.00016 1|2:167:44:59:168348388:100:0
chr2    168349309   .   T   C,* 21064.6 .   AC=2,4;AN=6;DP=42;MQ=60;NS=2    GT:GQ:DP:MQ:PS:PQ:MP    2|2:289:0:0:168348388:100:68.3  1|2:607:15:60:168348388:100:0.00016 2|1:819:27:60:168348388:100:0

Here are the realigned reads: