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](https://github.com/gringer/bioinfscripts/blob/master/sam2LongestBase.pl).