trouble with bedtools intersect and coverage

Question

[bio1234567]

Hello,

I encountered a problem using bedtools that I can't seem to solve myself. I am very new to bioinformatics and informatics in general, so I apologize if the answer is obvious.

My bedtools version is v2.28.0. I am trying to extract coverage data for a set of genomic regions (stored in a regions.bed-file) from a BAM file containing whole-exome data.

I started out trying the following command: bedtools coverage -a exome.bam -b regions.bed > output

However, the output always contains the coverage data for all regions of the BAM file. The regions.bed-file appears to be ignored. I have already made sure the chromosomes are named consistently in both files (chr1).

As I couldn't get this to work, I tried to create an intersection of the two files and use it as input for the coverage command. The command I used was bedtools intersect -wa -a exome.bam -b regions.bed > output The output, however, is illegible (like what you see when opening an archived file with cat).

Clearly, I am doing something wrong. Can you help me? Thanks a lot in advance!

Answer 1

[Vitis]

As of bedtools version 2.27.1, the information for "bedtools coverage" is the following:

Tool:    bedtools coverage (aka coverageBed)
Version: v2.27.1
Summary: Returns the depth and breadth of coverage of features from B
     on the intervals in A.

Usage:   bedtools coverage [OPTIONS] -a <bed/gff/vcf> -b <bed/gff/vcf>

So the region/feature file should come after "-a", while the BAM file should come after "-b". It seems that you got confused about these two parameters. I remember at some point, bedtools switched the "-a" and "-b" options and you might still remember the old scheme.

[bio1234567]

Hello Vitis,

thanks a lot for your reply! I'm afraid that is not what caused the problem though. I tried switching - a and -b but received an error message indicating that's just the wrong way around. The following note on the bedtools page confirms this, if I understand it correctly: "As of version 2.24.0, the coverage tool has changed such that the coverage is computed for the A file, not the B file. This changes the command line interface to be consistent with the other tools. Also, the coverage tool can accept multiple files for the -b option." Do you have any other ideas? :( My original command does give me a coverage output for my BAM file, it is just not restricted to the regions defined in my bed-file.

[Kevin Blighe]

How did you align the reads that are contained in the BAM? Have you checked that the aligner has formatted them correctly in SAM format?

[bio1234567]

Hi Kevin, I think I figured it out. Apparently, there was some problem with the bed-file, even though I'm not sure exactly what was wrong. I downloaded a new file, and now it seems to be working. Also, Vitis was right about the meaning of -a and -b. I still get an error message but a the output seems to be correct. Thanks you two for your help!