gseapy/stats.py
def calc_pvalues(query, gene_sets, background=20000, **kwargs):
    """ calculate pvalues for all categories in the graph

    :param set query: set of identifiers for which the p value is calculated
    :param dict gene_sets: gmt file dict after background was set
    :param set background: total number of genes in your annotated database.
    :returns: pvalues
              x: overlapped gene number
              n: length of gene_set which belongs to each terms
              hits: overlapped gene names.


    For 2*2 contingency table: 
    =============================================================================
                         |   in  query  |  not in query |    row total
    =>      in gene_set  |        a     |       b       |       a+b  
    =>  not in gene_set  |        c     |       d       |       c+d  
           column total                                 | a+b+c+d = anno database
    =============================================================================
    background genes number = a + b + c + d.

    Then, in R
        x=a     the number of white balls drawn without replacement 
                from an urn which contains both black and white balls.
        m=a+b   the number of white balls in the urn    
        n=c+d   the number of black balls in the urn    
        k=a+c   the number of balls drawn from the urn  
   
    In Scipy:
    for args in scipy.hypergeom.sf(k, M, n, N, loc=0):
        M: the total number of objects, 
        n: the total number of Type I objects. 
        k: the random variate represents the number of Type I objects in N drawn 
           without replacement from the total population.
    
    Therefore, these two functions are the same when using parameters from 2*2 table:
    R:     >   phyper(x-1, m, n, k, lower.tail=FALSE)
    Scipy: >>> hypergeom.sf(x-1, m+n, m, k)
     
    """

    # number of genes in your query data
    k = len(query) 
    query = set(query)
    vals = []
    # background should be all genes in annotated database
    # such as go, kegg et.al.
    if isinstance(background, set): 
        bg = len(background) # total number in your annotated database 
        # filter genes that not found in annotated database
        query = query.intersection(background)
    elif isinstance(background, int):
        bg = background
    else:
        raise ValueError("background should be set or int object")
    # pval
    subsets = sorted(gene_sets.keys())
    for s in subsets:
        category = gene_sets.get(s)
        m = len(category)
        hits = query.intersection(set(category))
        x = len(hits)
        if x < 1 : continue
        # pVal = hypergeom.sf(hitCount-1,popTotal,bgHits,queryTotal) 
        # p(X >= hitCounts)
        vals.append((s, hypergeom.sf(x-1, bg, m, k), x, m, hits))

    return zip(*vals)
def calc_pvalues(query, gene_sets, background=20000, **kwargs):
    """ calculate pvalues for all categories in the graph

    :param set query: set of identifiers for which the p value is calculated
    :param dict gene_sets: gmt file dict after background was set
    :param set background: total number of genes in your annotated database.
    :returns: pvalues
              x: overlapped gene number
              n: length of gene_set which belongs to each terms
              hits: overlapped gene names.


    For 2*2 contingency table: 
    =============================================================================
                         |   in  query  |  not in query |    row total
    =>      in gene_set  |        a     |       b       |       a+b  
    =>  not in gene_set  |        c     |       d       |       c+d  
           column total                                 | a+b+c+d = anno database
    =============================================================================
    background genes number = a + b + c + d.

    Then, in R
        x=a     the number of white balls drawn without replacement 
                from an urn which contains both black and white balls.
        m=a+b   the number of white balls in the urn    
        n=c+d   the number of black balls in the urn    
        k=a+c   the number of balls drawn from the urn  
   
    In Scipy:
    for args in scipy.hypergeom.sf(k, M, n, N, loc=0):
        M: the total number of objects, 
        n: the total number of Type I objects. 
        k: the random variate represents the number of Type I objects in N drawn 
           without replacement from the total population.
    
    Therefore, these two functions are the same when using parameters from 2*2 table:
    R:     >   phyper(x-1, m, n, k, lower.tail=FALSE)
    Scipy: >>> hypergeom.sf(x-1, m+n, m, k)
     
    """

    # number of genes in your query data
    k = len(query) 
    query = set(query)
    vals = []
    # background should be all genes in annotated database
    # such as go, kegg et.al.
    if isinstance(background, set): 
        bg = len(background) # total number in your annotated database 
        # filter genes that not found in annotated database
        query = query.intersection(background)
    elif isinstance(background, int):
        bg = background
    else:
        raise ValueError("background should be set or int object")
    # pval
    subsets = sorted(gene_sets.keys())
    for s in subsets:
        category = gene_sets.get(s)
        m = len(category)
        hits = query.intersection(set(category))
        x = len(hits)
        if x < 1 : continue
        # pVal = hypergeom.sf(hitCount-1,popTotal,bgHits,queryTotal) 
        # p(X >= hitCounts)
        vals.append((s, hypergeom.sf(x-1, bg, m, k), x, m, hits))

    return zip(*vals)
calculate pvalues for all categories in the graph

    :param set query: set of identifiers for which the p value is calculated
    :param dict gene_sets: gmt file dict after background was set
    :param set background: total number of genes in your annotated database.
    :returns: pvalues
              x: overlapped gene number
              n: length of gene_set which belongs to each terms
              hits: overlapped gene names.


    For 2*2 contingency table: 
    =============================================================================
                         |   in  query  |  not in query |    row total
    =>      in gene_set  |        a     |       b       |       a+b  
    =>  not in gene_set  |        c     |       d       |       c+d  
           column total                                 | a+b+c+d = anno database
    =============================================================================
    background genes number = a + b + c + d.

    Then, in R
        x=a     the number of white balls drawn without replacement 
                from an urn which contains both black and white balls.
        m=a+b   the number of white balls in the urn    
        n=c+d   the number of black balls in the urn    
        k=a+c   the number of balls drawn from the urn  
   
    In Scipy:
    for args in scipy.hypergeom.sf(k, M, n, N, loc=0):
        M: the total number of objects, 
        n: the total number of Type I objects. 
        k: the random variate represents the number of Type I objects in N drawn 
           without replacement from the total population.
    
    Therefore, these two functions are the same when using parameters from 2*2 table:
    R:     >   phyper(x-1, m, n, k, lower.tail=FALSE)
    Scipy: >>> hypergeom.sf(x-1, m+n, m, k)