Colleagues,

I have some BAM files with excessive coverage in certain regions that I want to cap to a given coverage. In regions with lower coverage, I do not want the coverage to change. Here's an illustration of what I want:

<img alt="" src="http://i.imgur.com/YQ7jsNR.jpg" style="height:450px; width:600px" />

I want read pairs to be intact in the output file.

Does anyone have a solution to this? I want to output a BAM file for other tools do use. Google has not been my friend.

(note: GATK PrintReads using `-dcov` looked promising until a realised that "For read-based traversals (ReadWalkers like BaseRecalibrator), it controls the maximum number of reads sharing the same alignment start position" from https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_engine_CommandLineGATK.php#--downsample_to_coverage. Close but no CIGAR (pun intended)).
I wrote two tools:

 - The first sorts a BAM on REF/read-name https://github.com/lindenb/jvarkit/wiki/SortSamRefName
 - The second reads the BAM above. For each set of reads having the same name, it checks if after adding those reads, the depth would be lower than the expected depth. If true, those reads are added to the BAM. **Edit**: doc is https://github.com/lindenb/jvarkit/wiki/Biostar154220

```
$ java -jar dist/sortsamrefname.jar -T bam /input.bam  |\
  java -jar dist/biostar154220.jar  -n 20 -T bam |\
  samtools sort - output
```

```
$ samtools mpileup output.bam  | cut -f 4 | sort | uniq -c
  12692 0
 596893 1
  94956 10
  56715 11
  76947 12
  57912 13
  66585 14
  51961 15
  63184 16
  47360 17
  65189 18
  65014 19
 364524 2
 169064 20
  72078 3
 118288 4
  54802 5
  82555 6
  53175 7
  78474 8
  54052 9
```

~~(what are those 21 and 22 ? I suppose it's because I don't consider supplementary, secondary, duplicate reads )~~ *(bug fixed)*. Anyway, can you tell me if that's what you need?

Pierre