# Double-counting coverage of overlapped read pairs

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  • 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 Nov 13, 2018 at 7:27 • 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! Nov 13, 2018 at 16:48 ## 4 Answers 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. • Does this method take bases that are not covered by a read into account? Nov 23, 2018 at 9:03 • Yes, that's what the -aa flag I mentioned does. Nov 23, 2018 at 9:45 • I tested this on data and it seems to handle the issue nicely! Thanks Devon. Nov 24, 2018 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. • added an edit, I do not want to remap reads or make any modifications to existing reads. Nov 20, 2018 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 \
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

• 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? Nov 24, 2018 at 14:23
• A sample.dict isn't needed. I edit my answer as @terdon asks for. Nov 24, 2018 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.

• 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. Nov 22, 2018 at 22:47
• 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 Nov 23, 2018 at 7:36