"""Enrichment Analysis Module. Contains different functions to perform enrichment
analysis.
Most things in this module are covered in https://www.youtube.com/watch?v=2NC1QOXmc5o
by Lars Juhl Jensen.
"""
from __future__ import annotations
import logging
import os
import re
import gseapy as gp
import pandas as pd
from pandera.typing.pandas import DataFrame
from acore.enrichment_analysis import permutations_enrichment
from acore.enrichment_analysis.annotate import annotate_features
from acore.enrichment_analysis.statistical_tests.fisher import run_fisher
from acore.enrichment_analysis.statistical_tests.kolmogorov_smirnov import (
run_kolmogorov_smirnov,
)
from acore.multiple_testing import apply_pvalue_correction
from acore.types.enrichment_analysis import EnrichmentAnalysisSchema
logger = logging.getLogger(__name__)
TYPE_COLS_MSG = """
columns: 'terms', 'identifiers', 'foreground',
'background', foreground_pop, background_pop, 'pvalue', 'padj' and 'rejected'.
"""
__all__ = [
"run_fisher",
"run_gsea",
"run_kolmogorov_smirnov",
"run_site_regulation_enrichment",
"run_up_down_regulation_enrichment",
]
# ! undocumented for now (find usage example)
[docs]
def run_site_regulation_enrichment(
regulation_data: pd.DataFrame,
annotation: pd.DataFrame,
identifier: str = "identifier",
groups: list[str] = ("group1", "group2"),
annotation_col: str = "annotation",
rejected_col: str = "rejected",
group_col: str = "group",
method: str = "fisher",
regex: str = "(\\w+~.+)_\\w\\d+\\-\\w+",
correction: str = "fdr_bh",
remove_duplicates: bool = False,
) -> DataFrame[EnrichmentAnalysisSchema]:
r"""
This function runs a simple enrichment analysis for significantly
regulated protein sites in a dataset.
:param regulation_data: pandas.DataFrame resulting from differential
regulation analysis.
:param annotation: pandas.DataFrame with annotations for features
(columns: 'annotation', 'identifier' (feature identifiers), and 'source').
:param str identifier: name of the column from annotation containing
feature identifiers.
:param list groups: column names from regulation_data containing
group identifiers.
:param str annotation_col: name of the column from annotation containing
annotation terms.
:param str rejected_col: name of the column from regulation_data containing
boolean for rejected null hypothesis.
:param str group_col: column name for new column in annotation dataframe
determining if feature belongs to foreground or background.
:param str method: method used to compute enrichment
(only 'fisher' is supported currently).
:param str regex: how to extract the annotated identifier from the site identifier
:return: pandas.DataFrame with columns: 'terms', 'identifiers', 'foreground',
'background', foreground_pop, background_pop, 'pvalue', 'padj' and 'rejected'.
:raises ValueError: if regulation_data is `None` or empty.
Example::
result = run_site_regulation_enrichment(regulation_data,
annotation,
identifier='identifier',
groups=['group1', 'group2'],
annotation_col='annotation',
rejected_col='rejected',
group_col='group',
method='fisher',
match="(\\w+~.+)_\\w\\d+\\-\\w+"
)
"""
result = pd.DataFrame()
if regulation_data is None or regulation_data.empty:
raise ValueError("regulation_data is empty")
new_ids = []
# find any identifiers with a PTM and save only prot+gene identifer
for ident in regulation_data[identifier].tolist():
match = re.search(regex, ident)
if match is not None:
new_ids.append(
match.group(1)
) # removes the PTM extension of the identifier of CKG
else:
new_ids.append(ident)
# so this is normalizing the identifiers to ignore the PTM extension
regulation_data[identifier] = new_ids # matches are used as identifiers
if remove_duplicates:
regulation_data = regulation_data.drop_duplicates(subset=[identifier])
result = run_regulation_enrichment(
regulation_data=regulation_data,
annotation=annotation,
identifier=identifier,
group_col=groups,
annotation_col=annotation_col,
rejected_col=rejected_col,
method=method,
correction=correction,
)
return result
[docs]
def run_up_down_regulation_enrichment(
regulation_data: pd.DataFrame,
annotation: pd.DataFrame,
identifier: str = "identifier",
groups: list[str] = ("group1", "group2"),
annotation_col: str = "annotation",
# rejected_col: str = "rejected", # could be passed
pval_col: str = "pval",
group_col: str = "group",
log2fc_col: str = "log2FC",
method: str = "fisher",
min_detected_in_set: int = 2,
correction: str = "fdr_bh",
correction_alpha: float = 0.05,
lfc_cutoff: float = 1,
permutations: int = 0,
seed: int | None = None,
) -> DataFrame[EnrichmentAnalysisSchema]:
"""
This function runs a simple enrichment analysis for significantly regulated proteins
distinguishing between upregulation per group defined by the `groups`.
Only overrepresentation is tested using one-sided Fisher's exact tests.
:param pandas.DataFrame regulation_data: pandas.DataFrame resulting from differential regulation
analysis (CKG's regulation table).
:param pandas.DataFrame annotation: pandas.DataFrame with annotations for features
(columns: 'annotation', 'identifier' (feature identifiers), and 'source').
:param str identifier: name of the column from annotation containing feature identifiers.
:param list[str] groups: column names from regulation_data containing group identifiers.
See `pandas.DataFrame.groupby`_ for more information.
.. _pandas.DataFrame.groupby: https://pandas.pydata.org/pandas-docs/stable/\
reference/api/pandas.DataFrame.groupby.html
:param str annotation_col: name of the column from annotation containing annotation terms.
:param str rejected_col: name of the column from regulation_data containing boolean for
rejected null hypothesis.
:param str group_col: column name for new column in annotation dataframe determining
if feature belongs to foreground or background.
:param str method: method used to compute enrichment
(only 'fisher' is supported currently).
:param str correction: method to be used for multiple-testing correction.
Use ``'efdr'`` to apply empirical FDR via permutation (requires
``permutations > 0``).
:param float alpha: adjusted p-value cutoff to define significance
:param float lfc_cutoff: log fold-change cutoff to define practical significance
:param int permutations: number of permutations used to estimate the empirical FDR
when ``correction='efdr'``. Ignored for other correction methods.
:param int | None seed: seed of the random number generator used for the
permutations when ``correction='efdr'``, for reproducible results.
:return: DataFrame adhering to EnrichmentAnalysisSchema
:rtype: DataFrame[EnrichmentAnalysisSchema]
Example::
result = run_up_down_regulation_enrichment(
regulation_data,
annotation,
identifier='identifier',
groups=['group1',
'group2'],
annotation_col='annotation',
rejected_col='rejected',
group_col='group',
method='fisher',
correction='fdr_bh',
alpha=0.05,
lfc_cutoff=1,
)
"""
if isinstance(groups, str):
groups = [groups]
if isinstance(groups, tuple):
groups = list(groups)
if len(groups) != 2:
raise ValueError("groups should contains exactly two columns.")
ret = []
# In case of multiple comparisons this is used to get all possible combinations
for g1, g2 in regulation_data.groupby(groups).groups:
df = regulation_data.groupby(groups).get_group((g1, g2))
df["up_pairwise_regulation"] = (df[pval_col] <= correction_alpha) & (
df[log2fc_col] >= lfc_cutoff
)
df["down_pairwise_regulation"] = (df[pval_col] <= correction_alpha) & (
df[log2fc_col] <= -lfc_cutoff
)
comparison_tag = str(g1) + "~" + str(g2)
if not regulation_data[identifier].is_unique:
logger.warning(
"Column '%s' in regulation_data contains duplicated values for comparison %s.",
identifier,
comparison_tag,
)
for rej_col, direction in zip(
("up_pairwise_regulation", "down_pairwise_regulation"),
(f"upregulated in {g1}", f"upregulated in {g2}"),
):
_enrichment = run_regulation_enrichment(
df,
annotation,
identifier=identifier,
annotation_col=annotation_col,
rejected_col=rej_col,
group_col=group_col,
method=method,
min_detected_in_set=min_detected_in_set,
correction=correction,
correction_alpha=correction_alpha,
permutations=permutations,
seed=seed,
)
if _enrichment.empty:
logger.warning(
"No significant enrichment found for comparison %s (%s).",
comparison_tag,
direction,
)
_enrichment["direction"] = direction
_enrichment["comparison"] = comparison_tag
ret.append(_enrichment)
ret = pd.concat(ret)
if not ret.empty:
ret["rejected"] = ret["rejected"].astype(bool)
# whole numbers only (e.g. eFDR of 0 or 1) are converted to int upstream
ret[["pvalue", "padj"]] = ret[["pvalue", "padj"]].astype("Float64")
ret = EnrichmentAnalysisSchema.validate(ret)
return ret
def run_regulation_enrichment(
regulation_data: pd.DataFrame,
annotation: pd.DataFrame,
identifier: str = "identifier",
annotation_col: str = "annotation",
rejected_col: str = "rejected",
group_col: str = "group",
method: str = "fisher",
min_detected_in_set: int = 2,
correction: str = "fdr_bh",
correction_alpha: float = 0.05,
permutations: int = 0,
seed: int | None = None,
) -> pd.DataFrame:
"""
This function runs a simple enrichment analysis for significantly regulated features
in a dataset. Only overrepresentation is tested using one-sided Fisher's exact tests.
:param regulation_data: pandas.DataFrame resulting from differential regulation analysis.
:param annotation: pandas.DataFrame with annotations for features
(columns: 'annotation', 'identifier' (feature identifiers), and 'source').
:param str identifier: name of the column from annotation containing feature identifiers.
It should be present both in `regulation_data` and `annotation`. In `regulation_data`
it should be unique, while in `annotation` it can contain duplicates as one
identifier can be part of multiple pathways.
:param str annotation_col: name of the column from annotation containing annotation terms.
:param str rejected_col: name of the column from `regulation_data` containing boolean for
rejected null hypothesis.
:param str group_col: column name for new column in annotation dataframe determining
if feature belongs to foreground or background.
:param str method: method used to compute enrichment (only 'fisher' is supported currently).
:param str correction: method to be used for multiple-testing correction.
Use ``'efdr'`` to apply empirical FDR via permutation (requires
``permutations > 0``).
:param int permutations: number of permutations used to estimate the empirical FDR
when ``correction='efdr'``. Ignored for other correction methods.
:param int | None seed: seed of the random number generator used for the
permutations when ``correction='efdr'``, for reproducible results.
:return: pandas.DataFrame with columns: 'terms', 'identifiers', 'foreground',
'background', 'foreground_pop', 'background_pop', 'pvalue', 'padj' and 'rejected'.
Example::
result = run_regulation_enrichment(
regulation_data,
annotation,
identifier='identifier',
annotation_col='annotation',
rejected_col='rejected',
group_col='group',
method='fisher',
min_detected_in_set=2,
correction='fdr_bh',
correction_alpha=0.05,
)
"""
# ? can we remove NA features in that column?
if regulation_data[rejected_col].isna().any():
raise ValueError(f"Rejected column '{rejected_col}' contains missing values.")
mask_rejected = regulation_data[rejected_col].astype(bool)
if not regulation_data[identifier].is_unique:
raise ValueError(f"Column '{identifier}' in regulation_data has to be unique.")
foreground_list = regulation_data.loc[mask_rejected, identifier]
background_list = regulation_data.loc[~mask_rejected, identifier]
foreground_pop = len(foreground_list)
background_pop = len(regulation_data[identifier])
# needs to allow for missing annotations
# ! this step needs unique identifiers in the regulation_data
# group_col contains either 'foreground', 'background' or NA
annotation[group_col] = annotate_features(
features=annotation[identifier],
in_foreground=foreground_list,
in_background=background_list,
)
annotation = annotation.dropna(subset=[group_col])
result = run_enrichment(
annotation,
foreground_id="foreground",
background_id="background",
foreground_pop=foreground_pop,
background_pop=background_pop,
annotation_col=annotation_col,
group_col=group_col,
identifier_col=identifier,
method=method,
correction=correction,
min_detected_in_set=min_detected_in_set,
correction_alpha=correction_alpha,
permutations=permutations,
seed=seed,
).convert_dtypes(convert_boolean=False)
return result
def run_enrichment(
data: pd.DataFrame,
foreground_id: str,
background_id: str,
foreground_pop: int,
background_pop: int,
min_detected_in_set: int = 2,
annotation_col: str = "annotation",
group_col: str = "group",
identifier_col: str = "identifier",
method: str = "fisher",
correction: str = "fdr_bh",
correction_alpha: float = 0.05,
permutations: int = 0,
seed: int | None = None,
) -> pd.DataFrame:
"""
Computes overrepresentation of the foreground relative to a given background,
using Fisher's exact test with ``alternative='greater'``, and corrects for
multiple hypothesis testing.
:param data: pandas.DataFrame with annotations for dataset features
(columns: 'annotation', 'identifier', 'group').
:param str foreground_id: group identifier of features that belong to the foreground.
:param str background_id: group identifier of features that belong to the background.
:param int foreground_pop: number of features in the foreground.
:param int background_pop: number of features in the background.
:param int min_detected_in_set: minimum number of features in the foreground
:param str annotation_col: name of the column containing annotation terms.
:param str group_col: name of column containing the group identifiers,
e.g. specifying belonging to 'foreground' or 'background'.
:param str identifier_col: name of column containing dependent variables identifiers.
:param str method: method used to compute enrichment
(only 'fisher' is supported currently).
:param str correction: method to be used for multiple-testing correction.
Use ``'efdr'`` to apply empirical FDR via permutation (requires
``permutations > 0``).
See `statsmodels.stats.multitest.multipletests`_ for other supported methods.
.. _statsmodels.stats.multitest.multipletests: https://www.statsmodels.org/dev/\
generated/statsmodels.stats.multitest.multipletests.html
:param float correction_alpha: adjusted p-value cutoff to define significance.
:param int permutations: number of permutations used to estimate the empirical FDR
when ``correction='efdr'``. Ignored for other correction methods.
:param int | None seed: seed of the random number generator used for the
permutations when ``correction='efdr'``, for reproducible results.
Must be greater than 0 when ``correction='efdr'``.
:return: pandas.DataFrame with columns: annotation terms, features,
number of foreground/background features in each term,
p-values and corrected p-values
Columns are: 'terms', 'identifiers',
'foreground', 'background', 'foreground_pop', 'background_pop',
'pvalue', 'padj' and 'rejected'.
:raises ValueError: if ``correction='efdr'`` and ``permutations`` is not > 0.
Example::
result = run_enrichment(
data,
foreground='foreground',
background='background',
foreground_pop=len(foreground_list),
background_pop=len(background_list),
annotation_col='annotation',
group_col='group',
identifier_col='identifier',
method='fisher',
)
Example with eFDR::
result = run_enrichment(
data,
foreground='foreground',
background='background',
foreground_pop=len(foreground_list),
background_pop=len(background_list),
annotation_col='annotation',
group_col='group',
identifier_col='identifier',
method='fisher',
correction='efdr',
permutations=100,
)
"""
if method != "fisher":
raise ValueError("Only Fisher's exact test is supported at the moment.")
if correction == "efdr" and permutations <= 0:
raise ValueError(
"permutations must be greater than 0 when correction='efdr'. "
"Left at zero so it is set explicitly to e.g. 10000."
)
terms = []
ids = []
pvalues = []
fnum = []
bnum = []
# Aggregate only the identifier column, without copying unrelated metadata.
counts = data.groupby([annotation_col, group_col], observed=True)[
identifier_col
].count()
counts = counts[
counts.index.get_level_values(group_col).isin([foreground_id, background_id])
]
foreground_ids = (
data.loc[data[group_col] == foreground_id, [annotation_col, identifier_col]]
.groupby(annotation_col, observed=True)[identifier_col]
.agg(",".join)
)
for annotation in foreground_ids.index:
num_foreground = int(counts.get((annotation, foreground_id), 0))
num_background = int(counts.get((annotation, background_id), 0))
if num_foreground < min_detected_in_set:
continue
_, pvalue = run_fisher(
[num_foreground, foreground_pop - num_foreground],
[num_background, background_pop - foreground_pop - num_background],
alternative="greater",
)
fnum.append(num_foreground)
bnum.append(num_background)
terms.append(annotation)
pvalues.append(pvalue)
ids.append(foreground_ids.get(annotation, ""))
if len(pvalues) >= 1:
if correction == "efdr":
rejected, padj = permutations_enrichment._compute_efdr_for_enrichment(
data=data.loc[
data[group_col].isin([foreground_id, background_id]),
[annotation_col, identifier_col],
],
observed_pvalues=pvalues,
foreground_pop=foreground_pop,
background_pop=background_pop,
min_detected_in_set=min_detected_in_set,
annotation_col=annotation_col,
identifier_col=identifier_col,
permutations=permutations,
alpha=correction_alpha,
seed=seed,
)
else:
rejected, padj = apply_pvalue_correction(
pvalues,
alpha=correction_alpha,
method=correction,
)
result = pd.DataFrame(
{
"terms": terms,
"identifiers": ids,
"foreground": fnum,
"background": bnum,
"foreground_pop": foreground_pop, # no. of foreground features, constant
"background_pop": background_pop, # no. of included features, constant
"pvalue": pvalues,
"padj": padj,
"rejected": rejected.astype(bool),
}
)
result = result.sort_values(by="padj", ascending=True)
else:
logger.warning(
"No significant enrichment found with the given parameters. "
"Returning an empty DataFrame."
)
# ToDo: Should we return an empty DataFrame with the expected columns?
result = pd.DataFrame()
return result
def run_ssgsea(
data: pd.DataFrame,
annotation: str,
set_index: list[str] | None = None,
annotation_col: str = "annotation",
identifier_col: str = "identifier",
outdir: str = "tmp",
min_size: int = 15,
max_size: int = 500,
scale: bool = False,
permutations: int = 0,
) -> pd.DataFrame:
"""
Project each sample within a data set onto a space of gene set enrichment scores using
the single sample gene set enrichment analysis (ssGSEA) projection methodology
described in Barbie et al., 2009:
https://www.nature.com/articles/nature08460#Sec3 (search "Single Sample" GSEA).
:param pd.DataFrame data: pandas.DataFrame with the quantified features
(i.e. subject x proteins)
:param str annotation: pandas.DataFrame with the annotation to be used in the
enrichment (i.e. CKG pathway annotation file)
:param list[str] set_index: column/s to be used as index. Enrichment will be
calculated for these values (i.e ["subject"] will return subjects x pathways
matrix of enrichment scores)
:param str annotation_col: name of the column containing annotation terms.
:param str identifier_col: name of column containing dependent variables identifiers.
:param str out_dir: directory path where results will be stored
(default None, tmp folder is used)
:param int min_size: minimum number of features (i.e. proteins) in enriched terms
(i.e. pathways)
:param int max_size: maximum number of features (i.e. proteins) in enriched terms
(i.e. pathways)
:param bool scale: whether or not to scale the data
:param int permutations: number of permutations used in the ssgsea analysis
:return: pandas.DataFrame containing unnormalized enrichment scores (`ES`)
for each sample, and normalized enrichment scores (`NES`)
with the enriched `Term` and sample `Name`.
:rtype: pandas.DataFrame
Example::
stproject = "P0000008"
p = project.Project(
stproject,
datasets={},
knowledge=None,
report={},
configuration_files=None,
)
p.build_project(False)
p.generate_report()
proteomics_dataset = p.get_dataset("proteomics")
annotations = proteomics_dataset.get_dataframe("pathway annotation")
processed = proteomics_dataset.get_dataframe('processed')
result = run_ssgsea(
processed,
annotations,
annotation_col='annotation',
identifier_col='identifier',
set_index=['group',
'sample',
'subject'],
outdir=None,
min_size=10,
scale=False,
permutations=0
)
"""
df = data.copy()
if not os.path.exists(outdir):
os.makedirs(outdir)
# Comine columns to create a unique name for each set (?)
name = []
if set_index is None:
index = data.index.to_frame()
set_index = index.columns.tolist()
else:
index = data[set_index]
df = df.drop(set_index, axis=1)
for _, row in index.iterrows():
name.append(
"_".join(row[set_index].tolist())
) # this assumes strings as identifiers
df["Name"] = name
index.index = name
df = df.set_index("Name").transpose()
if annotation_col not in annotation:
raise ValueError(
f"Missing Annotation Column: {annotation_col}"
" as specified by `annotation_col`"
)
if identifier_col not in annotation:
raise ValueError(
f"Missing Identifier Column: {identifier_col}"
" as specified by `identifier_col`"
)
# ! Legacy implementation that was missing the description column of gmt files:
# ?gmt-columns: term<TAB>description<TAB>gene1<TAB>gene2
# grouped_annotations = (
# annotation.groupby(annotation_col)[identifier_col].apply(list).reset_index()
# )
# gene_sets = os.path.join(outdir, "ss_gsea" + ".gmt")
# with open(gene_sets, "w", encoding="utf8") as out:
# for _, row in grouped_annotations.iterrows():
# out.write(
# row[annotation_col]
# + "\t"
# + "\t".join(list(filter(None, row[identifier_col])))
# + "\n"
# )
# In case annotations were not provided by acore itself:
annotation = annotation.dropna(subset=[identifier_col])
annotation = annotation[annotation[identifier_col] != ""]
gene_sets = annotation.groupby(annotation_col)[identifier_col].apply(list).to_dict()
enrichment = gp.ssgsea(
data=df,
gene_sets=gene_sets,
outdir=outdir,
min_size=min_size,
max_size=max_size,
scale=scale,
permutation_num=permutations,
no_plot=True,
processes=1,
seed=10,
format="png",
)
result = pd.DataFrame(enrichment.res2d).set_index("Name")
# potentially return wide format in separate format
# result = {"es": enrichment_es, "nes": enrichment_nes}
return result
[docs]
def run_gsea(
regulation_data: pd.DataFrame,
annotation: pd.DataFrame,
identifier: str = "identifier",
ranking_col: str = "log2FC",
annotation_col: str = "annotation",
identifier_col: str = "identifier",
min_size: int = 15,
max_size: int = 500,
permutation_num: int = 1000,
seed: int = 42,
weight: float = 1.0,
ascending: bool = False,
correction_alpha: float = 0.05,
) -> pd.DataFrame:
"""
Run Gene Set Enrichment Analysis (GSEA) on a pre-ranked gene list using
`gseapy.prerank <https://gseapy.readthedocs.io/en/latest/run.html#gseapy.prerank>`_.
Genes are ranked by `ranking_col` (e.g. log2 fold-change from differential
expression analysis) and enrichment is computed for each gene set defined in
`annotation`. Unlike overrepresentation analysis (ORA), GSEA does not require
a significance cutoff to split genes into foreground/background; it uses the
full ranked list instead.
:param pd.DataFrame regulation_data: DataFrame with differential regulation
results. Must contain the column specified by `identifier` and the
column specified by `ranking_col`. Each identifier must appear at most
once (duplicates are not allowed).
:param pd.DataFrame annotation: DataFrame with functional annotations.
Must contain `annotation_col` and `identifier_col`.
:param str identifier: Column in `regulation_data` containing feature
identifiers. These must match the identifiers in `annotation_col` of
the `annotation` DataFrame.
:param str ranking_col: Column in `regulation_data` used to rank the
features (e.g. ``'log2FC'``). Features are ranked in descending order
by default so that strongly up-regulated features appear at the top.
:param str annotation_col: Column in `annotation` containing annotation
term names (e.g. pathway names).
:param str identifier_col: Column in `annotation` containing feature
identifiers that map to `identifier` in `regulation_data`.
:param int min_size: Minimum number of features in a gene set that are also
present in the ranking. Gene sets smaller than this threshold are
excluded.
:param int max_size: Maximum number of features in a gene set that are also
present in the ranking. Gene sets larger than this threshold are
excluded.
:param int permutation_num: Number of permutations used to estimate nominal
p-values. Higher values give more precise p-values but increase runtime.
Default: 1000.
:param int seed: Random seed for reproducibility. Default: 42.
:param float weight: Exponent for the enrichment score weighting.
``weight=1`` (default) corresponds to weighted GSEA; ``weight=0``
gives classic (unweighted) GSEA.
:param bool ascending: Sort order for the ranking. ``False`` (default)
means descending (highest ranking metric first, suited for log2FC).
:param float correction_alpha: FDR q-value cutoff used to define the
``'rejected'`` column. Default: 0.05.
:return: DataFrame with one row per gene set. Columns are: ``'Term'``,
``'ES'`` (enrichment score), ``'NES'`` (normalised enrichment score),
``'NOM p-val'`` (nominal p-value), ``'FDR q-val'`` (gseapy's
permutation-based FDR q-value), ``'FWER p-val'``, ``'Tag %'``,
``'Gene %'``, ``'Lead_genes'`` (semi-colon-separated leading-edge genes)
and ``'rejected'`` (``'FDR q-val' <= correction_alpha``).
:rtype: pd.DataFrame
Example::
result = run_gsea(
regulation_data=diff_reg,
annotation=annotations,
identifier='identifier',
ranking_col='log2FC',
annotation_col='annotation',
identifier_col='identifier',
min_size=5,
permutation_num=100,
)
"""
if regulation_data is None or regulation_data.empty:
raise ValueError("regulation_data is empty.")
if identifier not in regulation_data.columns:
raise ValueError(
f"Column '{identifier}' not found in regulation_data. "
"Specify the correct column name via `identifier`."
)
if ranking_col not in regulation_data.columns:
raise ValueError(
f"Column '{ranking_col}' not found in regulation_data. "
"Specify the correct column name via `ranking_col`."
)
if not regulation_data[identifier].is_unique:
raise ValueError(
f"Column '{identifier}' in regulation_data contains duplicate values. "
"Each identifier must appear at most once."
)
if annotation_col not in annotation.columns:
raise ValueError(
f"Missing annotation column '{annotation_col}' in annotation DataFrame. "
"Specify the correct column name via `annotation_col`."
)
if identifier_col not in annotation.columns:
raise ValueError(
f"Missing identifier column '{identifier_col}' in annotation DataFrame. "
"Specify the correct column name via `identifier_col`."
)
# Build the pre-ranked Series: index = feature identifier, value = ranking metric
rnk = regulation_data.set_index(identifier)[ranking_col].dropna()
# Build gene sets dict from the annotation DataFrame
gene_sets = annotation.groupby(annotation_col)[identifier_col].apply(list).to_dict()
enrichment = gp.prerank(
rnk=rnk,
gene_sets=gene_sets,
min_size=min_size,
max_size=max_size,
permutation_num=permutation_num,
weight=weight,
ascending=ascending,
seed=seed,
no_plot=True,
outdir=None,
)
result = pd.DataFrame(enrichment.res2d).drop(columns=["Name"])
# 'Tag %' and 'Gene %' are strings like "5/20" and "10.00%", so they stay as is.
numeric_cols = ["ES", "NES", "NOM p-val", "FDR q-val", "FWER p-val"]
non_numeric_cols = result[numeric_cols].select_dtypes(exclude="number").columns
if len(non_numeric_cols):
result[non_numeric_cols] = result[non_numeric_cols].apply(
pd.to_numeric, errors="coerce"
)
fdr = result["FDR q-val"]
result["rejected"] = (fdr <= correction_alpha).where(fdr.notna()).astype("boolean")
return result
