bedtools intersect usage problem

Question

[ja4123]

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?

[Fatima]

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

Answer 1

[i.sudbery]

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.

[ja4123]

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?

[i.sudbery]

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.

[ja4123]

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

[i.sudbery]

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

[ja4123]

I will try that. Thanks for help!