Compare Super Enhancers in normal vs breast cancer cells

Question

[mropri]

I did ChIP seq using H3K27ac antibody to quantify super enhancer regions in normal vs breast cancer cells. I ran ROSE algorithm on normal mammary epithelial cells and breast cancer cells to get a list of super enhancer regions in each cell line. I have a list of super enhancer regions for both and I want to compare super enhancers between the two outputs I got. How would i compare the two lists and see if they have super enhancer regions that overlap or not to analyze acquired/lost super enhancers between the two cell lines. I know bedtools intersect might be able to do it, but am confused about the command I would input to compare the two files. Any help is appreciated. Thank you

Answer 1

[Friederike]

Not sure I grasp the full scope of the issue, but a typical bedtools intersect command would simply go like this:

bedtools intersect -a super_enhancer_list-A.bed -b super_enhancer_list-B.bed > super_enhancers_from_A_that_are_also_in_B.bed

[mropri]

Thank you, that is what I was looking for. Appreciate your help!

[mropri]

I want to compare the signal of the super enhancer regions that were classified by ROSE in normal and breast cancer cells and see calculate fold change of those regions to see if they are increasing or decreasing in intensity. I have a bed files of super enhancer regions for normal cells (MCF10A) and same for breast cancer cells (MCF10A-CA1). I know what I have to do. Briefly, collapse A.SE.bed and B.SE.bed into AllSEs.bed using bedtools merge. Quantify coverage of All SEs.bed using A.bam and B.bam and bedtools intersect. Normalize the read coverage to the sequencing depth, and calculate the fold-change between them. The problem I was having is how can I collapse or merge two or three bed files into one file that will contain all the super enhancer regions classified and quantify coverage? Any help is appreciated. Thank you so much!

Answer 2

[jared.andrews07]

If you are looking to compare signal for a common set of regions between groups, consider a method like csaw or diffBind rather than trying to come up with a home-grown solution.

From personal experience, I recommend that you compare signal for the constituent peaks rather than the entire SE, as differences for such large regions tend to get muddled due to the peak stitching inherent in the SE calling process.

[mropri]

I want to compare peaks and I think I can do that by collapse A.SE.bed and B.SE.bed into AllSEs.bed using bedtools merge. Quantify coverage of All SEs.bed using A.bam and B.bam and bedtools intersect. Normalize the read coverage to the sequencing depth, and calculate the fold-change between them. The problem I was having is how can I collapse or merge two or three bed files into one file that will contain all the super enhancer regions classified and quantify coverage? Hope that helps and appreciate your help as always

[ATpoint]

The merge itself can be done with bedtools merge, and the counts per interval you can get with featureCounts, but without replicates you will not have a reliable way to really make any statistics. Do you have replicates? I would also suggest to do a more sophisticated normalization than just per-million as you do (so by sequencing depth) as this might not be enough if library composition changes or data quality is quite different between samples, see A: ATAC-seq sample normalization (quantil normalization)

This all would be taken care of if using any of the tools suggested by jared.andrews07 above.

[jared.andrews07]

I agree with all of these comments and urge you to read the documentation of the linked tools. Without replicates, this is a difficult task and you really should view any results with a multiple grains of salt. bdgcmp from MACS may be worth a look as well.