Tutorial 06 — DESeq2: From Bulk Primer to Pseudobulk DE

First a clean bulk-RNA-seq walkthrough on airway, then pseudobulk DE on ifnb STIM vs CTRL

NoteRunning 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.

Note

Companion book chapter: Chapter 6 — DESeq2: Bulk Fundamentals + Pseudobulk DE — the long-form prose treatment of this tutorial’s material. · HPC version (Talapas): Talapas analysis pipeline — 06_pseudobulk_de.R

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).

TipHow to use this page

The rendered HTML shows the code but does not execute it. To run it:

1. Download the .qmd source: Tutorial_06_DESeq2_DE.qmd. 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 — it writes ifnb_integrated.rds.
4. Work through the chunks, flipping eval: true as you go.

TipSolutions / 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 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.

Tip

Can’t install the airway package? Download the preassembled airway_raw.rds and use airway <- readRDS("../data/airway_raw.rds") in place of data(airway) below — it’s the identical object. Details on the Datasets page.

WarningCommon 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 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

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

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

ImportantThink 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?

Show answers

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.

A.3 — Sanity-check the metadata

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

all(colnames(counts_data) %in% rownames(colData))
all(colnames(counts_data) == rownames(colData))

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

dds <- DESeqDataSetFromMatrix(
  countData = counts_data,
  colData   = colData,
  design    = ~ cell + dex
)
dds

ImportantThink about it

The design is ~ cell + dex. Why include cell if we only care about dex?

Show answers

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.

A.5 — Pre-filter low-count genes

keep <- rowSums(counts(dds)) >= 10
dds  <- dds[keep, ]
nrow(dds)

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

dds$dex <- relevel(dds$dex, ref = "untreated")

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

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

TipReading the output

The summary(res) printout (code) 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.

ImportantThink about it

summary(res) reports counts of “LFC > 0 (up)” and “outliers” and “low counts”. What do “outliers” and “low counts” mean here?

Show answers

• 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.

A.8 — Shrink the LFCs before ranking / plotting

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)")
})

TipReading 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.

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

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)

TipReading the output

Each panel (code) 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.

ImportantThink about it

Why is plotCounts() a more honest way to confirm a hit than just trusting the p-value?

Show answers

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.

A.10 — Sanity check by transformation + PCA

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

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)

TipReading 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.

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

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 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.

WarningCommon 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

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

set.seed(2026)

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

# 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

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

ImportantThink about it

Why is it not enough to just run FindMarkers(ifnb, ident.1 = "STIM", ident.2 = "CTRL") on the per-cell expression matrix?

Click for answer

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).

B.3 — Aggregate to pseudobulk

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:

# 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.

# 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"))

ImportantThink about it

Why is n_cells >= 10 reasonable? Why might you raise or lower it?

Click for answer

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.

B.5 — Run DESeq2 per cell type

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

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

TipReading 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.

ImportantThink 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?

Click for answer

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.

B.7 — Visualize

A simple per-cell-type volcano plot:

# 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)

TipReading 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.

ImportantThink 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?

Click for answer

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”.

B.8 — Save results

# 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.

TipGoing 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 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, 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.

See also

• Lecture 06 — DESeq2: Bulk Fundamentals + Pseudobulk DE
• Squair et al. 2021 — the paper that made pseudobulk the default (doi)
• Crowell et al. 2020 — muscat package paper (doi)

Credits

• Data: airway (Himes et al., 2014, PLoS ONE)