This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Problem with Bedtools

Hello everyone.

I'm having problems with bedtools coverage I have a .gff file with positions and a bam with the results from a RNA-seq experiment. I want to know how many reads has each of my positions however is not working.

I have used the following codes:

bedtools coverage -a.gff -b .bam
bedtools coverage -a.gff -b .bam -s

And it is not working :/

rna-seq
• 5,666 views • link
updated by Jeffin Rockey Jeffin Rockey • written by alvarocentron91 • alvarocentron91

Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time.
code_formatting

• 0 views • link
written by Ram Ram 
bedtools coverage -a.gff -b .bam
bedtools coverage -a.gff -b .bam -s

I run the code and nothing happens, however when I do:

bedtools coverage -abam .bam -b .gff -s

screenshot: https://ibb.co/jcpTtd

I get results but not the ones I want. It generates the coverage of features in the .gff on the features in the .bam and I want the opposite

• 0 views • link
written by alvarocentron91 • alvarocentron91

bam files do not have features, they have reads / sequences mapped.

I get results but not the ones I want. It generates the coverage of features in the .gff on the features in the .bam and I want the opposite

I am not sure what you want, could you explain? Maybe bedtools intersect can do what you want?

• 0 views • link
written by h.mon h.mon

Ok I want to know how many reads from the .bam overlap with each of the positions I have on my .gff.

I tried bedtools intersect as well, and it is the same.

May it work if I convert the bam to bed?

Thanks in advance

• 0 views • link
written by alvarocentron91 • alvarocentron91

Not working is a tremendously detailed description. Please elaborate.

• 0 views • link
written by ATpoint ATpoint

2 answers

And it is not working :/

"It is not working" is not really helpful, copy the error message, or, if no error message, explicitly say "there is no error message and bedtools produces an empty file", or something to this effect.

Anyway, I believe your command should be:

bedtools coverage -a file.bam -b file.gff

As per documentation:

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.

edit:

You are right, the command I gave above is the opposite of what you want. You really want 

bedtools coverage -a file.gff -b file.bam

Regarding the error you got (slurmstepd: error: task/cgroup: unable to add task[pid=XXXXXXX] to memory cg '(null)), it is related to excessive memory use. You can ask for more memory or try different queue. If you sort both bam and gff, you can pass the argument -sorted to bedtools, which is both faster and uses less memory.

May it work if I convert the bam to bed?

Yes, should work both with bam or bed.

• 0 views • link
written by h.mon h.mon

But that line generates just the opposite of what I want

• 0 views • link
written by alvarocentron91 • alvarocentron91

The problem with:

bedtools coverage -a file.gff -b file.bam

is that it is not running, it just keep "thinking" until the computer kills it

• 0 views • link
written by alvarocentron91 • alvarocentron91

  • You are using the Slurm Workload Manager, which controls how much resources (mainly memory, CPUs, time) you can use.

  • The error message you posted on that image relates to Slurm killing you job because it used too much resources.

  • I pointed to the fact if you sort the gff and bam file, memory use should be lower and running times should be faster, and you might successfully run bedtools coverage without Slurm killing it.

• 0 views • link
written by h.mon h.mon

Ahh ok, sorry either my english and bioinformatics skills are low. Ok I understand I will sort it and I will ask for more resources thank you!

• 0 views • link
written by alvarocentron91 • alvarocentron91

I finally got it, was a memory problem as u said, thank you!

• 0 views • link
written by alvarocentron91 • alvarocentron91

Use the -d option to get per base information using bedtools coverage

• 0 views • link
written by Jeffin Rockey Jeffin Rockey

Thank you, the output is closer but I think it is no exactly the one I need. Here it is an example of the output(in bold the new columns generated):

Um_chr06        frame_-2_peak_9719      noCDS   214798  215077  .       -       .       90      **197     11**
Um_chr06        frame_-2_peak_9719      noCDS   214798  215077  .       -       .       90      **198     11**
Um_chr06        frame_-2_peak_9719      noCDS   214798  215077  .       -       .       90      **199     11**
Um_chr06        frame_-2_peak_9719      noCDS   214798  215077  .       -       .       90      **200     11**

I need 1 result per peak not 10, I think the best option will be convert the bam to bed and do bedtools intersect

• 0 views • link
updated by Ram Ram • written by alvarocentron91 • alvarocentron91

Reminder: Use the code formatting option to better your post.

• 0 views • link
written by Ram Ram
  • Please give gff itself as -a in bedtools coverage. From the output shown seems like the command was executed the other way.
  • bedtools can work directly with bam for intersect in case it itself is required.
• 0 views • link
written by Jeffin Rockey Jeffin Rockey
  • Please give gff itself as -a in bedtools coverage. From the output shown seems like the command was executed the other way.
  • bedtools can work directly with bam for intersect in case it itself is required.
• 0 views • link
written by Jeffin Rockey Jeffin Rockey

Log in to answer this question.

More posts like this