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
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.

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.

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?

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.

0
Entering edit mode

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

1
Entering edit mode

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

0
Entering edit mode

I will try that. Thanks for help!