<p>This is a shot at it, first an example dataset:</p>
<pre><code>import matplotlib.pyplot as plt  
import numpy as np
import functools
from matplotlib_venn import venn3
    
# define universe
uni = [&quot;gene&quot;+str(i) for i in range(1000)]
# some overlap
gs1 = uni[250:300] + uni[900:950]
gs2 = uni[:300]
gs3 = uni[250:500]
</code></pre>
<p>The ever amazing venn diagram:</p>
<pre><code>venn3([set(gs1),set(gs2),set(gs3)],set_labels=[&quot;gs1&quot;,&quot;gs2&quot;,&quot;gs3&quot;])
</code></pre>
<p><a href="https://i.sstatic.net/ZV7Yl.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ZV7Yl.png" alt="enter image description here" /></a></p>
<p>Then a function to draw a set with length equivalent of each set, randomly from the universe and find length of intersection (all 3):</p>
<pre><code>def sim_intersect(uni,set_lengths):
    randomsets = [np.random.choice(uni,n) for n in set_lengths]
    return len(functools.reduce(np.intersect1d,randomsets))
</code></pre>
<p>We run this 1000 times:</p>
<pre><code>permuted_values = [sim_intersect(uni,[len(gs1),len(gs2),len(gs3)]) for i in range(1000)]

plt.hist(permuted_values,bins=range(50))
</code></pre>
<p><a href="https://i.sstatic.net/plpKC.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/plpKC.png" alt="enter image description here" /></a></p>
<p>The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this <a href="https://stats.stackexchange.com/questions/109207/p-values-equal-to-0-in-permutation-test">post</a>:</p>
<pre><code>(sum(np.array(permuted_values)&gt;obs_n)+1)/(1000+1)
0.000999000999000999
</code></pre>

CC BY-SA 4.0
2020-09-06T17:43:18.690Z
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?
Eliran Turgeman
CC BY-SA 4.0
2020-09-06T17:53:18.587Z
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
StupidWolf
CC BY-SA 4.0
2020-09-06T17:54:52.360Z
easiest thing to do is to re-run the above simulation with a small "universe"
StupidWolf
CC BY-SA 4.0
2020-09-03T14:54:56.783Z
2020-09-03T14:54:56.783Z
StupidWolf
<p>I was able to compute the significance of the overlap between 2 gene sets using the cdf function of <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.hypergeom.html" rel="nofollow noreferrer">scipy hypergeometric distribution</a>.</p>
<p>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?</p>
<p>Are there any websites that provides the same calculations over gene sets so I can validate my results?</p>

CC BY-SA 4.0
2020-09-03T13:32:39.140Z
Not so trivial, if you want to use a distribution. https://stats.stackexchange.com/questions/52004/overlapping-sets-3-significance-test
StupidWolf
CC BY-SA 4.0
2020-09-03T13:33:36.133Z
Easiest way is to simulate the data, and see how many times you end up with an intersection set larger than what you see
StupidWolf
<r>
Compute the significance of the overlap between 2 or more gene sets
# Compute the significance of the overlap between 2 or more gene sets

(Asked by: username_0 on 2020-09-02)

I was able to compute the significance of the overlap between 2 gene sets using the cdf function of [scipy hypergeometric distribution](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.hypergeom.html).

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?

**Tags:** r
> **Comments:**
> **username_1**: Not so trivial, if you want to use a distribution. https://stats.stackexchange.com/questions/52004/overlapping-sets-3-significance-test  
> **username_1**: Easiest way is to simulate the data, and see how many times you end up with an intersection set larger than what you see  

---

## 1 answers

### Answer by username_1 on 2020-09-03. Score: 3 | Accepted answer

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"])


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


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

> **Comments:**
> **username_0**: 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?  
> **username_1**: 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  
> **username_1**: easiest thing to do is to re-run the above simulation with a small "universe"  

---

