# Calculating the probability of gene list overlap between an RNA seq and a ChIP-chip data set

Hopefully someone on these forums can help me out with this basic problem in gene expression studies.

I did deep sequencing of an experimental and a control tissue. I then obtained fold enrichment values of genes in the experimental sample over control. The reference genome has ~15,000 genes. 3,000 out of 15,000 genes are enriched above a certain cut-off in my sample of interest compared to control.

So: A = total gene population = 15,000 B = RNA-Seq enriched subpopulation = 3,000.

In a previous ChIP-chip experiment, I found 400 genes that are enriched by ChIP-chip. Of the 400 ChIP-chip genes, 100 genes are in the group of 3,000 enriched RNA-Seq transcripts.

So: C= total # of ChIP-chip enriched genes = 400.

What is the probability that my 100 ChIP-chip genes would be enriched by RNA-Seq by chance alone? In other words, what is the most prudent way to calculate if my observed overlap between B and C (100 genes) is any better than that obtained by chance alone? From what I have read so far, the best way to test this is by using hypergeometric distribution.

I used an online calculator (stattrek.com) to set up a hypergeometric distribution test with the following parameters: - pop size=15,000 - # of successes in population=3,000 - sample size=400, -# of successes in sample=100. I get the following for Hypergeometric Probability P(x=100)= 0.00224050636447747

The actual # of genes overlapping between B and C = 100. Is this better than by chance alone? Doesn't look like it is if the chance of any one gene being enriched is 1:5 (3,000 out of 15,000). That's why I don't understand how come my P(x=100) I calculated above is 0.0022. That amounts to a 0.2% chance of the overlap occurring by chance. Shouldn't this be much higher?

If I sampled 400 random genes rom the big list of 15,000, then any 80 of these genes would be expected to be enriched by chance alone (1:5). The number of genes that are actually overlapping is 100, so this is just slightly better than by chance.

I also tried to come up with a solution using the dhyper or phyper functions in R (using what I saw in another post): A=all genes in the genome (15,000) B=RNA-Seq enriched genes (3,000) C=ChIP-chip enriched genes (400) Here's the R input/output (adapted from a previous stackexchange post):

> totalpop <- 15000
> sample1 <- 3000
> sample2 <- 400
> dhyper(0:2, sample1, totalpop-sample1, sample2)
[1] 4.431784e-40 4.584209e-38 2.364018e-36
> phyper(-1:2, sample1, totalpop-sample1, sample2)
[1] 0.000000e+00 4.431784e-40 4.628526e-38 2.410304e-36


I'm not sure how to interpret these numbers. I believe 2.36e-36 is the probability of getting a complete overlap between B and C by chance alone? But this makes no sense, since that probability is much closer to 1:5. If I start with 15,000 genes, 3,000 will be enriched. Similarly, if I start with 400 ChIP-chip genes, 80 of them should be enriched in the RNA-Seq alone due to the 1:5 chances of enrichment in that data set.

What is the proper way to calculate the p-value, according to the hypergeometric distribution, for the overlap of B and C?

-

You are close, with your use of dhyper and phyper, but I don't understand where 0:2 and -1:2 are coming from.

The p-value you want is the probability of getting 100 or more white balls in a sample of size 400 from an urn with 3000 white balls and 12000 black balls. Here are four ways to calculate it.

sum(dhyper(100:400, 3000, 12000, 400))
1 - sum(dhyper(0:99, 3000, 12000, 400))
phyper(99, 3000, 12000, 400, lower.tail=FALSE)
1-phyper(99, 3000, 12000, 400)


These give 0.0078.

dhyper(x, m, n, k) gives the probability of drawing exactly x. In the first line, we sum up the probabilities for 100 – 400; in the second line, we take 1 minus the sum of the probabilities of 0 – 99.

phyper(x, m, n, k) gives the probability of getting x or fewer, so phyper(x, m, n, k) is the same as sum(dhyper(0:x, m, n, k)).

The lower.tail=FALSE is a bit confusing. phyper(x, m, n, k, lower.tail=FALSE) is the same as 1-phyper(x, m, n, k), and so is the probability of x+1 or more. [I never remember this and so always have to double check.]

At that stattrek.com site, you want to look at the last row, "Cumulative Probability: P(X $\ge$ 100)," rather than the first row "Hypergeometric Probability: P(X = 100)."

Any particular number that you draw is going to have small probability (in fact, max(dhyper(0:400, 3000, 12000, 400)) gives $\sim$0.050), and getting 101 or 102 or any larger number is even more interesting that 100, and the p-value is the probability, if the null hypothesis were true, of getting a result as interesting or more so than what was observed.

Here's a picture of the hypergeometric distribution in this case. You can see that it's centered at 80 (20% of 400) and that 100 is pretty far out in the right tail.

-
Thanks a bunch for your help. I understand the logic behind your answer. But how do I explain to a bunch of biologists that this is greater than overlap observed due to chance alone? They will say that I have a 1:5 chance of overlap. Is my overlap significant because in a sample size of 400 balls (out of 15,000 total balls), my chance of getting a white ball is actually lower than 1:5 because I am sampling a smaller population (not the entire 15,000)? This doesn't make sense because even though 400 < 15,000, there is still a 1:5 ratio of white:black. Does this make sense? –  stlandroidfan Oct 1 '11 at 20:51
@stlandroidfan - I don't understand what you're finding confusing. I've added a figure; does this help? –  Karl Oct 2 '11 at 3:25

Look at it this way.. If you assumed it as a binomial, which may not be correct, but it should be fairly approximate.. your sigma^2 is .8*.2*400= 64, then sigma = 8. So from 80 to 100 you've gone 2.5 standard deviations.. This is pretty significant.. It should have a small p-value.

-
Thanks for your answer. Hypergeometric distributions tend to be used more frequently for gene list overlaps from what I have seen in literature. The question is what is the probability of getting 100 or more white balls in a sample of size 400 from an urn with 3000 white balls and 12000 black balls? I think I'm still perplexed by how to explain this to a bunch of biologists? The way they see it is 3000:12000 is a 1:5 chance of white:black. So in a sampling of 400, 80 should be white. So how come the probability of getting 100 or more is so much lower than 20% (1 in 5)? –  stlandroidfan Oct 1 '11 at 20:56

I hope I'm not being annoying with my follow-up questions. I know basic statistics, but I am a biologist with a background in genomics and neurobiology, and not biostatistics (or even statistics). I can grasp the basics of what a hypergeometric distribution is. I've had to 'defend' my low calculated probability of overlap against criticism. I have not been able to come up with a good, coherent answer as to why my result of 100 overlapping genes is MORE SIGNIFICANT than what I could have obtained by chance alone? To biologists, 80 vs. 100 is quite close. But I believe that the answers I received from both of you are quite sound statistics. So I guess my answer could be that when proper statistical tests are used, the overlap is significant at P(X ≥ 100)=0.0078. This relatively low p-value seems to indicate that this is significant.

When would this no longer be significant? If the p-value was greater than some commonly accepted p-value (let's say 0.05?).

I just want to make sure that I have a good answer backed by prudent statistics and not some hoaky answer. If critics are not convinced by my reasoning, there is not much I could do other than to refer them to some book on statistics.

I want to employ a test that would pass the criticism of an expert biostatistician. And I believe that the above calculation with the 0.0078 p-value is sound, valid and correct. Right?

Somewhere else, I saw the mentioning of 'block bootstrapping' as the method to use, but that sounds way too complicated.

Thanks a lot guys. I thought I would have had to consult a biostatistician at my med school.

-
I would run a simulation. If you need help with that, let me know. –  Adam Oct 2 '11 at 1:30
Thanks. How complicated would it be to run the simulation? Can I do it in R or would I have to use MATLAB? I'm somewhat familar with R, but never done MATLAB. –  stlandroidfan Oct 2 '11 at 1:37
Can def be done in R. –  Adam Oct 2 '11 at 2:51
Look at it this way.. a season of at bats for a baseball player is 400 at bats.. imagine a lifetime .200 hitter batting .250! Also note, the closer a hitter bats to 0.000, thee less likely is a. 050 improvement.. –  Adam Oct 2 '11 at 3:03
stlandroidfan. I've just reread the posts here, and I think your real problem is that you don't understand what a pvalue is. Pvalues which are close to zero mean they are unlikely to happen on chance alone. For example, if you flip a coin 100 times and get heads everytime..Calculate the pvalue = .000...000001 ((1/2)^100), right (basically impossible)? Now, what do you conclude.. "If the coin was fair, the chances of observing heads 100 times in a row is .0000..0001 " So what do you conclude.. These results are near impossible, the coin is double sided! –  Adam Oct 3 '11 at 3:13