samtools vs. bedtools

[bair]

Hi all,

I use samtools (depth) and bedtools (coverageBed -d) to calculate coverage for a given bedfile, the results are different

In the following dataset, the first three columns are generated by samtools, the rest are generated by bedtools. At most of the positions (>65%), cov1=cov2, at some positions the differences are huge. For example, at chr1:237433808, cov1 =522 (samtools), while cov2 = 649 (bedtools).

Ant explanation?

chrs posi cov1 start ends count cov2
chr1 237433798 647 237433797 237433916 1 649
chr1 237433799 647 237433797 237433916 2 649
chr1 237433800 647 237433797 237433916 3 648
chr1 237433808 522 237433797 237433916 11 649
chr1 237433826 653 237433797 237433916 29 654
chr1 237433828 653 237433797 237433916 31 654
chr1 237433829 654 237433797 237433916 32 655
chr1 237433830 654 237433797 237433916 33 655
chr1 237433831 654 237433797 237433916 34 655
chr1 237433833 654 237433797 237433916 36 655

[arvid]

While you mention "a given bedfile", I suppose you use a BAM as the original input (possibly converted for coverageBed) in both cases, right? Since you seem to have consistently higher counts with BEDTools, this could be due a different definition of coverage in the case of splices (or other non-matching CIGAR operations) of read alignments, not properly carried over into the bed file (covered with a bed entry although they should be excluded for a coverage calculation).

I had a similar issue with RNA-Seq splice-mapped reads and BEDTools, and could fix it by using the "-split" argument with the bamToBed converter (splice reads are then splitted into two bed entries) - but if it is due to other non-matching CIGAR operations you might need a closer look at the read alignments themselves at those differing positions.

[bair]

Thank you arvid,

I tried bamToBed with -split, no changes to the result, since my bedfile only contains CCDS regions.

coverageBed reports all the bases covered at the position, while samtools/depth may have some default CIGAR operation as filter, I need to look closer.

Thanks again