5

I have very noisy nanopore reads and am trying to call SNPs/indels.

I'm runinig into some trouble when truing to use the samtools mpileup | bcftools call combination.

It seems to incorrectly be calling very long indels where it looks like there is no support. For example:

samtools mpileup -f ref.fa sample.bam -r Chromosome:198940-198940 produces:

Chromosome    198940    C    37    .,.-1A.-1A...-1A.-1A.-1A.,-1a.-1A.,-1a,-1a.,.-1A,,.-1A,,+1a,+1a,.-1A,,-1a.,-1a,-1a,,-1a,-1a.-1A,-1a.    06?5/<:8>/4<?3691465<3354.08:23:11@3C

This looks fine... however when I try to call the genotype with bcftools I get a very long indel. samtools mpileup -ugf ref.fa sample.bam -r Chromosome:198940-198940 -t AD | bcftools call -mv produces:

Chromosome    198940    .    CAAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA    C    225    .    INDEL;IDV=32;IMF=0.5;DP=62;VDB=0.145385;SGB=-0.662043;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,7,2;MQ=60    GT:PL:AD    1/1:255,27,0:0,9

Using htsbox pileup -f ref.fa sample.bam | grep 198940 it shows there is one read with this indels together with many other reads with smaller indels.

Chromosome  198940  C   C,C-75AAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA,C-4AAAA,C-2AA,C-1A,T,C+1A,C+1T,C+3GGT 0/4:25,1,1,2,24,5,2,1,1

Is there any reason why it would be grouping these reads together into the big indel?

Any idea what is going wrong here?

gringer
  • 14,012
  • 5
  • 23
  • 79
user3715117
  • 53
  • 3

1 Answers1

1

It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021.

What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the pileup calls for the surrounding bases may help to diagnose this.

My general advice to people doing variant calling with nanopore is to be suspicious of deletions, especially if they include homopolymers. That usually means trying to confirm any important variants with another technology.

Newer, more accurate, base calling models are being created all the time by ONT, so if you still have the raw FAST5 data, it might be worth it to run through the most recent base caller.

gringer
  • 14,012
  • 5
  • 23
  • 79