---
title: "Tutorial 15 — scATAC-seq with Signac"
subtitle: "From peaks matrix to clustered cells (PBMC 10k v1, 10x)"
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 exercise

A parallel to the scRNA-seq tutorial series — but for **chromatin accessibility** (scATAC-seq). The workflow follows the standard **Signac** pipeline from the Stuart Lab on the 10x PBMC 10k scATAC-seq dataset.

-   Reference vignette: [Signac PBMC scATAC-seq vignette](https://stuartlab.org/signac/articles/pbmc_vignette)
-   Companion lecture: [Lecture 15 — Bulk + scATAC-seq](../Lecture_Folder/Lecture_15_ATACseq.html)
-   Alternative toolkit walkthrough: [scNotebooks Module 11](https://integrativebioinformatics.github.io/scNotebooks/modules/en/module11/module11.html) does the same kind of analysis with **ArchR** (peak calling, motif enrichment, TF regulators, scATAC↔scRNA integration, trajectory) on Kumegawa *et al.* 2022 breast-cancer data — a useful cross-toolkit reference

::: {.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_15_scATACseq_Signac.qmd.txt" download="Tutorial_15_scATACseq_Signac.qmd">Tutorial_15_scATACseq_Signac.qmd</a>. If your browser saves the file as `Tutorial_15_scATACseq_Signac.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  Download the data files (block below) into `../data/`.
3.  Work through the chunks, flipping `eval: true` as you go.
:::

::: {.callout-tip title="What's different about ATAC-seq?"}
ATAC-seq measures **open chromatin** (regions accessible to Tn5), not transcription. The raw unit is a **fragment** (a pair of Tn5 cut sites). The "features" are **peaks** — genomic intervals of accessibility — rather than genes. Most of the Seurat idioms carry over, but the **normalization**, **dimensionality reduction**, and **QC metrics** are ATAC-specific.
:::



::: callout-note
**Companion book chapter:** [Chapter 15 — scATAC-seq](../Resources_Folder/Chapter_15_ATAC.html) — the long-form prose treatment of this tutorial's material, with cross-references to the prerequisite appendices.
:::


## Dataset — PBMC 10k scATAC-seq v1 (10x Genomics, hg19)

**What it is.** Approximately ten-thousand human **peripheral blood mononuclear cells (PBMCs)** from a single healthy donor, run on the 10x Chromium **Single Cell ATAC v1** assay. This is the canonical tutorial dataset used across Signac, ArchR, and EpiScanpy vignettes.

-   Landing page: [10x Genomics — PBMC 10k scATAC-seq v1 (hg19)](https://www.10xgenomics.com/datasets/10-k-peripheral-blood-mononuclear-cells-pbm-cs-from-a-healthy-donor-1-standard-1-1-0)
-   Reference genome: **hg19 / GRCh37** (hence `EnsDb.Hsapiens.v75` below)
-   Chemistry: 10x Chromium Single Cell ATAC v1
-   Cell count: \~8,700 post-filter cells

::: callout-warning
Only the **v1** release is pinned to hg19. The *Next GEM v1.1* re-release of the same sample is mapped to **hg38** and uses file names prefixed with `atac_v1_pbmc_10k_nextgem_*`. This exercise assumes the original hg19 release.
:::

### Download the data

The four v1 (hg19) files load directly into Signac — there is no R-package download for them, so fetch them once from the 10x CDN into your project's `data/` folder (`../data/` when you run the commands from `scripts/`):

```{bash}
#| label: M15-mkdir_p_data_cd
#| eval: false
mkdir -p ../data && cd ../data

wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbi
```

::: callout-tip
The optional `pre-processed scRNA-seq` companion (used for label transfer in the Signac vignette) lives at <https://signac-objects.s3.amazonaws.com/pbmc_10k_v3.rds> — pull it the same way if you want to reproduce the RNA-anchored annotation step.
:::

## Learning goals

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

-   Load peak-count and fragment data into a Signac `ChromatinAssay` inside a Seurat object
-   Compute and interpret ATAC-specific QC metrics (nucleosome signal, TSS enrichment, blacklist ratio, peak-read fraction)
-   Run TF-IDF normalization and SVD (LSI) on accessibility data
-   Produce a UMAP and cluster scATAC-seq cells


::: {.callout-warning title="Common errors / things that bite"}
**`EnsDb.Hsapiens.v75` (hg19) doesn't match the data** — the workshop's PBMC 10k v1 dataset is on hg19. The newer "Next GEM v1.1" release of the same sample is on hg38 and has filenames prefixed `atac_v1_pbmc_10k_nextgem_*`. Mixing them produces bizarre coordinate errors. Use `v75` (hg19) with the v1 hg19 dataset, or `v86` (hg38) with the Next GEM hg38 dataset.

**`fragments.tsv.gz` is missing its `.tbi` index** — Signac uses `Rsamtools` to query the fragments file by genomic region; without `.tbi` it can't seek. Both files must be present and the `.tbi` must be newer than the `.tsv.gz`. Run `tabix fragments.tsv.gz` to regenerate.

**`CreateChromatinAssay` warns about "duplicate fragments"** — usually safe to ignore; happens when the same fragment is recorded under different cell barcodes. If the warning rate is high (> 5%), suspect a sample-mixing issue upstream.
:::

::: {.callout-tip title="When you need to liftOver — hg19 ↔ hg38"}
If you ever need to bring a third-party peak set or annotation onto the build you're working on, do a coordinate **liftOver** with `rtracklayer` and a UCSC chain file. Don't try to remap by gene name — coordinate offsets between hg19 and hg38 are nontrivial.

```{r}
#| label: M15-load_hg19_bed_style
#| eval: false
library(rtracklayer)

# Load hg19 BED-style coordinates (chr, start, end)
coords <- read.delim("../data/peaks_hg19.bed", header = FALSE)
names(coords) <- c("chromosome", "start", "end")
gr <- GRanges(coords)

# UCSC chain file: download once from
# https://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz
chain   <- import.chain("chain_files/hg19ToHg38.over.chain")
results <- as.data.frame(liftOver(gr, chain))
```

A few coordinates may map to multiple intervals (or none) — `liftOver()` returns a `GRangesList`, so always check `lengths(results) == 1` before downstream filtering. The reverse direction uses `hg38ToHg19.over.chain`. Mouse equivalents: `mm9ToMm10`, `mm10ToMm39`.
:::


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

## Setup

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

# Packages are installed in Tutorial 00 (Setup → bonus modules) — load them here.
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v75)
library(tidyverse)
library(patchwork)   # combine + annotate multi-panel Signac/Seurat plots

set.seed(2026)

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

## Data files and data dictionary

You will need three files from the 10x PBMC 10k scATAC-seq release, placed under `../data/`. All three are available under "Output and supplemental files" on the [10x landing page](https://www.10xgenomics.com/datasets/10-k-peripheral-blood-mononuclear-cells-pbm-cs-from-a-healthy-donor-next-gem-v-1-1-1-1-standard-2-0-0).

| File | Size | Purpose |
|---|---|---|
| `atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5` | \~30 MB | peak × barcode counts (sparse, filtered to called cells) |
| `atac_v1_pbmc_10k_fragments.tsv.gz` + `.tbi` | \~1.8 GB + \~1 MB | per-fragment record (tabix-indexed) |
| `atac_v1_pbmc_10k_singlecell.csv` | \~20 MB | per-barcode Cell Ranger ATAC summary |

### `filtered_peak_bc_matrix.h5`

| Dimension | Type | Description |
|---|---|---|
| rows | character | **Peak coordinates** as `chrom:start-end` strings (e.g. `chr1:10279-10779`). Peaks are fixed-width genomic intervals called by Cell Ranger. |
| cols | character | **Cell barcodes** (16 bp + `-1` GEM-well suffix). |
| values | integer | Number of **Tn5 fragments** from that cell whose midpoint falls inside that peak. |

### `fragments.tsv.gz` (bgzipped, tabix-indexed)

Tab-separated, one row per unique fragment. Five columns:

| # | Column | Description |
|---|---|---|
| 1 | `chrom` | Chromosome (`chr1` … `chrY`, plus contigs). |
| 2 | `start` | 0-based start of the fragment (= position of one Tn5 cut). |
| 3 | `end` | 0-based exclusive end of the fragment (= position of the other Tn5 cut). |
| 4 | `cell_barcode` | 10x cell barcode (matches column names of the peak matrix). |
| 5 | `duplicate_count` | Number of PCR-duplicate reads collapsed into this fragment (> 1 means the same fragment was sequenced multiple times). |

### `singlecell.csv` (per-barcode Cell Ranger summary)

One row per barcode. Relevant columns (not exhaustive — 10x produces \~20 columns):

| Column | Description |
|---|---|
| `barcode` | 10x barcode string (row names in `read.csv(..., row.names = 1)`). |
| `is__cell_barcode` | Boolean flag, Cell Ranger's call of "this barcode is a real cell". |
| `total` | Total read pairs assigned to this barcode. |
| `passed_filters` | Reads passing all Cell Ranger ATAC QC filters. |
| `peak_region_fragments` | Fragments falling within a called peak. Numerator for `pct_reads_in_peaks`. |
| `blacklist_region_fragments` | Fragments falling inside the ENCODE blacklist. Numerator for `blacklist_ratio`. |
| `TSS_fragments` | Fragments within ± flank of a TSS. |
| `mitochondrial` | Fragments mapping to chrM. |
| `duplicate` | PCR-duplicate fragments. |

These columns feed directly into the ATAC QC metrics computed in Step 4.

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

## Step 1 — Peek at the fragment file

```{r}
#| label: M15-peek_fragments

frag.file <- read.delim(
  "../data/atac_v1_pbmc_10k_fragments.tsv.gz",
  header = FALSE,
  nrows  = 10
)
head(frag.file)

# --- Table out: first rows of the raw fragment file ------------------------
frag_preview <- frag.file
colnames(frag_preview) <- c("chrom", "start", "end",
                            "cell_barcode", "duplicate_count")
readr::write_csv(tibble::as_tibble(frag_preview),
                 file.path(out_dir, "Mod15_C2_fragment_preview.csv"))
```

::: {.callout-important title="Think about it"}
1.  A fragment file has one line per Tn5 fragment. What are the five columns typically recorded?
2.  Why do we keep the fragment file around instead of just using the peak-counts matrix?

<details>

<summary>Show answers</summary>

1.  `chrom`, `start`, `end`, `cell_barcode`, `duplicate_count` (number of reads collapsed to that fragment).
2.  Many ATAC analyses (TSS enrichment, nucleosome signal, coverage tracks, motif activity) work at the fragment level, not the peak level. The fragment file is the raw, peak-independent record; peaks can be recomputed later.

</details>
:::

## Step 2 — Load the peaks matrix and create the object

```{r}
#| label: M15-load_matrix

counts <- Read10X_h5("../data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
counts[1:5, 1:5]

chrom_assay <- CreateChromatinAssay(
  counts       = counts,
  sep          = c(":", "-"),
  fragments    = "../data/atac_v1_pbmc_10k_fragments.tsv.gz",
  min.cells    = 10,
  min.features = 200
)

metadata <- read.csv(
  file      = "../data/atac_v1_pbmc_10k_singlecell.csv",
  header    = TRUE,
  row.names = 1
)

pbmc <- CreateSeuratObject(
  counts    = chrom_assay,
  meta.data = metadata,
  assay     = "ATAC"
)

pbmc
```

::: {.callout-important title="Think about it"}
1.  The row names of an RNA assay are gene symbols. What do row names of a `ChromatinAssay` look like, and why?
2.  We pass both `counts` and `fragments` when building the ChromatinAssay. Why is that redundant-looking pairing actually useful?

<details>

<summary>Show answers</summary>

1.  Peak coordinates — strings like `chr1:100234-100734`. Peaks are genomic intervals, not annotated features, so the natural identifier is the interval itself.
2.  `counts` is the fast aggregated view (peak × cell). `fragments` retains per-fragment, per-base precision for downstream peak recomputation, coverage plots, TSS enrichment, motif accessibility, etc. You want both.

</details>
:::

## Step 3 — Add gene annotation

```{r}
#| label: M15-annotation

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
seqlevels(annotations) <- paste0("chr", seqlevels(annotations))  # UCSC naming
Annotation(pbmc) <- annotations
```

::: {.callout-important title="Think about it"}
1.  The data is mapped to **hg19** (GRCh37). Why `EnsDb.Hsapiens.v75` specifically?
2.  What would you change for a dataset mapped to hg38 or to mm10?

<details>

<summary>Show answers</summary>

1.  Ensembl v75 is the last annotation release pinned to GRCh37 / hg19. Using a matching annotation is essential — otherwise TSSs and gene bodies will be in the wrong places.
2.  hg38: use a current human EnsDb (e.g. `EnsDb.Hsapiens.v86` or newer). mm10: use `EnsDb.Mmusculus.v79` or equivalent. Chromosome-name style also differs (`1` vs `chr1`); convert to match the counts.

</details>
:::

## Step 4 — Compute ATAC-specific QC metrics

```{r}
#| label: M15-qc_compute

pbmc <- NucleosomeSignal(pbmc)
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

pbmc$blacklist_ratio    <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100

colnames(pbmc@meta.data)
```

::: {.callout-important title="Think about it"}
1.  What does a **high nucleosome signal** indicate?
2.  Why is **TSS enrichment** a quality metric, not a biology metric?
3.  What is the ENCODE blacklist, and why exclude it?

<details>

<summary>Show answers</summary>

1.  High nucleosome signal means an excess of multi-nucleosome-length fragments (vs sub-nucleosomal). Good ATAC libraries show a clear bias toward sub-nucleosomal fragments in open chromatin. High signal flags over-digested or low-quality cells.
2.  Because active TSSs are *reliably* open in essentially all cells. Good cells show a sharp accessibility peak at TSSs; low enrichment at TSSs signals poor Tn5 activity or barcode contamination regardless of cell type.
3.  A curated list of regions with anomalously high signal in *any* experiment (repeats, centromeres, mappability artefacts). Reads falling there are technical noise; a high blacklist ratio indicates a problematic cell or library.

</details>
:::

## Step 5 — Visualize QC

```{r}
#| label: M15-qc_viz
#| fig-cap: "ATAC QC density scatters (depth/nucleosome signal vs TSS enrichment) and per-metric violins for the PBMC 10k cells."

a1 <- DensityScatter(pbmc, x = "nCount_ATAC",
                     y = "TSS.enrichment", log_x = TRUE, quantiles = TRUE) +
  labs(
    title = "Depth vs TSS enrichment",
    x     = "nCount_ATAC (fragments per cell, log scale)",
    y     = "TSS enrichment score"
  )
a2 <- DensityScatter(pbmc, x = "nucleosome_signal",
                     y = "TSS.enrichment", log_x = TRUE, quantiles = TRUE) +
  labs(
    title = "Nucleosome signal vs TSS enrichment",
    x     = "Nucleosome signal (log scale)",
    y     = "TSS enrichment score"
  )
p_qc_density <- (a1 | a2) +
  patchwork::plot_annotation(
    title    = "scATAC-seq QC density scatters — PBMC 10k v1",
    subtitle = "Good cells: sufficient depth, high TSS enrichment, low nucleosome signal",
    caption  = "Module 15 · scATAC-seq with Signac"
  )
p_qc_density

p_qc_vln <- VlnPlot(
  object   = pbmc,
  features = c("nCount_ATAC", "nFeature_ATAC", "TSS.enrichment",
               "nucleosome_signal", "blacklist_ratio", "pct_reads_in_peaks"),
  pt.size  = 0.1,
  ncol     = 6
) &
  theme(plot.title = element_text(size = 10))
p_qc_vln <- p_qc_vln +
  patchwork::plot_annotation(
    title    = "Per-cell ATAC QC metric distributions — PBMC 10k v1",
    subtitle = "Depth, complexity, TSS enrichment, nucleosome signal, blacklist ratio and peak-read fraction",
    caption  = "Module 15 · scATAC-seq with Signac"
  )
p_qc_vln

# --- Figures out -----------------------------------------------------------
save_fig(file.path(out_dir, "Mod15_C6_qc_density_scatter.png"), p_qc_density,
       width = 10, height = 5, dpi = 300)
save_fig(file.path(out_dir, "Mod15_C6_qc_violins.png"),         p_qc_vln,
       width = 14, height = 4, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The **density scatters** (left pair) show each cell as a point; color indicates local density. Good cells cluster in the upper-middle of each panel — moderate-to-high `nCount_ATAC` and high `TSS.enrichment` (above ~2), and low `nucleosome_signal` (below ~2). A cloud of cells at high depth but low TSS enrichment signals a poorly enriched library or debris. The **violin panel** gives the overall distribution of all six QC metrics: look for well-separated, narrow violins without extreme upper tails — long tails on `blacklist_ratio` or `nucleosome_signal` point to cells you'll remove in Step 6. An alternative to these fixed thresholds is `scDblFinder`-style MAD-based outlier detection via `scater::isOutlier()`.
:::

::: {.callout-important title="Think about it"}
1.  In the `nCount_ATAC` × `TSS.enrichment` density plot, where should *good* cells sit?
2.  What would a blob at high `nCount_ATAC` and **low** TSS enrichment suggest?

<details>

<summary>Show answers</summary>

1.  Middle-to-high `nCount_ATAC` with high TSS enrichment — enough reads, concentrated at bona-fide open regions.
2.  Many reads but they're **not** at TSSs. Classic signature of a bad library, ambient reads, or reads from dead cells — not useful signal.

</details>
:::

## Step 6 — Filter low-quality cells

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

pbmc <- subset(
  x = pbmc,
  subset = nCount_ATAC       > 3000  &
           nCount_ATAC       < 30000 &
           pct_reads_in_peaks > 15   &
           blacklist_ratio   < 0.05  &
           nucleosome_signal < 4     &
           TSS.enrichment    > 3
)
pbmc
```

::: {.callout-important title="Think about it"}
1.  These thresholds come from the Signac PBMC vignette. Would you use them as-is on another dataset?
2.  What's the analogue of `percent.mt` in ATAC, and is it relevant?

<details>

<summary>Show answers</summary>

1.  No — ATAC metric distributions shift with tissue, depth, and protocol version. Use the vignette thresholds as a **starting point** and always inspect the violins and density plots first.
2.  Direct analogue doesn't exist (ATAC doesn't measure RNA). The quality proxies are **TSS enrichment**, **nucleosome signal**, **blacklist ratio**, and **pct_reads_in_peaks** — each captures a different failure mode.

</details>
:::

## Step 7 — TF-IDF normalization, top features, SVD

[]{#lst-M15-lsi}

```{r}
#| label: M15-lsi
#| fig-cap: "Correlation between each LSI component and sequencing depth; component 1 is the depth-confounded one to drop."

pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = "q0")
pbmc <- RunSVD(pbmc)

p_depthcor <- DepthCor(pbmc) +
  labs(
    title    = "LSI component correlation with sequencing depth",
    subtitle = "Component 1 typically tracks depth, not biology — exclude it downstream",
    x        = "LSI component",
    y        = "Correlation with depth"
  )
p_depthcor

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

::: {.callout-tip title="Reading the output"}
The `DepthCor` bar chart shows the Pearson correlation of each LSI component with total fragment count (`nCount_ATAC`). Component 1 almost always has a correlation close to −1 or +1 — meaning it tracks library size rather than biology. When you run this chunk, look for that first bar to stand clearly apart from the rest; components 2 onward should hover near zero, meaning they capture chromatin structure. If multiple early components show strong depth correlation, the library preparation may be uneven and worth inspecting.
:::

::: {.callout-important title="Think about it"}
1.  RNA uses `LogNormalize`. ATAC uses **TF-IDF**. Why the different normalization?
2.  Why is the **first LSI component often discarded** (see [code](#lst-M15-lsi))?

<details>

<summary>Show answers</summary>

1.  ATAC data is **binary-ish** — most peaks are either accessible (1) or not (0) per cell. TF-IDF (term-frequency × inverse document-frequency) weighs peaks by rarity, which is well-suited to sparse binary data. LogNormalize is tuned for count RNA data.
2.  The first SVD component typically tracks **sequencing depth / library size** rather than biology. `DepthCor` shows correlation of each LSI component with depth; if component 1 is highly correlated, skip it and use dims 2:30 downstream.

</details>
:::

## Step 8 — Non-linear embedding and clustering

```{r}
#| label: M15-umap_cluster
#| fig-cap: "UMAP of the PBMC 10k scATAC-seq cells coloured by SLM cluster (LSI dims 2:30)."

pbmc <- RunUMAP(pbmc, reduction = "lsi", dims = 2:30)
pbmc <- FindNeighbors(pbmc, reduction = "lsi", dims = 2:30)
pbmc <- FindClusters(pbmc, algorithm = 3)   # SLM

p_umap <- DimPlot(pbmc, label = TRUE) +
  NoLegend() +
  labs(
    title    = "scATAC-seq clusters — PBMC 10k v1",
    subtitle = "UMAP on LSI components 2:30, clustered with SLM (algorithm 3)",
    x        = "UMAP 1",
    y        = "UMAP 2"
  )
p_umap

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod15_C9_umap_clusters.png"), p_umap,
       width = 7, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each point on the UMAP is a cell, colored by SLM cluster. For PBMC scATAC-seq you should see roughly 10–15 well-separated clusters corresponding to the major blood-cell lineages (T cells, B cells, monocytes, NK cells, dendritic cells). Clusters that form tight, round blobs with clear gaps between them indicate good separation; a continuous smear without distinct islands suggests the resolution parameter may need adjustment or that a depth-confounded LSI component was included. The cluster labels are arbitrary integers at this stage — Step 8's stretch goals show how to assign cell-type identities using `GeneActivity` scores or label transfer from a matched scRNA-seq object.
:::

::: {.callout-important title="Think about it"}
1.  Why `dims = 2:30` instead of `1:30`?
2.  `algorithm = 3` selects the SLM (Smart Local Moving) algorithm. Name one case when you'd use it instead of the default Louvain.

<details>

<summary>Show answers</summary>

1.  Because LSI component 1 is the depth-confounded one we saw in `DepthCor`. Starting at `2` keeps only biologically-meaningful dimensions.
2.  SLM often produces more stable partitions than basic Louvain on large, noisy graphs (e.g. scATAC) because it allows within-cluster re-optimization. Leiden (algorithm 4) is usually the modern default; SLM remains a reasonable alternative.

</details>
:::

## Save the object

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

saveRDS(pbmc, file = "../data/pbmc_atac_clustered.rds")
```

## Stretch goals

-   **Gene activity** — convert peaks to per-gene accessibility scores with `GeneActivity()`, then run a Seurat-style RNA-like analysis on top
-   **Motif enrichment** — use `chromVAR` or Signac's `RunChromVAR` for TF-motif activity per cell
-   **Multi-modal integration** — integrate scRNA-seq + scATAC-seq from the same cells (multiome) using `FindTransferAnchors`

### Stretch — motif enrichment with HOMER (alternative to chromVAR)

Signac's `FindMotifs()` and `chromVAR` are the in-R toolchain. **HOMER** (`findMotifsGenome.pl`) is the classic command-line alternative — same question, very different stack. Useful when you want to share results with a collaborator who lives in HOMER, or when you want a known-motif `p`-value table without standing up `chromVAR`.

The input HOMER expects is a tab-separated peak file with `chr / start / end / strand` columns. From your Signac peak set:

```{r}
#| label: M15-peaks_granges_pbmc_peaks
#| eval: false
peaks <- granges(pbmc[["peaks"]])
df <- data.frame(
  chr    = as.character(seqnames(peaks)),
  start  = start(peaks),
  end    = end(peaks),
  strand = "+"
)
write.table(df, "results/peaks_homer.bed", quote = FALSE,
            sep = "\t", row.names = FALSE, col.names = FALSE)
```

Then run HOMER from the shell:

```{bash}
#| label: M15-background_all_open_regions
#| eval: false
# Background = all open regions; foreground = your peaks of interest
findMotifsGenome.pl results/peaks_homer.bed hg19 results/homer_out -size 200

# Locate instances of one motif inside the same peaks
findMotifsGenome.pl results/peaks_homer.bed hg19 results/homer_out \
    -find results/homer_out/knownResults/known1.motif \
  > results/homer_motif_instances.txt

# Or: annotate peaks with known motifs in a single pass
annotatePeaks.pl results/peaks_homer.bed hg19 \
    -m results/homer_out/knownResults/known1.motif \
  > results/homer_motif_annotated.txt
```

::: callout-warning
**HOMER's `-size` flag matters.** `-size 200` centres a 200 bp window on each peak — appropriate for ATAC-seq peaks (typically 500–1000 bp wide). For wider regions or full-peak coverage use `-size given`. See the [HOMER peakMotifs notes](http://homer.ucsd.edu/homer/ngs/peakMotifs.html).
:::

::: callout-tip
**Genome assembly must match your peak coordinates.** This workshop's PBMC 10k v1 data is **hg19**, so use `hg19` here. Run `configureHomer.pl -list` to see what's installed locally, or `configureHomer.pl -install hg19` to fetch it. If your peaks are on a different build (e.g. you re-aligned to hg38), liftOver them first using the call-out earlier in this tutorial — don't pass mismatched coordinates to HOMER.
:::

## Credits

-   Reference vignette: [Signac PBMC scATAC-seq vignette](https://stuartlab.org/signac/articles/pbmc_vignette)
-   Data: 10x Genomics PBMC 10k scATAC-seq