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):
    """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)
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.