<p><a href="https://broadinstitute.github.io/picard/command-line-overview.html#CollectHsMetrics" rel="nofollow noreferrer">picard CollectHsMetrics</a> has the parameter <code>PER_BASE_COVERAGE</code> which take this into account.</p>

<p>You can run it like this:</p>

<pre><code> $ java -jar picard.jar CollectHsMetrics \
      I=input.bam \
      O=hs_metrics.txt \
      R=reference_sequence.fasta \
      BAIT_INTERVALS=bait.interval_list \
      TARGET_INTERVALS=target.interval_list \
      PER_BASE_COVERAGE=base_coverage.csv
</code></pre>

<p><code>BAIT_INTERVALS</code> and <code>TARGET_INTERVALS</code> could be the same file. You will need to convert your <code>bed</code> file to it using <a href="http://broadinstitute.github.io/picard/command-line-overview.html#BedToIntervalList" rel="nofollow noreferrer">picard BedToIntervalList</a>.</p>

<p>Another useful tool is <a href="https://github.com/brentp/mosdepth" rel="nofollow noreferrer">mosdepth</a>. </p>

<pre><code>$ mosdepth --by regions.bed SampleName input.bam
</code></pre>

<p>fin swimmer </p>

CC BY-SA 4.0
2018-11-24T14:23:54.463Z
Forgive me because I’be only looked into the tool, never used it. Do I have to generate some kind of sample.dict file for every sample I want to do this on or is it straightforward. It’s no big deal doing this on a reference but it would have too much overhead running it on every sample. Could you provide the steps involved in running this?
d_kennetz
CC BY-SA 4.0
2018-11-24T17:44:12.757Z
A `sample.dict` isn't needed. I edit my answer as @terdon asks for.
finswimmer
CC BY-SA 4.0
2018-11-24T12:59:09.717Z
2018-11-24T17:40:35.280Z
2018-11-24T17:40:35.280Z
finswimmer
<p>One way to deal with this would be to first merge paired-end reads based upon their overlapping regions, and then map them and calculate the coverage. This way you're only counting once per unique sequence.</p>

<p>Programs like <a href="https://github.com/jstjohn/SeqPrep" rel="nofollow noreferrer">SeqPrep</a>, <a href="https://sco.h-its.org/exelixis/web/software/pear/doc.html#cl-usage" rel="nofollow noreferrer">PEAR</a> (Paired-End reAd mergeR), and <a href="https://github.com/ExpressionAnalysis/ea-utils/blob/wiki/FastqJoin.md" rel="nofollow noreferrer">fastq-join</a> can do this fairly quickly. However, this isn't a perfect solution as the reads will not be joined unless they meet some minimum overlap criteria.</p>

CC BY-SA 4.0
2018-11-20T19:50:58.353Z
added an edit, I do not want to remap reads or make any modifications to existing reads.
d_kennetz
CC BY-SA 4.0
2018-11-20T19:48:20.707Z
2018-11-20T19:48:20.707Z
LucasBoatwright
<p>Bases from overlapping reads of the same fragment provide extra information, as their base qualities come from different sequencing cycles on the sequencer. This information can be used by downstream tools. For example, the merging tool of <a href="https://drive5.com/usearch/manual/merge_pair.html" rel="nofollow noreferrer">USEARCH</a> increases base quality scores of bases that match between overlapping reads. Depending on the context under which coverage information will be used in the future, it may be a good idea to "double-count" overlapping bases.</p>

CC BY-SA 4.0
2018-11-22T22:47:04.420Z
The reads aren't from different clusters in the traditional sense, if that were the case they wouldn't be assigned as mates. The traditional method to exploit overlapping reads is to increase/decrease phred scores on one of the mates according to whether the base calls in overlapping regions are the same or different, respectively. That's performed internally in things like `samtools mpileup`.
Devon Ryan
CC BY-SA 4.0
2018-11-23T07:36:06.447Z
I agree. The reads come from the same cluster.  However, you could match up mates from different clusters. The Illumina sequencer uses a fragment index to match up mates from the two mate sequencing runs. It also appears from your comment that samtools incorporates the extra information gained from having two reads cover the same base. Just so we're all on the same page, here is a video describing Illumina paired-end sequencing: https://youtu.be/fCd6B5HRaZ8
winni2k
CC BY-SA 4.0
2018-11-21T16:22:17.793Z
2018-12-11T18:16:36.957Z
2018-12-11T18:16:36.957Z
winni2k
<p><code>samtools depth</code> is a simplified version of <code>samtools mpileup</code>, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use <code>samtools mpileup</code> instead:</p>

<pre><code>$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  &lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  &lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;&lt;  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...
</code></pre>

<p>I haven't shown the entirety of the output and you'd want to use the <code>-aa</code> flag to ensure that all bases are returned.</p>

<p>Note that <code>samtools mpileup</code> is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to <code>-Q 1</code> (the default is <code>-Q 13</code>, which you'd want to change).</p>

<p>As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use <code>bamPEFragmentSize</code> from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.</p>

CC BY-SA 4.0
2018-11-23T09:03:57.810Z
Does this method take bases that are not covered by a read into account?
winni2k
CC BY-SA 4.0
2018-11-23T09:45:58.603Z
Yes, that's what the `-aa` flag I mentioned does.
Devon Ryan
CC BY-SA 4.0
2018-11-24T04:57:15.430Z
I tested this on data and it seems to handle the issue nicely! Thanks Devon.
d_kennetz
CC BY-SA 4.0
2018-11-22T23:09:03.470Z
2018-11-22T23:09:03.470Z
Devon Ryan
<p>EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region.</p>

<p>I used samtools depth to calculate the depth of coverage for samples in the whole Exome region using a GRCh37_ref.bed. These samples are sorted and duplicate marked. I ran this calculation on a few hundred samples to determine how much more sequencing needed to be done and found something interesting. Some samples had 100X coverage at 96,000,000 mapped reads and some samples had 100X coverage at 83,000,000 mapped reads. The samples that had 100X coverage at 83,000,000 reads had read pairs overlapping certain regions of the bedfile (read 1 was covering the same coordinates as read2 to some extent). </p>

<p>I then found out that samtools depth double counts these overlapped regions even though they are technically from the same molecule in sequencing and would be a duplicate bases in the read. See the image for an example:</p>

<p><a href="https://i.sstatic.net/7dkOV.png" rel="noreferrer"><img src="https://i.sstatic.net/7dkOV.png" alt="Duplicate Region"></a></p>

<p>Each base in this Duplicate Overlap region is double counted when it shouldn't be. How do you guys deal with this duplicate overlap region when calculating coverage across a bam file?</p>

<p>samtools (v1.7):</p>

<pre><code>samtools depth -a -b exome_targets_GRCh37.bed sample1.bam &gt; sample1_depth.txt
awk '{s+=$3}END{print s/NR}' sample1_depth.txt
</code></pre>

CC BY-SA 4.0
2018-11-13T07:27:46.277Z
There is already code in htslib for resolving overlapping paired-end reads, but it only seems to be used in samtools mpileup and not in samtools depth. I guess you could modify bam2depth.c and add the [bam_mplp_init_overlaps() function call](https://github.com/samtools/samtools/blob/02d93a154a3c27f1acea6021e9e4d75898c2586c/bam_plcmd.c#L566). I used it in my own [program](https://github.com/pmenzel/stranded-coverage) for strand-specific coverage calculation, see [here](https://github.com/pmenzel/stranded-coverage/blob/master/strand_cov.c#L69)
Peter Menzel
CC BY-SA 4.0
2018-11-13T16:48:42.607Z
That is interesting indeed and I will definitely look into the bam2depth.c modifications in order to drop overlap from the function call. Was it difficult to implement this into your program? I see your program and the overlapping ends counting feature appears to be what I am looking for. I will need to locally install samtools as I work on a cluster. Thanks for the recommendation!
d_kennetz
<samtools><coverage>
Double-counting coverage of overlapped read pairs
# Double-counting coverage of overlapped read pairs

(Asked by: username_0 on 2018-11-12)

EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region.
I used samtools depth to calculate the depth of coverage for samples in the whole Exome region using a GRCh37_ref.bed. These samples are sorted and duplicate marked. I ran this calculation on a few hundred samples to determine how much more sequencing needed to be done and found something interesting. Some samples had 100X coverage at 96,000,000 mapped reads and some samples had 100X coverage at 83,000,000 mapped reads. The samples that had 100X coverage at 83,000,000 reads had read pairs overlapping certain regions of the bedfile (read 1 was covering the same coordinates as read2 to some extent). 
I then found out that samtools depth double counts these overlapped regions even though they are technically from the same molecule in sequencing and would be a duplicate bases in the read. See the image for an example:

Each base in this Duplicate Overlap region is double counted when it shouldn't be. How do you guys deal with this duplicate overlap region when calculating coverage across a bam file?
samtools (v1.7):
samtools depth -a -b exome_targets_GRCh37.bed sample1.bam > sample1_depth.txt
awk '{s+=$3}END{print s/NR}' sample1_depth.txt

**Tags:** samtools, coverage
> **Comments:**
> **username_1**: There is already code in htslib for resolving overlapping paired-end reads, but it only seems to be used in samtools mpileup and not in samtools depth. I guess you could modify bam2depth.c and add the [bam_mplp_init_overlaps() function call](https://github.com/samtools/samtools/blob/02d93a154a3c27f1acea6021e9e4d75898c2586c/bam_plcmd.c#L566). I used it in my own [program](https://github.com/pmenzel/stranded-coverage) for strand-specific coverage calculation, see [here](https://github.com/pmenzel/stranded-coverage/blob/master/strand_cov.c#L69)  
> **username_0**: That is interesting indeed and I will definitely look into the bam2depth.c modifications in order to drop overlap from the function call. Was it difficult to implement this into your program? I see your program and the overlapping ends counting feature appears to be what I am looking for. I will need to locally install samtools as I work on a cluster. Thanks for the recommendation!  

---

## 4 answers

### Answer by username_5 on 2018-11-22. Score: 5 | Accepted answer

samtools depth is a simplified version of samtools mpileup, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use samtools mpileup instead:
$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...

I haven't shown the entirety of the output and you'd want to use the -aa flag to ensure that all bases are returned.
Note that samtools mpileup is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to -Q 1 (the default is -Q 13, which you'd want to change).
As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use bamPEFragmentSize from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.

> **Comments:**
> **username_4**: Does this method take bases that are not covered by a read into account?  
> **username_5**: Yes, that's what the `-aa` flag I mentioned does.  
> **username_0**: I tested this on data and it seems to handle the issue nicely! Thanks Devon.  

---

### Answer by username_2 on 2018-11-24. Score: 2

picard CollectHsMetrics has the parameter PER_BASE_COVERAGE which take this into account.
You can run it like this:
 $ java -jar picard.jar CollectHsMetrics \
      I=input.bam \
      O=hs_metrics.txt \
      R=reference_sequence.fasta \
      BAIT_INTERVALS=bait.interval_list \
      TARGET_INTERVALS=target.interval_list \
      PER_BASE_COVERAGE=base_coverage.csv

BAIT_INTERVALS and TARGET_INTERVALS could be the same file. You will need to convert your bed file to it using picard BedToIntervalList.
Another useful tool is mosdepth. 
$ mosdepth --by regions.bed SampleName input.bam

fin swimmer

> **Comments:**
> **username_0**: Forgive me because I’be only looked into the tool, never used it. Do I have to generate some kind of sample.dict file for every sample I want to do this on or is it straightforward. It’s no big deal doing this on a reference but it would have too much overhead running it on every sample. Could you provide the steps involved in running this?  
> **username_2**: A `sample.dict` isn't needed. I edit my answer as @terdon asks for.  

---

### Answer by username_3 on 2018-11-20. Score: 2

One way to deal with this would be to first merge paired-end reads based upon their overlapping regions, and then map them and calculate the coverage. This way you're only counting once per unique sequence.

Programs like [SeqPrep](https://github.com/jstjohn/SeqPrep), [PEAR](https://sco.h-its.org/exelixis/web/software/pear/doc.html#cl-usage) (Paired-End reAd mergeR), and [fastq-join](https://github.com/ExpressionAnalysis/ea-utils/blob/wiki/FastqJoin.md) can do this fairly quickly. However, this isn't a perfect solution as the reads will not be joined unless they meet some minimum overlap criteria.

> **Comments:**
> **username_0**: added an edit, I do not want to remap reads or make any modifications to existing reads.  

---

### Answer by username_4 on 2018-11-21. Score: 1

Bases from overlapping reads of the same fragment provide extra information, as their base qualities come from different sequencing cycles on the sequencer. This information can be used by downstream tools. For example, the merging tool of [USEARCH](https://drive5.com/usearch/manual/merge_pair.html) increases base quality scores of bases that match between overlapping reads. Depending on the context under which coverage information will be used in the future, it may be a good idea to "double-count" overlapping bases.

> **Comments:**
> **username_5**: The reads aren't from different clusters in the traditional sense, if that were the case they wouldn't be assigned as mates. The traditional method to exploit overlapping reads is to increase/decrease phred scores on one of the mates according to whether the base calls in overlapping regions are the same or different, respectively. That's performed internally in things like `samtools mpileup`.  
> **username_4**: I agree. The reads come from the same cluster.  However, you could match up mates from different clusters. The Illumina sequencer uses a fragment index to match up mates from the two mate sequencing runs. It also appears from your comment that samtools incorporates the extra information gained from having two reads cover the same base. Just so we're all on the same page, here is a video describing Illumina paired-end sequencing: https://youtu.be/fCd6B5HRaZ8  

---

