gseapy/algorithm.py
def gsea_compute_tensor(data, gmt, n, weighted_score_type, permutation_type,
                 method, pheno_pos, pheno_neg, classes, ascending,
                 processes=1, seed=None, single=False, scale=False):
    """compute enrichment scores and enrichment nulls.

        :param data: preprocessed expression dataframe or a pre-ranked file if prerank=True.
        :param dict gmt: all gene sets in .gmt file. need to call load_gmt() to get results.
        :param int n: permutation number. default: 1000.
        :param str method: ranking_metric method. see above.
        :param str pheno_pos: one of labels of phenotype's names.
        :param str pheno_neg: one of labels of phenotype's names.
        :param list classes: a list of phenotype labels, to specify which column of dataframe belongs to what category of phenotype.
        :param float weighted_score_type: default:1
        :param bool ascending: sorting order of rankings. Default: False.
        :param seed: random seed. Default: np.random.RandomState()
        :param bool scale: if true, scale es by gene number.

        :return: a tuple contains::

                | zipped results of es, nes, pval, fdr.
                | nested list of hit indices of input gene_list.
                | nested list of ranked enrichment score of each input gene_sets.
                | list of enriched terms

    """
    w = weighted_score_type
    subsets = sorted(gmt.keys())
    rs = np.random.RandomState(seed)
    genes_mat, cor_mat = data.index.values, data.values
    base = 5 if data.shape[0] >= 5000 else 10
    block = ceil(len(subsets) / base)

    if permutation_type == "phenotype":
        # shuffling classes and generate random correlation rankings
        logging.debug("Start to permutate classes..............................")
        genes_ind = []
        cor_mat = []
        temp_rnk = []
        pool_rnk = Pool(processes=processes)
        # split large array into smaller blocks to avoid memory overflow
        i=1
        while i <= block:
            rs = np.random.RandomState(seed)
            temp_rnk.append(pool_rnk.apply_async(ranking_metric_tensor,
                                         args=(data, method, base, pheno_pos, pheno_neg, classes,
                                                ascending, rs)))
            i +=1
        pool_rnk.close()
        pool_rnk.join()

        for k, temp in enumerate(temp_rnk):
            gi, cor = temp.get()
            if k+1 == block:
               genes_ind.append(gi)
               cor_mat.append(cor)
            else:
                genes_ind.append(gi[:-1])
                cor_mat.append(cor[:-1])

        genes_ind, cor_mat = np.vstack(genes_ind), np.vstack(cor_mat)
        # convert to tuple
        genes_mat = (data.index.values, genes_ind)

    logging.debug("Start to compute es and esnulls........................")
    # Prerank, ssGSEA, GSEA
    es = []
    RES = []
    hit_ind = []
    esnull = []
    temp_esnu = []
    pool_esnu = Pool(processes=processes)
    # split large array into smaller blocks to avoid memory overflow
    i, m = 1, 0
    while i <= block:
        # you have to reseed, or all your processes are sharing the same seed value
        rs = np.random.RandomState(seed)
        gmtrim = {k: gmt.get(k) for k in subsets[m:base * i]}
        temp_esnu.append(pool_esnu.apply_async(enrichment_score_tensor,
                                               args=(genes_mat, cor_mat,
                                                     gmtrim, w, n, rs,
                                                     single, scale)))
        m = base * i
        i += 1
    pool_esnu.close()
    pool_esnu.join()
    # esn is a list, don't need to use append method.
    for si, temp in enumerate(temp_esnu):
        e, enu, hit, rune = temp.get()
        esnull.append(enu)
        es.append(e)
        RES.append(rune)
        hit_ind += hit
    # concate results
    es, esnull, RES = np.hstack(es), np.vstack(esnull), np.vstack(RES)

    return gsea_significance(es, esnull), hit_ind, RES, subsets
def gsea_compute_tensor(data, gmt, n, weighted_score_type, permutation_type,
                 method, pheno_pos, pheno_neg, classes, ascending,
                 processes=1, seed=None, single=False, scale=False):
    """compute enrichment scores and enrichment nulls.

        :param data: preprocessed expression dataframe or a pre-ranked file if prerank=True.
        :param dict gmt: all gene sets in .gmt file. need to call load_gmt() to get results.
        :param int n: permutation number. default: 1000.
        :param str method: ranking_metric method. see above.
        :param str pheno_pos: one of labels of phenotype's names.
        :param str pheno_neg: one of labels of phenotype's names.
        :param list classes: a list of phenotype labels, to specify which column of dataframe belongs to what category of phenotype.
        :param float weighted_score_type: default:1
        :param bool ascending: sorting order of rankings. Default: False.
        :param seed: random seed. Default: np.random.RandomState()
        :param bool scale: if true, scale es by gene number.

        :return: a tuple contains::

                | zipped results of es, nes, pval, fdr.
                | nested list of hit indices of input gene_list.
                | nested list of ranked enrichment score of each input gene_sets.
                | list of enriched terms

    """
    w = weighted_score_type
    subsets = sorted(gmt.keys())
    rs = np.random.RandomState(seed)
    genes_mat, cor_mat = data.index.values, data.values
    base = 5 if data.shape[0] >= 5000 else 10
    block = ceil(len(subsets) / base)

    if permutation_type == "phenotype":
        # shuffling classes and generate random correlation rankings
        logging.debug("Start to permutate classes..............................")
        genes_ind = []
        cor_mat = []
        temp_rnk = []
        pool_rnk = Pool(processes=processes)
        # split large array into smaller blocks to avoid memory overflow
        i=1
        while i <= block:
            rs = np.random.RandomState(seed)
            temp_rnk.append(pool_rnk.apply_async(ranking_metric_tensor,
                                         args=(data, method, base, pheno_pos, pheno_neg, classes,
                                                ascending, rs)))
            i +=1
        pool_rnk.close()
        pool_rnk.join()

        for k, temp in enumerate(temp_rnk):
            gi, cor = temp.get()
            if k+1 == block:
               genes_ind.append(gi)
               cor_mat.append(cor)
            else:
                genes_ind.append(gi[:-1])
                cor_mat.append(cor[:-1])

        genes_ind, cor_mat = np.vstack(genes_ind), np.vstack(cor_mat)
        # convert to tuple
        genes_mat = (data.index.values, genes_ind)

    logging.debug("Start to compute es and esnulls........................")
    # Prerank, ssGSEA, GSEA
    es = []
    RES = []
    hit_ind = []
    esnull = []
    temp_esnu = []
    pool_esnu = Pool(processes=processes)
    # split large array into smaller blocks to avoid memory overflow
    i, m = 1, 0
    while i <= block:
        # you have to reseed, or all your processes are sharing the same seed value
        rs = np.random.RandomState(seed)
        gmtrim = {k: gmt.get(k) for k in subsets[m:base * i]}
        temp_esnu.append(pool_esnu.apply_async(enrichment_score_tensor,
                                               args=(genes_mat, cor_mat,
                                                     gmtrim, w, n, rs,
                                                     single, scale)))
        m = base * i
        i += 1
    pool_esnu.close()
    pool_esnu.join()
    # esn is a list, don't need to use append method.
    for si, temp in enumerate(temp_esnu):
        e, enu, hit, rune = temp.get()
        esnull.append(enu)
        es.append(e)
        RES.append(rune)
        hit_ind += hit
    # concate results
    es, esnull, RES = np.hstack(es), np.vstack(esnull), np.vstack(RES)

    return gsea_significance(es, esnull), hit_ind, RES, subsets
def gsea_compute_tensor ( data , gmt , n , weighted_score_type , permutation_type , method , pheno_pos , pheno_neg , classes , ascending , processes = 1 , seed = None , single = False , scale = False ) : w = weighted_score_type subsets = sorted ( gmt . keys ( ) ) rs = np . random . RandomState ( seed ) genes_mat , cor_mat = data . index . values , data . values base = 5 if data . shape [ 0 ] >= 5000 else 10 block = ceil ( len ( subsets ) / base ) if permutation_type == "phenotype" : # shuffling classes and generate random correlation rankings logging . debug ( "Start to permutate classes.............................." ) genes_ind = [ ] cor_mat = [ ] temp_rnk = [ ] pool_rnk = Pool ( processes = processes ) # split large array into smaller blocks to avoid memory overflow i = 1 while i <= block : rs = np . random . RandomState ( seed ) temp_rnk . append ( pool_rnk . apply_async ( ranking_metric_tensor , args = ( data , method , base , pheno_pos , pheno_neg , classes , ascending , rs ) ) ) i += 1 pool_rnk . close ( ) pool_rnk . join ( ) for k , temp in enumerate ( temp_rnk ) : gi , cor = temp . get ( ) if k + 1 == block : genes_ind . append ( gi ) cor_mat . append ( cor ) else : genes_ind . append ( gi [ : - 1 ] ) cor_mat . append ( cor [ : - 1 ] ) genes_ind , cor_mat = np . vstack ( genes_ind ) , np . vstack ( cor_mat ) # convert to tuple genes_mat = ( data . index . values , genes_ind ) logging . debug ( "Start to compute es and esnulls........................" ) # Prerank, ssGSEA, GSEA es = [ ] RES = [ ] hit_ind = [ ] esnull = [ ] temp_esnu = [ ] pool_esnu = Pool ( processes = processes ) # split large array into smaller blocks to avoid memory overflow i , m = 1 , 0 while i <= block : # you have to reseed, or all your processes are sharing the same seed value rs = np . random . RandomState ( seed ) gmtrim = { k : gmt . get ( k ) for k in subsets [ m : base * i ] } temp_esnu . append ( pool_esnu . apply_async ( enrichment_score_tensor , args = ( genes_mat , cor_mat , gmtrim , w , n , rs , single , scale ) ) ) m = base * i i += 1 pool_esnu . close ( ) pool_esnu . join ( ) # esn is a list, don't need to use append method. for si , temp in enumerate ( temp_esnu ) : e , enu , hit , rune = temp . get ( ) esnull . append ( enu ) es . append ( e ) RES . append ( rune ) hit_ind += hit # concate results es , esnull , RES = np . hstack ( es ) , np . vstack ( esnull ) , np . vstack ( RES ) return gsea_significance ( es , esnull ) , hit_ind , RES , subsets
compute enrichment scores and enrichment nulls.

        :param data: preprocessed expression dataframe or a pre-ranked file if prerank=True.
        :param dict gmt: all gene sets in .gmt file. need to call load_gmt() to get results.
        :param int n: permutation number. default: 1000.
        :param str method: ranking_metric method. see above.
        :param str pheno_pos: one of labels of phenotype's names.
        :param str pheno_neg: one of labels of phenotype's names.
        :param list classes: a list of phenotype labels, to specify which column of dataframe belongs to what category of phenotype.
        :param float weighted_score_type: default:1
        :param bool ascending: sorting order of rankings. Default: False.
        :param seed: random seed. Default: np.random.RandomState()
        :param bool scale: if true, scale es by gene number.

        :return: a tuple contains::

                | zipped results of es, nes, pval, fdr.
                | nested list of hit indices of input gene_list.
                | nested list of ranked enrichment score of each input gene_sets.
                | list of enriched terms
compute enrichment scores and enrichment nulls .