Question: bedtools intersect usage problem
0
gravatar for ja4123
26 days ago by
ja41230
ja41230 wrote:

Hello! I need to filter peaks which are withwin gene regions. I have .bam file from CHIP-seq experiment and .bam file from RNA-seq. I am trying to do this with this command:

bedtools intersect -a A.bam -b B.bam > AB.bam

The problem is that the output AB.bam is exactly the same as A.bam (it is visuable in IGV). Can someone give some advices?

bam intersect bedtools • 144 views
ADD COMMENTlink modified 26 days ago by i.sudbery11k • written 26 days ago by ja41230
1

Can you add a few lines of A.bam and B.bam? especially lines that show some of the intersections? You can modify the data if it's confidential, but if you add more details you'll get better answers. Also, check out similar questions that appears on the right section of your post.

Try adding one of these options:

-wa Write the original entry in A for each overlap.

-wb Write the original entry in B for each overlap. Useful for knowing what A overlaps. Restricted by -f and -r.

https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

ADD REPLYlink modified 26 days ago • written 26 days ago by Fatima930
2
gravatar for i.sudbery
26 days ago by
i.sudbery11k
Sheffield, UK
i.sudbery11k wrote:

I don't know what is going wrong with the command you are using. But to filter peaks in that are in gene regions, I would normally do

 bedtools intersect -u -a peaks.bed -b geneset.gtf

Where peaks.bed is the output of a peak caller (and if it was MACS then the file would be .narrowPeak or .broadPeak not .bed) and geneset.gtf is either annotation downloaded from an annotation database (e.g. refseq or ENSEMBL) or is created from the RNAseq using a transcript assembly program like stringtie. The -u means Write original A entry once if any overlaps found in B.

ADD COMMENTlink written 26 days ago by i.sudbery11k

Thanks for the answer! Could you explain to me what is the output from the command you wrote and if it is possible to store it in .bam file or any other file?

ADD REPLYlink modified 26 days ago • written 26 days ago by ja41230
1

The output from this command is a bed file of peaks. No, if would be difficult to store this in a BAM file. However, you said you wanted peaks that overlapped genes, and BAM files do not store peaks.

ADD REPLYlink written 26 days ago by i.sudbery11k

It is possible to get list of genes (gene id) in which peaks fall into?

ADD REPLYlink written 26 days ago by ja41230
1

You probably want bedtools closest -a peaks.bed -b geneset.gtf for that.

ADD REPLYlink written 26 days ago by i.sudbery11k

I will try that. Thanks for help!

ADD REPLYlink written 26 days ago by ja41230
Please log in to add an answer.

Help
Access

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