---
title: "Tutorial 01 — scRNA-seq: QC & Preprocessing"
subtitle: "Load a two-sample dataset (ifnb), filter cells, normalize, pick HVGs, scale"
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

The first hands-on tutorial walking the standard **Seurat** scRNA-seq workflow on a **laptop-friendly** dataset. We use the **`ifnb`** dataset (Kang *et al.* 2017) — peripheral blood mononuclear cells (PBMCs) from eight lupus patients, split into a **control** sample and an **interferon-β-stimulated** sample. Two samples, \~24,000 cells, no 10x click-through.

::: callout-note
**Companion lecture:** [Lecture 01 — Pre-processing: QC, doublets, normalization, HVGs](../Lecture_Folder/Lecture_01_Preprocessing.html) · **Companion book chapter:** [Chapter 1 — Preprocessing](../Resources_Folder/Chapter_01_Preprocessing.html) — the long-form prose treatment of this tutorial's material, with cross-references to the prerequisite appendices. · **HPC version (Talapas):** [Talapas analysis pipeline — `01_qc_preprocessing.R`](Tutorial_10_Talapas_Pipeline.html)
:::

Why two samples from the start? Because Tutorial 05 will teach **integration** — and integration only makes sense if you've already loaded, QC-filtered, and clustered samples that come from different conditions / batches. We carry both `STIM` and `CTRL` cells through every step.

::: callout-tip
**Want to run this on the cluster, or scale up later?** The [**Talapas HPC track**](Tutorial_10_Talapas_Pipeline.html) runs this *same* `ifnb` workflow as R + SLURM scripts — see [Tutorial — SLURM Basics](Tutorial_09_SLURM_Basics.html). To go all the way from **raw FASTQs** on a larger real dataset, see the optional [Full Talapas run from raw FASTQs](../Resources_Folder/Talapas_FullRun_FromFASTQs.html) (a self-paced bonus track).
:::

::: {.callout-tip title="How to use this page"}
The rendered HTML shows the code but does **not** execute it (chunks default to `eval: false`). To actually run the analysis:

1.  Download the `.qmd` source: <a href="Tutorial_01_QC_Preprocessing.qmd.txt" download="Tutorial_01_QC_Preprocessing.qmd">Tutorial_01_QC_Preprocessing.qmd</a>. If your browser saves the file as `Tutorial_01_QC_Preprocessing.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  Install the `ifnb` dataset once (one R command — see [Setup](#setup)).
3.  Work through the chunks, flipping `eval: true` as you go.
:::

## Dataset — `ifnb` (Kang *et al.* 2017)

**What it is.** Peripheral blood mononuclear cells (PBMCs) from eight lupus patients. Each donor's PBMCs were split into two aliquots:

- **`CTRL`** — control (no stimulation), \~12,000 cells
- **`STIM`** — treated with interferon-β for 6 hours, \~12,000 cells

Together \~24,000 annotated singlets across two samples on the **10x Chromium 3' v1** assay (the `muscData` build of `ifnb` is larger than the older SeuratData `ifnb` you may have seen quoted at \~14k). **Cell-type labels** (e.g. `CD14+ Monocytes`, `CD4 T cells`, `B cells`, `NK cells`, `Dendritic cells`) are bundled in the metadata so you can sanity-check clustering and integration without doing your own annotation.

- **Reference:** Kang *et al.* (2017) *Nat Biotechnol* 36: 89–94. [doi:10.1038/nbt.4042](https://doi.org/10.1038/nbt.4042)
- **Origin:** GEO accession [GSE96583](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96583)
- **Distributed via:** the Bioconductor [`muscData`](https://bioconductor.org/packages/muscData/) package as `Kang18_8vs8()` — cached locally through `ExperimentHub` on first download
- **Used in:** Tutorials 01 → 04 of this workshop, and the official [Seurat integration vignette](https://satijalab.org/seurat/articles/integration_introduction.html)

::: callout-note
**Why this dataset for teaching.** It's small enough to run on a laptop, it has a clean two-condition structure (so integration is meaningful), and the cell types are PBMCs — well-known, well-marked, and easy to reason about. The interferon response is a strong, biologically real signal that drives both batch-style separation (use case for integration) **and** condition-level differential expression (use case for pseudobulk DE in Tutorial 06).
:::

### Disk-space budget

| Item                                     | Size          |
|------------------------------------------|---------------|
| `ifnb` raw counts (compressed cache)     | \~25 MB       |
| `ifnb` Seurat object in memory           | \~600 MB      |
| `.rds` snapshots saved between tutorials | \~500 MB each |

Plan for **\~1.5 GB total** including intermediate `.rds` files for Tutorials 01–08.

### Data dictionary — what's in the loaded `ifnb` Seurat object

| Slot | Class | Contents |
|----|----|----|
| `ifnb@assays$RNA@counts` | `dgCMatrix` | Raw UMI counts. \~35,000 genes × \~24,000 cells (sparse). |
| `ifnb@meta.data$stim` | factor | `"CTRL"` or `"STIM"` — the **sample-of-origin** label. This becomes the integration covariate. |
| `ifnb@meta.data$seurat_annotations` | factor | Author-curated cell-type labels from the original paper. Useful as ground truth in Tutorial 03. |
| `ifnb@meta.data$nCount_RNA`, `nFeature_RNA` | numeric | Auto-populated by `CreateSeuratObject` — UMIs and unique genes per cell. |

## Learning goals

By the end of this tutorial you will be able to:

- Fetch a published scRNA-seq dataset from Bioconductor's `muscData` package and assemble it as a Seurat object
- Inspect a multi-sample Seurat object and split it by sample
- Add per-cell QC metrics (`percent.mt`, `percent.rb`) and visualize them per sample
- Choose filtering thresholds that are defensible for *your* dataset
- Run `NormalizeData` → `FindVariableFeatures` → `ScaleData` and understand what each does
- Save the QC'd object so Tutorial 02 can pick it up

::: {.callout-warning title="Common errors / things that bite"}
**`library(muscData)` errors with "there is no package called 'muscData'"** — it should already be installed from the [Tutorial 00 setup](Tutorial_00_Setup_RStudio_Packages.html); if it's missing, re-run that setup block.

**`scDblFinder` returns a different result than the workshop** — `scDblFinder` is stochastic. The workshop uses `set.seed(2026)` *immediately before* calling it; if you set the seed earlier or not at all, the doublet calls will differ. The 5–10% range is what to expect; exact cell IDs will vary.

**SoupX errors with "no raw matrix found"** (or `subscript out of bounds` from `load10X`) — this is correct behaviour. The `Kang18_8vs8()` download only ships the *filtered* counts. SoupX needs both. The SoupX block ([code](#lst-M1-soupx)) is an **illustrative template** for your own 10x data and now stops with a clear message if you run it with the placeholder path. Run it for real on data that has a Cell Ranger `raw_feature_bc_matrix/` folder.

**Warnings (not errors) during doublet detection** — when Step 6 runs you'll see `Layer 'data'/'scale.data' is empty` and `'normalizeCounts'/'librarySizeFactors' is deprecated`. Both are harmless: doublet detection runs on **raw counts before normalization** (so those layers are legitimately empty), and the deprecation notes come from `scran`/`scuttle` *inside* `scDblFinder`, not your code. As long as you see `Creating … artificial doublets`, it's working.

**Plotting warnings before normalization (harmless)** — the QC plots ([code](#lst-M1-qc_plots)) and doublet plots ([code](#lst-M1-dbl)) draw per-cell QC *metadata* (`nFeature_RNA`, `nCount_RNA`, …) before normalization, so Seurat prints `Default search for "data" layer … utilizing "counts" layer instead`. The QC values are metadata, so the layer it uses is irrelevant. Likewise the HVG plot ([code](#lst-M1-hvgs)) uses a log10 axis and prints `log-10 transformation introduced infinite values` for genes with zero mean expression (they're not variable features). The tutorial wraps these plots in `suppressWarnings()` so you shouldn't see them, but they're safe either way.

**Gene names look like `SYMBOL_ENSG…` / `SYMBOL-ENSG…`** — `muscData`'s `Kang18_8vs8()` stores each feature as `SYMBOL_ENSEMBLID` (e.g. `ISG15_ENSG00000187608`). The loader strips the trailing Ensembl-ID suffix so only the gene symbol remains (`ISG15`), which is what canonical markers, Azimuth, and `org.Hs.eg.db` match on. It strips *only* the Ensembl suffix — not everything after the first dash — so real dashed symbols (`HLA-A`, `MT-CO1`) are preserved, and it reports how many names it simplified.
:::

::: {.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_01_QC_Preprocessing.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 {#setup}

::: {.callout-important title="Set up your project folder first — read this once"}
These tutorials assume a simple **project layout**. Make a project folder, put a **`scripts/`** subfolder inside it, and save every downloaded tutorial `.qmd` into `scripts/`:

```         
my_scrnaseq_project/
├── scripts/      ← put the downloaded Tutorial_*.qmd files here, and run them from here
├── data/         ← auto-created on first run: the .rds objects passed between tutorials
└── output/       ← auto-created on first run: each module's figures (Mod1/, Mod2/, …)
```

Run each tutorial **from `scripts/`** (open the `.qmd` in RStudio from there). The setup chunk below then creates `data/` and `output/` as **siblings of `scripts/`** — i.e. at the project root — using the paths `../data` and `../output/Mod1`. (The `..` is "one level up from `scripts/`".) Two practical rules:

1.  **`data/` is the hand-off.** Tutorial 01 writes `../data/ifnb_preprocessed.rds`; Tutorial 02 reads it back from the same `../data/`. Keep all the `.qmd` files together in `scripts/` so they share one project `data/`.
2.  **Can't find a file?** The setup chunk prints the exact absolute path (`This module writes its figures & tables to: …`). If it's not where you expect, check `getwd()` — it should be your `scripts/` folder. See the downloadable **README** (linked on the [Materials](../Materials.html) page) for the full layout.
:::

```{r}
#| label: M1-setup

# Packages are installed once in Tutorial 00 (Setup) — here we just load them.
library(Seurat)
library(muscData)
library(SingleCellExperiment)
library(tidyverse)
library(patchwork)   # combine + annotate multi-panel Seurat plots

set.seed(2026)

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

# ---------------------------------------------------------------------------
# Output directory for this module's figures and tables.
# Every figure/table chunk below writes a file named  Mod1_C<chunk>_<name>
# into ../output/Mod1/ so it can be cross-referenced from the rest of the site.
#   Mod1 = Module 1 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod1"
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))
```

::: callout-note
`muscData` (installed in the [Tutorial 00 setup](Tutorial_00_Setup_RStudio_Packages.html)) ships the data through Bioconductor's `ExperimentHub`: the first call to `Kang18_8vs8()` downloads \~25 MB into your local cache, and subsequent calls are instant.
:::

::: callout-tip
**Can't install `muscData`?** Download the preassembled <a href="../rds/ifnb_raw.rds" download>`ifnb_raw.rds`</a> instead and replace the load block below with `ifnb <- readRDS("../data/ifnb_raw.rds")`. It is the exact object Step 1 produces, so you pick up right at QC. Details on the [Datasets](../Datasets.html#core-data-load) page.
:::

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

## Step 1 — Load the `ifnb` dataset

```{r}
#| label: M1-load

# Fetch the Kang et al. 2017 PBMC dataset (CTRL + IFN-beta STIM) from
# Bioconductor's ExperimentHub. The first call downloads ~25 MB and caches it;
# every subsequent call is a no-op.
sce <- Kang18_8vs8()

# Keep only true singlets that have an author-assigned cell-type label (drop the
# multiplexed multiplets and unlabeled droplets). This leaves ~24,000 cells — the
# muscData build is larger than the historical SeuratData `ifnb` (~14k).
sce <- sce[, sce$multiplets == "singlet" & !is.na(sce$cell)]

# Gene (feature) name cleanup. muscData's Kang18_8vs8 stores each feature as
# SYMBOL_ENSEMBLID (e.g. "ISG15_ENSG00000187608"). We keep only the readable
# SYMBOL so canonical markers ("CD14", "MS4A1", ...) and reference tools (Azimuth,
# org.Hs.eg.db) match by name. Without this, Seurat would also rewrite the '_' to
# '-' and warn ("Feature names cannot have underscores ...").
#
# We strip ONLY a trailing "_ENSEMBLID" suffix, on the original underscore names.
# We deliberately do NOT cut at the first dash — real symbols contain dashes
# (HLA-A, MT-CO1, HLA-DRB1) and must be preserved. If a feature has no symbol
# (blank / pure Ensembl ID) we keep its original name so nothing is dropped, and
# make.unique() de-duplicates any symbols that collapse to the same name.
cm   <- counts(sce)
orig <- rownames(cm)
sym  <- sub("_ENS[A-Z]*[0-9]+(\\.[0-9]+)?$", "", orig)   # SYMBOL_ENSG... -> SYMBOL
blank <- is.na(sym) | sym == "" | grepl("^ENS[A-Z]*[0-9]+", sym)
sym[blank] <- orig[blank]
rownames(cm) <- make.unique(gsub("_", "-", sym))         # Seurat-safe, unique
n_changed <- sum(rownames(cm) != orig)
if (n_changed > 0) {
  i <- which(rownames(cm) != orig)[1]
  message(n_changed, " of ", length(orig), " feature names simplified to gene ",
          "symbols (e.g. '", orig[i], "' -> '", rownames(cm)[i], "')")
}

# Convert to a Seurat object, mirroring the metadata column names the rest of
# this tutorial series expects:
#   - `stim` as a factor with levels c("CTRL", "STIM")
#   - `seurat_annotations` for the author-curated cell-type labels
ifnb <- CreateSeuratObject(
  counts    = cm,
  meta.data = as.data.frame(colData(sce))
)
ifnb$stim <- factor(toupper(ifnb$stim), levels = c("CTRL", "STIM"))
ifnb$seurat_annotations <- factor(ifnb$cell)

ifnb

# What samples are in here?
table(ifnb$stim)

# Author-curated cell-type labels (we'll use these as ground truth in Tut 03)
table(ifnb$seurat_annotations)

# --- Tables out: per-sample and per-cell-type cell counts ------------------
tibble::enframe(table(ifnb$stim),
                name = "sample", value = "n_cells") |>
  dplyr::mutate(n_cells = as.integer(n_cells)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C2_cells_per_sample.csv"))

tibble::enframe(table(ifnb$seurat_annotations),
                name = "cell_type", value = "n_cells") |>
  dplyr::mutate(n_cells = as.integer(n_cells)) |>
  dplyr::arrange(dplyr::desc(n_cells)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C2_cells_per_celltype.csv"))
```

::: {.callout-important title="Think about it"}
1.  Why do real integration tutorials almost always have at least two samples?
2.  What's the difference between a "sample" (`stim`) and a "cell type" (`seurat_annotations`) in this dataset?

<details>

<summary>Show answers</summary>

1.  Integration aligns the **shared biological signal** between samples while preserving any **real differences** between them. With one sample there's nothing to align — clustering on PCs alone is enough. The `ifnb` dataset has CTRL and STIM as two samples drawn from the same donors, which is the perfect minimal setup to teach the alignment step.
2.  `stim` is the **technical / experimental covariate** — which condition the cell came from. `seurat_annotations` is the **biological identity** — what kind of cell it is. Integration aims to remove `stim`-driven structure from the embedding *without* erasing `seurat_annotations`-driven structure.

</details>
:::

## Step 2 — Inspect dimensions, split by sample

The `ifnb` object already arrives as a single Seurat object containing both conditions. For a two-sample QC workflow, it's useful to look at QC metrics **per sample** rather than globally — that way you'll catch a sample-specific quality issue (e.g. one sample with much lower depth) before it confounds the analysis.

```{r}
#| label: M1-split

# Record the starting matrix dimensions
dim_raw <- dim(ifnb)
dim_raw

# --- Table out: starting matrix dimensions (features x barcodes) -----------
tibble::tibble(
  stage    = "Step 1 — loaded (singlets, annotated)",
  features = dim_raw[1],
  barcodes = dim_raw[2]
) |>
  readr::write_csv(file.path(out_dir, "Mod1_C3_starting_dimensions.csv"))
```

::: {.callout-important title="Think about it"}
1.  The `ifnb` object holds both CTRL and STIM cells in a single counts matrix. What information distinguishes a CTRL barcode from a STIM barcode?
2.  Why might it matter to split visualizations by `stim` even though clustering happens on the joint matrix?

<details>

<summary>Show answers</summary>

1.  The `stim` column in `meta.data`. Each barcode also carries the `-1` / `-2` GEM-well suffix from Cell Ranger, but you should never rely on barcode suffixes for sample membership — the metadata is the source of truth.
2.  Per-sample plots reveal **batch-level QC issues** (lower median depth in one sample, different mitochondrial fractions, different cell counts per cluster) that can be hidden by a global histogram. If one sample's distribution is shifted, you may need different filter thresholds for it — or you may need to investigate before going further.

</details>
:::

## Step 3 — Add QC metrics

```{r}
#| label: M1-qc_metrics

ifnb[["percent.mt"]] <- PercentageFeatureSet(ifnb, pattern = "^MT-")
ifnb[["percent.rb"]] <- PercentageFeatureSet(ifnb, pattern = "^RP[SL]")

head(ifnb@meta.data)

# --- Table out: per-sample summary of the QC metrics just computed ---------
ifnb@meta.data |>
  dplyr::group_by(stim) |>
  dplyr::summarise(
    n_cells           = dplyr::n(),
    median_nFeature   = median(nFeature_RNA),
    median_nCount     = median(nCount_RNA),
    median_percent_mt = round(median(percent.mt), 2),
    median_percent_rb = round(median(percent.rb), 2),
    .groups = "drop"
  ) |>
  readr::write_csv(file.path(out_dir, "Mod1_C4_qc_metric_summary.csv"))
```

::: callout-note
The `ifnb` raw matrix uses HUGO gene symbols, so `^MT-` matches the 13 protein-coding mitochondrial genes (MT-ND1, MT-CO1, …) and `^RP[SL]` matches the ribosomal protein genes (RPS\*, RPL\*). For mouse data the conventions are lowercase: `^mt-` and `^Rp[sl]`.
:::

::: {.callout-important title="Think about it"}
1.  What biological signal does `percent.mt` capture?
2.  Why might `percent.rb` (ribosomal-protein fraction) be useful?

<details>

<summary>Show answers</summary>

1.  Dying or stressed cells leak cytoplasmic RNA while retaining mitochondrial RNA, so their mitochondrial fraction rises. High `percent.mt` flags apoptotic / compromised cells.
2.  Ribosomal-protein genes are **highly expressed** and **broadly variable** across cell states. A very low or very high `percent.rb` per cluster is sometimes a tell-tale sign of a metabolic / activation difference rather than a cell-type difference — useful to separate "this cluster is biologically distinct" from "this cluster is just translationally quiet."

</details>
:::

## Step 4 — Visualize QC, per sample

```{r}
#| label: M1-qc_plots

# Violin panel: four QC metrics, split by sample, fully labelled.
# These features are per-cell QC METADATA plotted before normalization, so Seurat
# warns that the "data" layer is empty and falls back to "counts". That's
# irrelevant here (metadata isn't read from an expression layer), so we silence
# the expected note.
p_qc_vln <- suppressWarnings(
  VlnPlot(ifnb,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
        group.by = "stim",
        ncol = 4,
        pt.size = 0) &
  xlab("Sample (stim condition)") &
  theme(plot.title = element_text(size = 11)))
p_qc_vln <- p_qc_vln +
  patchwork::plot_annotation(
    title    = "Per-cell QC metrics by sample — ifnb (CTRL vs IFN-β STIM)",
    subtitle = "Genes/cell, UMIs/cell, mitochondrial % and ribosomal % distributions",
    caption  = "Module 1 · QC & Preprocessing"
  )
p_qc_vln

# Scatter: UMIs vs unique genes, coloured by sample, with a linear trend
p_qc_scatter <- suppressWarnings(FeatureScatter(ifnb,
               feature1 = "nCount_RNA",
               feature2 = "nFeature_RNA",
               group.by = "stim")) +
  geom_smooth(method = "lm") +
  labs(
    title    = "Sequencing depth vs library complexity",
    subtitle = "Each point is a cell; the line is a per-group linear fit",
    x        = "nCount_RNA (total UMIs per cell)",
    y        = "nFeature_RNA (unique genes per cell)",
    colour   = "Sample"
  )
p_qc_scatter

# --- Figures out -----------------------------------------------------------
save_fig(file.path(out_dir, "Mod1_C5_qc_violins.png"),  p_qc_vln,
       width = 12, height = 4, dpi = 300)
save_fig(file.path(out_dir, "Mod1_C5_qc_scatter.png"),  p_qc_scatter,
       width = 6,  height = 5, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
**Violin panel** — one violin per sample (x-axis = `CTRL` vs `STIM`) for each of the four QC metrics; the width of a violin at a given height is how many cells sit at that value, so you're reading the *shape* of each per-cell distribution. Healthy cells form the fat central body; the thin upper tail on `nFeature_RNA`/`nCount_RNA` is where doublets hide, and the upper tail on `percent.mt` is where dying cells sit. **Scatter** — each point is one cell: x = total UMIs (`nCount_RNA`), y = unique genes (`nFeature_RNA`). The two scale together along a tight curve; points far *below* the trend are low-complexity/dying cells, points far *above* at the top-right are doublet-like. The per-sample fit lines let you spot a sample with systematically different depth. (`scater::isOutlier()` is a data-driven alternative to eyeballing these thresholds.)
:::

::: {.callout-important title="Think about it"}
1.  Compare the violins between CTRL and STIM. Are the medians similar? Does either sample have a much fatter upper tail on `nFeature_RNA`?
2.  What would a doublet look like on the `nCount_RNA` vs `nFeature_RNA` scatter? A dying cell?

<details>

<summary>Show answers</summary>

1.  In `ifnb` the two samples are very similar in median QC metrics — that's a feature, not a bug, because the donors are matched. A noticeably fatter upper tail in one sample would suggest more **doublets** there (perhaps a higher loading concentration). A noticeably *lower* median `nFeature_RNA` would suggest sequencing depth or cell-quality differences.
2.  Doublets sit **high on both axes**, often above the main trend line. Dying cells are **low** on both and tend to have **high `percent.mt`** (color the scatter by `percent.mt` or facet by it).

</details>
:::

## Step 5 — Optional: ambient-RNA correction with SoupX

In a 10x droplet experiment, lysed-cell debris in the bulk solution leaks **ambient mRNA** into every droplet. After Cell Ranger, this shows up as a low background of "wrong" transcripts in every cell — most visibly when a cell expresses a marker for a *different* cell type at low level. **`SoupX`** estimates the ambient profile and corrects each cell's counts. It runs *before* QC filtering and replaces the counts matrix with a cleaned one.

::: callout-warning
**SoupX needs both the raw and filtered Cell Ranger matrices** (`raw_feature_bc_matrix/` and `filtered_feature_bc_matrix/`) because it learns the soup profile from empty droplets. `Kang18_8vs8()` only ships the **filtered** counts, so you cannot run SoupX on `ifnb` directly. The code below is the canonical pattern — use it on **your own** 10x data, or in the optional [Full Talapas run from raw FASTQs](../Resources_Folder/Talapas_FullRun_FromFASTQs.html), where the raw matrix is available.
:::

*Optional reference pattern — left `#| eval: false` because it needs a raw Cell Ranger matrix (not shipped with `ifnb`). You don't need to change anything here; use it on your own 10x data.*

```{r}
#| label: M1-soupx
#| eval: false

# Pattern for any 10x dataset where you have the raw Cell Ranger output.
# This block does NOT run on ifnb (which ships only *filtered* counts) and will
# NOT run as-is: edit CR_OUTS to point at a real Cell Ranger 'outs/' directory
# that contains BOTH raw_feature_bc_matrix/ and filtered_feature_bc_matrix/.
library(SoupX)
library(DropletUtils)

# 1. Point at the Cell Ranger 'outs/' directory (containing raw_feature_bc_matrix/
#    and filtered_feature_bc_matrix/)
CR_OUTS <- "path/to/cellranger/outs"   # <-- EDIT to your own Cell Ranger outs/

# Fail early with a clear message rather than a cryptic "subscript out of bounds"
# if the path is still the placeholder or is missing the two required matrices.
if (!dir.exists(file.path(CR_OUTS, "raw_feature_bc_matrix")) ||
    !dir.exists(file.path(CR_OUTS, "filtered_feature_bc_matrix"))) {
  stop("SoupX needs a real Cell Ranger 'outs/' folder containing BOTH ",
       "raw_feature_bc_matrix/ and filtered_feature_bc_matrix/.\n",
       "  '", CR_OUTS, "' is not valid — ifnb ships only the filtered counts, so ",
       "this block is a template for your OWN 10x data, not a step you run on ifnb.")
}

sc <- load10X(CR_OUTS)

# 2. (Recommended) provide preliminary clusters so SoupX can refine the soup
#    profile per cluster — compute them with a quick Seurat workflow first:
# sc <- setClusters(sc, setNames(seu$seurat_clusters, colnames(seu)))
# sc <- setDR(sc, Embeddings(seu, "umap"))

# 3. Estimate the ambient contamination fraction
sc  <- autoEstCont(sc)
rho <- sc$fit$rhoEst
message("Estimated ambient fraction: ", round(rho, 3))

# 4. Adjust the counts (round to integers so downstream UMI tools are happy)
adj <- adjustCounts(sc, roundToInt = TRUE)

# 5. Use `adj` as the input to CreateSeuratObject() in Step 1 instead of the raw
#    filtered matrix; everything else in this tutorial then runs on the cleaned counts.
```

::: {.callout-important title="Think about it"}
1.  Why does SoupX need the **raw** (pre-filter) matrix? What information does it pull from there?
2.  When is SoupX overkill? When is it essential?

<details>

<summary>Show answers</summary>

1.  The raw matrix contains droplets that captured **no real cell** — they hold only ambient RNA. SoupX uses those empty droplets to estimate the soup's expression profile, then subtracts a fraction of that profile from every called cell.
2.  **Overkill:** small contamination (rho ≲ 0.05), well-prepared samples, no marker bleed-through. **Essential:** complex tissues with high cell death / large cell-type differences (tumor microenvironment, mucosal tissues), or any time you see a marker like `HBB` or `IGKC` showing up faintly in cells that shouldn't express it.

</details>
:::

## Step 6 — Doublet detection with `scDblFinder`

A **doublet** is one barcode's worth of counts that came from two cells getting trapped in the same droplet. They show up as cells with anomalously high `nCount_RNA` and `nFeature_RNA` and a mixed expression profile (e.g. T-cell markers AND B-cell markers in the same cell). The simple `nFeature_RNA < 2500` filter applied later ([code](#lst-M1-filter)) catches some, but doublets that pair two *different* cell types of similar size will pass that filter — you need a model.

`scDblFinder` simulates artificial doublets, trains a classifier to distinguish them from singletons, and returns a per-cell score and class. Critically, **doublets form within a sample, not across samples**, so we run it per sample.

```{r}
#| label: M1-dbl

library(scDblFinder)

set.seed(2026)

# Run scDblFinder per sample, then merge the results back into the Seurat object.
# NOTE on the messages you'll see here (all harmless, NOT errors):
#   - "Layer 'data'/'scale.data' is empty" — doublet detection runs BEFORE
#     normalization (Step 8), so those layers aren't populated yet; scDblFinder
#     works on raw counts and normalizes internally. We wrap the conversion in
#     suppressWarnings() to quiet those expected notes.
#   - "'normalizeCounts'/'librarySizeFactors' is deprecated" — these come from
#     scran/scuttle *inside* scDblFinder, not from this code; safe to ignore.
ifnb_list <- SplitObject(ifnb, split.by = "stim")
ifnb_list <- lapply(ifnb_list, function(x) {
  sce <- suppressWarnings(as.SingleCellExperiment(x))
  sce <- scDblFinder(sce, returnType = "sce")
  x$scDblFinder.score <- colData(sce)$scDblFinder.score
  x$scDblFinder.class <- colData(sce)$scDblFinder.class    # "singlet" or "doublet"
  x
})
ifnb <- merge(ifnb_list[[1]], y = ifnb_list[-1])

# `merge()` in Seurat v5 leaves the RNA assay with per-sample split layers
# (counts.CTRL / counts.STIM, etc.). Downstream tools like FindMarkers expect
# a single joined layer, so rejoin them now.
ifnb <- JoinLayers(ifnb)

# How many doublets were called per sample?
table(sample = ifnb$stim, class = ifnb$scDblFinder.class)

# --- Table out: doublet calls per sample -----------------------------------
as.data.frame(table(sample = ifnb$stim, class = ifnb$scDblFinder.class)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C7_doublet_counts.csv"))

# Inspect: doublets should sit at the high end of nCount/nFeature.
# (Same pre-normalization metadata plot as Step 4 — silence the expected
# "data layer empty, using counts" note.)
p_dbl <- suppressWarnings(
  VlnPlot(ifnb,
        features = c("nFeature_RNA", "nCount_RNA"),
        group.by = "scDblFinder.class", pt.size = 0) &
  xlab("scDblFinder class") &
  theme(plot.title = element_text(size = 11)))
p_dbl <- p_dbl +
  patchwork::plot_annotation(
    title    = "Library size and complexity for called singlets vs doublets",
    subtitle = "Doublets are expected to sit at the high end of both metrics",
    caption  = "Module 1 · QC & Preprocessing"
  )
p_dbl

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod1_C7_doublet_violins.png"), p_dbl,
       width = 7, height = 4, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The `table()` is a 2×2 of sample × class: the key number is the **doublet fraction** per sample (`doublet / (singlet + doublet)`), which should land around **5–10%** for a normal 10x run — much higher points to chip overloading. In the violins (x-axis = `singlet` vs `doublet`), the `doublet` group should sit visibly **higher** on both `nFeature_RNA` and `nCount_RNA`; that separation is the sanity check that the classifier keyed on the expected signal. `scDblFinder.score` (0–1) is the per-cell confidence if you want to threshold more or less aggressively than the default class call.
:::

::: {.callout-important title="Think about it"}
1.  Why run `scDblFinder` **per sample** instead of on the merged object?
2.  scDblFinder marks \~5–10% of cells as doublets in a typical 10x run. If your run shows 25%, what likely happened?
3.  Some pipelines run doublet detection *after* clustering and remove doublet-enriched clusters wholesale. Why might that be more conservative than removing per-cell doublet calls?

<details>

<summary>Click for answer</summary>

1.  Doublets are made when **two cells from the same droplet load** are co-encapsulated. Cells from different samples (CTRL and STIM here) were prepared in separate channels, so a STIM and a CTRL cell can never be a doublet. Running per-sample also lets each sample's doublet-rate be modeled correctly.
2.  Either you over-loaded the chip (the loading concentration was too high — typical for "nominal 10k cells, recovered 18k") or there's a quality problem (debris, large cell aggregates) making the simulated-doublet classifier confused. Inspect `nCount_RNA`/`nFeature_RNA` distributions and the ambient fraction.
3.  Per-cell doublet calls have false positives. If you trust them and remove individual cells, you may *also* be removing legitimate cycling cells, large activated cells, or megakaryocytes (which look doublet-like by sheer transcriptional output). Removing doublet-enriched **clusters** is more conservative — only act on a cluster where the doublet fraction is clearly elevated and the markers don't make biological sense.

</details>
:::

## Step 7 — Filter low-quality cells

```{r}
#| label: M1-filter

ifnb <- subset(
  ifnb,
  subset = nFeature_RNA > 200 &
           nFeature_RNA < 2500 &
           percent.mt   < 5 &
           scDblFinder.class == "singlet"
)
ifnb
table(ifnb$stim)

dim_after_filter <- dim(ifnb)

# --- Table out: cells surviving the QC filter, per sample ------------------
tibble::enframe(table(ifnb$stim),
                name = "sample", value = "n_cells_after_filter") |>
  dplyr::mutate(n_cells_after_filter = as.integer(n_cells_after_filter)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C8_cells_after_filter.csv"))
```

::: {.callout-important title="Think about it"}
1.  Did the CTRL or STIM sample lose more cells in this filter? If they lost noticeably different fractions, what would that suggest?
2.  How would you set thresholds *empirically* rather than by default?

<details>

<summary>Show answers</summary>

1.  If one sample loses a much larger fraction at `percent.mt < 5` it usually means its cells were under more stress at dissociation or stimulation. For `ifnb`, IFN-β stimulation can mildly elevate `percent.mt` in some donors — flag the sample, but don't necessarily relax the threshold without thinking.
2.  Look at the per-sample violin / histogram shape; choose thresholds that trim the clear **tails** (not the mode). For an automated rule, use `scater::isOutlier()` or `median ± k·MAD` — and consider running it **per sample** so a globally lopsided distribution doesn't push the cutoff in the wrong direction for one sample.

</details>
:::

## Step 8 — Normalize

```{r}
#| label: M1-normalize

ifnb <- NormalizeData(ifnb)
# Equivalent explicit form:
# NormalizeData(ifnb,
#               normalization.method = "LogNormalize",
#               scale.factor = 10000)

dim_after_normalize <- dim(ifnb)
```

::: {.callout-important title="Think about it"}
1.  What does `LogNormalize` do algebraically?
2.  When would you prefer `SCTransform()` instead?

<details>

<summary>Show answers</summary>

1.  For each cell: divide each gene's count by the cell's total count, multiply by `scale.factor` (default 10 000), then take `log1p`. In formula form: `log(1 + (count_gc / total_c) * 10000)`.
2.  `SCTransform` models the mean–variance relationship per gene with a regularized negative binomial, which is usually better for **low-count** datasets and more faithful downstream variance. It also replaces the NormalizeData → FindVariableFeatures → ScaleData chain in a single call. For `ifnb` either approach works fine.

</details>
:::

## Step 9 — Highly variable features

```{r}
#| label: M1-hvgs

ifnb <- FindVariableFeatures(
  ifnb,
  selection.method = "vst",
  nfeatures = 2000
)

top10 <- head(VariableFeatures(ifnb), 10)
plot1 <- VariableFeaturePlot(ifnb)
p_hvg <- LabelPoints(plot = plot1, points = top10, repel = TRUE) +
  labs(
    title    = "Highly variable genes (vst) — top 2,000 selected",
    subtitle = "Top 10 most variable genes labelled; many are interferon-stimulated genes",
    x        = "Average expression (log scale)",
    y        = "Standardized variance",
    colour   = "Selected as HVG"
  )
# VariableFeaturePlot draws mean expression on a log10 x-axis. Genes with zero
# mean (undetected across all retained cells) become -Inf and ggplot drops them
# with "log-10 transformation introduced infinite values" — harmless, since those
# genes are not variable features. Silence it at draw and save time.
suppressWarnings(print(p_hvg))

# --- Outputs: HVG plot + the labelled top-10 genes as a table --------------
suppressWarnings(
  save_fig(file.path(out_dir, "Mod1_C10_variable_features.png"), p_hvg,
         width = 7, height = 5, dpi = 300))

tibble::tibble(rank = seq_along(top10), gene = top10) |>
  readr::write_csv(file.path(out_dir, "Mod1_C10_top10_hvgs.csv"))
```

::: {.callout-tip title="Reading the output"}
Each dot is a gene: x-axis = **average expression** (log scale), y-axis = **standardized variance** (how much more variable the gene is than expected for its mean). The \~2,000 genes coloured/highlighted as variable are the ones that rise above the fitted mean–variance curve — these are the only genes PCA will use in Tutorial 02. The labelled top-10 are the most variable of all; in `ifnb` they're dominated by **interferon-stimulated genes** (`ISG15`, `IFI6`, `IFIT1`…) because the STIM half of the data switches them on and the CTRL half doesn't. Seeing ISGs at the top here is your first preview of the batch-style split you'll visualize in Tutorial 02 and correct in Tutorial 05.
:::

::: {.callout-important title="Think about it"}
1.  Which genes are at the top of `top10`? Many will be **interferon-stimulated genes** (ISGs) like `ISG15`, `IFI6`, `IFIT1` — why?
2.  What happens if you set `nfeatures = 500`? Or `nfeatures = 10000`?

<details>

<summary>Show answers</summary>

1.  Because half the dataset (the STIM sample) has IFN-β–induced ISG expression while the other half (CTRL) does not, ISGs have **enormous variance across cells** and rise to the top of the HVG list. This is exactly the signal we'll see drive batch-style separation in Tutorial 02 — and it's exactly what integration in Tutorial 05 will need to handle without erasing.
2.  Too few (500) risks missing subtle cell-type signals and is brittle to noisy gene choices. Too many (10 000) pulls in low-information genes and dilutes the signal. 1 500–3 000 is a common sweet spot; 2 000 is a good default.

</details>
:::

## Step 10 — Scale

```{r}
#| label: M1-scale

ifnb <- ScaleData(ifnb)
```

By default `ScaleData` only scales the variable features identified in Step 9 (\~2,000 HVGs) — that's all PCA needs, and it keeps the `scale.data` slot (and the saved `.rds`) small. If you later need scaled values for a specific non-HVG gene (e.g. to draw it on a heatmap), re-run `ScaleData(ifnb, features = <those genes>)` at that point.

::: {.callout-important title="Think about it"}
1.  What does `ScaleData` do to each gene numerically?
2.  Why is scaling necessary *before* PCA?
3.  Could you `vars.to.regress = "stim"` here and skip integration in Tutorial 05?

<details>

<summary>Show answers</summary>

1.  Centers each gene's expression to mean 0 and scales to unit variance across cells (z-scoring).
2.  PCA finds directions of maximum variance. Without scaling, highly-expressed genes (which inherently have large raw variance) would dominate PCs, regardless of whether they are biologically informative.
3.  You *could*, but regressing on a binary covariate is a **brittle** form of batch correction — it removes a linear effect uniformly across all cells, but real batch effects are usually non-linear and cell-type-specific. Modern integration methods (Harmony, anchors, scVI) are much better tools. We use `ScaleData` here as the standard pre-PCA step and reserve the integration question for Tutorial 05.

</details>
:::

## Counts matrix progression — summary

```{r}
#| label: M1-dim_summary

dim_summary <- tibble::tribble(
  ~step,                                                    ~features_after,        ~barcodes_after,
  "Step 1 — Kang18_8vs8 (singlets, raw)",                    dim_raw[1],             dim_raw[2],
  "Step 7 — subset() QC filter (nFeature_RNA, percent.mt)",  dim_after_filter[1],    dim_after_filter[2],
  "Step 8 — NormalizeData (no dimension change)",            dim_after_normalize[1], dim_after_normalize[2]
)
dim_summary

# --- Table out: counts-matrix progression summary --------------------------
readr::write_csv(dim_summary, file.path(out_dir, "Mod1_C12_dimension_progression.csv"))
```

You should see roughly:

| Step | Features after | Barcodes after |
|----|---:|---:|
| Step 1 — `Kang18_8vs8` (singlets, annotated) | \~35,635 | \~24,500 |
| Step 7 — `subset()` QC filter (`nFeature_RNA`, `percent.mt`) | \~35,635 | \~23,000 |
| Step 8 — `NormalizeData` (no dimension change) | \~35,635 | \~23,000 |

(Numbers will vary by a few hundred depending on package versions.)

## Save and continue to Tutorial 02

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

# Make sure data/ exists
dir.create("../data", showWarnings = FALSE)

saveRDS(ifnb, file = "../data/ifnb_preprocessed.rds")
```

::: callout-tip
The saved `.rds` is the starting point for **Tutorial 02 — Dimensionality reduction & clustering**. We deliberately keep both `CTRL` and `STIM` cells in the same object — Tutorial 02 will run PCA + clustering on the unintegrated joint data so you can **see the batch effect** before fixing it in Tutorial 05.
:::

## Credits

- Dataset: Kang *et al.* (2017) *Nat Biotechnol* 36: 89–94, [doi:10.1038/nbt.4042](https://doi.org/10.1038/nbt.4042) — fetched via Bioconductor's [`muscData::Kang18_8vs8()`](https://bioconductor.org/packages/muscData/)
- Reference vignettes: [Seurat integration tutorial](https://satijalab.org/seurat/articles/integration_introduction.html) · [Seurat PBMC 3k tutorial](https://satijalab.org/seurat/articles/pbmc3k_tutorial.html)
- Inspiration for the multi-step QC pacing: [scNotebooks](https://integrativebioinformatics.github.io/scNotebooks/) (Sengupta lab)