gseapy/algorithm.py
def gsea_significance(enrichment_scores, enrichment_nulls):
    """Compute nominal pvals, normalized ES, and FDR q value.

        For a given NES(S) = NES* >= 0. The FDR is the ratio of the percentage of all (S,pi) with
        NES(S,pi) >= 0, whose NES(S,pi) >= NES*, divided by the percentage of
        observed S wih NES(S) >= 0, whose NES(S) >= NES*, and similarly if NES(S) = NES* <= 0.
    """
    # For a zero by zero division (undetermined, results in a NaN),
    np.seterr(divide='ignore', invalid='ignore')
    # import warnings
    # warnings.simplefilter("ignore")
    es = np.array(enrichment_scores)
    esnull = np.array(enrichment_nulls)
    logging.debug("Start to compute pvals..................................")
    # compute pvals.
    enrichmentPVals = gsea_pval(es, esnull).tolist()

    logging.debug("Compute nes and nesnull.................................")
    # nEnrichmentScores, nEnrichmentNulls = normalize(es, esnull)

    # new normalized enrichment score implementation.
    # this could speed up significantly.
    esnull_pos = (esnull*(esnull>=0)).mean(axis=1)
    esnull_neg = (esnull*(esnull<0)).mean(axis=1)
    nEnrichmentScores  = np.where(es>=0, es/esnull_pos, -es/esnull_neg)
    nEnrichmentNulls = np.where(esnull>=0, esnull/esnull_pos[:,np.newaxis],
                                          -esnull/esnull_neg[:,np.newaxis])

    logging.debug("start to compute fdrs..................................")

    # FDR null distribution histogram
    # create a histogram of all NES(S,pi) over all S and pi
    # Use this null distribution to compute an FDR q value,
    # vals = reduce(lambda x,y: x+y, nEnrichmentNulls, [])
    # nvals = np.array(sorted(vals))
    # or
    nvals = np.sort(nEnrichmentNulls.flatten())
    nnes = np.sort(nEnrichmentScores)
    fdrs = []
    # FDR computation
    for i in range(len(enrichment_scores)):
        nes = nEnrichmentScores[i]
        # use the same pval method to calculate fdr
        if nes >= 0:
            allPos = int(len(nvals) - np.searchsorted(nvals, 0, side="left"))
            allHigherAndPos = int(len(nvals) - np.searchsorted(nvals, nes, side="left"))
            nesPos = len(nnes) - int(np.searchsorted(nnes, 0, side="left"))
            nesHigherAndPos = len(nnes) - int(np.searchsorted(nnes, nes, side="left"))
            # allPos = (nvals >= 0).sum()
            # allHigherAndPos = (nvals >= nes).sum()
            # nesPos = (nnes >=0).sum()
            # nesHigherAndPos = (nnes >= nes).sum()
        else:
            allPos = int(np.searchsorted(nvals, 0, side="left"))
            allHigherAndPos = int(np.searchsorted(nvals, nes, side="right"))
            nesPos = int(np.searchsorted(nnes, 0, side="left"))
            nesHigherAndPos = int(np.searchsorted(nnes, nes, side="right"))
            # allPos = (nvals < 0).sum()
            # allHigherAndPos = (nvals < nes).sum()
            # nesPos = (nnes < 0).sum()
            # nesHigherAndPos = (nnes < nes).sum()
        
        try:
            pi_norm = allHigherAndPos/float(allPos) 
            pi_obs = nesHigherAndPos/float(nesPos)
            fdr =  pi_norm / pi_obs
            fdrs.append(fdr if fdr < 1 else 1.0)
        except:
            fdrs.append(1000000000.0)

    logging.debug("Statistical testing finished.............................")

    return zip(enrichment_scores, nEnrichmentScores, enrichmentPVals, fdrs)
def gsea_significance(enrichment_scores, enrichment_nulls):
    # For a zero by zero division (undetermined, results in a NaN),
    np.seterr(divide='ignore', invalid='ignore')
    # import warnings
    # warnings.simplefilter("ignore")
    es = np.array(enrichment_scores)
    esnull = np.array(enrichment_nulls)
    logging.debug("Start to compute pvals..................................")
    # compute pvals.
    enrichmentPVals = gsea_pval(es, esnull).tolist()

    logging.debug("Compute nes and nesnull.................................")
    # nEnrichmentScores, nEnrichmentNulls = normalize(es, esnull)

    # new normalized enrichment score implementation.
    # this could speed up significantly.
    esnull_pos = (esnull*(esnull>=0)).mean(axis=1)
    esnull_neg = (esnull*(esnull<0)).mean(axis=1)
    nEnrichmentScores  = np.where(es>=0, es/esnull_pos, -es/esnull_neg)
    nEnrichmentNulls = np.where(esnull>=0, esnull/esnull_pos[:,np.newaxis],
                                          -esnull/esnull_neg[:,np.newaxis])

    logging.debug("start to compute fdrs..................................")

    # FDR null distribution histogram
    # create a histogram of all NES(S,pi) over all S and pi
    # Use this null distribution to compute an FDR q value,
    # vals = reduce(lambda x,y: x+y, nEnrichmentNulls, [])
    # nvals = np.array(sorted(vals))
    # or
    nvals = np.sort(nEnrichmentNulls.flatten())
    nnes = np.sort(nEnrichmentScores)
    fdrs = []
    # FDR computation
    for i in range(len(enrichment_scores)):
        nes = nEnrichmentScores[i]
        # use the same pval method to calculate fdr
        if nes >= 0:
            allPos = int(len(nvals) - np.searchsorted(nvals, 0, side="left"))
            allHigherAndPos = int(len(nvals) - np.searchsorted(nvals, nes, side="left"))
            nesPos = len(nnes) - int(np.searchsorted(nnes, 0, side="left"))
            nesHigherAndPos = len(nnes) - int(np.searchsorted(nnes, nes, side="left"))
            # allPos = (nvals >= 0).sum()
            # allHigherAndPos = (nvals >= nes).sum()
            # nesPos = (nnes >=0).sum()
            # nesHigherAndPos = (nnes >= nes).sum()
        else:
            allPos = int(np.searchsorted(nvals, 0, side="left"))
            allHigherAndPos = int(np.searchsorted(nvals, nes, side="right"))
            nesPos = int(np.searchsorted(nnes, 0, side="left"))
            nesHigherAndPos = int(np.searchsorted(nnes, nes, side="right"))
            # allPos = (nvals < 0).sum()
            # allHigherAndPos = (nvals < nes).sum()
            # nesPos = (nnes < 0).sum()
            # nesHigherAndPos = (nnes < nes).sum()
        
        try:
            pi_norm = allHigherAndPos/float(allPos) 
            pi_obs = nesHigherAndPos/float(nesPos)
            fdr =  pi_norm / pi_obs
            fdrs.append(fdr if fdr < 1 else 1.0)
        except:
            fdrs.append(1000000000.0)

    logging.debug("Statistical testing finished.............................")

    return zip(enrichment_scores, nEnrichmentScores, enrichmentPVals, fdrs)
def gsea_significance ( enrichment_scores , enrichment_nulls ) : # For a zero by zero division (undetermined, results in a NaN), np . seterr ( divide = 'ignore' , invalid = 'ignore' ) # import warnings # warnings.simplefilter("ignore") es = np . array ( enrichment_scores ) esnull = np . array ( enrichment_nulls ) logging . debug ( "Start to compute pvals.................................." ) # compute pvals. enrichmentPVals = gsea_pval ( es , esnull ) . tolist ( ) logging . debug ( "Compute nes and nesnull................................." ) # nEnrichmentScores, nEnrichmentNulls = normalize(es, esnull) # new normalized enrichment score implementation. # this could speed up significantly. esnull_pos = ( esnull * ( esnull >= 0 ) ) . mean ( axis = 1 ) esnull_neg = ( esnull * ( esnull < 0 ) ) . mean ( axis = 1 ) nEnrichmentScores = np . where ( es >= 0 , es / esnull_pos , - es / esnull_neg ) nEnrichmentNulls = np . where ( esnull >= 0 , esnull / esnull_pos [ : , np . newaxis ] , - esnull / esnull_neg [ : , np . newaxis ] ) logging . debug ( "start to compute fdrs.................................." ) # FDR null distribution histogram # create a histogram of all NES(S,pi) over all S and pi # Use this null distribution to compute an FDR q value, # vals = reduce(lambda x,y: x+y, nEnrichmentNulls, []) # nvals = np.array(sorted(vals)) # or nvals = np . sort ( nEnrichmentNulls . flatten ( ) ) nnes = np . sort ( nEnrichmentScores ) fdrs = [ ] # FDR computation for i in range ( len ( enrichment_scores ) ) : nes = nEnrichmentScores [ i ] # use the same pval method to calculate fdr if nes >= 0 : allPos = int ( len ( nvals ) - np . searchsorted ( nvals , 0 , side = "left" ) ) allHigherAndPos = int ( len ( nvals ) - np . searchsorted ( nvals , nes , side = "left" ) ) nesPos = len ( nnes ) - int ( np . searchsorted ( nnes , 0 , side = "left" ) ) nesHigherAndPos = len ( nnes ) - int ( np . searchsorted ( nnes , nes , side = "left" ) ) # allPos = (nvals >= 0).sum() # allHigherAndPos = (nvals >= nes).sum() # nesPos = (nnes >=0).sum() # nesHigherAndPos = (nnes >= nes).sum() else : allPos = int ( np . searchsorted ( nvals , 0 , side = "left" ) ) allHigherAndPos = int ( np . searchsorted ( nvals , nes , side = "right" ) ) nesPos = int ( np . searchsorted ( nnes , 0 , side = "left" ) ) nesHigherAndPos = int ( np . searchsorted ( nnes , nes , side = "right" ) ) # allPos = (nvals < 0).sum() # allHigherAndPos = (nvals < nes).sum() # nesPos = (nnes < 0).sum() # nesHigherAndPos = (nnes < nes).sum() try : pi_norm = allHigherAndPos / float ( allPos ) pi_obs = nesHigherAndPos / float ( nesPos ) fdr = pi_norm / pi_obs fdrs . append ( fdr if fdr < 1 else 1.0 ) except : fdrs . append ( 1000000000.0 ) logging . debug ( "Statistical testing finished............................." ) return zip ( enrichment_scores , nEnrichmentScores , enrichmentPVals , fdrs )
Compute nominal pvals, normalized ES, and FDR q value.

        For a given NES(S) = NES* >= 0. The FDR is the ratio of the percentage of all (S,pi) with
        NES(S,pi) >= 0, whose NES(S,pi) >= NES*, divided by the percentage of
        observed S wih NES(S) >= 0, whose NES(S) >= NES*, and similarly if NES(S) = NES* <= 0.
Compute nominal pvals normalized ES and FDR q value .