Hello,

I need to calculate for my BED file with ChIP-seq peaks fold-enrichment of genomic regions containing or intersecting with promoters (versus genome average). How to do this?

I have tried GENOMATIX, which has this option, but it seems that it calculates fold-enrichment incorrectly (does not take into account different lengths of my genomic regions). I have also tried GREAT and CHIP-ENRICH, but what they both calculate is something like the distribution from each peak's midpoint to the nearest gene, which is not what I need.

Do you know any software to do this task?

PS. I do not need to calculate the read enrichment. I need to calculate the enrichment of peaks contained in my BED file, which intersect with promoters, versus the amount of peaks of this size which would intersect with promoters by chance based on genome-average probability to encounter a promoter.

Like this: My reads contain XX% of reads intersecting with promoters. It is expected that one would get YY% of reads of this length intersecting with promoters by chance, and therefore, the enrichment of promoters is equal to XX/YY-fold in comparison with what one would get by chance. So, I am looking for the solution of this problem. I am pretty sure it is realized already in some software, please advise where!

Thanks

There are several existing options:

• bedtools jaccard (command line; bp-based)
• bedtools fisher (command line; interval-based)
• IntervalStats (command line)
• permutation testing with pybedtools (Python)
• Genometricorr (R package)
• poverlap (Python/command line)
• GAT (Python/command line)
• regioneR (R)

Also check out The dilemma of choosing the ideal permutation strategy while estimating statistical significance of genome-wide enrichment which gives a lot of good insight and has a lot of references to read up on.

For what it's worth, I've tried many of these (haven't tried GAT or regioneR) and have settled on IntervalStats. Since it's an asymmetric metric, it has the potential to give more biological insight. And hopefully it avoids bias from permutation strategy choice.

(edit: added GAT, and clarified that I haven't tried it yet)

(edit 3: added regioneR as per @bernatgel's answer)

EDIT: After the PO edited his/her question my answer is no longer relevant, anyway in case someone is interested here it is...:

I have a script localEnrichmentBed.py which might be useful for you. It calculates the enrichment (as log2 fold change) and significance (p value from Fisher test) in target regions, promoters in your case, relative to the regions flanking the targets. So it doesn't use the genomic background but a local background.

I don't know what you want to do downstream but using the whole genome as background (effectively normalizing for library size) is going to give bias results read coverage can be quite uneven.

This is from the program's help:

localEnrichmentBed.py -h
usage: localEnrichmentBed.py [-h] --target TARGET --bam BAM --genome GENOME
                             [--slop SLOP] [--blacklist BLACKLIST]
                             [--tmpdir TMPDIR] [--keeptmp] [--verbose]
                             [--version]

DESCRIPTION
    Compute the read enrichment in target intervals relative to local background.
    
    Typical use case: A ChIP-Seq experiment on a sample returns a number of regions
    of enrichment. We want to know how enriched these regions are in a *different*
    sample. Note that enrichment is quantified relative to the local background
    not relative to an input control.
    
    See also localEnrichmentScore.R to combine replicates and compare
    treatment vs control.
    
OUTPUT:
    bed file with header and columns:
1. chrom
2. start
3. end
4. targetID
5. flank_cnt
6. target_cnt
7. flank_len
8. target_len
9. log10_pval
10. log2fc

EXAMPLE
    localEnrichmentBed.py -b rhh047.bam -t rhh047.macs_peaks.bed -g genome.fa.fai -bl blacklist.bed > out.bed
    
    ## Using pipes:
    samtools view -u rhh047.bam chr18 \
    | localEnrichmentBed.py -b - -t rhh047.macs_peaks.bed -g genome.fa.fai > out.bed

    
Useful tip: Get genome file from bam file:

    samtools view -H rhh047.bam \
    | grep -P "@SQ\tSN:" \
    | sed 's/@SQ\tSN://' \
    | sed 's/\tLN:/\t/' > genome.txt
    
REQUIRES:
    - bedtools suite
    - numpy, scipy

NOTES:
For PE reads, the second read in pair is excluded (by samtools view -F 128)
since coverageBed double counts pairs.

optional arguments:
  -h, --help            show this help message and exit
  --target TARGET, -t TARGET
                        Target bed file where enrichment is to be computed.
                        Use - to read from stdin.
                                        
  --bam BAM, -b BAM     Bam file of the library for which enrichment is to be
                        computed. Use - to read from stdin.
                                           
  --genome GENOME, -g GENOME
                        A genome file giving the length of the chromosomes.
                        A tab separated file with columns <chrom> <chrom lenght>.
                        NB: It can be created from the header of the bam file (see tip above).
                                           
  --slop SLOP, -S SLOP  Option passed to slopBed to define the flanking region (aka background).
                        If `int` each target will be extended left and right this many bases.
                        If `float` each target is extended left and right this many times its size.
                        E.g. 5.0 (default) extends each target regions 5 times its length left and right.
                                           
  --blacklist BLACKLIST, -bl BLACKLIST
                        An optional bed file of regions to ignore to compute
                        the local background. These might be unmappable regions with 0-counts which would
                        inflate the target enrichment.
                                           
  --tmpdir TMPDIR       Temp dir to use for intermediate files. If not set
                        python will get one. A subdir will be created here.
                                           
  --keeptmp             If set, the tmp dir is not deleted at the end of the
                        job (useful for debugging).
                                           
  --verbose, -V         Print to stderr the commands that are executed.
                                           
  --version             show program's version number and exit

Although this is an old thread, since it includes a thorough list of tools, I'll include another one we recently released.

If you are using R, you can use regioneR, that will perform the process sketched by @TriS in another answer. It uses a permutation test approach, so it's not the fastest way to do that, but does not need to rely in modeling and since it implements different randomization strategies and accepts arbitrary genome masks, the test can be tailored to your specific needs. In addition to testing for the overlap between regions, it can test for the distance, the % of overlap or actually any custom function you feed to it.

For example, if you have your peaks and promoters in data.frame's you can do something like:

library(regioneR)
pt <- overlapPermTest(A=peaks, B=promoters)
plot(pt)

or, if you want to evaluate the size of your peaks, you can define a function returning the mean height of the peaks overlapping the promoters and use it in your test:

meanHeight <- function(A, B, ...) {
  overlapping.A <- subsetByOverlaps(A, B)
  return(mean(overlapping.A$height))
}
pt <- permTest(A=peaks, B=promoters, evaluate.function=meanHeight, randomize.function=randomizeRegions, ntimes=100)
plot(pt)

a pretty precise way of doing it is to bootstrap random genomic regions and get the % that overlaps with TSS regions, (edit - it's pretty similar to the approach the Ryan wrote, if you have fun coding, go for it)...to do this you:

- take the size of your set i.e. 100 regions with your xx% intersection

- randomly sample from the genome groups of 100 regions of the same size (bedtools allows you to do so) (with replacement)

- calculate the % of those that overlap and save value in array

- repeat 1000 or 10.000 times

- plot the distribution of this (expected) vs the xx% intersection that you have from your initial regions (this will be a histogram for the expected and 1 vertical line for the observed values)

- calculate the empirical pvalue as: ((number elements in simulated data  >= xx%) +1 / ((number sampling) + 1) (the +1 is to avoid dividing by zero)

to get your "fold enrichment" take the average of the simulated data vs xx%

or you could just use hypergeometric distribution, but you won't have the distribution of the simulated data

I did not mean to calculate the read enrichment. I meant to calculate the enrichment of promoters in my reads. Like this:

My reads contain XX% of reads intersecting with promoters. It is expected that one would get YY% of reads of this length intersecting with promoters by chance, and therefore, the enrichment of promoters is equal to XX/YY-fold in comparison with what one would get by chance. I am looking for the solution of this problem.

I have tried bedtools jaccard -- it is great, except that it does not do exactly what I need. It calculates the enrichment of basepairs intersecting in my BED file and in the BED file with promoters.

However, what I need the enrichment of regions in my BED file intersecting (having at least one bp overlap) with promoters. Could you please advise with respect to this? Thanks!