Trouble with bedtools groupby - not summarizing by columns
3
0
Entering edit mode
7.3 years ago
yarmda ▴ 40

I'm trying to use bedtools intersect and group by to get a count of reads by region they overlap. This should be straightforward, but the output I'm getting is just a number (presumably, the count of all reads that overlapped any region).

Process:

bamToBed mybam.bam > mybam.bed
sort -k1,1 -k2,2n -k3,3n mybam.bed > sort.mybam.bed
bedtools intersect -a sort.mybam.bed -b regions.bed -wa -wb | groupBy -g 1-3 -c 8 -o count > groups.bed

As far as I can tell from the documentation, that should do just what I need. I receive no errors when running this, but no columns print in the output. I understand that the input files must be sorted on the columns to be output, but I think I've done that.

Any thoughts?

bedtools groupby sort bam • 3.9k views
ADD COMMENT
1
Entering edit mode

I'm trying to use bedtools intersect and group by to get a count of reads by region they overlap

use GATK DepthOfCoverage

ADD REPLY
2
Entering edit mode
7.3 years ago

What version of bedtools are you using? In v2.26.0 there are some issues with groupby and the one you describe may be similar to #435. If this is the case try downgrading to 2.25.

By the way, for what you need you could use intersectbed with -c option: For each entry in A, report the number of hits in B. I.e. like: intersectBed -c -a regions.bed -b sort.mybam.bed. Also consider using the -sorted option for speed.

ADD COMMENT
0
Entering edit mode

Thanks! I didn't even see the -c option for intersectBed. And I do have v2.26 installed, which may have been the issue.

ADD REPLY
0
Entering edit mode
7.0 years ago

I had same problem. By using v2.25, I got expected result

ADD COMMENT
0
Entering edit mode
7.0 years ago
Eijy Nagai ▴ 90

Yes, I got exactly the same problem today. Seems that Bedtools version v.2.26.0 still have the problem with groupby function. I tested with v.2.25 and worked.

ADD COMMENT

Login before adding your answer.

Traffic: 1440 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6