---
title: "Tutorial 07 — Functional Analysis (GO + GSEA) with clusterProfiler"
subtitle: "Take a per-cell-type pseudobulk DE list and turn it into biology"
author: "Single Cell RNA-seq Workshop"
format:
  html:
    toc: true
    toc-depth: 3
    code-fold: false
    code-overflow: wrap
    highlight-style: github
    embed-resources: true
execute:
  eval: false
  echo: true
  warning: false
  message: false
editor: visual
---

::: {.callout-note title="Running the code in this tutorial"}
Every code chunk here is tagged with `#| eval: false`, so the published page shows the code **without running it**. To run it yourself, open the downloaded `.qmd` in **RStudio**:

- **Run one chunk:** click the green ▶ (*Run Current Chunk*) at the chunk's top-right, or press *Ctrl/Cmd + Enter*. This runs the chunk **regardless of its `eval` setting** — the easiest way to work through the tutorial interactively.
- **Run a chunk on render:** change that chunk's `#| eval: false` to `#| eval: true` (or delete the line) so it executes when the document is rendered.
- **Render the whole tutorial:** click the **Render** button in RStudio (or run `quarto render` in a terminal) to execute the chunks top-to-bottom and knit a finished HTML. To run everything on render, set `eval: true` once in the YAML header at the top.
:::


## About this tutorial

A short follow-on to **[Tutorial 06 — DESeq2: Bulk Primer + Pseudobulk DE](Tutorial_06_DESeq2_DE.html)**. We take the per-cell-type DE table you wrote at the end of Tutorial 06 (`ifnb_pseudobulk_de.csv`) and run two of the most common downstream functional-enrichment analyses with `clusterProfiler`:


::: callout-note
**Companion lecture:** [Lecture 07 — Functional Analysis: Turning DE Lists into Biology](../Lecture_Folder/Lecture_07_FunctionalAnalysis.html) · **Companion book chapter:** [Chapter 7 — Functional Analysis](../Resources_Folder/Chapter_07_FunctionalAnalysis.html) — the long-form prose treatment of this tutorial's material, with cross-references to the prerequisite appendices. · **HPC version (Talapas):** [Talapas analysis pipeline — `07_functional_analysis.R`](Tutorial_10_Talapas_Pipeline.html)
:::


-   **GO over-representation analysis (ORA)** — given a *list* of significant genes, which Gene Ontology terms are enriched compared to the genome background?
-   **GSEA — gene-set enrichment analysis** — given a *ranked* list of all genes (signed by direction of change), which terms are enriched at the top or bottom?

We'll do both per cell type and then **compare cell types** with `compareCluster()` — a single dot plot that shows which biological themes are shared between cell types vs cell-type–specific.

::: callout-note
**Why two methods?** ORA is fast and easy to read but throws away the gene ranking — every "significant" gene contributes equally. GSEA uses the *whole* ranked list and is much more sensitive when the effect is broad-but-modest (a coordinated shift across many ISGs, for example). For STIM-vs-CTRL on `ifnb`, both will work; in practice you should run both and look for agreement.
:::

::: {.callout-tip title="How to use this page"}
1.  Download the `.qmd` source: <a href="Tutorial_07_FunctionalAnalysis.qmd.txt" download="Tutorial_07_FunctionalAnalysis.qmd">Tutorial_07_FunctionalAnalysis.qmd</a>. If your browser saves the file as `Tutorial_07_FunctionalAnalysis.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  Make sure you have completed **[Tutorial 06](Tutorial_06_DESeq2_DE.html)** — it writes `ifnb_pseudobulk_de.csv`.
3.  Work through the chunks, flipping `eval: true` as you go.
:::

## Dataset — resuming from Tutorial 06

| File | Produced by | What it is |
|---|---|---|
| `../data/ifnb_pseudobulk_de.csv` | Tutorial 06 | Long table of per-cell-type DE results: columns `gene`, `baseMean`, `log2FoldChange`, `lfcSE`, `pvalue`, `padj`, `celltype`. The contrast was `STIM` vs `CTRL`. |


::: {.callout-warning title="Common errors / things that bite"}
**Many gene symbols don't map via `bitr`** — gene-symbol-to-ENTREZ-ID mapping is incomplete (deprecated symbols, lncRNAs, pseudogenes). A 5–10% unmapped rate is normal. Check the unmapped list with `setdiff(unique(de$gene), de_mapped$symbol)` — if it's mostly pseudogenes / lncRNAs, ignore. If it's all gene symbols, your `org.Hs.eg.db` annotation may simply be out of date (update it through Bioconductor when you next refresh your packages).

**`enrichGO` returns 0 enriched terms** — most often you forgot the `universe = uni` argument and `enrichGO` is using "all annotated genes in the genome" as the background. That over-inflates the p-value of common terms and makes everything non-significant. Always pass the genes you actually tested.

**GSEA produces "ties detected" warning** — `log2FoldChange` ties from multiple genes with identical effect sizes break the rank-based test. `arrange(desc(log2FoldChange))` followed by `distinct(ENTREZID, .keep_all = TRUE)` handles it. The tutorial does this in the helper.

**`object '.' not found` inside a `bitr(...)` / pipe** — the **native pipe** `|>` has no `.` placeholder (that's magrittr `%>%`). Refer to the data frame by name (e.g. `unique(de$gene)`) instead of `.$symbol`.

**`object 'de_mapped' not found` (or similar) further down** — that's a *consequence* of an earlier chunk failing (e.g. the `bitr` step), so the object was never created. Fix the first error and re-run from there; don't debug the downstream chunk.

**Bare `select`/`filter`/`desc` misbehave** — `org.Hs.eg.db`/`clusterProfiler` (and their `AnnotationDbi`/`IRanges` deps) mask those dplyr verbs. The setup chunk re-binds them to dplyr (`select <- dplyr::select`, …) so the pipelines work; keep that block if you adapt the tutorial.

**The GSEA highlight figure (Mod7_C6) is missing** — it's only drawn for one cell type. The tutorial picks it data-driven (a monocyte type if present, else the first with GSEA results); hardcoding a name that isn't in `names(gse_list)` silently skips the plot.
:::


::: {.callout-tip title="Solutions / instructor copy"}
The *Think about it* prompts on this page have inline answers (click "Show answer" / "Click for answer" to expand). For the full **runnable instructor copy** with every chunk pre-evaluated and outputs shown:

- **For students:** the same `.qmd` you downloaded already contains every solution inline — just expand the collapsed answers
- **For instructors:** render with `quarto render --profile solutions` from the `Exercise_Folder/` directory. The fully-evaluated HTML is written to `docs/Exercise_Folder/_solutions/Tutorial_07_FunctionalAnalysis.html`. See [Exercise_Folder/_quarto-solutions.yml](https://github.com/wcresko/scRNAseq_tutorial/blob/main/Exercise_Folder/_quarto-solutions.yml) for the build profile
:::

## Setup

```{r}
#| label: M7-setup
library(tidyverse)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ReactomePA)
library(patchwork)   # combine + annotate multi-panel plots

set.seed(2026)

# ---------------------------------------------------------------------------
# Output directory for this module's figures and tables.
# Every figure/table chunk below writes a file named  Mod7_C<chunk>_<name>
# into ../output/Mod7/ so it can be cross-referenced from the rest of the site.
#   Mod7 = Module 7 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod7"
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
dir.create("../data",  showWarnings = FALSE)   # shared .rds hand-off folder (sibling of scripts/, created one level up)
message("This module writes its figures & tables to: ", normalizePath(out_dir))

# Prefer dplyr's verbs over identically-named functions from the Bioconductor
# packages loaded above: S4Vectors/IRanges/AnnotationDbi/matrixStats (pulled in by
# org.Hs.eg.db, DESeq2, miloR, SingleCellExperiment, ...) mask select, filter,
# count, desc, slice, rename, first. Re-bind the ones used below so the dplyr
# pipelines work regardless of package load order.
select <- dplyr::select; filter <- dplyr::filter; count <- dplyr::count
desc   <- dplyr::desc;   slice  <- dplyr::slice
rename <- dplyr::rename; first  <- dplyr::first

# --- Figure saving --------------------------------------------------------
# save_fig() writes every figure as a .png (for viewing) AND a .svg (vector,
# editable in Illustrator / Inkscape for a manuscript). Pass the .png path; the
# .svg is written alongside with the same basename. (.svg uses the 'svglite'
# package, installed in Tutorial 00.)
save_fig <- function(filename, plot, width, height, dpi = 300, ...) {
  ggplot2::ggsave(filename, plot, width = width, height = height, dpi = dpi, ...)
  svg_path <- paste0(tools::file_path_sans_ext(filename), ".svg")
  tryCatch(ggplot2::ggsave(svg_path, plot, width = width, height = height, ...),
           error = function(e) message("  (could not write ", basename(svg_path),
                                       " - install 'svglite'? ", conditionMessage(e), ")"))
}

de <- read_csv("../data/ifnb_pseudobulk_de.csv")
glimpse(de)
table(de$celltype)
```

## Step 1 — Map gene symbols to ENTREZ IDs

Most enrichment tools work in ENTREZ-ID space because that's how `OrgDb` packages key their annotations. Convert once with `bitr()`:

```{r}
#| label: M7-bitr

de_mapped <- de |>
  mutate(symbol = gene) |>
  inner_join(
    # NB: the native pipe |> has no `.` placeholder (that's magrittr %>%), so refer
    # to `de` directly here rather than `.$symbol`.
    bitr(unique(de$gene),
         fromType = "SYMBOL",
         toType   = "ENTREZID",
         OrgDb    = org.Hs.eg.db),
    by = c("symbol" = "SYMBOL")
  )

# How many genes failed to map?
cat("Mapped:   ", n_distinct(de_mapped$symbol),   "\n")
cat("Unmapped: ", n_distinct(de$gene) - n_distinct(de_mapped$symbol), "\n")
```

::: {.callout-important title="Think about it"}
Why do some symbols fail to map? What kinds of genes typically don't have ENTREZ IDs?

<details>

<summary>Show answers</summary>

`bitr` won't map: deprecated symbols (HGNC renamed the gene since `org.Hs.eg.db` was last built), pseudogenes, lncRNAs that lack an ENTREZ entry, and any symbol with a typo. The unmapped fraction is usually small (<5%) and rarely matters for ORA/GSEA — but inspect the unmapped list if it's larger than that, since you may have an `org.Hs.eg.db` version mismatch.

</details>
:::

## Step 2 — GO over-representation per cell type

For each cell type, take the genes that are **significantly DE** (we'll require `padj < 0.05` and `|log2FC| > 0.5`) and ask: against the *universe* of tested genes for that cell type, which GO Biological Process terms are over-represented?

[]{#lst-M7-enrichGO}

```{r}
#| label: M7-enrichGO

run_enrichGO <- function(ct, df) {
  sub <- df |> filter(celltype == ct)
  sig <- sub |> filter(padj < 0.05, abs(log2FoldChange) > 0.5) |> pull(ENTREZID)
  uni <- sub |> pull(ENTREZID)
  if (length(sig) < 10) return(NULL)
  ego <- enrichGO(
    gene          = sig,
    universe      = uni,
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",                  # Biological Process
    keyType       = "ENTREZID",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.10,
    readable      = TRUE
  )
  ego@result$celltype <- ct
  ego
}

ego_list <- map(unique(de_mapped$celltype),
                run_enrichGO,
                df = de_mapped) |>
  set_names(unique(de_mapped$celltype)) |>
  compact()

# Top 5 BP terms per cell type
ego_top5 <- map_dfr(names(ego_list),
        \(ct) head(ego_list[[ct]]@result, 5) |>
              mutate(celltype = ct) |>
              select(celltype, ID, Description, Count, p.adjust))
ego_top5

# --- Table out: top-5 GO BP terms per cell type ----------------------------
readr::write_csv(ego_top5, file.path(out_dir, "Mod7_C3_enrichgo_top5.csv"))
```

::: {.callout-tip title="Reading the output"}
The `ego_top5` table ([code](#lst-M7-enrichGO)) shows the top 5 GO Biological Process terms per cell type, columns `Description` (the term name), `Count` (how many of your significant genes hit this term), and `p.adjust` (BH-corrected p-value). For `ifnb` STIM vs CTRL, expect `type I interferon signaling pathway`, `response to virus`, and related immune-defense terms to dominate — that is the biology of the contrast. The interesting signal is in the *cell-type-specific* rows that appear for one cell type only, such as T-cell-receptor signaling in T cells or complement activation in monocytes. A `Count` of 2–3 with a very small `p.adjust` often indicates a small but tightly defined term — useful, but verify the gene list. `enrichplot::dotplot(ego_list[[ct]])` renders this as a ranked dot plot for any single cell type.
:::

::: {.callout-important title="Think about it"}
1.  Why pass `universe = uni` rather than letting `enrichGO` use the whole genome?
2.  Some GO terms (`response to virus`, `defense response`, `type I interferon signaling pathway`) will dominate every cell type. Why?

<details>

<summary>Show answers</summary>

1.  The "universe" should be the set of genes you actually **tested**. If a cell type has very few expressed genes (because, say, granulocytes have shallow libraries), using the whole genome as the background inflates the apparent enrichment of common terms. Always use the tested-gene universe.
2.  Because IFN-β stimulation **literally** activates the type I interferon response — those terms are the right answer for the contrast you asked about. The interesting thing is *which* cell types show the strongest enrichment (typically monocytes) and which show the weakest (sometimes B cells), and what *additional* terms come up cell-type-specifically (e.g. `T cell receptor signaling` only in T cells).

</details>
:::

## Step 3 — Cross-cell-type comparison with `compareCluster`

A single dot plot showing which BP terms are enriched in which cell type:

```{r}
#| label: M7-compareCluster

# Build the input list: named list of significant ENTREZIDs per cell type
sig_by_ct <- de_mapped |>
  filter(padj < 0.05, abs(log2FoldChange) > 0.5) |>
  split(~ celltype) |>
  map(~ unique(.x$ENTREZID))

# Drop empty
sig_by_ct <- compact(sig_by_ct[lengths(sig_by_ct) >= 10])

cc <- compareCluster(
  geneClusters  = sig_by_ct,
  fun           = "enrichGO",
  OrgDb         = org.Hs.eg.db,
  ont           = "BP",
  keyType       = "ENTREZID",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05
)

p_cc <- dotplot(cc, showCategory = 8) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    title    = "GO BP enrichment per cell type — STIM vs CTRL",
    subtitle = "Top 8 over-represented Biological Process terms per cell type",
    x        = "Cell type",
    y        = "GO Biological Process term",
    colour   = "Adjusted p-value",
    size     = "Gene ratio"
  )
p_cc

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod7_C4_comparecluster_dotplot.png"), p_cc,
       width = 10, height = 8, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The x-axis is cell type and the y-axis is GO Biological Process term. Each dot marks one enriched term in one cell type: **dot size** encodes the gene ratio (fraction of significant genes in that term), and **dot colour** encodes the BH-adjusted p-value (darker = more significant). Terms shared across many cell types appear as a horizontal row of dots — in this dataset, `type I interferon signaling pathway` will run across almost every column because IFN-β stimulation is the contrast. Cell-type-specific terms (a dot in only one column) are the biologically interesting findings. If every column shows exactly the same terms with equal significance, raise the `|log2FC|` threshold in `sig_by_ct` to filter to only the strongest responders — that often surfaces the cell-type-specific signal beneath the shared core response.
:::

::: {.callout-important title="Think about it"}
What do you do when *all* cell types show the same top terms?

<details>

<summary>Show answers</summary>

Step out of GO BP and into a more focused ontology — try `ont = "MF"` (Molecular Function) or `ont = "CC"` (Cellular Component), or use a tighter database like **Reactome pathways** (Step 5 below). You can also raise the LFC threshold so only the *strongest* responders contribute — that often surfaces cell-type-specific themes that get washed out by the shared core response.

</details>
:::

## Step 4 — GSEA per cell type

Where ORA throws away the ranking, GSEA uses it. Build a **named, decreasing-sorted vector** of `log2FoldChange`s for each cell type and run `gseGO`:

[]{#lst-M7-gseGO}

```{r}
#| label: M7-gseGO

run_gseGO <- function(ct, df) {
  sub <- df |>
    filter(celltype == ct, !is.na(ENTREZID)) |>
    arrange(desc(log2FoldChange)) |>
    distinct(ENTREZID, .keep_all = TRUE)

  rk <- setNames(sub$log2FoldChange, sub$ENTREZID)

  gse <- gseGO(
    geneList      = rk,
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",
    keyType       = "ENTREZID",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    eps           = 0,
    verbose       = FALSE
  )
  gse@result$celltype <- ct
  gse
}

# This may take 1-3 minutes per cell type
gse_list <- map(unique(de_mapped$celltype), run_gseGO, df = de_mapped) |>
  set_names(unique(de_mapped$celltype)) |>
  compact()

# Top GSEA hits per cell type, ranked by NES (Normalized Enrichment Score)
gse_top5 <- map_dfr(names(gse_list),
        \(ct) gse_list[[ct]]@result |>
              arrange(desc(NES)) |>
              head(5) |>
              mutate(celltype = ct) |>
              select(celltype, Description, NES, p.adjust))
gse_top5

# --- Table out: top-5 GSEA BP terms per cell type (by NES) -----------------
readr::write_csv(gse_top5, file.path(out_dir, "Mod7_C5_gsego_top5.csv"))
```

::: {.callout-tip title="Reading the output"}
The `gse_top5` table ([code](#lst-M7-gseGO)) shows the top 5 GSEA GO BP hits per cell type ordered by **NES** (Normalized Enrichment Score). A positive NES means the gene set is enriched at the top of the ranked list (genes most up in STIM), while a negative NES means enrichment at the bottom (genes most up in CTRL). Unlike ORA, GSEA uses the whole ranked list, so modest but coordinated shifts across many ISGs produce a strongly positive NES even if few individual genes pass an FDR threshold. For `ifnb`, expect interferon-signaling terms with NES > 2 in monocytes; a term with NES close to 1 or –1 and a borderline `p.adjust` is likely noise. Use `enrichplot::gseaplot2()` to draw the running enrichment curve for any term and visually confirm the hit.
:::

A classic GSEA "running enrichment" plot for a single term in a single cell type:

[]{#lst-M7-gsea_plot}

```{r}
#| label: M7-gsea_plot
#| fig-cap: "GSEA running-enrichment plot for the top 3 GO BP terms in CD14 monocytes."

# Pick a cell type to highlight — data-driven so it matches your labels: prefer a
# monocyte type (strongest IFN response), else the first with GSEA results. (The
# muscData ifnb labels are "CD14+ Monocytes" etc., not "CD14 Mono".)
ct_pick <- grep("Mono", names(gse_list), value = TRUE)[1]
if (is.na(ct_pick)) ct_pick <- if (length(gse_list)) names(gse_list)[1] else NA_character_
if (!is.na(ct_pick) && nrow(gse_list[[ct_pick]]@result) > 0) {
  p_gsea <- gseaplot2(gse_list[[ct_pick]],
            geneSetID = head(gse_list[[ct_pick]]@result$ID, 3),
            title     = paste0(ct_pick, " — top 3 BP enrichments"))
  print(p_gsea)

  # --- Figure out ----------------------------------------------------------
  save_fig(file.path(out_dir, "Mod7_C6_gsea_running_enrichment.png"), p_gsea,
         width = 8, height = 6, dpi = 300)
}
```

::: {.callout-tip title="Reading the output"}
The GSEA running-enrichment plot ([code](#lst-M7-gsea_plot)) has three panels per gene set. The top panel shows the **running enrichment score** (y-axis) as you walk from the most up-regulated genes (left) to the most down-regulated (right) in your ranked list — the score climbs steeply for a positively enriched set, reaches a peak, then falls. The middle tick-mark panel marks where each gene from the set falls in the ranked list. The bottom bar shows the ranked list metric (`log2FoldChange`) as a gradient from positive (red) to negative (blue). A healthy enrichment shows an early, steep rise in the top panel; a flat, noisy curve with an ambiguous peak suggests the gene set is not coherently regulated. The **NES** and `p.adjust` values from the `gse_top5` table tell you whether the peak is statistically reliable.
:::

## Step 5 — Reactome pathways (cleaner than GO for signaling biology)

Reactome's curated pathways are typically more interpretable than GO BP for immune signaling. Same input, different database:

[]{#lst-M7-reactome}

```{r}
#| label: M7-reactome

reactome_list <- map(unique(de_mapped$celltype), function(ct) {
  sub <- de_mapped |> filter(celltype == ct)
  sig <- sub |> filter(padj < 0.05, abs(log2FoldChange) > 0.5) |> pull(ENTREZID)
  if (length(sig) < 10) return(NULL)
  enrichPathway(gene = sig,
                organism = "human",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                readable      = TRUE)
}) |> set_names(unique(de_mapped$celltype)) |> compact()

# Top hits per cell type
reactome_top5 <- map_dfr(names(reactome_list),
        \(ct) head(reactome_list[[ct]]@result, 5) |>
              mutate(celltype = ct) |>
              select(celltype, Description, p.adjust))
reactome_top5

# --- Table out: top-5 Reactome pathways per cell type ----------------------
readr::write_csv(reactome_top5, file.path(out_dir, "Mod7_C7_reactome_top5.csv"))
```

::: {.callout-tip title="Reading the output"}
The `reactome_top5` table ([code](#lst-M7-reactome)) lists the top 5 Reactome pathways per cell type, with `Description` (the pathway name) and `p.adjust` (BH-corrected). Reactome pathways are more specific and less redundant than GO BP: you'll see named signaling cascades (`Interferon alpha/beta signaling`, `Cytokine Signaling in Immune system`, `MHC class II antigen presentation`) rather than the vague meta-terms GO often produces. Compare these results to your GO BP hits — concordance between the two databases is strong evidence for a real biological signal. Reactome divergence (a term appearing in Reactome but not GO, or vice versa) often reflects how each database drew its boundaries around overlapping biology, not a discrepancy in your data.
:::

::: {.callout-important title="Think about it"}
1.  GO BP, GO MF, Reactome, KEGG, MSigDB Hallmark, WikiPathways — there are a lot of databases. How do you pick?
2.  Hallmark vs Reactome vs GO — what's the trade-off?

<details>

<summary>Show answers</summary>

1.  **Pick by question, then sanity-check with one alternative.** If you're asking about *signaling*, start with Reactome or KEGG. If you're asking about *cell process* (apoptosis, proliferation, autophagy), start with GO BP. If you want a *small, well-curated* hallmark of broad biological states, start with MSigDB Hallmark (50 well-known gene sets). Then run a second database to confirm the result isn't a database artifact.
2.  **Hallmark** (50 sets, 200 genes each, MSigDB v7+) — coarse but highly curated; fewer false positives, less granularity. **Reactome** — pathway-centric, well-organized hierarchy, excellent for signaling. **GO BP** — very fine-grained but high redundancy; you'll see the same response named ten different ways. Use Hallmark for headline figures, Reactome / GO for the mechanistic story.

</details>
:::

## Step 6 — Save the enrichment tables

```{r}
#| label: M7-save

# GO BP per cell type
ego_long <- map_dfr(names(ego_list),
                    \(ct) ego_list[[ct]]@result |>
                          mutate(celltype = ct))
readr::write_csv(ego_long, file.path(out_dir, "Mod7_C8_enrichgo_bp.csv"))

# Reactome per cell type
reactome_long <- map_dfr(names(reactome_list),
                         \(ct) reactome_list[[ct]]@result |>
                               mutate(celltype = ct))
readr::write_csv(reactome_long, file.path(out_dir, "Mod7_C8_reactome.csv"))

# GSEA per cell type
gse_long <- map_dfr(names(gse_list),
                    \(ct) gse_list[[ct]]@result |>
                          mutate(celltype = ct))
readr::write_csv(gse_long, file.path(out_dir, "Mod7_C8_gsego_bp.csv"))
```

## Wrap-up

You now have:

-   GO BP, GO MF (rerun the helper with `ont = "MF"`), and Reactome ORA tables per cell type
-   A GSEA table per cell type with NES and `p.adjust`
-   A `compareCluster` dot plot showing cell-type–specific vs shared biology
-   A workflow you can re-point at *any* DE result — drop in your own table from DESeq2 / edgeR / limma and re-run

::: {.callout-tip title="Going further"}
-   **MSigDB hallmark** sets via `msigdbr::msigdbr(category = "H")` plug straight into `enricher()` and `GSEA()`.
-   **Network plots:** `cnetplot()` and `emapplot()` in `enrichplot` show how enriched terms relate to each other and to the contributing genes — much more interpretable than long bar plots when there are >20 hits.
-   **Cross-condition trajectories:** if you have multiple conditions, run pseudobulk DE per pairwise contrast, then `compareCluster` with `fun = "GSEA"` to see which pathways shift along the gradient.
-   The same `clusterProfiler` framework supports KEGG (`enrichKEGG`/`gseKEGG`), DO disease ontology, and custom MSigDB collections.
:::

## See also

-   [Tutorial 06 — DESeq2: Bulk Primer + Pseudobulk DE](Tutorial_06_DESeq2_DE.html) — produces the input table
-   [Tutorial 08 — Differential Abundance](Tutorial_08_DifferentialAbundance.html) — the *complement*: instead of "what genes change", "do cell-type proportions change?"
-   [scNotebooks Module 04](https://integrativebioinformatics.github.io/scNotebooks/modules/en/module04/module04.html) — covers the same `clusterProfiler` walkthrough on COVID-19 BAL data
-   `clusterProfiler` book: <https://yulab-smu.top/biomedical-knowledge-mining-book/>