Add GSEA (pre-ranked) enrichment analysis to acore

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
 

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",

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
 #

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.

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)

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)