2
$\begingroup$

I was able to compute the significance of the overlap between 2 gene sets using the cdf function of scipy hypergeometric distribution.

I wish to be able to perform the same calculation for more than 2 gene sets; should I use the multivariate hypergeometric distribution cdf function for that?

Are there any websites that provides the same calculations over gene sets so I can validate my results?

$\endgroup$
2
  • $\begingroup$ Not so trivial, if you want to use a distribution. stats.stackexchange.com/questions/52004/… $\endgroup$
    – StupidWolf
    Sep 3 '20 at 13:32
  • $\begingroup$ Easiest way is to simulate the data, and see how many times you end up with an intersection set larger than what you see $\endgroup$
    – StupidWolf
    Sep 3 '20 at 13:33
3
$\begingroup$

This is a shot at it, first an example dataset:

import matplotlib.pyplot as plt  
import numpy as np
import functools
from matplotlib_venn import venn3
    
# define universe
uni = ["gene"+str(i) for i in range(1000)]
# some overlap
gs1 = uni[250:300] + uni[900:950]
gs2 = uni[:300]
gs3 = uni[250:500]

The ever amazing venn diagram:

venn3([set(gs1),set(gs2),set(gs3)],set_labels=["gs1","gs2","gs3"])

enter image description here

Then a function to draw a set with length equivalent of each set, randomly from the universe and find length of intersection (all 3):

def sim_intersect(uni,set_lengths):
    randomsets = [np.random.choice(uni,n) for n in set_lengths]
    return len(functools.reduce(np.intersect1d,randomsets))

We run this 1000 times:

permuted_values = [sim_intersect(uni,[len(gs1),len(gs2),len(gs3)]) for i in range(1000)]

plt.hist(permuted_values,bins=range(50))

enter image description here

The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this post:

(sum(np.array(permuted_values)>obs_n)+1)/(1000+1)
0.000999000999000999
$\endgroup$
3
  • $\begingroup$ How would a population size parameter fit in this approach? Doesn't the p value should be affected by the total 'available genes' ? For example,. if I have two gene sets A, B. |A| = 250, |B| = 220, |A&B| = 65 As far as I understand, the p value should change dramatically if the population size is 300 or 30,000. Could you please explain where/how does it fit in this approach? $\endgroup$ Sep 6 '20 at 17:43
  • 1
    $\begingroup$ Yes you are right, the total available genes would matter. In the code above, you would define the set of available genes as uni . It is not like a hypergeometric where you input a number of genes parameter. You simply simulate the results and see how likely you end up with a result as extreme as yours $\endgroup$
    – StupidWolf
    Sep 6 '20 at 17:53
  • 1
    $\begingroup$ easiest thing to do is to re-run the above simulation with a small "universe" $\endgroup$
    – StupidWolf
    Sep 6 '20 at 17:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.