---
title: "Tutorial 06 — DESeq2: From Bulk Primer to Pseudobulk DE"
subtitle: "First a clean bulk-RNA-seq walkthrough on `airway`, then pseudobulk DE on `ifnb` STIM vs CTRL"
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

This tutorial covers the **full DESeq2 story** in two parts:

- **Part A — DESeq2 primer on `airway`** (the Bioconductor-shipped 8-sample bulk RNA-seq dataset). Gets the mechanics out of the way on a clean, small bulk dataset before you confront the single-cell aggregation step.
- **Part B — Pseudobulk DE on `ifnb`** (STIM vs CTRL within each cell type). The full single-cell DE workflow: aggregate to *(donor × condition × cell type)* counts and run a proper bulk-style test per cell type.

Paired with **[Lecture 06 — DESeq2: Bulk Fundamentals + Pseudobulk DE](../Lecture_Folder/Lecture_06_DESeq2_DE.html)**.

::: callout-note
**Companion book chapter:** [Chapter 6 — DESeq2: Bulk Fundamentals + Pseudobulk DE](../Resources_Folder/Chapter_06_DESeq2_DE.html) — the long-form prose treatment of this tutorial's material. · **HPC version (Talapas):** [Talapas analysis pipeline — `06_pseudobulk_de.R`](Tutorial_10_Talapas_Pipeline.html)
:::

Part A is a standard bulk DESeq2 walkthrough on the `airway` data. Part B starts from the integrated, annotated object saved at the end of Tutorial 05 (`ifnb_integrated.rds`).

::: {.callout-tip title="How to use this page"}
The rendered HTML shows the code but does **not** execute it. To run it:

1.  Download the `.qmd` source: <a href="Tutorial_06_DESeq2_DE.qmd.txt" download="Tutorial_06_DESeq2_DE.qmd">Tutorial_06_DESeq2_DE.qmd</a>. If your browser saves the file as `Tutorial_06_DESeq2_DE.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  For Part A: install the `airway` Bioconductor package — no other download needed.
3.  For Part B: make sure you've completed **[Tutorial 05](Tutorial_05_Integration.html)** — it writes `ifnb_integrated.rds`.
4.  Work through the chunks, flipping `eval: true` as you go.
:::

::: {.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_06_DESeq2_DE.html`. See [Exercise_Folder/_quarto-solutions.yml](https://github.com/wcresko/scRNAseq_tutorial/blob/main/Exercise_Folder/_quarto-solutions.yml) for the build profile
:::

------------------------------------------------------------------------

# Part A — DESeq2 Primer on `airway`

## Dataset — `airway`

**What it is.** RNA-seq of four primary human airway smooth muscle cell lines, each **untreated** vs. **dexamethasone**-treated → **8 samples**. Standard reference dataset for teaching DESeq2.

-   Bioconductor landing page: <https://bioconductor.org/packages/airway/>
-   Original paper: Himes et al. (2014) *PLoS ONE* 9: e99625

There's no separate download — the counts and metadata ride along with the package, which is installed in the [Tutorial 00 setup](Tutorial_00_Setup_RStudio_Packages.html).

::: callout-tip
**Can't install the `airway` package?** Download the preassembled <a href="../rds/airway_raw.rds" download>`airway_raw.rds`</a> and use `airway <- readRDS("../data/airway_raw.rds")` in place of `data(airway)` below — it's the identical object. Details on the [Datasets](../Datasets.html#core-data-load) page.
:::

::: {.callout-warning title="Common errors / things that bite"}
**Wrong sign on `log2FoldChange`** — `relevel()` matters. The reference level becomes the *denominator*. `dds$dex <- relevel(dds$dex, ref = "untrt")` makes `trt vs untrt` the comparison, with positive LFCs meaning higher in treated. Skip this and your "up" / "down" labels are reversed.

**`lfcShrink(... type = "apeglm")` errors with "requires installing the Bioconductor package 'apeglm'"** — `apeglm` is only a *Suggests* of DESeq2, so it isn't pulled in automatically. It's in the [Tutorial 00 setup](Tutorial_00_Setup_RStudio_Packages.html) list; if you set up earlier, install it once with `BiocManager::install("apeglm")`.

**`apeglm` errors with "coef ... not in resultsNames" or a `fitType` message** — `apeglm` needs the design **coefficient name**, not a contrast vector. This tutorial renames the `dex` levels to `treated`/`untreated`, so the coefficient is `dex_treated_vs_untreated` (not `dex_trt_vs_untrt`). Run `resultsNames(dds)` to see the exact names available.
:::

## A.1 — Setup

```{r}
#| label: M6-setup_airway
library(DESeq2)
library(airway)
library(tidyverse)
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  Mod6_C<chunk>_<name>
# into ../output/Mod6/ so it can be cross-referenced from the rest of the site.
#   Mod6 = Module 6 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod6"
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))

# --- 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), ")"))
}
# save_base_fig() does the same for BASE-graphics figures (plot()/heatmap()/etc.):
# pass a no-argument function that draws the figure; it is rendered to .png and .svg.
save_base_fig <- function(filename, draw, width = 7, height = 5, res = 300) {
  grDevices::png(filename, width = width, height = height, units = "in", res = res)
  draw(); grDevices::dev.off()
  svg_path <- paste0(tools::file_path_sans_ext(filename), ".svg")
  tryCatch({ grDevices::svg(svg_path, width = width, height = height); draw(); grDevices::dev.off() },
           error = function(e) message("  (could not write ", basename(svg_path), ")"))
}
```

## A.2 — Load the data and pull out counts + metadata

```{r}
#| label: M6-load_airway
data(airway)
airway

counts_data <- assay(airway)         # genes × samples integer matrix
colData     <- as.data.frame(colData(airway))

# Tidy the metadata: keep cell line + treatment,
# rename "trt"/"untrt" to something readable.
colData <- colData[, c("cell", "dex")]
colData$dex <- factor(gsub("trt",   "treated",   colData$dex))
colData$dex <- factor(gsub("untrt", "untreated", colData$dex))

head(counts_data[, 1:4])
colData

# --- Tables out: sample metadata + a peek at the counts matrix --------------
colData |>
  tibble::rownames_to_column("sample") |>
  readr::write_csv(file.path(out_dir, "Mod6_C3_airway_coldata.csv"))

as.data.frame(counts_data) |>
  tibble::rownames_to_column("gene") |>
  readr::write_csv(file.path(out_dir, "Mod6_C3_airway_counts.csv"))
```

::: {.callout-important title="Think about it"}
1.  What's in the rows of `counts_data`? In the columns?
2.  Why do we want `dex` as a `factor` rather than a `character`?

<details>
<summary>Show answers</summary>

1.  Rows are **Ensembl gene IDs**; columns are **sample IDs** (`SRR…` accessions). Values are integer read counts.
2.  DESeq2 builds a model matrix from factors and uses **factor levels** to decide which group is the *reference* (denominator) in `log2FoldChange`. A character column would be coerced silently and might pick the wrong reference.
</details>
:::

## A.3 — Sanity-check the metadata

DESeq2 fails fast and unhelpfully if `colData` rows don't line up with `countData` columns. Always check.

```{r}
#| label: M6-sanity_airway
all(colnames(counts_data) %in% rownames(colData))
all(colnames(counts_data) == rownames(colData))
```

::: callout-warning
Both should return `TRUE`. If the second is `FALSE`, **reorder `colData`** to match — don't reorder the counts (you'll forget which orientation you fixed). `colData <- colData[colnames(counts_data), ]`.
:::

## A.4 — Build the `DESeqDataSet`

```{r}
#| label: M6-build_dds_airway
dds <- DESeqDataSetFromMatrix(
  countData = counts_data,
  colData   = colData,
  design    = ~ cell + dex
)
dds
```

::: {.callout-important title="Think about it"}
The design is `~ cell + dex`. Why include `cell` if we only care about `dex`?

<details>
<summary>Show answers</summary>

The eight samples come from **four different cell lines**. Each cell line is a paired (treated, untreated) block. Adding `cell` as a covariate **soaks up the cell-line effect** so the test for `dex` is within-line — the bulk equivalent of a paired analysis.
</details>
:::

## A.5 — Pre-filter low-count genes

```{r}
#| label: M6-prefilter_airway
keep <- rowSums(counts(dds)) >= 10
dds  <- dds[keep, ]
nrow(dds)
```

::: callout-tip
**Pre-filtering is for speed only**, not for FDR control. DESeq2 already does *independent filtering* at the `results()` step, which protects multiple testing. Dropping rows with `< 10` total counts just makes `DESeq()` faster.
:::

## A.6 — Set the reference level

```{r}
#| label: M6-relevel_airway
dds$dex <- relevel(dds$dex, ref = "untreated")
```

::: callout-warning
**This step is easy to forget and easy to misread.** With `ref = "untreated"`, a positive `log2FoldChange` means the gene is **higher in treated**. Without `relevel()`, R's alphabetical default would flip the sign on you.
:::

## A.7 — Run `DESeq()` and pull results

[]{#lst-M6-run_deseq_airway}

```{r}
#| label: M6-run_deseq_airway
dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "treated", "untreated"))
summary(res)

# --- Table out: full DESeq2 results (treated vs untreated) ------------------
as.data.frame(res) |>
  tibble::rownames_to_column("gene") |>
  readr::write_csv(file.path(out_dir, "Mod6_C8_airway_deseq_results.csv"))
```

::: {.callout-tip title="Reading the output"}
The `summary(res)` printout ([code](#lst-M6-run_deseq_airway)) reports the total number of genes tested and how many have `LFC > 0` (up) or `LFC < 0` (down) at the chosen `alpha` threshold (default 0.1). The "outliers" and "low counts" lines tell you how many genes were **removed from testing** before FDR adjustment: outliers have a Cook's-distance flag (one sample is wildly different) and are set to `padj = NA`; low-count genes fail DESeq2's independent filtering step because they lack statistical power at the chosen FDR. A healthy summary shows a handful of outliers (often 0–50), a moderate low-count fraction, and a biologically plausible number of significant genes — for `airway` treated vs untreated, expect several hundred significant genes.
:::

::: {.callout-important title="Think about it"}
`summary(res)` reports counts of "LFC > 0 (up)" and "outliers" and "low counts". What do "outliers" and "low counts" mean here?

<details>
<summary>Show answers</summary>

-   **Outliers** are genes flagged by Cook's distance — a single sample with a wildly different count from the rest. DESeq2 sets `padj = NA` for these so they don't drive results.
-   **Low counts** are genes filtered by *independent filtering*: their mean expression is so low they have no power to be detected at the chosen FDR. DESeq2 estimates the optimal cut adaptively.
</details>
:::

## A.8 — Shrink the LFCs before ranking / plotting

```{r}
#| label: M6-lfcshrink_airway
res_shrunk <- lfcShrink(dds,
                        coef = "dex_treated_vs_untreated",
                        type = "apeglm")
plotMA(res_shrunk, ylim = c(-3, 3),
       main = "airway: shrunken LFC (treated vs untreated)",
       xlab = "Mean of normalized counts",
       ylab = "log2 fold change (apeglm-shrunken)")

# --- Figure out: base-graphics MA plot via png device ----------------------
save_base_fig(file.path(out_dir, "Mod6_C9_airway_ma_plot.png"), width = 7, height = 5, draw = function() {
  plotMA(res_shrunk, ylim = c(-3, 3),
         main = "airway: shrunken LFC (treated vs untreated)",
         xlab = "Mean of normalized counts",
         ylab = "log2 fold change (apeglm-shrunken)")
})
```

::: {.callout-tip title="Reading the output"}
The x-axis is average normalized expression (`baseMean`); the y-axis is the apeglm-shrunken `log2FoldChange` (positive = higher in treated). Each point is a gene; red points are significant at the chosen FDR. Notice that the cloud of points is narrow and centered near zero for low-`baseMean` genes — that is shrinkage working: raw fold changes for lowly-expressed genes are pulled back toward zero because they carry little information. Well-expressed, genuinely regulated genes (ISGs in `ifnb`, corticosteroid-response genes in `airway`) will sit away from the band and remain red even after shrinkage. A healthy MA plot looks like a horizontal band with a handful of off-axis red points; a funnel that opens up at low baseMean without shrinkage signals that you should always apply `lfcShrink` before ranking.
:::

::: callout-tip
**Always rank by shrunk LFCs.** Genes with `baseMean = 2` can show wild raw LFCs purely by sampling noise; `apeglm` shrinks those toward zero and de-clutters your top-hits list.
:::

## A.9 — Inspect a top hit

[]{#lst-M6-top_hit_airway}

```{r}
#| label: M6-top_hit_airway
#| fig-cap: "Per-sample normalized counts for the top 5 airway hits (by padj), each panel labelled with its gene, split by treatment."
top_genes <- as.data.frame(res_shrunk) |>
  rownames_to_column("gene") |>
  arrange(padj) |>
  head(5)
top_genes

# Gather per-sample normalized counts for the top 5 genes (plotCounts can return
# the data instead of drawing), so we can plot them together with gene labels.
counts_top <- purrr::map_dfr(top_genes$gene, function(g) {
  d <- plotCounts(dds, gene = g, intgroup = "dex", returnData = TRUE)
  d$gene <- g
  d
})
# Keep the panels ordered by significance (padj rank), not alphabetically.
counts_top$gene <- factor(counts_top$gene, levels = top_genes$gene)

p_top <- ggplot(counts_top, aes(x = dex, y = count, colour = dex)) +
  geom_jitter(width = 0.12, height = 0, size = 2, alpha = 0.85) +
  facet_wrap(~ gene, nrow = 1, scales = "free_y") +   # facet strip = the gene label
  scale_y_log10() +
  labs(
    title    = "Top 5 airway hits (by padj) — per-sample normalized counts",
    subtitle = "Each panel is a gene; points are samples, split by treatment",
    x        = "Treatment (dex)",
    y        = "Normalized count (log scale)",
    colour   = "Treatment"
  ) +
  theme(legend.position = "none")
p_top

# --- Table out: top 5 hits by padj -----------------------------------------
readr::write_csv(top_genes, file.path(out_dir, "Mod6_C10_airway_top_hits.csv"))

# --- Figure out: per-sample counts for the top 5 hits (png + svg) -----------
save_fig(file.path(out_dir, "Mod6_C10_airway_top_hit_counts.png"), p_top,
         width = 11, height = 4, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each panel ([code](#lst-M6-top_hit_airway)) shows one gene (ordered left to right by significance); the x-axis splits by treatment (`untreated` vs `treated`) and each point is one sample on a log-scale y-axis of normalized counts. What to look for: a clear, consistent shift in the median across all four cell lines — if any panel has one outlier sample dragging the effect, that gene deserves closer scrutiny before reporting it. For `airway`, the top hits should show tight, concordant up-regulation in all four treated replicates. A gene where two treated replicates are high and two are near the untreated level is a borderline call, regardless of the p-value.
:::

::: {.callout-important title="Think about it"}
Why is `plotCounts()` a more honest way to confirm a hit than just trusting the p-value?

<details>
<summary>Show answers</summary>

It shows the **per-sample** values. A single outlier sample can drag the model into "significance" that disappears under inspection. Eyeballing the per-sample dots beats a single p-value every time.
</details>
:::

## A.10 — Sanity check by transformation + PCA

Quality-control your samples *post-hoc* with a variance-stabilizing transform (`vst`) and a PCA:

```{r}
#| label: M6-vst_pca_airway
vsd <- vst(dds, blind = FALSE)
p_pca <- plotPCA(vsd, intgroup = c("dex", "cell")) +
  labs(
    title    = "airway: sample PCA on variance-stabilized counts",
    subtitle = "Look for treated vs untreated separation within each cell line",
    colour   = "Treatment · cell line"
  )
p_pca

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod6_C11_airway_vst_pca.png"), p_pca,
       width = 7, height = 5, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each point is one of the eight samples; the axes are PC1 and PC2 from the variance-stabilized count matrix. Look for two things: (1) samples of the same treatment (`treated` vs `untreated`) should cluster together, ideally along the same PC; (2) the four cell lines should form paired clusters, with treated and untreated replicates from each line sitting near each other. In `airway`, PC1 typically captures the treatment effect and PC2 the cell-line differences — exactly the structure the `~ cell + dex` design exploits. If samples cluster purely by cell line with no treatment separation, the treatment effect is too small relative to the between-line variance to trust the DE calls.
:::

::: callout-tip
**What you want to see.** Treated and untreated separated on PC1 or PC2 *within* each cell line. If samples cluster by cell line only and not by `dex`, your treatment effect is small relative to between-line variation — interpret the DE results accordingly.
:::

## Part A — What you've practiced

-   Building a `DESeqDataSet` from a counts matrix + sample table
-   Adding a covariate to a design
-   Reading DESeq2 output (LFC, lfcSE, padj, baseMean)
-   Shrinking LFCs for honest ranking
-   QCing samples with `vst` + PCA

You're now ready for **Part B**, where the same DESeq2 mechanics apply to a per-cell-type aggregated count matrix.

------------------------------------------------------------------------

# Part B — Pseudobulk DE on `ifnb` (STIM vs CTRL)

We start from the **integrated, annotated object saved at the end of Tutorial 05** (`ifnb_integrated.rds`) and learn to:

-   Aggregate single-cell counts to a *(donor × condition × cell type)* pseudobulk matrix
-   Run a proper, well-calibrated DE test with `DESeq2`, one cell type at a time
-   Apply LFC shrinkage for ranking and visualize results
-   Avoid the most common pseudobulk pitfalls

::: callout-note
**The contrast is real this time.** The `ifnb` dataset has a built-in STIM (IFN-β stimulated) vs CTRL (control) contrast. In Tutorial 05 we **integrated** to remove the STIM-vs-CTRL axis from the embedding so cell types co-cluster. Here we **measure** the STIM-vs-CTRL effect *within each cell type* — the biology that integration was careful **not** to erase.
:::

## Dataset — resuming from Tutorial 05

| File | Produced by | What it is |
|---|---|---|
| `../data/ifnb_integrated.rds` | Tutorial 05 | Seurat object with integrated PCA / UMAP / clusters, the original `RNA` assay restored as default, and `seurat_annotations` cell-type labels per cell. |

### About donor information

The Kang *et al.* (2017) study profiled PBMCs from **eight lupus donors**, each split into a control and an IFN-β–stimulated aliquot. The `muscData::Kang18_8vs8()` loader we use in [Tutorial 01](Tutorial_01_QC_Preprocessing.html) **carries the real per-cell donor IDs** in the `ind` column, so we have genuine biological replicates: eight donors, each measured in both conditions.

Pseudobulk DE needs **biological replicates within each condition**, and here we have them: the eight real donors become the replicates, one pseudobulk profile per *(donor × condition × cell type)*. The per-cell-type test below is `~ condition`; because each donor appears in both conditions you can upgrade to the paired `~ donor + condition` (see *Going further*). The code detects the donor column automatically:

-   **`ind` is present (the expected case)** — we use the eight real donors as biological replicates. This is what you'll see when you load `ifnb` as in Tutorial 01.
-   **Fallback only — donor column absent** — if an object somehow reaches this step without donor metadata (e.g. an older SeuratData-sourced `ifnb`), the code creates **synthetic pseudo-donors** by random cell assignment. This is a teaching shortcut, **not** something you'd ever do in a paper: the contrast still lets DESeq2 fit, but the resulting p-values aren't biologically meaningful.

::: {.callout-warning title="Common errors / things that bite"}
**`AggregateExpression` returns a different structure than expected** — the API changed across Seurat versions. In v5+ it returns an Assay5 object; subset with `$RNA$counts`. In v4, it returns a sparse matrix list — use `$RNA`. The tutorial's parsing handles both, but if you're on a custom branch this can bite.

**Underscores in donor or cell-type names break the column-name parsing** — if your sample is `donor_03_treated`, the `sub("^[^_]+_", "", colnames(pb))` parse will get confused. Defensive parsing approach: pull the metadata directly from the Seurat object's `@meta.data` rather than parsing column names. The tutorial shows this.

**`colData` comes back all `NA` / numeric donor IDs get a `g` prefix** — `ifnb`'s donor column (`ind`) is **numeric** (101, 107, …), and `AggregateExpression()` prepends `g` to any `group.by` value that starts with a digit (it prints *"…starts with a number, appending 'g'…"*). The aggregated columns become `g101_…`, which no longer match a reconstructed `101_…` key, so the join yields all `NA` and every group is dropped. Fix: make donor a valid name up front — the tutorial uses `ifnb$donor <- paste0("d", ind)`.

**`count(...)` errors with "Argument 'x' is not a vector: list"** — loading `DESeq2`/`SummarizedExperiment` attaches `matrixStats`, whose `count()` **masks** `dplyr::count()` (along with `desc`, `slice`, `rename`, `first` from S4Vectors/IRanges). Call the dplyr ones explicitly, e.g. `dplyr::count(...)`. Watch the package-startup messages for the list of masked names.

**Volcano plot errors with "Faceting variables must have at least one value"** — the cell types you asked to plot don't exist in `de_all$celltype`. The `muscData` labels are `CD14+ Monocytes`, `CD4 T cells`, `B cells`, … — not the older `CD14 Mono`/`B`. The tutorial now selects the cell types **data-driven** (the 3 with the most significant genes), so it adapts to your labels; if you hardcode names, match them to `unique(de_all$celltype)`.

**Benign warnings during the per-cell-type loop** — `converting counts to integer mode` (pseudobulk sums stored as doubles; the loop now `round()`s them), `some variables in design formula are characters, converting to factors` (silenced by making `condition` an explicit factor), and `nbinomGLM … line search routine failed` (an `apeglm` optimizer note on a few noisy genes — harmless, not fixable by factor coercion). The DE table is correct regardless.

**Synthetic pseudo-donors are a fallback, not the default** — objects loaded the way this series does (`muscData::Kang18_8vs8()`, Tutorial 01) carry real donor IDs in `ind`, so you'll take the real-donor path automatically. The synthetic `set.seed(2026)` 4-pseudo-donor split only triggers if no donor column is present (e.g. an older SeuratData `ifnb`); it's a teaching device whose p-values are not biologically meaningful — say so in any writeup.
:::

## B.1 — Setup

```{r}
#| label: M6-setup_pb

library(Seurat)
library(DESeq2)
library(tidyverse)
library(patchwork)

set.seed(2026)

ifnb <- readRDS("../data/ifnb_integrated.rds")
DefaultAssay(ifnb) <- "RNA"
ifnb
```

```{r}
#| label: M6-peek_meta

# Look for a real donor column
candidate_cols <- intersect(c("donor_id","donor","ind"), colnames(ifnb@meta.data))
candidate_cols
```

## B.2 — Define donor (real or synthetic) and condition

```{r}
#| label: M6-donor

if (length(candidate_cols) > 0) {
  # Prefix with a letter: the ifnb donor IDs (`ind`) are numbers (101, 107, ...),
  # and AggregateExpression() prepends "g" to any group.by value that starts with
  # a digit. Making donor a valid name up front ("d101") keeps the aggregated
  # column names stable so the colData join below matches.
  ifnb$donor <- paste0("d", ifnb@meta.data[[ candidate_cols[1] ]])
  message("Using real donor column: ", candidate_cols[1])
} else {
  # Synthetic pseudo-donors: 4 per stim, assigned by random cell index
  set.seed(2026)
  ifnb$donor <- paste0(
    ifnb$stim, "_d",
    sample(1:4, size = ncol(ifnb), replace = TRUE)
  )
  message("Created synthetic pseudo-donors (teaching shortcut). ",
          "DO NOT do this for a real publication.")
}

ifnb$condition <- ifnb$stim   # CTRL or STIM
ifnb$celltype  <- ifnb$seurat_annotations

table(ifnb$donor, ifnb$condition)

# --- Table out: cells per donor x condition --------------------------------
as.data.frame(table(donor = ifnb$donor, condition = ifnb$condition)) |>
  readr::write_csv(file.path(out_dir, "Mod6_C14_donor_by_condition.csv"))
```

::: {.callout-important title="Think about it"}
Why is it not enough to just run `FindMarkers(ifnb, ident.1 = "STIM", ident.2 = "CTRL")` on the per-cell expression matrix?

<details>

<summary>Click for answer</summary>

Because cells from one donor are **not statistically independent**. Running a per-cell test treats every cell as a separate replicate, which inflates the effective sample size and produces **anti-conservative** p-values — you'll get hundreds of "significant" genes that wouldn't replicate. The number of independent biological replicates is the **number of donors** (or pseudo-donors), not the number of cells. Pseudobulk fixes this by aggregating to one observation per (donor × condition × cell type).

For a thorough demonstration of how badly per-cell DE fails on this dataset specifically, see Squair *et al.* 2021 (Nat Commun, [doi](https://doi.org/10.1038/s41467-021-25960-2)).

</details>
:::

## B.3 — Aggregate to pseudobulk

```{r}
#| label: M6-aggregate

pb <- AggregateExpression(
  ifnb,
  assays    = "RNA",
  slot      = "counts",
  group.by  = c("donor", "condition", "celltype"),
  return.seurat = FALSE
)$RNA

dim(pb)               # genes × (donor × condition × celltype) replicates
head(colnames(pb), 6)
```

The column names look like `<donor>_<condition>_<celltype>`. Build the matching `colData`:

```{r}
#| label: M6-colData

# Seurat's AggregateExpression replaces "_" with "-" by default in some versions.
# Be defensive and parse from the metadata directly instead.
group_meta <- ifnb@meta.data |>
  distinct(donor, condition, celltype) |>
  mutate(group_id = paste(donor, condition, celltype, sep = "_"))

# Match the colnames that AggregateExpression actually used. Different
# Seurat versions differ; this should align them.
fix_id <- function(x) gsub("[ /]", "-", x)
group_meta$group_id_clean <- fix_id(group_meta$group_id)
colnames(pb) <- fix_id(colnames(pb))

meta_pb <- tibble(group_id_clean = colnames(pb)) |>
  left_join(group_meta, by = "group_id_clean")

head(meta_pb)

# --- Table out: pseudobulk column metadata (donor x condition x celltype) ---
readr::write_csv(meta_pb, file.path(out_dir, "Mod6_C16_pseudobulk_coldata.csv"))
```

## B.4 — Filter low-count groups

A pseudobulk column built from very few cells is noise. Drop them.

```{r}
#| label: M6-filter_pb

# NOTE: use dplyr::count explicitly — loading DESeq2/SummarizedExperiment attaches
# matrixStats, whose count() masks dplyr::count() (you'd get "Argument 'x' is not
# a vector: list").
cells_per_group <- ifnb@meta.data |>
  dplyr::count(donor, condition, celltype) |>
  mutate(group_id_clean = fix_id(paste(donor, condition, celltype, sep = "_")))

meta_pb <- meta_pb |>
  left_join(cells_per_group |> select(group_id_clean, n_cells = n),
            by = "group_id_clean")

# Keep groups with at least 10 cells
keep    <- !is.na(meta_pb$n_cells) & meta_pb$n_cells >= 10
pb      <- pb[, keep]
meta_pb <- meta_pb[keep, ]

table(meta_pb$celltype, meta_pb$condition)

# --- Table out: surviving pseudobulk groups per celltype x condition --------
as.data.frame(table(celltype = meta_pb$celltype, condition = meta_pb$condition)) |>
  readr::write_csv(file.path(out_dir, "Mod6_C17_groups_per_celltype.csv"))
```

::: {.callout-important title="Think about it"}
Why is `n_cells >= 10` reasonable? Why might you raise or lower it?

<details>

<summary>Click for answer</summary>

The threshold balances signal against noise: pseudobulk built from 1–2 cells is dominated by stochastic dropout and is unlikely to represent the cell type. A floor of \~10 cells is the de-facto default in the field (see Squair *et al.* 2021 and Crowell *et al.* 2020). Raise it (e.g. 20–50) for very sparse cell types or shallow libraries; lower it only if you also lower your statistical-significance expectations.

</details>
:::

## B.5 — Run DESeq2 per cell type

```{r}
#| label: M6-run_deseq_pb

run_de <- function(ct) {
  keep <- meta_pb$celltype == ct
  if (sum(keep) < 4 || length(unique(meta_pb$condition[keep])) < 2) return(NULL)

  cd <- as.data.frame(meta_pb[keep, ])
  # Make condition an explicit factor with CTRL as the reference. This silences
  # DESeq2's "variables in design formula are characters, converting to factors"
  # note and guarantees the STIM-vs-CTRL direction.
  cd$condition <- factor(cd$condition, levels = c("CTRL", "STIM"))

  dds <- DESeqDataSetFromMatrix(
    countData = round(pb[, keep]),   # round() avoids the "converting counts to integer mode" note
    colData   = cd,
    design    = ~ condition
  )
  dds <- DESeq(dds, quiet = TRUE)

  res <- results(dds, contrast = c("condition","STIM","CTRL"))
  res <- lfcShrink(dds, coef = "condition_STIM_vs_CTRL",
                   res = res, type = "apeglm")

  as.data.frame(res) |>
    rownames_to_column("gene") |>
    mutate(celltype = ct)
}

de_all <- map_dfr(unique(meta_pb$celltype), run_de) |>
  filter(!is.na(padj))
glimpse(de_all)

# --- Table out: full per-cell-type pseudobulk DE results --------------------
readr::write_csv(de_all, file.path(out_dir, "Mod6_C18_pseudobulk_de_all.csv"))
```

## B.6 — Inspect the top hits

```{r}
#| label: M6-top_hits_pb

de_top <- de_all |>
  group_by(celltype) |>
  slice_min(padj, n = 10) |>
  ungroup()
de_top

# --- Table out: top 10 hits per cell type by padj --------------------------
readr::write_csv(de_top, file.path(out_dir, "Mod6_C19_pseudobulk_de_top.csv"))
```

::: {.callout-tip title="Reading the output"}
The printed tibble is grouped by `celltype`; within each group the rows are ordered by `padj` (smallest first). The `log2FoldChange` column uses apeglm-shrunken estimates: a value of 2 means the gene is four-fold higher in STIM than CTRL on average across donors. `baseMean` is the average normalized count across all pseudobulk samples for that gene — a very low `baseMean` with a large `log2FoldChange` is a warning sign (little evidence, large noise). The key sanity check: are the top rows dominated by well-known interferon-stimulated genes (`ISG15`, `MX1`, `OAS1`, `IFIT1`)? If so, the contrast is working exactly as expected.
:::

::: {.callout-important title="Think about it"}
1.  Many of the top hits per cell type will be **interferon-stimulated genes** (`ISG15`, `IFI6`, `IFIT1/2/3`, `MX1`, `OAS1`, `RSAD2`). Are they up in STIM or CTRL?
2.  Are **the same** ISGs at the top of every cell type, or do you see cell-type-specific responses?

<details>

<summary>Click for answer</summary>

1.  Up in STIM — the IFN-β response is what STIM is, mechanistically. `log2FoldChange > 0` means higher in STIM (because we set `contrast = c("condition","STIM","CTRL")` with STIM as the numerator).
2.  A core ISG cluster (`ISG15`, `IFI6`, `MX1`) tops every cell type, **but the magnitude varies**. Monocytes (CD14, CD16) typically show the **strongest** ISG induction; T and B cells respond more modestly; some rare populations may be near-flat. This cell-type-specific quantitative difference is exactly the kind of biology pseudobulk DE is designed to find — and the kind of biology that the per-cell `FindMarkers` test cannot resolve confidently.

</details>
:::

## B.7 — Visualize

A simple per-cell-type volcano plot:

```{r}
#| label: M6-volcano

# Pick the 3 cell types with the most significant genes — data-driven, so it works
# whatever your cell-type labels are (the muscData ifnb uses "CD14+ Monocytes",
# "CD4 T cells", "B cells", ... — not the older "CD14 Mono"/"B" names).
plot_types <- de_all |>
  dplyr::filter(padj < 0.05) |>
  dplyr::count(celltype, sort = TRUE) |>
  head(3) |>
  dplyr::pull(celltype)
if (length(plot_types) == 0) plot_types <- head(unique(de_all$celltype), 3)

p_volcano <- de_all |>
  filter(celltype %in% plot_types) |>
  ggplot(aes(log2FoldChange, -log10(padj))) +
  geom_point(alpha = 0.4, size = 1) +
  geom_hline(yintercept = -log10(0.05), linetype = 2) +
  geom_vline(xintercept = c(-1, 1), linetype = 3, color = "grey50") +
  facet_wrap(~ celltype) +
  labs(title    = "Pseudobulk DE: STIM vs CTRL (log2FC > 0 = up in STIM)",
       subtitle = "Dashed line: padj = 0.05; dotted lines: |log2FC| = 1",
       x        = "log2 fold change (apeglm-shrunken)",
       y        = "-log10(padj)",
       caption  = "Module 6 · DESeq2 pseudobulk DE")
p_volcano

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod6_C20_pseudobulk_volcano.png"), p_volcano,
       width = 10, height = 4, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each panel is one cell type; each point is one gene. The x-axis is the apeglm-shrunken `log2FoldChange` (positive = higher in STIM), and the y-axis is `-log10(padj)` — so genes in the upper-right corner are strongly induced in STIM with high statistical confidence. The horizontal dashed line marks `padj = 0.05`; the vertical dotted lines mark `|log2FC| = 1` (twofold change). For `ifnb`, you should see a dense cloud of significant points in the upper-right of monocyte panels (strong ISG response) and a sparser, flatter cloud in T and B cell panels. Points that are statistically significant but close to the x-axis origin have large sample sizes and small effects — biologically real but possibly not the most actionable hits. An enrichment tool (Tutorial 07) will help you interpret them collectively.
:::

::: {.callout-important title="Think about it"}
Some cell types will show many DE genes; others almost none. What can drive this pattern *that has nothing to do with biology*?

<details>

<summary>Click for answer</summary>

Three usual suspects:

1.  **Cell counts per group.** Cell types that are abundant in every donor have well-estimated dispersions; rare types do not.
2.  **Sample size.** With 4 pseudo-donors split 4-vs-4 (or 8-vs-8 with real donors), you have moderate power on a clean signal — but very little if a cell type is missing in some donors.
3.  **Library composition.** A cell type whose library is dominated by a few hyper-expressed genes will have noisier dispersion estimates and fewer detected hits.

Always inspect `n_cells` per group and the dispersion plot before reporting "X genes are significantly DE".

</details>
:::

## B.8 — Save results

```{r}
#| label: M6-save_de

# Pipeline hand-off: Tutorial 07 reads these from data/ — keep them there.
write_csv(de_all,  "../data/ifnb_pseudobulk_de.csv")
write_csv(de_top,  "../data/ifnb_pseudobulk_de_top.csv")

# --- Tables out: student-facing copies in the module output dir -------------
readr::write_csv(de_all, file.path(out_dir, "Mod6_C21_pseudobulk_de.csv"))
readr::write_csv(de_top, file.path(out_dir, "Mod6_C21_pseudobulk_de_top.csv"))
```

## Wrap-up

You now have:

-   A working understanding of the DESeq2 mechanics from a clean bulk dataset (`airway`).
-   A pseudobulk matrix per *(donor, condition, cell type)* and per-cell-type DE results with shrunken LFCs.
-   A workflow you can re-point at *any* condition contrast in your own data — replace `condition = stim` with your own covariate and you're done.

::: {.callout-tip title="Going further"}
-   Add covariates (donor, sex, age) to the design: `~ donor + condition` (when donors are real).
-   Use `IHW` or `ashr` for adaptive multiple-testing.
-   **Interaction terms.** Try `design = ~ cell + dex + cell:dex` on `airway` and compare results.
-   **Volcano plot.** Plot `-log10(padj)` vs. `log2FoldChange` from `res_shrunk` and label your top 20 hits.
-   Hand the gene lists to `clusterProfiler::compareCluster` for cross-cell-type enrichment — see the GO-enrichment section of [scNotebooks Module 04](https://integrativebioinformatics.github.io/scNotebooks/modules/en/module04/module04.html) for a worked example.
-   **Differential abundance** with `miloR` — the *complement* to differential expression. DE asks "for cells of the same type, what genes change between conditions?" Differential abundance asks "do the *cell-type proportions themselves* change between conditions?"
-   Compare against `muscat::pbDS` (formal pseudobulk wrapper) and `glmGamPoi`.
:::

## Continue to Tutorial 07 — functional interpretation

The DE table you wrote (`ifnb_pseudobulk_de.csv`) is the input to **[Tutorial 07 — Functional Analysis with `clusterProfiler`](Tutorial_07_FunctionalAnalysis.html)**, which runs GO over-representation, GSEA, and Reactome enrichment per cell type and shows you how to compare the biological themes across cell types.

For the *complementary* condition-level question — "do cell-type **proportions** change?" rather than "which genes change?" — see **[Tutorial 08 — Differential Abundance with miloR](Tutorial_08_DifferentialAbundance.html)**.

## See also

-   [Lecture 06 — DESeq2: Bulk Fundamentals + Pseudobulk DE](../Lecture_Folder/Lecture_06_DESeq2_DE.html)
-   Squair *et al.* 2021 — the paper that made pseudobulk the default ([doi](https://doi.org/10.1038/s41467-021-25960-2))
-   Crowell *et al.* 2020 — `muscat` package paper ([doi](https://doi.org/10.1038/s41467-020-19894-4))

## Credits

-   Data: [`airway`](https://bioconductor.org/packages/airway/) (Himes et al., 2014, *PLoS ONE*)