This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How Do You Calculate If Two Sets Of Genomic Regions Overlap Significantly?

Given two sets of genomic locations, such as two bed files of ChIP-Seq peaks, how would you calculate if they overlap more than would be expected by chance?

I have been trying a couple of methods.

1) Using the BEDTools shuffleBed tool to randomise one or both of the files repeatedly (say 1000 times) then running intersectBed to determine a distribution of expected overlaps. I then work out an empirical p-value based on the number of the random samples that overlaps to a greater degree than my real overlap.

2) Using the GSC (genome structural correction) tools from Encode.

I have been thinking that method 1 limited to specific regions, rather than whole genome, might be a better approach. I think the assumption that the regions could be anywhere on the genome when shuffled is not accurate. This is something that method two tries to address, but I just can not seem to get it work reliably.

So how should I calculate the significance of my overlaps? My aim is to look at the overlap of peaks between biological replicates of CHiP-Seq and between different peaks callers etc.

chip-seq coverage statistics bed
• 22,400 views • link
updated by bwd • bwd • written by Stew Stew

To anyone arriving to this question, at this other question there's a useful list of tools and some more detailed answers

• 0 views • link
written by bernatgel bernatgel

3 answers

When comparing set A to set B, that are from ChIP-seq data, I have used a method where I obtain random regions to resemble B (the same length) from the mapable region of, in my case, the human genome. This data can be obtained from the test UCSC browser 'CRG mapability'. The overlap between A and B, to A and random can be tested by a Fisher Exact Chi-square. This way random regions cannot come from areas of the genome that could never be binding regions.

I am also looking at alternative methods.

You might be interested in the The Genome HyperBrowser http://hyperbrowser.uio.no/hb/. Where you can ask questions just like this. It offers are range of null hypotheses, some complex like maintaining inter-feature distance, and randomising the feature calls, etc. I quite like the look of this and intend to test it out myself.

David Huen (Cambridge, UK) gave a seminar and discussed this problem. He has an R package called 'cooccur', see http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2910723/ which does a similar sort of comparison.

• 1 views • link
written by Ian Ian

Thanks for the answer, I am trying Cooccur now.

• 0 views • link
written by Stew Stew

Hi! I am also on my way to use Cooccur package, but I am experiencing som problem like the package works pretty well with the example data. but when I load in my data

library(Cooccur) a<-read.table("path") b<-read.table("path") (it works fine till here, but on next step it gives error) explore_pairs(c("a","b")) ERROR: Error in sum(sapply(tf1_s, score_sample, tf2_hits = tf2_s, hit_list = hit_l)) : invalid 'type' (list) of argument could you please tell me what is going wrong. Thank you for your time and your article on 21 TF in drosophila is quite informative.abhisheksinghnl@gmail.com

• 0 views • link
written by Dataminer Dataminer

Hi! [?] I am also on my way to use Cooccur package, but I am experiencing som problem like the package works pretty well with the example data. but when I load in my data

a<-read.table("path") [?] b<-read.table("path") (it works fine till here, but on next step it gives error) [?] explore_pairs(c("a","b")) [?] ERROR: Error in sum(sapply(tf1_s, score_sample, tf2_hits = tf2_s, hit_list = hit_l)) : invalid 'type' (list) of argument could you please tell me what is going wrong. Thank you for your time and your article on 21 TF in drosophila is quite informative.abhisheksinghnl@gmail.com

• 0 views • link
written by Dataminer Dataminer

Hi! I am also on my way to use Cooccur package, but I am experiencing some problem like the package works pretty well with the example data. but when I load in my data

>library(Cooccur)
a<-read.table("path")
b<-read.table("path")

everything goes is ok till here.

>explore_pairs(c("a","b"))

This is where I get stuck.
ERROR:

Error in sum(sapply(tf1_s, score_sample, tf2_hits = tf2_s, hit_list = hit_l)) :
  invalid 'type' (list) of argument

Thank you, It will be very kind of you if you could email me on my id abhisheksinghnl@gmail.com

• 0 views • link
updated by Ram Ram • written by Dataminer Dataminer

make sure your data have exactly three columns named 'seq' 'start' 'end' in that order

• 0 views • link
written by Malcolm.Cook Malcolm.Cook

This is a trickier problem than it first seems, and there appears to be no consensus yet on the best method to do this. For a summary of most methods available, see the review by Fu and Adryan on this topic: http://www.ncbi.nlm.nih.gov/pubmed/19763325

For interval-based data, I second Ian's suggestion about Huen's method, which follows on from work by Haiminen et al.: http://www.ncbi.nlm.nih.gov/pubmed/18691400.

• 1 views • link
written by Casey Bergman Casey Bergman

That review looks interesting, but unfortunately it is not open access, and we do not have a subscription. Peter Bickel's group at Berkeley did a lot of work on this for ENCODE, that we used on the last CHiP-chip paper, but his methods do not seem to be taking off in the CHiP-Seq world.

• 0 views • link
written by Stew Stew

MSPC is designed specifically for this purpose!

Download it from this page, and use it as the following: 

MSPC -i rep1.bed -i rep2.bed -i rep3.bed -r biological -s 1E-8 -w 1E-4 -g 1E-9 -c 2 -a 0.005

where

  • -s sets stringent threshold: peaks in input samples with p-value more significant than this threshold, are considered "stringent" peaks;
  • -w sets weak threshold: peaks in input samples with p-values less significant that this threshold, are considered "weak" peaks;
  • -g how strong should be p-value of overlapping peaks so to consider them "significant";
  • -c a peak should overlap with how many peaks at minimum to be considered for analysis;
  • -a BH multiple testing correction.

See this page for a complete list of arguments.

• 0 views • link
written by bwd • bwd

Log in to answer this question.

More posts like this