diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 93e7b34..05c8d6b 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -120,8 +120,7 @@ Before you submit a pull request, check that it meets these guidelines:
    into a function with a docstring, and add the feature to the list in README.rst.
 3. The pull request should pass the GitHub workflows.
 
-See the PR template example:
-[Add module PR template](https://github.com/Multiomics-Analytics-Group/acore/blob/main/.github/workflows/PULL_REQUEST_TEMPLATE/module.md)
+See the checklist if you open a PR!
 
 ## Deploying
 
diff --git a/docs/examples_single_analysis/enrichment_analysis.ipynb b/docs/examples_single_analysis/enrichment_analysis.ipynb
index cdbb032..ede67f7 100644
--- a/docs/examples_single_analysis/enrichment_analysis.ipynb
+++ b/docs/examples_single_analysis/enrichment_analysis.ipynb
@@ -47,6 +47,8 @@
     "\n",
     "import dsp_pandas\n",
     "import pandas as pd\n",
+    "import plotly.graph_objects as go\n",
+    "from vuecore import viz\n",
     "\n",
     "import acore\n",
     "import acore.differential_regulation\n",
@@ -544,7 +546,7 @@
     "enrichtments = acore.enrichment_analysis.run_ssgsea(\n",
     "    data=df_omics,\n",
     "    annotation=annotations,\n",
-    "    min_size=1,\n",
+    "    min_size=2,\n",
     ")\n",
     "enrichtments"
    ]
@@ -646,16 +648,97 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import plotly.graph_objects as go\n",
-    "from vuecore import viz\n",
-    "\n",
-    "args = {\"factor\": 2, \"loadings\": 1}  # increase number of loadings or scaling factor\n",
+    "args = {\"factor\": 2, \"loadings\": 2}  # increase number of loadings or scaling factor\n",
     "#! pca_results has three items, but docstring requests only two -> double check\n",
     "figure = viz.get_pca_plot(data=pca_result, identifier=\"PCA enrichment\", args=args)\n",
     "figure = go.Figure(data=figure[\"data\"], layout=figure[\"layout\"])\n",
     "figure"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "eb6bd6cd",
+   "metadata": {},
+   "source": [
+    "# GSEA (Gene Set Enrichment Analysis)\n",
+    "GSEA ranks all measured features (e.g. proteins) by a continuous metric such as the\n",
+    "log2 fold-change from differential expression analysis and tests whether the members of\n",
+    "a gene set are enriched at the top or bottom of this ranked list.\n",
+    "\n",
+    "Unlike overrepresentation analysis (ORA), GSEA does **not** require a significance\n",
+    "cutoff to define a foreground set: it uses the full ranked list and is therefore less\n",
+    "sensitive to the choice of an arbitrary threshold.\n",
+    "\n",
+    "Here we use `gseapy.prerank` under the hood, which implements the classic weighted\n",
+    "Kolmogorov–Smirnov enrichment statistic described in\n",
+    "[Subramanian et al. (2005)](https://doi.org/10.1073/pnas.0506580102).\n",
+    "\n",
+    "See [`gseapy.prerank`](https://gseapy.readthedocs.io/en/latest/run.html#gseapy.prerank)\n",
+    "for more details."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ef98e5b1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gsea_results = acore.enrichment_analysis.run_gsea(\n",
+    "    regulation_data=diff_reg,\n",
+    "    annotation=annotations,\n",
+    "    identifier=\"identifier\",\n",
+    "    ranking_col=\"log2FC\",\n",
+    "    annotation_col=\"annotation\",\n",
+    "    identifier_col=\"identifier\",\n",
+    "    min_size=1,\n",
+    "    permutation_num=100,\n",
+    "    seed=42,\n",
+    "    correction_alpha=0.05,\n",
+    ")\n",
+    "gsea_results"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f5f96579",
+   "metadata": {},
+   "source": [
+    "The normalised enrichment score (NES) summarises the enrichment: positive values\n",
+    "indicate that the gene set is enriched among the **up**-regulated features and negative\n",
+    "values indicate enrichment among the **down**-regulated features."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fc0abd44",
+   "metadata": {},
+   "source": [
+    "gseapy reports a permutation-based `FDR q-val` for each gene set (derived from the NES\n",
+    "null distribution), so no further correction is needed. The `rejected` column flags\n",
+    "gene sets with `FDR q-val <= correction_alpha`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fdfaa460",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gsea_results[[\"Term\", \"NES\", \"NOM p-val\", \"FDR q-val\", \"rejected\"]]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fd02b00c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ax = gsea_results[\"NES\"].plot.hist(title=\"Distribution of NES\", xlabel=\"NES\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "4be5a8c8",
diff --git a/docs/examples_single_analysis/enrichment_analysis.py b/docs/examples_single_analysis/enrichment_analysis.py
index b05df4e..8913ebe 100644
--- a/docs/examples_single_analysis/enrichment_analysis.py
+++ b/docs/examples_single_analysis/enrichment_analysis.py
@@ -22,6 +22,8 @@
 
 import dsp_pandas
 import pandas as pd
+import plotly.graph_objects as go
+from vuecore import viz
 
 import acore
 import acore.differential_regulation
@@ -305,7 +307,7 @@
 enrichtments = acore.enrichment_analysis.run_ssgsea(
     data=df_omics,
     annotation=annotations,
-    min_size=1,
+    min_size=2,
 )
 enrichtments
 
@@ -350,15 +352,60 @@
 # for this, which is also developed by the Multiomics Analytics Group.
 
 # %%
-import plotly.graph_objects as go
-from vuecore import viz
-
-args = {"factor": 2, "loadings": 1}  # increase number of loadings or scaling factor
+args = {"factor": 2, "loadings": 2}  # increase number of loadings or scaling factor
 # #! pca_results has three items, but docstring requests only two -> double check
 figure = viz.get_pca_plot(data=pca_result, identifier="PCA enrichment", args=args)
 figure = go.Figure(data=figure["data"], layout=figure["layout"])
 figure
 
+# %% [markdown]
+# # GSEA (Gene Set Enrichment Analysis)
+# GSEA ranks all measured features (e.g. proteins) by a continuous metric such as the
+# log2 fold-change from differential expression analysis and tests whether the members of
+# a gene set are enriched at the top or bottom of this ranked list.
+#
+# Unlike overrepresentation analysis (ORA), GSEA does **not** require a significance
+# cutoff to define a foreground set: it uses the full ranked list and is therefore less
+# sensitive to the choice of an arbitrary threshold.
+#
+# Here we use `gseapy.prerank` under the hood, which implements the classic weighted
+# Kolmogorov–Smirnov enrichment statistic described in
+# [Subramanian et al. (2005)](https://doi.org/10.1073/pnas.0506580102).
+#
+# See [`gseapy.prerank`](https://gseapy.readthedocs.io/en/latest/run.html#gseapy.prerank)
+# for more details.
+
+# %%
+gsea_results = acore.enrichment_analysis.run_gsea(
+    regulation_data=diff_reg,
+    annotation=annotations,
+    identifier="identifier",
+    ranking_col="log2FC",
+    annotation_col="annotation",
+    identifier_col="identifier",
+    min_size=1,
+    permutation_num=100,
+    seed=42,
+    correction_alpha=0.05,
+)
+gsea_results
+
+# %% [markdown]
+# The normalised enrichment score (NES) summarises the enrichment: positive values
+# indicate that the gene set is enriched among the **up**-regulated features and negative
+# values indicate enrichment among the **down**-regulated features.
+
+# %% [markdown]
+# gseapy reports a permutation-based `FDR q-val` for each gene set (derived from the NES
+# null distribution), so no further correction is needed. The `rejected` column flags
+# gene sets with `FDR q-val <= correction_alpha`.
+
+# %%
+gsea_results[["Term", "NES", "NOM p-val", "FDR q-val", "rejected"]]
+
+# %%
+ax = gsea_results["NES"].plot.hist(title="Distribution of NES", xlabel="NES")
+
 # %% [markdown]
 # # Compare two distributions - KS test
 #
diff --git a/docs/markdown_ref/acore.enrichment_analysis.md b/docs/markdown_ref/acore.enrichment_analysis.md
index 245cfc9..a3c2158 100644
--- a/docs/markdown_ref/acore.enrichment_analysis.md
+++ b/docs/markdown_ref/acore.enrichment_analysis.md
@@ -23,6 +23,76 @@ odds, pvalue = stats.fisher_exact(group1=[a, b],
                 )
 ```
 
+### run_gsea(regulation_data: DataFrame, annotation: DataFrame, identifier: [str](https://docs.python.org/3/builtins/stdtypes.html#str) = 'identifier', ranking_col: [str](https://docs.python.org/3/builtins/stdtypes.html#str) = 'log2FC', annotation_col: [str](https://docs.python.org/3/builtins/stdtypes.html#str) = 'annotation', identifier_col: [str](https://docs.python.org/3/builtins/stdtypes.html#str) = 'identifier', min_size: [int](https://docs.python.org/3/builtins/functions.html#int) = 15, max_size: [int](https://docs.python.org/3/builtins/functions.html#int) = 500, permutation_num: [int](https://docs.python.org/3/builtins/functions.html#int) = 1000, seed: [int](https://docs.python.org/3/builtins/functions.html#int) = 42, weight: [float](https://docs.python.org/3/builtins/functions.html#float) = 1.0, ascending: [bool](https://docs.python.org/3/builtins/functions.html#bool) = False, correction_alpha: [float](https://docs.python.org/3/builtins/functions.html#float) = 0.05) → 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.
+
+* **Parameters:**
+  * **regulation_data** (*pd.DataFrame*) – 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).
+  * **annotation** (*pd.DataFrame*) – DataFrame with functional annotations.
+    Must contain annotation_col and identifier_col.
+  * **identifier** ([*str*](https://docs.python.org/3/builtins/stdtypes.html#str)) – Column in regulation_data containing feature
+    identifiers. These must match the identifiers in annotation_col of
+    the annotation DataFrame.
+  * **ranking_col** ([*str*](https://docs.python.org/3/builtins/stdtypes.html#str)) – 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.
+  * **annotation_col** ([*str*](https://docs.python.org/3/builtins/stdtypes.html#str)) – Column in annotation containing annotation
+    term names (e.g. pathway names).
+  * **identifier_col** ([*str*](https://docs.python.org/3/builtins/stdtypes.html#str)) – Column in annotation containing feature
+    identifiers that map to identifier in regulation_data.
+  * **min_size** ([*int*](https://docs.python.org/3/builtins/functions.html#int)) – Minimum number of features in a gene set that are also
+    present in the ranking. Gene sets smaller than this threshold are
+    excluded.
+  * **max_size** ([*int*](https://docs.python.org/3/builtins/functions.html#int)) – Maximum number of features in a gene set that are also
+    present in the ranking. Gene sets larger than this threshold are
+    excluded.
+  * **permutation_num** ([*int*](https://docs.python.org/3/builtins/functions.html#int)) – Number of permutations used to estimate nominal
+    p-values. Higher values give more precise p-values but increase runtime.
+    Default: 1000.
+  * **seed** ([*int*](https://docs.python.org/3/builtins/functions.html#int)) – Random seed for reproducibility. Default: 42.
+  * **weight** ([*float*](https://docs.python.org/3/builtins/functions.html#float)) – Exponent for the enrichment score weighting.
+    `weight=1` (default) corresponds to weighted GSEA; `weight=0`
+    gives classic (unweighted) GSEA.
+  * **ascending** ([*bool*](https://docs.python.org/3/builtins/functions.html#bool)) – Sort order for the ranking. `False` (default)
+    means descending (highest ranking metric first, suited for log2FC).
+  * **correction_alpha** ([*float*](https://docs.python.org/3/builtins/functions.html#float)) – FDR q-value cutoff used to define the
+    `'rejected'` column. Default: 0.05.
+* **Returns:**
+  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`).
+* **Return type:**
+  pd.DataFrame
+
+Example:
+
+```default
+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,
+)
+```
+
 ### run_kolmogorov_smirnov(dist1: [list](https://docs.python.org/3/builtins/stdtypes.html#list)[[float](https://docs.python.org/3/builtins/functions.html#float)], dist2: [list](https://docs.python.org/3/builtins/stdtypes.html#list)[[float](https://docs.python.org/3/builtins/functions.html#float)], alternative: [str](https://docs.python.org/3/builtins/stdtypes.html#str) = 'two-sided') → [tuple](https://docs.python.org/3/builtins/stdtypes.html#tuple)[[float](https://docs.python.org/3/builtins/functions.html#float), [float](https://docs.python.org/3/builtins/functions.html#float)]
 
 Compute the Kolmogorov-Smirnov statistic on 2 samples.
diff --git a/docs/markdown_ref/acore.md b/docs/markdown_ref/acore.md
index 6627268..1e0ab9b 100644
--- a/docs/markdown_ref/acore.md
+++ b/docs/markdown_ref/acore.md
@@ -68,6 +68,7 @@
     * [`run_loess_drift_correction()`](acore.drift_correction.md#acore.drift_correction.loess_drift_correction.run_loess_drift_correction)
 * [acore.enrichment_analysis package](acore.enrichment_analysis.md)
   * [`run_fisher()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_fisher)
+  * [`run_gsea()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_gsea)
   * [`run_kolmogorov_smirnov()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_kolmogorov_smirnov)
   * [`run_site_regulation_enrichment()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_site_regulation_enrichment)
   * [`run_up_down_regulation_enrichment()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_up_down_regulation_enrichment)
diff --git a/docs/markdown_ref/modules.md b/docs/markdown_ref/modules.md
index e128c15..d8ef499 100644
--- a/docs/markdown_ref/modules.md
+++ b/docs/markdown_ref/modules.md
@@ -41,6 +41,7 @@
       * [acore.drift_correction.loess_drift_correction module](acore.drift_correction.md#module-acore.drift_correction.loess_drift_correction)
     * [acore.enrichment_analysis package](acore.enrichment_analysis.md)
       * [`run_fisher()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_fisher)
+      * [`run_gsea()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_gsea)
       * [`run_kolmogorov_smirnov()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_kolmogorov_smirnov)
       * [`run_site_regulation_enrichment()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_site_regulation_enrichment)
       * [`run_up_down_regulation_enrichment()`](acore.enrichment_analysis.md#acore.enrichment_analysis.run_up_down_regulation_enrichment)
diff --git a/src/acore/enrichment_analysis/__init__.py b/src/acore/enrichment_analysis/__init__.py
index 3538e0e..6314c8a 100644
--- a/src/acore/enrichment_analysis/__init__.py
+++ b/src/acore/enrichment_analysis/__init__.py
@@ -10,7 +10,6 @@
 import logging
 import os
 import re
-import uuid
 
 import gseapy as gp
 import pandas as pd
@@ -33,6 +32,7 @@
 
 __all__ = [
     "run_fisher",
+    "run_gsea",
     "run_kolmogorov_smirnov",
     "run_site_regulation_enrichment",
     "run_up_down_regulation_enrichment",
@@ -228,6 +228,12 @@ def run_up_down_regulation_enrichment(
                 correction=correction,
                 correction_alpha=correction_alpha,
             )
+            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)
@@ -563,23 +569,28 @@ def run_ssgsea(
             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()
 
-    grouped_annotations = (
-        annotation.groupby(annotation_col)[identifier_col].apply(list).reset_index()
-    )
-    fid = uuid.uuid4()
-    file_path = os.path.join(outdir, str(fid) + ".gmt")
-    with open(file_path, "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"
-            )
     enrichment = gp.ssgsea(
         data=df,
-        gene_sets=str(file_path),
+        gene_sets=gene_sets,
         outdir=outdir,
         min_size=min_size,
         max_size=max_size,
@@ -594,3 +605,142 @@ def run_ssgsea(
     # potentially return wide format in separate format
     # result = {"es": enrichment_es, "nes": enrichment_nes}
     return result
+
+
+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
diff --git a/tests/enrichment_analysis/test_enrichment.py b/tests/enrichment_analysis/test_enrichment.py
index a72df6c..9018e1d 100644
--- a/tests/enrichment_analysis/test_enrichment.py
+++ b/tests/enrichment_analysis/test_enrichment.py
@@ -600,3 +600,80 @@ def test_run_regulation_enrichment_with_duplicates():
             annotation=annotation,
             min_detected_in_set=1,
         )
+
+
+def test_run_gsea_basic():
+    """Integration test for run_gsea: checks that the function returns
+    a DataFrame with the expected columns and one row per gene set."""
+    regulation_data = pd.DataFrame(
+        {
+            "identifier": ["gene1", "gene2", "gene3", "gene4", "gene5", "gene6"],
+            "log2FC": [3.0, 2.5, 1.0, -1.0, -2.5, -3.0],
+        }
+    )
+    annotation = pd.DataFrame(
+        {
+            "annotation": ["path1", "path1", "path1", "path2", "path2", "path2"],
+            "identifier": ["gene1", "gene2", "gene3", "gene4", "gene5", "gene6"],
+        }
+    )
+
+    result = ea.run_gsea(
+        regulation_data=regulation_data,
+        annotation=annotation,
+        identifier="identifier",
+        ranking_col="log2FC",
+        annotation_col="annotation",
+        identifier_col="identifier",
+        min_size=2,
+        permutation_num=10,
+        seed=42,
+    )
+
+    assert isinstance(result, pd.DataFrame)
+    assert not result.empty
+    expected_cols = {"Term", "ES", "NES", "NOM p-val", "FDR q-val", "Lead_genes"}
+    assert expected_cols.issubset(result.columns)
+    assert "rejected" in result.columns
+    fdr = result["FDR q-val"].astype(float)
+    assert (result["rejected"] == (fdr <= 0.05)).all()
+    assert set(result["Term"]) == {"path1", "path2"}
+    # path1 genes are up-regulated -> positive NES
+    nes_path1 = result.loc[result["Term"] == "path1", "NES"].iloc[0]
+    nes_path2 = result.loc[result["Term"] == "path2", "NES"].iloc[0]
+    assert nes_path1 > 0, "path1 (up-regulated genes) should have a positive NES"
+    assert nes_path2 < 0, "path2 (down-regulated genes) should have a negative NES"
+
+
+def test_run_gsea_raises_on_empty_data():
+    """run_gsea should raise ValueError when regulation_data is empty."""
+    with pytest.raises(ValueError, match="regulation_data is empty"):
+        ea.run_gsea(
+            regulation_data=pd.DataFrame(),
+            annotation=pd.DataFrame({"annotation": [], "identifier": []}),
+        )
+
+
+def test_run_gsea_raises_on_missing_identifier_col():
+    """run_gsea should raise ValueError when identifier column is missing."""
+    regulation_data = pd.DataFrame({"gene": ["g1"], "log2FC": [1.0]})
+    annotation = pd.DataFrame({"annotation": ["path1"], "identifier": ["g1"]})
+    with pytest.raises(ValueError, match="identifier"):
+        ea.run_gsea(
+            regulation_data=regulation_data,
+            annotation=annotation,
+            identifier="identifier",  # not a column in regulation_data
+        )
+
+
+def test_run_gsea_raises_on_duplicate_identifiers():
+    """run_gsea should raise ValueError when identifier column has duplicates."""
+    regulation_data = pd.DataFrame(
+        {"identifier": ["gene1", "gene1"], "log2FC": [1.0, 2.0]}
+    )
+    annotation = pd.DataFrame({"annotation": ["path1"], "identifier": ["gene1"]})
+    with pytest.raises(ValueError, match="duplicate"):
+        ea.run_gsea(
+            regulation_data=regulation_data,
+            annotation=annotation,
+        )