Tutorial 15 — scATAC-seq with Signac

From peaks matrix to clustered cells (PBMC 10k v1, 10x)

Author

Single Cell RNA-seq Workshop

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

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_15_scATACseq_Signac.qmd. 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.
TipWhat’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.

Note

Companion book chapter: Chapter 15 — scATAC-seq — 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.

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/):

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

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

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.

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_15_scATACseq_Signac.html. See Exercise_Folder/_quarto-solutions.yml for the build profile

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.

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

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"))
ImportantThink 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?
Show answers
  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.

Step 2 — Load the peaks matrix and create the object

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
ImportantThink 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?
Show answers
  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.

Step 3 — Add gene annotation

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
seqlevels(annotations) <- paste0("chr", seqlevels(annotations))  # UCSC naming
Annotation(pbmc) <- annotations
ImportantThink 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?
Show answers
  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.

Step 4 — Compute ATAC-specific QC metrics

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)
ImportantThink 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?
Show answers
  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.

Step 5 — Visualize QC

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

ImportantThink 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?
Show answers
  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.

Step 6 — Filter low-quality cells

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
ImportantThink 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?
Show answers
  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.

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

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

ImportantThink 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)?
Show answers
  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.

Step 8 — Non-linear embedding and clustering

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

ImportantThink 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.
Show answers
  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.

Save the object

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:

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:

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

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