bedtools intersect usage problem
1
0
Entering edit mode
3 months ago
ja4123 • 0

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?

bedtools intersect bam • 214 views
ADD COMMENT
1
Entering edit mode

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 REPLY
2
Entering edit mode
3 months ago

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 COMMENT
0
Entering edit mode

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 REPLY
1
Entering edit mode

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 REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

I will try that. Thanks for help!

ADD REPLY

Login before adding your answer.

Traffic: 1486 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