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 :-)

-
It is difficult to understand what is being asked here, because "track" lacks a clear definition and no examples are provided. Could you perhaps edit your question to rectify this? –  whuber Dec 7 '12 at 14:47
@whuber: Done, has it become a bit more clear now? –  datamole Dec 7 '12 at 15:41
Great! Let's get to the question, then: you mention "chance." What is your probability model for these intervals? (Without a model, we cannot quantify what "chance" means, so we would have no hope of finding a p-value.) What do these intervals represent? How do they arise? –  whuber Dec 7 '12 at 15:43
@whuber: They are empirical data (chromosome positions) that can come from different types of experiments. I don't think that you can presume any common distribution of the interval-positions per se.Having two sets of intervals of a certain density with random distribution there should be some random overlap as a background –  datamole Dec 7 '12 at 15:54
I'm afraid you have lost me. If you cannot make any assumptions about distributions, you cannot perform any probabilistic analysis, so obtaining a p-value is out of the question. We have to know something about your "random distribution" in order to proceed. –  whuber Dec 7 '12 at 15:58

Absent a probability model, you can consider a permutation test. This requires two things:

1. A measure of overlap.

2. 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")


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.

-

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".

-
Thanks! The paper looks very interesting indeed - though I still have to digest it. I actually had had a look into the Genomic Hyperbrowser discussed in the other blog that you've linked and haven't really got useful results using that one ... where can I find the associated software of the Encode-related paper and have you tried that one? –  datamole Dec 17 '12 at 11:19
Here you go: encodestatistics.org with direct link to the software here: encodestatistics.org/releases/block_bootstrap-0.8.1.zip –  Henry B. Dec 17 '12 at 11:23