<p>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 <em>mapped</em> reads. This method reduces the genome coverage, but preserves reads where coverage is low.</p>

<p>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:</p>

<pre><code>  if(($F[2] ne $seqName) || ($F[3] != $pos) || (length($bestSeq) &lt;= 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];
</code></pre>

<p>The full code (which works as a filter on uncompressed SAM files) can be found <a href="https://github.com/gringer/bioinfscripts/blob/master/sam2LongestBase.pl" rel="nofollow noreferrer">here</a>.</p>

CC BY-SA 3.0
2017-06-05T09:44:15.057Z
To clarify: this is Perl code?
Konrad Rudolph
CC BY-SA 3.0
2017-06-05T10:27:19.797Z
yes, Perl code. I can tell by the dollar symbols everywhere, but I guess it's not obvious to people who are not familiar with Perl.
gringer
CC BY-SA 3.0
2017-06-05T10:37:09.437Z
Well I do speak Perl but it’s hardly the only language with dollars before identifiers. ;-)
Konrad Rudolph
CC BY-SA 3.0
2017-06-05T10:39:38.757Z
Indeed. My mind was blank, but my wife reminded me of bash/sh.
gringer
CC BY-SA 3.0
2017-06-04T23:34:21.633Z
2017-06-05T10:26:20.287Z
2017-06-05T10:26:20.287Z
gringer
<p>The <code>samtools mpileup</code> command has quite a neat feature that it is able to correct mapping errors associated with misalignment of INDELs. By default, the <code>mpileup</code> 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.</p>

<p>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).</p>

<p>I see that @rightskewed has mentioned the downsampling ability of samtools with <code>samtools view -s &lt;float&gt;</code> (see <a href="https://bioinformatics.stackexchange.com/a/406/73">here</a>), which seems like it might work as a solution for this if used prior to the mpileup operation.</p>

<sam><samtools><variant-calling><indel>
How can I speed up INDEL calling/correction on BAM files?
# How can I speed up INDEL calling/correction on BAM files?

(Asked by: username_0 on 2017-06-04)

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`  (see [here](https://bioinformatics.stackexchange.com/a/406/73)), which seems like it might work as a solution for this if used prior to the mpileup operation.

**Tags:** sam, samtools, variant-calling, indel

---

## 1 answers

### Answer by username_0 on 2017-06-04. Score: 4 | Accepted answer

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.

> **Comments:**
> **username_1**: To clarify: this is Perl code?  
> **username_0**: yes, Perl code. I can tell by the dollar symbols everywhere, but I guess it's not obvious to people who are not familiar with Perl.  
> **username_1**: Well I do speak Perl but it’s hardly the only language with dollars before identifiers. ;-)  
> **username_0**: Indeed. My mind was blank, but my wife reminded me of bash/sh.  

---

