Clarification on Bedtools Intersect Usage

Question

[Noam Teyssier]

So I am recently getting my foot in the door in using bedtools and I had a question about correct usage for my end goal

I am looking to essentially filter a BAM file for reads that only have 100% intersection on regions specified with a bed file.

  [====================BAM===================] (A)

                (-----------------)                          (-----------------------)      (B)

I would only want to pull down reads in A that have complete overlap of reads in B.

The command I am using is:

bedtools intersect -F 1.0 -a file.bam -b regions.bed

However when I run

samtools mpileup -l regions.bed file.bam

I see that there are inconsistencies in the depth at each of the positions in the regions. This doesn't quite make sense to me intuitively because I would expect an even depth across each position in the region because the intersect should filter out reads with incomplete coverage across the region.

Perhaps I am thinking of the usages of these tools wrong though and could definitely use some help in clarifying either inconsistencies in my intuition or commands that don't match my intentions.

Cheers!

Answer 1

[Kevin Blighe]

samtools mpileup by default only counts bases with a base quality (BQ) ≥ 13. It also sets a 'ceiling' of maximum depth to 250 at each position. There are even more default thresholds set and you should take a look at them. bedtools intersect will not make any such assumption about BQ or maximum depth.

A better comparison would be bedtools intersect (or bedtools coverage -d) with samtools depth.

Kevin