gseapy/algorithm.py
def enrichment_score_tensor(gene_mat, cor_mat, gene_sets, weighted_score_type, nperm=1000,
                            rs=np.random.RandomState(), single=False, scale=False):
    """Next generation algorithm of GSEA and ssGSEA.

        :param gene_mat:        the ordered gene list(vector) with or without gene indices matrix.
        :param cor_mat:         correlation vector or matrix  (e.g. signal to noise scores)
                                corresponding to the genes in the gene list or matrix.
        :param dict gene_sets:  gmt file dict.
        :param float weighted_score_type:     weighting by the correlation.
                                options: 0(classic), 1, 1.5, 2. default:1 for GSEA and 0.25 for ssGSEA.
        :param int nperm:       permutation times.
        :param bool scale:      If True, normalize the scores by number of genes_mat.
        :param bool single:     If True, use ssGSEA algorithm, otherwise use GSEA.
        :param rs:              Random state for initialize gene list shuffling.
                                Default: np.random.RandomState(seed=None)
        :return: a tuple contains::

                 | ES: Enrichment score (real number between -1 and +1), for ssGSEA, set scale eq to True.
                 | ESNULL: Enrichment score calculated from random permutation.
                 | Hits_Indices: Indices of genes if genes are included in gene_set.
                 | RES: The running enrichment score for all locations in the gene list.

    """
    # gene_mat -> 1d: prerank, ssSSEA or 2d: GSEA
    keys = sorted(gene_sets.keys())

    if weighted_score_type == 0:
        # don't bother doing calcuation, just set to 1
        cor_mat = np.ones(cor_mat.shape)
    elif weighted_score_type > 0:
        pass
    else:
        logging.error("Using negative values of weighted_score_type, not allowed")
        sys.exit(0)

    cor_mat = np.abs(cor_mat)
    if cor_mat.ndim ==1:
        # ssGSEA or Prerank
        # genestes->M, genes->N, perm-> axis=2
        N, M = len(gene_mat), len(keys)
        # generate gene hits matrix
        # for 1d ndarray of gene_mat, set assume_unique=True,
        # means the input arrays are both assumed to be unique,
        # which can speed up the calculation.
        tag_indicator = np.vstack([np.in1d(gene_mat, gene_sets[key], assume_unique=True) for key in keys])
        tag_indicator = tag_indicator.astype(int)
        # index of hits
        hit_ind = [ np.flatnonzero(tag).tolist() for tag in tag_indicator ]
        # generate permutated hits matrix
        perm_tag_tensor = np.repeat(tag_indicator, nperm+1).reshape((M,N,nperm+1))
        # shuffle matrix, last matrix is not shuffled when nperm > 0
        if nperm: np.apply_along_axis(lambda x: np.apply_along_axis(rs.shuffle,0,x),1, perm_tag_tensor[:,:,:-1])
        # missing hits
        no_tag_tensor = 1 - perm_tag_tensor
        # calculate numerator, denominator of each gene hits
        rank_alpha = (perm_tag_tensor*cor_mat[np.newaxis,:,np.newaxis])** weighted_score_type

    elif cor_mat.ndim == 2:
        # GSEA
        # 2d ndarray, gene_mat and cor_mat are shuffled already
        # reshape matrix
        cor_mat = cor_mat.T
        # gene_mat is a tuple contains (gene_name, permuate_gene_name_indices)
        genes, genes_ind = gene_mat
        # genestes->M, genes->N, perm-> axis=2
        # don't use assume_unique=True in 2d array when use np.isin().
        # elements in gene_mat are not unique, or will cause unwanted results
        tag_indicator = np.vstack([np.in1d(genes, gene_sets[key], assume_unique=True) for key in keys])
        tag_indicator = tag_indicator.astype(int)
        perm_tag_tensor = np.stack([tag.take(genes_ind).T for tag in tag_indicator], axis=0)
        #index of hits
        hit_ind = [ np.flatnonzero(tag).tolist() for tag in perm_tag_tensor[:,:,-1] ]
        # nohits
        no_tag_tensor = 1 - perm_tag_tensor
        # calculate numerator, denominator of each gene hits
        rank_alpha = (perm_tag_tensor*cor_mat[np.newaxis,:,:])** weighted_score_type
    else:
        logging.error("Program die because of unsupported input")
        sys.exit(0)

    # Nhint = tag_indicator.sum(1)
    # Nmiss =  N - Nhint
    axis=1
    P_GW_denominator = np.sum(rank_alpha, axis=axis, keepdims=True)
    P_NG_denominator = np.sum(no_tag_tensor, axis=axis, keepdims=True)
    REStensor = np.cumsum(rank_alpha / P_GW_denominator - no_tag_tensor / P_NG_denominator, axis=axis)
    # ssGSEA: scale es by gene numbers ?
    # https://gist.github.com/gaoce/39e0907146c752c127728ad74e123b33
    if scale: REStensor = REStensor / len(gene_mat)
    if single:
        #ssGSEA
        esmatrix = REStensor.sum(axis=axis)
    else:
        #GSEA
        esmax, esmin = REStensor.max(axis=axis), REStensor.min(axis=axis)
        esmatrix = np.where(np.abs(esmax)>np.abs(esmin), esmax, esmin)

    es, esnull, RES = esmatrix[:,-1], esmatrix[:,:-1], REStensor[:,:,-1]

    return es, esnull, hit_ind, RES
def enrichment_score_tensor(gene_mat, cor_mat, gene_sets, weighted_score_type, nperm=1000,
                            rs=np.random.RandomState(), single=False, scale=False):
    """Next generation algorithm of GSEA and ssGSEA.

        :param gene_mat:        the ordered gene list(vector) with or without gene indices matrix.
        :param cor_mat:         correlation vector or matrix  (e.g. signal to noise scores)
                                corresponding to the genes in the gene list or matrix.
        :param dict gene_sets:  gmt file dict.
        :param float weighted_score_type:     weighting by the correlation.
                                options: 0(classic), 1, 1.5, 2. default:1 for GSEA and 0.25 for ssGSEA.
        :param int nperm:       permutation times.
        :param bool scale:      If True, normalize the scores by number of genes_mat.
        :param bool single:     If True, use ssGSEA algorithm, otherwise use GSEA.
        :param rs:              Random state for initialize gene list shuffling.
                                Default: np.random.RandomState(seed=None)
        :return: a tuple contains::

                 | ES: Enrichment score (real number between -1 and +1), for ssGSEA, set scale eq to True.
                 | ESNULL: Enrichment score calculated from random permutation.
                 | Hits_Indices: Indices of genes if genes are included in gene_set.
                 | RES: The running enrichment score for all locations in the gene list.

    """
    # gene_mat -> 1d: prerank, ssSSEA or 2d: GSEA
    keys = sorted(gene_sets.keys())

    if weighted_score_type == 0:
        # don't bother doing calcuation, just set to 1
        cor_mat = np.ones(cor_mat.shape)
    elif weighted_score_type > 0:
        pass
    else:
        logging.error("Using negative values of weighted_score_type, not allowed")
        sys.exit(0)

    cor_mat = np.abs(cor_mat)
    if cor_mat.ndim ==1:
        # ssGSEA or Prerank
        # genestes->M, genes->N, perm-> axis=2
        N, M = len(gene_mat), len(keys)
        # generate gene hits matrix
        # for 1d ndarray of gene_mat, set assume_unique=True,
        # means the input arrays are both assumed to be unique,
        # which can speed up the calculation.
        tag_indicator = np.vstack([np.in1d(gene_mat, gene_sets[key], assume_unique=True) for key in keys])
        tag_indicator = tag_indicator.astype(int)
        # index of hits
        hit_ind = [ np.flatnonzero(tag).tolist() for tag in tag_indicator ]
        # generate permutated hits matrix
        perm_tag_tensor = np.repeat(tag_indicator, nperm+1).reshape((M,N,nperm+1))
        # shuffle matrix, last matrix is not shuffled when nperm > 0
        if nperm: np.apply_along_axis(lambda x: np.apply_along_axis(rs.shuffle,0,x),1, perm_tag_tensor[:,:,:-1])
        # missing hits
        no_tag_tensor = 1 - perm_tag_tensor
        # calculate numerator, denominator of each gene hits
        rank_alpha = (perm_tag_tensor*cor_mat[np.newaxis,:,np.newaxis])** weighted_score_type

    elif cor_mat.ndim == 2:
        # GSEA
        # 2d ndarray, gene_mat and cor_mat are shuffled already
        # reshape matrix
        cor_mat = cor_mat.T
        # gene_mat is a tuple contains (gene_name, permuate_gene_name_indices)
        genes, genes_ind = gene_mat
        # genestes->M, genes->N, perm-> axis=2
        # don't use assume_unique=True in 2d array when use np.isin().
        # elements in gene_mat are not unique, or will cause unwanted results
        tag_indicator = np.vstack([np.in1d(genes, gene_sets[key], assume_unique=True) for key in keys])
        tag_indicator = tag_indicator.astype(int)
        perm_tag_tensor = np.stack([tag.take(genes_ind).T for tag in tag_indicator], axis=0)
        #index of hits
        hit_ind = [ np.flatnonzero(tag).tolist() for tag in perm_tag_tensor[:,:,-1] ]
        # nohits
        no_tag_tensor = 1 - perm_tag_tensor
        # calculate numerator, denominator of each gene hits
        rank_alpha = (perm_tag_tensor*cor_mat[np.newaxis,:,:])** weighted_score_type
    else:
        logging.error("Program die because of unsupported input")
        sys.exit(0)

    # Nhint = tag_indicator.sum(1)
    # Nmiss =  N - Nhint
    axis=1
    P_GW_denominator = np.sum(rank_alpha, axis=axis, keepdims=True)
    P_NG_denominator = np.sum(no_tag_tensor, axis=axis, keepdims=True)
    REStensor = np.cumsum(rank_alpha / P_GW_denominator - no_tag_tensor / P_NG_denominator, axis=axis)
    # ssGSEA: scale es by gene numbers ?
    # https://gist.github.com/gaoce/39e0907146c752c127728ad74e123b33
    if scale: REStensor = REStensor / len(gene_mat)
    if single:
        #ssGSEA
        esmatrix = REStensor.sum(axis=axis)
    else:
        #GSEA
        esmax, esmin = REStensor.max(axis=axis), REStensor.min(axis=axis)
        esmatrix = np.where(np.abs(esmax)>np.abs(esmin), esmax, esmin)

    es, esnull, RES = esmatrix[:,-1], esmatrix[:,:-1], REStensor[:,:,-1]

    return es, esnull, hit_ind, RES
def enrichment_score_tensor ( gene_mat , cor_mat , gene_sets , weighted_score_type , nperm = 1000 , rs = np . random . RandomState ( ) , single = False , scale = False ) : # gene_mat -> 1d: prerank, ssSSEA or 2d: GSEA keys = sorted ( gene_sets . keys ( ) ) if weighted_score_type == 0 : # don't bother doing calcuation, just set to 1 cor_mat = np . ones ( cor_mat . shape ) elif weighted_score_type > 0 : pass else : logging . error ( "Using negative values of weighted_score_type, not allowed" ) sys . exit ( 0 ) cor_mat = np . abs ( cor_mat ) if cor_mat . ndim == 1 : # ssGSEA or Prerank # genestes->M, genes->N, perm-> axis=2 N , M = len ( gene_mat ) , len ( keys ) # generate gene hits matrix # for 1d ndarray of gene_mat, set assume_unique=True, # means the input arrays are both assumed to be unique, # which can speed up the calculation. tag_indicator = np . vstack ( [ np . in1d ( gene_mat , gene_sets [ key ] , assume_unique = True ) for key in keys ] ) tag_indicator = tag_indicator . astype ( int ) # index of hits hit_ind = [ np . flatnonzero ( tag ) . tolist ( ) for tag in tag_indicator ] # generate permutated hits matrix perm_tag_tensor = np . repeat ( tag_indicator , nperm + 1 ) . reshape ( ( M , N , nperm + 1 ) ) # shuffle matrix, last matrix is not shuffled when nperm > 0 if nperm : np . apply_along_axis ( lambda x : np . apply_along_axis ( rs . shuffle , 0 , x ) , 1 , perm_tag_tensor [ : , : , : - 1 ] ) # missing hits no_tag_tensor = 1 - perm_tag_tensor # calculate numerator, denominator of each gene hits rank_alpha = ( perm_tag_tensor * cor_mat [ np . newaxis , : , np . newaxis ] ) ** weighted_score_type elif cor_mat . ndim == 2 : # GSEA # 2d ndarray, gene_mat and cor_mat are shuffled already # reshape matrix cor_mat = cor_mat . T # gene_mat is a tuple contains (gene_name, permuate_gene_name_indices) genes , genes_ind = gene_mat # genestes->M, genes->N, perm-> axis=2 # don't use assume_unique=True in 2d array when use np.isin(). # elements in gene_mat are not unique, or will cause unwanted results tag_indicator = np . vstack ( [ np . in1d ( genes , gene_sets [ key ] , assume_unique = True ) for key in keys ] ) tag_indicator = tag_indicator . astype ( int ) perm_tag_tensor = np . stack ( [ tag . take ( genes_ind ) . T for tag in tag_indicator ] , axis = 0 ) #index of hits hit_ind = [ np . flatnonzero ( tag ) . tolist ( ) for tag in perm_tag_tensor [ : , : , - 1 ] ] # nohits no_tag_tensor = 1 - perm_tag_tensor # calculate numerator, denominator of each gene hits rank_alpha = ( perm_tag_tensor * cor_mat [ np . newaxis , : , : ] ) ** weighted_score_type else : logging . error ( "Program die because of unsupported input" ) sys . exit ( 0 ) # Nhint = tag_indicator.sum(1) # Nmiss =  N - Nhint axis = 1 P_GW_denominator = np . sum ( rank_alpha , axis = axis , keepdims = True ) P_NG_denominator = np . sum ( no_tag_tensor , axis = axis , keepdims = True ) REStensor = np . cumsum ( rank_alpha / P_GW_denominator - no_tag_tensor / P_NG_denominator , axis = axis ) # ssGSEA: scale es by gene numbers ? # https://gist.github.com/gaoce/39e0907146c752c127728ad74e123b33 if scale : REStensor = REStensor / len ( gene_mat ) if single : #ssGSEA esmatrix = REStensor . sum ( axis = axis ) else : #GSEA esmax , esmin = REStensor . max ( axis = axis ) , REStensor . min ( axis = axis ) esmatrix = np . where ( np . abs ( esmax ) > np . abs ( esmin ) , esmax , esmin ) es , esnull , RES = esmatrix [ : , - 1 ] , esmatrix [ : , : - 1 ] , REStensor [ : , : , - 1 ] return es , esnull , hit_ind , RES
Next generation algorithm of GSEA and ssGSEA.

        :param gene_mat:        the ordered gene list(vector) with or without gene indices matrix.
        :param cor_mat:         correlation vector or matrix  (e.g. signal to noise scores)
                                corresponding to the genes in the gene list or matrix.
        :param dict gene_sets:  gmt file dict.
        :param float weighted_score_type:     weighting by the correlation.
                                options: 0(classic), 1, 1.5, 2. default:1 for GSEA and 0.25 for ssGSEA.
        :param int nperm:       permutation times.
        :param bool scale:      If True, normalize the scores by number of genes_mat.
        :param bool single:     If True, use ssGSEA algorithm, otherwise use GSEA.
        :param rs:              Random state for initialize gene list shuffling.
                                Default: np.random.RandomState(seed=None)
        :return: a tuple contains::

                 | ES: Enrichment score (real number between -1 and +1), for ssGSEA, set scale eq to True.
                 | ESNULL: Enrichment score calculated from random permutation.
                 | Hits_Indices: Indices of genes if genes are included in gene_set.
                 | RES: The running enrichment score for all locations in the gene list.
Next generation algorithm of GSEA and ssGSEA .