How can I speed up INDEL calling/correction on BAM files?

Question

The samtools mpileup command has quite a neat feature that it is able to correct mapping errors associated with misalignment of INDELs. By default, the mpileup command will not work for reads that have more than 250X coverage of the reference genome. While this limit can be increased, very high coverage causes the mpileup program to grind to a halt, so it'd be nice to know if there's some easy way to make that faster.

To add a bit more context, I've been doing this with mitochondrial genome reads that were extracted from both Illumina whole-genome sequencing (coverage ~1000X), and from targeted amplicon sequencing done on the IonTorrent (coverage up to ~4000X).

I see that @rightskewed has mentioned the downsampling ability of samtools with samtools view -s <float> (see here), which seems like it might work as a solution for this if used prior to the mpileup operation.

Answer 1

I wasn't aware of the samtools subsampling when I had this problem a couple of years ago, so ended up writing my own digital normalisation method to deal with mapped reads. This method reduces the genome coverage, but preserves reads where coverage is low.

Because I was working with IonTorrent reads (which have variable length), I came up with the idea of selecting the longest read that mapped to each location in the genome (assuming such a read existed). This meant that the highly variable coverage for different samples (sometimes as low as 200X, sometimes as high as 4000X) was flattened out to a much more consistent coverage of about 100-200X. Here's the core of the Perl code that I wrote:

  if(($F[2] ne $seqName) || ($F[3] != $pos) || (length($bestSeq) <= length($F[9]))){
    if(length($bestSeq) == length($F[9])){
      ## reservoir sampling with a reservoir size of 1
      ## See https://en.wikipedia.org/wiki/Reservoir_sampling
      ## * with probability 1/i, keep the new item instead of the current item
      $seenCount++;
      if(!rand($seenCount)){
        ## i.e. if rand($seenCount) == 0, then continue with replacement
        next;
      }
    } else {
      $seenCount = 1;
    }
    if(($F[2] ne $seqName) || ($F[3] != $pos)){
      if($output eq "fastq"){
        printSeq($bestID, $bestSeq, $bestQual);
      } elsif($output eq "sam"){
        print($bestLine);
      }
    }
    $seqName = $F[2];
    $pos = $F[3];
    $bestLine = $line;
    $bestID = $F[0];
    $bestFlags = $F[1];
    $bestSeq = $F[9];
    $bestQual = $F[10];

The full code (which works as a filter on uncompressed SAM files) can be found here.