p-value for correlation/co-occurrence in sets of intervals I'd like to have a probability measure for the co-occurrence of intervals in two datasets of intervals.
dataset1: [125,500],[900,1300],[2220,2500], ...
dataset2: [600,800],[1200,1400],[3020,3500], ...

(Two exemplary datasets with an overlap in the second intervals)
The intervals within a dataset are represented by the coordinates of their start and end (just two integer values). Within a dataset intervals are non-overlapping.
The number, average lengths and distribution of intervals may be very different between two datasets, but the length of the dataset that is the maximal integer a coordinate can have is identical and known for both.
I am looking for a kind of p-value for the probability that there is a "correlation" between two datasets, that is that intervals overlap more frequently (or on the whole more widely in terms of intersection length) than can be expected by chance.
I would be grateful for suggestions on how to address the problem in practical as well as statistical terms (perhaps using R).
Perhaps there is even a solution for the "correlation" between multiple datasets ... 
Thanks a lot :-)
 A: Absent a probability model, you can consider a permutation test.  This requires two things:


*

*A measure of overlap.

*A way to re-use the data to create alternative hypothetical datasets for comparison.
A test statistic
One possible measure of overlap is obtained by considering each dataset to represent a union of intervals--that is, a subset of the number line--and to compute the total amount of their set-theoretic intersection.  Let's start our programming by implementing this.  It will be convenient to represent the data as $2\times n$ arrays, with the first row giving the start and the second row giving the end of each interval, thus:
d.1 <- matrix(c(125,500, 900,1300, 2220,2500), nrow=2)
d.2 <- matrix(c(600,800, 1200,1400, 3020,3500), nrow=2)

It should not matter that the intervals are sorted within each matrix, because conceptually these are just collections of the intervals that were observed.  (I won't bother to program any checks to enforce the requirement that intervals within a dataset have no mutual overlaps: the user will be responsible for assuring this.  The presence of overlaps will not break any of the code that is later written.)
The amount of overlap can be computed as
overlap <- function(x,y) {
  o <- function(i,j) max(0, min(c(i[2]-i[1],i[2]-j[1],j[2]-i[1],j[2]-j[1])))
  sum(apply(x, 2, function(u) apply(y, 2, function(v) o(u,v))))
}
stat <- overlap(d.1, d.2)

This is a little crude--it compares all possible ordered pairs of intervals--but it's fast enough to illustrate the ideas.  If the permutation test turns out to take too long, this is the function to optimize (by sorting the intervals within each dataset, its timing can be reduced to a linear function of the number of intervals rather than a quadratic function).
I have saved the overlap of the actual data in stat for future reference.  For these data, it equals $100$, which is the length of $[1200, 1300]$, the intersection of the two datasets.
The permutation test
One way to modify the data, while preserving whatever structure is evident--I don't know whether this is appropriate for such an experiment but I have nothing else to go on--is to view the interval lengths as independent realizations of one random variable $X$ and the gaps between the intervals as independent realizations of another variable $Y$.  The null hypothesis (that the two datasets do not really differ) is that each is obtained by a suitable number of draws from $(X,Y)$.  The permutation test draws without replacement from the combined realizations $(x_i,y_i)$ obtained from pooling the data.
For instance, the within-interval gaps, $x_i$, in the first dataset are $375=500-125$, $400=1300-900$, and $280=2500-2220$.  The between-interval gaps, $y_i$, in the first dataset are $125=125-0$, $400=900-500$, and $920=2220-1300$.  This function computes them:
gaps <- function(x) {
  y <- as.vector(x[, order(x[1,])])
  matrix(diff(c(0, y)), nrow=2, dimnames=list(c("Between", "Within"), NULL))
}

(Notice how it protects itself against unordered input by sorting by the starts of each interval.)
Later, after permuting the between-interval gaps among themselves and the within-interval gaps among themselves, we will need to reassemble the gaps into a facsimile of a dataset:
assemble <- function(x) matrix(cumsum(x), nrow=2)

Using these preliminaries, it is now easy to construct a permutation test or even for bootstrapping the distribution (where sampling is done with replacement): an optional parameter, mysteriously indicated by ... here, will determine what happens.  The parameter n requests the number of iterations.  The two datasets are matrices x and y:
boot <- function(x, y, n=1, ...) {
  m.x <- dim(x)[2]; m.y <- dim(y)[2]
  g <- cbind(gaps(x), gaps(y))

  trial <- function() {
    z <- t(apply(g, 1, sample, ...))
    overlap(assemble(z[, 1:m.x]), assemble(z[, 1:m.y + m.x]))
  }
  replicate(n, trial())
}

The optional parameter ... is used by sample.  To sample with replacement, set that parameter to replace=TRUE.
Example
Let's perform both a bootstrap and a permutation test, since we're capable of doing both with the same code. To make the results directly comparable, they each use exactly the same sequence of random numbers:
set.seed(17); sim.1 <- boot(d.1, d.2, 1000)               # Permutation test
set.seed(17); sim.2 <- boot(d.1, d.2, 1000, replace=TRUE) # Bootstrap

The p-value depends on the alternative hypothesis, but in any event will be found in terms of the number of trial results having an overlap less than or equal to what was observed:
weight <- function(x) ifelse(x < 0, 1, ifelse(x == 0, 1/2, 0))
p.1 <- mean(weight(sim.1)); p.2 <- mean(weight(sim.2))

For this calculation, these values turn out to be $0.047$ and $0.0465$, respectively.  That's just small enough to give a whiff of "significance" to a one-sided test of whether there is too little overlap.
A single number, like p.1 or p.2, does not fully reveal the simulation results. Let's plot them:
breaks <- hist(c(sim.1,sim.2), plot=FALSE)$breaks # Get common bins for the plots

par(mfrow=c(1,2))
hist(sim.1, probability=TRUE, breaks=breaks, main="Permutation distribution")
abline(v = stat, lwd=2, col="Red")
hist(sim.2, probability=TRUE, breaks=breaks, main="Bootstrap distribution")
abline(v = stat, lwd=2, col="Red")


Comments
To use a different measure of overlap, recode overlap and proceed as above.
Other probability models can be accommodated within this framework, but require more coding.  For instance, we might postulate (perhaps based on scientific grounds) that the width of an interval might be correlated with the gaps to its neighboring intervals.  If that is the case, we would not draw the between-interval widths and within-interval widths independently of each other.  But now we're starting to fit an actual probability model to the data and we might want to start thinking about a parametric test or a parametric bootstrap.
A: Have you considered this paper "Subsampling methods for genomic inference" ( http://arxiv.org/pdf/1101.0947.pdf ) from the ENCODE project.  They have associated software.  There are a few approaches such as Cooccur package for R.  A more complete discussion can be found here: http://www.biostars.org/p/5484/
EDIT:  I got a message saying my answer was a little too short, so I'm elaborating on why I provided it. The reason I'm mentioning this solution is that overlap of genomic features is a well-researched topic in bioinformatics with a lot of associated literature.  There are lots of aspects of genomes that might not be anticipated by an approach from first principles.  ENCODE is a large research project for aggregating genomic data in humans, and the method I mentioned (which was developed within the ENCODE project) is thus quite widely used.
I have actually spent a lot of my time over the last few weeks reading and rereading this 'Subsampling' paper.  It's very involved. They model the genome as a "piecewise stationary ergodic random process".
