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