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.