Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-06T17:43:18.690Z
Answers.Text: 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?
Answers.UserDisplayName: Eliran Turgeman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-06T17:53:18.587Z
Answers.Text: 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
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-06T17:54:52.360Z
Answers.Text: easiest thing to do is to re-run the above simulation with a small "universe"
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-03T14:54:56.783Z
Answers.LastActivityDate: 2020-09-03T14:54:56.783Z
Answers.UserDisplayName: StupidWolf
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-03T13:32:39.140Z
Comments.Text: Not so trivial, if you want to use a distribution. https://stats.stackexchange.com/questions/52004/overlapping-sets-3-significance-test
Comments.UserDisplayName: StupidWolf
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-03T13:33:36.133Z
Comments.Text: Easiest way is to simulate the data, and see how many times you end up with an intersection set larger than what you see
Comments.UserDisplayName: StupidWolf
Tags: <r>
Title: Compute the significance of the overlap between 2 or more gene sets
ThreadText: # 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"])
  
  ```
  
  [![enter image description here](https://i.sstatic.net/ZV7Yl.png)](https://i.sstatic.net/ZV7Yl.png)
  
  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](https://i.sstatic.net/plpKC.png)](https://i.sstatic.net/plpKC.png)
  
  The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this [post](https://stats.stackexchange.com/questions/109207/p-values-equal-to-0-in-permutation-test):
  
  ```
  (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"  
  
  ---
  
  