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

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
  • 1
    $\begingroup$ 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. I used it in my own program for strand-specific coverage calculation, see here $\endgroup$ – Peter Menzel Nov 13 '18 at 7:27
  • $\begingroup$ 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! $\endgroup$ – d_kennetz Nov 13 '18 at 16:48

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.

  • $\begingroup$ Does this method take bases that are not covered by a read into account? $\endgroup$ – winni2k Nov 23 '18 at 9:03
  • 1
    $\begingroup$ Yes, that's what the -aa flag I mentioned does. $\endgroup$ – Devon Ryan Nov 23 '18 at 9:45
  • $\begingroup$ I tested this on data and it seems to handle the issue nicely! Thanks Devon. $\endgroup$ – d_kennetz Nov 24 '18 at 4:57

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, PEAR (Paired-End reAd mergeR), and fastq-join 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.

  • 1
    $\begingroup$ added an edit, I do not want to remap reads or make any modifications to existing reads. $\endgroup$ – d_kennetz Nov 20 '18 at 19:50

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 \

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

  • $\begingroup$ 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? $\endgroup$ – d_kennetz Nov 24 '18 at 14:23
  • 1
    $\begingroup$ A sample.dict isn't needed. I edit my answer as @terdon asks for. $\endgroup$ – finswimmer Nov 24 '18 at 17:44

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 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.

  • $\begingroup$ 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. $\endgroup$ – Devon Ryan Nov 22 '18 at 22:47
  • $\begingroup$ 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: youtu.be/fCd6B5HRaZ8 $\endgroup$ – winni2k Nov 23 '18 at 7:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.