Question: trouble with bedtools intersect and coverage
gravatar for bio1234567
17 months ago by
bio123456710 wrote:


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!

ADD COMMENTlink modified 17 months ago by Vitis2.4k • written 17 months ago by bio123456710
gravatar for Vitis
17 months ago by
New York
Vitis2.4k wrote:

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.

ADD COMMENTlink written 17 months ago by Vitis2.4k

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.

ADD REPLYlink written 17 months ago by bio123456710

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?

ADD REPLYlink written 17 months ago by Kevin Blighe67k

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!

ADD REPLYlink written 17 months ago by bio123456710
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1387 users visited in the last hour