Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-24T14:23:54.463Z
Answers.Text: 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?
Answers.UserDisplayName: d_kennetz
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-24T17:44:12.757Z
Answers.Text: A `sample.dict` isn't needed. I edit my answer as @terdon asks for.
Answers.UserDisplayName: finswimmer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-24T12:59:09.717Z
Answers.LastActivityDate: 2018-11-24T17:40:35.280Z
Answers.LastEditDate: 2018-11-24T17:40:35.280Z
Answers.UserDisplayName: finswimmer
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-20T19:50:58.353Z
Answers.Text: added an edit, I do not want to remap reads or make any modifications to existing reads.
Answers.UserDisplayName: d_kennetz
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-20T19:48:20.707Z
Answers.LastActivityDate: 2018-11-20T19:48:20.707Z
Answers.UserDisplayName: LucasBoatwright
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-22T22:47:04.420Z
Answers.Text: 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`.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-23T07:36:06.447Z
Answers.Text: 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
Answers.UserDisplayName: winni2k
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-21T16:22:17.793Z
Answers.LastActivityDate: 2018-12-11T18:16:36.957Z
Answers.LastEditDate: 2018-12-11T18:16:36.957Z
Answers.UserDisplayName: winni2k
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-23T09:03:57.810Z
Answers.Text: Does this method take bases that are not covered by a read into account?
Answers.UserDisplayName: winni2k
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-23T09:45:58.603Z
Answers.Text: Yes, that's what the `-aa` flag I mentioned does.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-24T04:57:15.430Z
Answers.Text: I tested this on data and it seems to handle the issue nicely! Thanks Devon.
Answers.UserDisplayName: d_kennetz
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-22T23:09:03.470Z
Answers.LastActivityDate: 2018-11-22T23:09:03.470Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-13T07:27:46.277Z
Comments.Text: 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)
Comments.UserDisplayName: Peter Menzel
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-13T16:48:42.607Z
Comments.Text: 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!
Comments.UserDisplayName: d_kennetz
Tags: <samtools><coverage>
Title: Double-counting coverage of overlapped read pairs
ThreadText: # 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:
  
  [![Duplicate Region](https://i.sstatic.net/7dkOV.png)](https://i.sstatic.net/7dkOV.png)
  
  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 username_5.  
  
  ---
  
  ### Answer by username_2 on 2018-11-24. Score: 2
  
  [picard CollectHsMetrics](https://broadinstitute.github.io/picard/command-line-overview.html#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](http://broadinstitute.github.io/picard/command-line-overview.html#BedToIntervalList).
  
  Another useful tool is [mosdepth](https://github.com/brentp/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  
  
  ---
  
  