How to Analyze Single-Cell ATAC-seq Data — A Complete Beginner’s Guide Part 8: Cell-Type-Specific Peak Recalling with MACS3 and Co-Accessibility with Cicero

Meta description: Recall scATAC-seq peaks per cell type with MACS3, compare them with Cell Ranger consensus peaks, and map co-accessible regulatory elements with Cicero — full R code.

Pay off a promise from Part 1: call peaks separately in every annotated cell type with MACS3, merge them into one unified peak set, measure what it recovers that the Cell Ranger consensus missed, and use the new peaks to find co-accessible regulatory elements with Cicero.

Part 1 of this series made a promise and then deferred it. In Peak Calling: Cell Ranger ATAC vs MACS3, we argued that the decisive reason to use MACS3 is not that it is a better whole-sample peak caller. It is that MACS3 lets you call peaks separately within each cell type, so regulatory elements open only in a small population are not diluted into the pooled signal. That requires knowing your cell types, and in Part 1 we did not. So the series carried Cell Ranger’s peaks forward: Part 3 merged them into a 185,581-peak consensus set, and Parts 4 to 7 all ran on it.

We now have cell types. Part 4 assigned a label to every one of the 30,571 nuclei, and Parts 5 to 7 put those labels to work. It is time to pay the debt.

This tutorial does two things with the same new peak set:

1. Cell-type-specific peak recalling. Run MACS3 once per annotated cell type, merge the results into a unified peak set, and compare it with the Cell Ranger consensus: how many peaks, how wide, how much overlap, and which cell types gained the most.
2. Co-accessibility. Ask which of those regulatory elements open and close together across cells, using Cicero. This links distal elements to each other and to promoters without any RNA data, and cell-type-specific peaks are exactly the right input for it.

Prerequisites: This tutorial assumes you have completed Part 1: From FASTQ to Peaks, Part 2: Thorough Quality Control with Signac, Part 3: Integration and Clustering and Part 4: Cell Type Identification. You need the atac_annotated.rds object from Part 4, the four fragments files it points to (Part 1), the consensus_peaks.rds file saved at the end of Part 3, the ENCODE blacklist from Part 2, and the same Pixi workspace used through Part 7. Parts 5 to 7 are not required, but the co-accessibility example ties back to their CD14 monocyte findings. Everything in this tutorial runs in R.

🧭 Introduction: Why Recall Peaks After Cell Type Annotation?

A Quick Refresher: How Cell Ranger and MACS3 Called Peaks in Part 1

A peak is a genomic region with more Tn5 insertions than background — a putative regulatory element. Peaks are the features of scATAC-seq data, the way genes are the features of scRNA-seq. Part 1 met two ways of calling them.

Cell Ranger ATAC calls peaks automatically inside cellranger-atac count. Since version 2.0 it applies a wavelet-based caller to the aggregate signal: all barcodes of a sample are pooled into one pseudobulk profile, and enriched regions are reported in peaks.bed. The peaks are broad and fairly uniform in width. Part 3 then built a consensus across the four samples:

Part 3 also noticed that the number of peaks each sample called tracked how many cells it captured. That is the first hint of the problem this tutorial addresses.

MACS3 is the community-standard peak caller from the MACS project (Zhang et al. 2008), the same tool used in our bulk ATAC-seq guide. In Part 1 we ran it on one whole sample with -f FRAG, its parser for 10x fragments files, and noted that MACS3 peaks are narrower than Cell Ranger’s. We also said that, on a whole sample, that run was mainly a sanity check. MACS3’s real value comes later.

The Dilution Problem: Why Aggregate Peak Calling Misses Rare Cell Types

Both Cell Ranger and our Part 1 MACS3 run pooled every nucleus of a sample before calling peaks. Pooling works well for regulatory elements open in common cell types, because those contribute strong signal. An element open only in a rare population contributes almost nothing to the pool, and its peak can fall below the detection threshold.

Our dataset has a concrete example. Part 4 labelled 87 nuclei as dendritic cells, 0.28 percent of the 30,571-cell cohort. Cell Ranger called peaks one sample at a time, and in no sample did dendritic cells make up even half a percent of the nuclei:

A regulatory element open only in dendritic cells has to be detected from fewer than 50 nuclei mixed into thousands of others. Taking the union of four samples in Part 3 cannot rescue it, because a union only keeps what at least one sample called.

The fix is to call peaks after annotation, once per cell type, on that cell type’s own pooled fragments. Inside the dendritic-cell pseudobulk, a dendritic-cell element is no longer a whisper among 15,000 nuclei; it is part of the dominant signal. Calling peaks per group of cells and then combining them is what Signac’s CallPeaks() does when you give it a group.by variable, and what ArchR does with its cluster pseudobulk replicates (Stuart et al. 2021; Granja et al. 2021).

Beginner’s Tip: Think of peak calling as listening for voices in a room. Cell Ranger records the whole room at once, so a small group talking quietly in one corner is lost in the crowd. Calling peaks per cell type is like recording each group at its own table. The quiet group was always there; now you can hear it.

When Do You Need Cell-Type-Specific Peak Calling?

Recalling peaks costs time and disk space, and it changes the feature set every downstream analysis uses. It is worth doing when the question depends on what the aggregate peaks cannot see.

The common thread: once you know your cell types, every analysis that asks “what does this cell type do?” should be run on peaks that this cell type helped define.

From Peaks to Regulatory Networks: What Co-Accessibility Adds

A peak tells you where chromatin is open. It does not tell you what that open region does, or which other regions it works with. Enhancers often sit tens to hundreds of kilobases from the genes they regulate, and nothing in a peak list connects the two.

With matched RNA in the same nuclei (multiome data), you could correlate a peak’s accessibility with a nearby gene’s expression — that is what Signac’s LinkPeaks() does. Our dataset has no same-cell RNA. Co-accessibility is the RNA-free alternative: it asks which pairs of peaks tend to be open in the same cells and closed in the same cells. Regions that switch together are candidates for working together.

Cicero (Pliner et al. 2018) is the established method for this, and Signac maintains a vignette for using it. Four ideas from its design shape how you will read the results:

1. It pools similar cells first. Single-cell ATAC data are near-binary and very sparse, so correlations between two peaks across individual cells are mostly noise. Cicero samples groups of cells that sit close together in a low-dimensional embedding (50 cells per group by default) and sums their counts, which turns thousands of sparse cells into a smaller number of denser “aggregate cells”.
2. It works in overlapping windows of the genome. Cicero estimates co-accessibility for all pairs of peaks within overlapping 500 kb windows.
3. It penalizes distance. Within each window it fits a graphical lasso — a regularized estimate of the correlation structure (Friedman et al. 2008) — with a penalty that grows with genomic distance, so that few pairs far apart (beyond 250 kb by default) end up connected.
4. It returns a score per pair, and networks of pairs. The co-accessibility score ranges from -1 to 1; the authors describe positive scores as indicating a connection, and raising the threshold leaves fewer connections with higher confidence. Groups of peaks that are highly co-accessible with each other form cis-co-accessibility networks (CCANs), which the authors propose as candidate chromatin hubs.

Two caveats matter from the start. Cicero links are predictions from correlated accessibility, not measured physical contacts, and the authors caution that they need experimental validation. And co-accessibility can only connect elements that are in the peak set, which is why it belongs in the same tutorial as peak recalling.

Beginner’s Tip: Co-accessibility is guilt by association. If two regions are open in the same cells and closed in the same cells, often enough and more than their distance would predict, Cicero calls them co-accessible. It is evidence that two regions may belong to the same regulatory program, not proof that one controls the other.

The Part 8 Workflow

atac_annotated.rds (Part 4)              consensus_peaks.rds (Part 3)
        |                                         |
        v                                         |
[1] Choose cell types for peak calling     Step 4 |
        |                                         |
        v                                         |
[2] CallPeaks(): split fragments by cell   Step 5 |
    type, run MACS3 once per cell type            |
        |                                         |
        v                                         |
[3] Merge into one unified peak set,       Steps 6-7
    record which cell types called each peak      |
    width, chromosome, blacklist filters          |
        |                                         |
        v                                         v
[4] Compare with the Cell Ranger consensus: size, width,    Steps 8-9
    overlap, new peaks per cell type
        |
        v
[5] Count fragments in the unified peaks (new assay), FRiP  Step 10
        |
        v
[6] Visualization                                           Steps 11-16
        |
        v
[7] Cicero co-accessibility in CD14 monocytes, chromosome 20  Steps 17-20
    -> links, CCANs, links at CEBPB

⚙️ Setting Up the Analysis Environment

What Is New in Part 8

Most of what this tutorial needs is already in the Pixi workspace. MACS3 was installed in Part 1, rtracklayer in Part 2, Signac and Seurat throughout, BSgenome.Hsapiens.UCSC.hg38 in Part 4 and ggrepel in Part 7. Three packages are new:

#-----------------------------------------------
# STEP 1: Add Cicero and ggforce to the Pixi workspace
#-----------------------------------------------

# Run from a LOGIN node (package downloads need internet access)
cd /projects/mylab/shared/scatac-analysis

# HSMMSingleCell is a Bioconductor data package, which Bioconda ships as a stub (see Part 4).
# Let Pixi run Bioconda post-link scripts, so the stub is completed on install.
pixi config set --local run-post-link-scripts insecure

pixi add bioconductor-cicero bioconductor-hsmmsinglecell r-ggforce

# Confirm MACS3 and the R packages this tutorial loads
pixi run macs3 --version
pixi run Rscript -e 'sapply(c("Signac", "cicero", "HSMMSingleCell", "ggforce", "rtracklayer", "BSgenome.Hsapiens.UCSC.hg38"), requireNamespace, quietly = TRUE)'

Expected output: a macs3 version line, then all six package names followed by TRUE.

Loading Libraries and Creating the Part 8 Directory Structure

This tutorial runs in one R session. Launch it from the workspace on a compute node with several CPUs and generous memory — Step 10 recounts fragments for all 29,537 cells we analyze, the same kind of job as Part 3’s requantification:

cd /projects/mylab/shared/scatac-analysis
pixi run R

Launching R through pixi run matters more than usual here: CallPeaks() runs MACS3 as a command-line program, so macs3 must be on the PATH of the R session, and pixi run puts it there.

#-----------------------------------------------
# STEP 2: Load packages and set options
#-----------------------------------------------

library(Signac)
library(Seurat)
library(GenomicRanges)
library(GenomeInfoDb)                 # keepStandardChromosomes()
library(rtracklayer)                  # import() for the ENCODE blacklist
library(BSgenome.Hsapiens.UCSC.hg38)  # chromosome lengths for Cicero
library(ggplot2)
library(ggrepel)                      # non-overlapping labels in Step 14
library(patchwork)
library(future)                       # parallel backend for FeatureMatrix()

set.seed(1234)

# FeatureMatrix() in Step 10 counts fragments in parallel, as in Part 3
plan("multicore", workers = 4)
options(future.globals.maxSize = 50 * 1024^3)

We do not load Cicero yet. It attaches the original Monocle and several of its dependencies, so we load it only when the co-accessibility section begins (Step 17).

#-----------------------------------------------
# STEP 3: Define and create the Part 8 directories
#-----------------------------------------------

base_dir   <- "/projects/mylab/shared/scatac_tutorial"
integ_dir  <- file.path(base_dir, "integrated")      # Part 3: consensus peaks
annot_dir  <- file.path(base_dir, "annotation")      # Part 4: annotated object
recall_dir <- file.path(base_dir, "peak_recalling")  # this tutorial

obj_dir   <- file.path(recall_dir, "objects")
plot_dir  <- file.path(recall_dir, "plots")
table_dir <- file.path(recall_dir, "tables")
macs_dir  <- file.path(recall_dir, "macs3")          # MACS3 working files
split_dir <- file.path(recall_dir, "split_fragments") # temporary per-cell-type fragments

invisible(lapply(c(obj_dir, plot_dir, table_dir, macs_dir, split_dir),
                 dir.create, recursive = TRUE, showWarnings = FALSE))

# ENCODE unified hg38 blacklist (ENCFF356LFX), downloaded in Part 2
blacklist_bed <- "/projects/mylab/shared/reference/blacklist/ENCFF356LFX.bed"

Your project tree now looks like this:

/projects/mylab/shared/scatac_tutorial/
├── results/               # Cell Ranger ATAC output (Part 1)    <- fragments files
├── downstream_qc/         # per-sample QC objects (Part 2)
├── integrated/            # consensus peaks, integrated object (Part 3)  <- comparison
├── annotation/            # annotated object (Part 4)                    <- input
├── differential/          # Part 5
├── footprinting/          # Part 6
├── tf_activity/           # Part 7
└── peak_recalling/        # this tutorial
    ├── objects/
    ├── plots/
    ├── tables/
    ├── macs3/
    └── split_fragments/

Why split_fragments/ lives on project storage. Step 5 writes one uncompressed fragment file per cell type before MACS3 reads them. By default Signac writes these to R’s temporary directory, which on many HPC compute nodes is a small local /tmp. The split files together hold every fragment of every cell included, uncompressed, so budget up to the uncompressed size of your four fragments files. Signac deletes each split file after MACS3 has used it.

📦 Example Data: The Annotated Object From Part 4

Why Start From the Part 4 Object

Peak recalling needs exactly two things: cell-type labels and the fragments files. The Part 4 object has both. It is also the common ancestor of Parts 5, 6 and 7, so anything this tutorial finds describes the same cells and the same labels the rest of the series used.

Loading the Object and Selecting the Ten Cell Types

#-----------------------------------------------
# STEP 4: Load the annotated object from Part 4
#-----------------------------------------------

atac <- readRDS(file.path(annot_dir, "objects", "atac_annotated.rds"))
DefaultAssay(atac) <- "peaks"

# The ten cell types analyzed in this tutorial
cell_types <- c("CD8 T effector", "CD14 monocyte", "NK/T cell unresolved", "CD4 T naive",
                "CD4 T memory", "B cell", "NK cell", "CD8 T naive", "CD16 monocyte",
                "Dendritic cell")
atac <- subset(atac, subset = cell_type %in% cell_types)
atac$cell_type <- droplevels(as.factor(atac$cell_type))  # forget labels with no cells left

# How many cells carry each label, and how they split across donors
sort(table(atac$cell_type), decreasing = TRUE)
table(atac$cell_type, atac$sample_id)

Output:

      CD8 T effector        CD14 monocyte NK/T cell unresolved
                5326                 5213                 4013
         CD4 T naive         CD4 T memory               B cell
                3828                 3148                 2716
             NK cell          CD8 T naive        CD16 monocyte
                1931                 1690                 1585
      Dendritic cell
                  87

                       healthy1 healthy2 severe1 severe2
  CD4 T naive               843     1815      73    1097
  CD4 T memory             1077     1418     170     483
  CD8 T naive               208     1370      27      85
  CD8 T effector           2100     2508     306     412
  NK cell                   308     1425      99      99
  NK/T cell unresolved      729     1864     677     743
  B cell                    269      913     695     839
  CD14 monocyte             807     2575    1323     508
  CD16 monocyte             260     1038      27     260
  Dendritic cell              8       46      16      17

Reading this output. Ten cell types and 29,537 cells, spanning a 60-fold range: from 5,326 CD8 effector T cells down to 87 dendritic cells. Each cell type will become one MACS3 run, so this table is also a preview of how much signal each run has to work with.

Pseudobulks pool donors. Each cell type’s peak call pools that type’s cells from all four donors. That is what gives rare types enough signal, but it also means a donor with many cells of a type has more say in which peaks are called. Keep that in mind when you interpret a peak called in only one cell type.

🔬 Calling Peaks per Cell Type With MACS3

How Signac’s CallPeaks() Runs MACS3 for Each Cell Type

Signac’s CallPeaks() does three things when you give it a group.by variable. We read its source code (Signac 1.17.1) to describe each one exactly:

1. Split the fragments by cell type. It writes one fragment file per group, combining all four fragments files and keeping only the cells assigned to that group.
2. Run MACS3 on each file. For every group, it runs macs3 callpeak and reads the resulting .narrowPeak file back into R as a GRanges object.
3. Optionally merge. With combine.peaks = TRUE (the default), it merges overlapping peaks from all groups with GenomicRanges::reduce() and records which groups called each merged peak in a column called peak_called_in. The documentation notes that merging drops each peak’s statistics (p-value, q-value, fold change).

Step 1 quietly solves a problem that would trip up a hand-written MACS3 command. All four libraries were made with the same 10x barcode whitelist, so the same barcode string occurs in all four fragments files and refers to a different nucleus in each — Part 3 added sample prefixes to cell names for exactly this reason. Signac keeps, for each fragments file, a map from the prefixed cell names in the object back to that file’s barcodes, so a barcode is assigned to the right cell type file by file. Selecting cells by barcode string across all four files at once would mix cells from different donors and cell types.

For each cell type, the command Signac builds is:

macs3 callpeak -t <cell_type>.bed -g 2.7e+09 -f BED --nomodel --extsize 200 --shift -100 \
  -n <cell_type> --outdir <outdir> <additional.args>

Here is what the less obvious flags do:

How this differs from Part 1. Part 1 used -f FRAG, MACS3’s single-cell fragments mode, which piles up whole fragments. CallPeaks() only offers the BED and BEDPE formats, and it uses BED by default, building the signal from insertion sites instead. Both are reasonable ways to model ATAC-seq signal; what matters for a comparison is to use the same settings for every cell type, which CallPeaks() guarantees.

Three settings deserve a closer look before we run anything:

• macs2.path. The argument still carries MACS2’s name, and if you leave it empty Signac looks for a program called macs2. Pass the path to macs3 explicitly. Signac’s own peak-calling vignette does the same.
• --keep-dup all. In BED mode MACS3 keeps, by default, at most one tag at the same position and strand (--keep-dup 1). A cell type’s pseudobulk contains fragments from thousands of nuclei, and two nuclei can carry an insertion at the same base without either being a PCR duplicate — Cell Ranger already removed duplicates within each cell. MACS3’s own FRAG mode always behaves like --keep-dup all for this reason, and Part 1 passed --keep-dup all too. We pass it through additional.args.
• Effective genome size. Signac passes 2.7e9 by default. MACS3’s own hs shortcut, which Part 1’s -g hs used, is 2,913,022,398, a value MACS3 takes from deepTools. We keep Signac’s default. It affects only the genome-wide background estimate, and it is the same for every cell type here.

The q-value threshold is MACS3’s default, 0.05, since we do not override it.

Step 5: Running MACS3 on Each Cell Type

#-----------------------------------------------
# STEP 5: Call peaks with MACS3, once per cell type
#-----------------------------------------------

peaks_by_type <- CallPeaks(
  object           = atac,
  group.by         = "cell_type",
  macs2.path       = unname(Sys.which("macs3")),  # the argument is called macs2.path, but takes MACS3
  outdir           = macs_dir,
  fragment.tempdir = split_dir,
  combine.peaks    = FALSE,                       # keep one peak set per cell type
  additional.args  = "--keep-dup all",            # fragments are already deduplicated
  verbose          = FALSE
)

# Name each peak set by its cell type (spaces and "/" become "_")
peaks_by_type <- Filter(Negate(is.null), peaks_by_type)
names(peaks_by_type) <- vapply(peaks_by_type, function(gr) gr$ident[1], character(1))

saveRDS(peaks_by_type, file.path(obj_dir, "macs3_peaks_by_celltype.rds"))

sapply(peaks_by_type, length)

Output:

       CD14_monocyte NK_T_cell_unresolved               B_cell
               99577                48032                68496
      CD8_T_effector          CD4_T_naive              NK_cell
               64507                58329                55471
        CD4_T_memory        CD16_monocyte          CD8_T_naive
               58525                25566                54568
      Dendritic_cell
               30962

Why combine.peaks = FALSE? With the default TRUE, CallPeaks() would return one merged set and drop every peak’s MACS3 statistics. We keep the ten cell-type sets separate, save them, and merge them ourselves in Step 6. That keeps the summit q-values, which Step 16 uses to pick a peak to plot, and lets you rebuild the unified set later with a different merging rule without rerunning MACS3.

Cell type names change slightly. Signac replaces spaces and / in group names with _ before using them as file names, and those cleaned names end up in the ident column. So CD14 monocyte becomes CD14_monocyte, and NK/T cell unresolved becomes NK_T_cell_unresolved. The rest of this tutorial uses the cleaned names when it refers to peak sets.

Reading this output. One number per cell type: how many peaks MACS3 called in that cell type’s pseudobulk at q < 0.05. The order is the order in which the labels first appear in the object, not size.

• All ten cell types returned peaks, from 25,566 (CD16 monocytes) to 99,577 (CD14 monocytes).
• Peak counts do not simply follow cell numbers. CD14 monocytes (5,213 cells) called 99,577 peaks, while CD8 effector T cells (5,326 cells) called 64,507. B cells called 68,496 peaks from only 2,716 cells.
• Two numbers look backwards. Dendritic cells called 30,962 peaks from 87 nuclei — more than CD16 monocytes called from 1,585. Before reading anything into that, look at how strong the peaks are.

Peak counts alone hide how confident each call is. CallPeaks() returned each peak’s MACS3 statistics, so one more summary shows the quality behind each count:

# Width and MACS3 statistics of each cell type's peaks
peak_stats <- t(sapply(peaks_by_type, function(gr) c(
  median_width  = median(width(gr)),
  pct_200bp     = round(100 * mean(width(gr) == 200), 1),  # peaks at the narrowest width possible
  median_fold   = round(median(gr$fold_change), 2),
  median_log10q = round(median(gr$neg_log10qvalue_summit), 2)
)))
peak_stats

Output:

                     median_width pct_200bp median_fold median_log10q
CD14_monocyte                 378       3.4        9.11         40.06
NK_T_cell_unresolved          418       4.1        8.30         31.11
B_cell                        371       4.5       10.00         33.55
CD8_T_effector                430       2.8       12.01         69.01
CD4_T_naive                   392       3.7       11.58         46.71
NK_cell                       358       5.1       11.21         30.31
CD4_T_memory                  389       3.6       12.03         47.90
CD16_monocyte                 458       8.4        7.44         16.58
CD8_T_naive                   355       5.8       10.96         27.58
Dendritic_cell                232      39.8        5.45          4.95

Reading this table. The dendritic-cell row stands apart from every other.

• Its peaks are the weakest by a wide margin. The median dendritic-cell peak has a -log10 q-value of 4.95 (a q-value of about 1e-5). The next weakest cell type, CD16 monocytes, has 16.58, and CD8 effector T cells have 69.01. Its median fold enrichment, 5.45, is also the lowest.
• 40 percent of its peaks sit at the 200 bp floor (39.8 percent), against 2.8 to 8.4 percent for every other cell type. A peak this narrow is a single pile of overlapping 200 bp windows: the smallest signal MACS3 can call.
• So “more peaks from fewer cells” is not a paradox. With only 87 nuclei, the dendritic-cell pseudobulk is thin, and MACS3 estimates its background from that same thin signal. Many small piles of insertions cross the q < 0.05 threshold. Those calls are statistically supported, but they are the least certain in the whole peak set. Keep this in mind for everything that follows about dendritic cells.
• CD16 monocytes are the low-signal population of the dataset. They called the fewest peaks, and their second-weakest q-values agree with what Part 3 found for their cluster (cluster 7): a median of 32.1 percent of reads in peaks, a TSS enrichment of 5.3 and a nucleosome signal of 0.82, all among the poorest of the 17 clusters.
• CD14 monocytes called many peaks with mid-range strength. These outputs do not say why they called so many more than CD8 effector T cells. One fact worth remembering when you interpret them: the label pools three Part 3 clusters, including clusters 12 and 13, which came almost entirely from severe1 (Step 4).

Runtime and resources. Splitting four fragments files and running ten MACS3 jobs takes a while, and it runs on one CPU. If the session is interrupted, the saved macs3_peaks_by_celltype.rds means you never have to repeat it, and any split files left behind in split_fragments/ can be deleted.

Step 6: Merging the Cell-Type Peak Sets Into One Unified Peak Set

Ten peak sets are useful for questions about one cell type, but every downstream analysis needs one shared feature space — the lesson of Part 3. A region open in both B cells and monocytes will have been called twice, with slightly different boundaries, so we merge overlapping peaks with reduce(), exactly as CallPeaks() would with combine.peaks = TRUE. The revmap it returns records which input peaks went into each merged peak, which tells us which cell types called it.

#-----------------------------------------------
# STEP 6: Merge the cell-type peak sets into one unified peak set
#-----------------------------------------------

all_peaks <- unlist(GRangesList(peaks_by_type), use.names = FALSE)

# reduce() merges overlapping peaks; revmap records which input peaks went into each
unified <- reduce(all_peaks, with.revmap = TRUE)

# Record the cell types each merged peak was called in
called_in <- lapply(unified$revmap, function(i) unique(all_peaks$ident[i]))
unified$called_in      <- CharacterList(called_in)
unified$peak_called_in <- vapply(called_in, paste, character(1), collapse = ",")
unified$n_cell_types   <- lengths(called_in)
unified$revmap <- NULL

length(unified)

Output:

[1] 156037

The ten sets held 564,033 peaks between them (the sum of Step 5’s counts), and reduce() merged them into 156,037 regions.

Each merged peak now carries three new columns:

A note on merging. reduce() merges any peaks that overlap, and Part 3 showed that this can chain staggered peaks into one long region. ArchR avoids this with an iterative overlap procedure (Corces et al. 2018; Granja et al. 2021). We keep reduce() because it is what Signac does and because it lets Step 8 compare like with like — Part 3 merged the Cell Ranger peaks the same way. Step 8 checks whether chaining became a problem.

Step 7: Applying the Same Filters as the Part 3 Consensus

A fair comparison needs both peak sets to pass the same filters. Part 3 applied three to the Cell Ranger union; we apply the same three here, plus one more.

#-----------------------------------------------
# STEP 7: Apply the same filters as the Part 3 consensus set
#-----------------------------------------------

# Filter 1: peak width, the bounds used in Part 3
unified <- unified[width(unified) > 20 & width(unified) < 10000]

# Filter 2: standard chromosomes, then drop chrM (mitochondrial DNA is not chromatin)
unified <- keepStandardChromosomes(unified, pruning.mode = "coarse")
unified <- dropSeqlevels(unified, "chrM", pruning.mode = "coarse")

# Filter 3: ENCODE blacklist
blacklist_hg38 <- import(blacklist_bed, format = "bed")
unified <- subsetByOverlaps(unified, blacklist_hg38, invert = TRUE)

length(unified)
table(unified$n_cell_types)

Output:

[1] 155853

    1     2     3     4     5     6     7     8     9    10
68021 21130 12349  7778  5772  5707  5787 10422  6911 11976

Why drop chrM? keepStandardChromosomes() keeps the mitochondrial genome, because it is a standard sequence of the human assembly. But mitochondrial DNA is not packaged into nuclear chromatin, so accessibility peaks on it do not describe regulatory elements. If MACS3 called none there, this line changes nothing.

Reading this output. 155,853 peaks go forward. The filters removed 184 of 156,037 (0.12 percent), about the same small fraction as in Part 3 (186,150 down to 185,581).

The n_cell_types table is the most informative output of the whole peak-calling section. It counts how many merged peaks were called in exactly 1, 2, 3 … 10 cell types.

• 68,021 peaks (43.6 percent) were called in exactly one cell type. These are the peaks only one pseudobulk had the signal to call, and they are where per-cell-type calling can add something new. Step 9 breaks them down by cell type.
• 11,976 peaks (7.7 percent) were called in all ten. This is the shared chromatin landscape of blood cells, such as promoters of genes used everywhere. Cell Ranger finds these easily.
• The bump at 8 cell types has a simple cause. Counts fall from 1 to 6 cell types, then jump to 10,422 at 8. Listing which cell types are missing from those peaks (from the peak_called_in column) shows that 9,201 of the 10,422 were called in every cell type except CD16 monocytes and dendritic cells — the two weakest peak sets in Step 5. These are mostly shared peaks that the two thinnest pseudobulks did not reach significance on, not elements specific to eight cell types. The same holds at 9: the missing cell type is dendritic cells for 4,140 of 6,911 peaks and CD16 monocytes for 2,385.

📏 Comparing the MACS3 Peak Set With the Cell Ranger Consensus

Step 8: Size, Width and Overlap

Three questions frame the comparison: is the new set larger or smaller, are its peaks wider or narrower, and how much of each set the other one contains?

#-----------------------------------------------
# STEP 8: Compare the unified MACS3 peaks with the Part 3 consensus peaks
#-----------------------------------------------

cr_peaks <- readRDS(file.path(integ_dir, "objects", "consensus_peaks.rds"))

# Size of each peak set: number of peaks, and genome covered (Mb)
data.frame(peak_set   = c("Cell Ranger consensus (Part 3)", "MACS3 per cell type (Part 8)"),
           n_peaks    = c(length(cr_peaks), length(unified)),
           genome_Mb  = round(c(sum(width(reduce(cr_peaks))), sum(width(unified))) / 1e6, 1))

summary(width(cr_peaks))
summary(width(unified))

# Overlap in both directions
unified$in_cellranger <- overlapsAny(unified, cr_peaks)
cr_peaks$recalled     <- overlapsAny(cr_peaks, unified)

table(unified$in_cellranger)
table(cr_peaks$recalled)

# Cell Ranger peaks that contain two or more separate MACS3 peaks
table(pmin(countOverlaps(cr_peaks, unified), 3))

# Width and novelty by the number of cell types that called each peak
round(tapply(width(unified), unified$n_cell_types, median))
round(100 * tapply(!unified$in_cellranger, unified$n_cell_types, mean), 1)

# Peaks at the narrowest width this MACS3 setup can produce (the 200 bp extension)
table(width_200 = width(unified) == 200, new = !unified$in_cellranger)

saveRDS(unified, file.path(obj_dir, "unified_macs3_peaks.rds"))

Output:

                        peak_set n_peaks genome_Mb
1 Cell Ranger consensus (Part 3)  185581     164.7
2   MACS3 per cell type (Part 8)  155853      75.6

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  105.0   836.0   901.0   887.6   944.0  2598.0
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    200     277     389     485     594    3804

 FALSE   TRUE
 43063 112790

 FALSE   TRUE
 73626 111955

     0      1      2      3
 73626 108830   3101     24

  1   2   3   4   5   6   7   8   9  10
280 383 447 484 501 532 545 548 738 980

   1    2    3    4    5    6    7    8    9   10
54.2 21.7  8.3  4.1  2.4  1.1  0.4  0.1  0.0  0.0

         new
width_200  FALSE   TRUE
    FALSE 110729  34107
    TRUE    2061   8956

The Cell Ranger width summary reproduces Part 3 exactly (median 901 bp, quartiles 836 and 944 bp, maximum 2,598 bp), which confirms we loaded the same consensus set.

Reading this output, line by line.

• The MACS3 set is smaller in every sense. It has 155,853 peaks against 185,581 (16 percent fewer), and covers 75.6 Mb of the genome against 164.7 Mb — less than half. Keep that ratio in mind for Step 10.
• Its peaks are much narrower. The median is 389 bp against 901 bp, and half of the peaks fall between 277 and 594 bp. That is the pattern Part 1 described for MACS3 peaks.
• Merging did not chain peaks badly. The widest unified peak is 3,804 bp, wider than Cell Ranger’s maximum (2,598 bp) but far below the 10 kb filter. Width grows steadily with the number of cell types that called a peak: a median of 280 bp for peaks called in one cell type, 980 bp for peaks called in all ten. That is what reduce() does to one element called ten times with slightly different boundaries — it keeps the union of all ten.
• 43,063 unified peaks (27.6 percent) are new: they overlap no peak in the Cell Ranger consensus. The second-to-last table shows where they come from. 54.2 percent of the single-cell-type peaks are new, 21.7 percent of the two-cell-type peaks, and almost none of the peaks called in six or more cell types. New peaks are, overwhelmingly, peaks that only one or two cell types called.
• 73,626 Cell Ranger peaks (39.7 percent) were not recalled by MACS3 in any of the ten cell types. That is a large share, and these outputs cannot say which explanation dominates. Candidates include regions with diffuse signal that Cell Ranger’s wavelet caller accepts and MACS3 does not, and the broad flanks of Cell Ranger peaks, which MACS3 peaks inside them do not need to reach.
• Cell Ranger peaks rarely split into several MACS3 peaks. Of the 111,955 recalled Cell Ranger peaks, 108,830 contain exactly one MACS3 peak, 3,101 contain two and 24 contain three or more. The resolution gain is mostly a narrower peak inside the same element, not two elements separated.
• The last table counts the peaks at the 200 bp floor: 11,017 of them (2,061 + 8,956), 7.1 percent of the set. 8,956 of those (81.3 percent) are new. These are the smallest calls MACS3 can make, and Step 5 showed where most of them come from: dendritic cells.

Step 9: Which Cell Types Gained the Most New Peaks?

The dilution argument makes a testable prediction: if Cell Ranger misses rare-population elements, the peaks private to small cell types should be absent from the Cell Ranger consensus more often than the peaks private to large ones.

#-----------------------------------------------
# STEP 9: Which cell types gained the most new peaks?
#-----------------------------------------------

n_cells <- table(gsub("[ /]", "_", atac$cell_type))

type_summary <- do.call(rbind, lapply(names(peaks_by_type), function(ct) {
  in_type <- any(unified$called_in == ct)          # unified peaks called in this cell type
  private <- in_type & unified$n_cell_types == 1   # ... and in no other cell type
  data.frame(cell_type        = ct,
             n_cells          = as.integer(n_cells[ct]),
             n_peaks          = sum(in_type),
             n_private        = sum(private),
             n_private_new    = sum(private & !unified$in_cellranger),
             pct_private_new  = round(100 * sum(private & !unified$in_cellranger) / max(sum(private), 1), 1))
}))
type_summary <- type_summary[order(type_summary$n_cells), ]
rownames(type_summary) <- NULL
type_summary

write.csv(type_summary, file.path(table_dir, "peaks_by_celltype_summary.csv"), row.names = FALSE)

# Do smaller cell types have a larger share of new private peaks?
cor.test(type_summary$n_cells, type_summary$pct_private_new, method = "spearman")

Output:

              cell_type n_cells n_peaks n_private n_private_new pct_private_new
1        Dendritic_cell      87   30332      4992          4515            90.4
2         CD16_monocyte    1585   24525       735           488            66.4
3           CD8_T_naive    1690   52217      1618          1062            65.6
4               NK_cell    1931   52397      2363          1470            62.2
5                B_cell    2716   66410     11571          5472            47.3
6          CD4_T_memory    3148   56361      1614           940            58.2
7           CD4_T_naive    3828   56733      2505          1553            62.0
8  NK_T_cell_unresolved    4013   46503       544           267            49.1
9         CD14_monocyte    5213   98141     37276         18712            50.2
10       CD8_T_effector    5326   63767      4803          2414            50.3

	Spearman's rank correlation rho

data:  type_summary$n_cells and type_summary$pct_private_new
S = 290, p-value = 0.01592
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho
-0.7575758

Reading this table. Read it from the top, where the smallest cell types are.

• The dilution argument holds, and dendritic cells show it most strongly. 4,515 of the 4,992 dendritic-cell private peaks (90.4 percent) are absent from the Cell Ranger consensus. For the next three smallest cell types the share is 62 to 66 percent, and for the three largest it is 49 to 50 percent. The Spearman correlation between cell number and that share is -0.76 (p = 0.016): the fewer cells a type has, the more of its private peaks Cell Ranger missed. With ten cell types this is a trend, not a law: B cells (47.3 percent) sit lower than their size predicts, and CD4 naive T cells (62.0 percent) higher.
• Weigh the dendritic-cell result against Step 5. Those 4,515 peaks are the most novel in the set and also the weakest. Both facts follow from the same 87 nuclei: aggregate calling could not see this population at all, and per-cell-type calling sees it only faintly.
• CD14 monocytes contribute most of the private peaks in absolute terms: 37,276, which is 55 percent of all 68,021 private peaks, and 18,712 of them are new. B cells come second with 11,571. Neither has a close relative that shares its peaks: B cells have none in this dataset, and the only close relative of CD14 monocytes, CD16 monocytes, called few peaks of its own.
• T-cell subsets have few private peaks (1,614 to 4,803). CD4 naive, CD4 memory, CD8 naive and CD8 effector T cells share most of their regulatory landscape, so most of their peaks are called in several T-cell types at once. That is biology, not a failure of the method.
• NK_T_cell_unresolved has the fewest private peaks of all (544). That fits Part 4’s difficulty in naming it: its open chromatin is mostly shared with the NK and T cells it sits between.
• Private does not mean exclusive. A peak called only in dendritic cells may still be weakly accessible in other cell types. “Private” here means only one pseudobulk crossed MACS3’s threshold. Step 16 shows what that looks like in the raw signal.

Step 10: Quantifying the Unified Peaks in Every Cell

New peaks need new counts — Part 3’s central lesson. FeatureMatrix() rereads the four fragments files and counts, for every cell, the fragments in every unified peak. We store the result as a second assay called macs3, next to the original peaks assay, so both feature sets live in the same object.

#-----------------------------------------------
# STEP 10: Quantify the unified peaks in every cell
#-----------------------------------------------

macs3_counts <- FeatureMatrix(
  fragments = Fragments(atac),
  features  = unified,
  cells     = colnames(atac),
  verbose   = FALSE
)

# A second assay next to "peaks"; same fragments, same gene annotation
atac[["macs3"]] <- CreateChromatinAssay(
  counts     = macs3_counts,
  fragments  = Fragments(atac),
  annotation = Annotation(atac[["peaks"]])
)

# Fraction of each cell's fragments that fall in each peak set
atac <- FRiP(atac, assay = "peaks", total.fragments = "passed_filters",
             col.name = "FRiP_cellranger", verbose = FALSE)
atac <- FRiP(atac, assay = "macs3", total.fragments = "passed_filters",
             col.name = "FRiP_macs3", verbose = FALSE)

frip_by_type <- aggregate(cbind(FRiP_cellranger, FRiP_macs3) ~ cell_type,
                          data = atac[[]], FUN = median)
frip_by_type$difference <- frip_by_type$FRiP_macs3 - frip_by_type$FRiP_cellranger
frip_by_type[order(frip_by_type$difference), ]

saveRDS(atac, file.path(obj_dir, "atac_macs3_peaks.rds"))

Output:

              cell_type FRiP_cellranger FRiP_macs3    difference
9         CD16 monocyte       0.3609660  0.3376516 -0.0233143341
6  NK/T cell unresolved       0.4953984  0.4834559 -0.0119425287
1           CD4 T naive       0.6721615  0.6676047 -0.0045568352
2          CD4 T memory       0.7304714  0.7289859 -0.0014854255
3           CD8 T naive       0.7188150  0.7176845 -0.0011304925
5               NK cell       0.7429083  0.7427256 -0.0001827709
4        CD8 T effector       0.7631408  0.7634949  0.0003540776
7                B cell       0.6858093  0.6903836  0.0045742566
8         CD14 monocyte       0.6565171  0.6679814  0.0114642976
10       Dendritic cell       0.6367603  0.6702114  0.0334511246

FRiP — the fraction of reads in peaks — is each cell’s peak counts (the column sum of the assay) divided by its total fragments. We use passed_filters, Cell Ranger’s per-barcode count of fragments passing its filters, which has been in the metadata since Part 2. Because both FRiP columns use the same numerator rule and the same denominator, they are directly comparable; only the peak set differs.

Reading this output.

• The MACS3 set captures almost the same share of fragments in less than half the genome. For T cells, NK cells and B cells, the two median FRiP values differ by less than 0.005, even though the MACS3 peaks cover 75.6 Mb against Cell Ranger’s 164.7 Mb (Step 8). Fragments are packed about twice as densely into the MACS3 peaks. The extra width of the Cell Ranger peaks holds few additional fragments — which is why narrower peaks are an advantage for motif analysis.
• The three cell types that gain are the three with the most new private peaks. Dendritic cells gain the most (+0.033), then CD14 monocytes (+0.011) and B cells (+0.005); Step 9 credited them with 4,515, 18,712 and 5,472 new private peaks. More of their fragments now fall in peaks that describe them — the per-cell-type benefit in a second, independent form.
• The cell types that lose had little say in the new peak set. CD16 monocytes (-0.023) and NK_T_cell_unresolved (-0.012) had the fewest private peaks (735 and 544). They lean on peaks called mainly in other cell types.
• CD16 monocytes have the lowest FRiP of any cell type under either set (0.36 and 0.34), consistent with the low signal-to-noise Part 3 saw for their cluster and Step 5 saw in their peak statistics.

📊 Visualizing the Peak Set Comparison

The same R session continues. Six figures go from the overall picture down to a single locus.

Number of Peaks Called per Cell Type

#-----------------------------------------------
# STEP 11: Number of peaks called per cell type
#-----------------------------------------------

p_counts <- ggplot(type_summary,
                   aes(x = reorder(cell_type, n_cells), y = n_peaks)) +
  geom_col(fill = "#2E86AB") +
  geom_text(aes(label = paste0(n_cells, " cells")), hjust = -0.1, size = 3) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.25))) +
  labs(x = NULL, y = "Unified peaks called in this cell type",
       title = "MACS3 peaks per cell type") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "01_peaks_per_celltype.png"), p_counts,
       width = 7, height = 5, dpi = 300)

Reading this figure: one bar per cell type, ordered by cell number with the largest at the top, and labelled with that number. The bar is the number of unified peaks the cell type called.

CD14 monocytes stand out with 98,141 peaks, half again as many as the next cell type (B cells, 66,410). Bar length only loosely follows cell number: B cells, with half as many cells as CD8 effector T cells, called more peaks.

The two shortest bars are the two weakest peak sets from Step 5. CD16 monocytes have the shortest bar (24,525 peaks from 1,585 cells). Dendritic cells come next (30,332 from 87 cells), and Step 5 showed that many of those peaks are minimal calls. A bar this long from so few cells is a reason to look at peak strength, not a measure of how open the chromatin is.

Peak Width: Cell Ranger Consensus vs MACS3

#-----------------------------------------------
# STEP 12: Peak width distributions
#-----------------------------------------------

width_df <- data.frame(
  peak_set = rep(c("Cell Ranger consensus", "MACS3 per cell type"),
                 times = c(length(cr_peaks), length(unified))),
  width    = c(width(cr_peaks), width(unified))
)

p_width <- ggplot(width_df, aes(x = width, fill = peak_set)) +
  geom_histogram(bins = 100, position = "identity", alpha = 0.5) +
  scale_x_log10(breaks = c(100, 250, 500, 1000, 2500, 5000, 10000)) +
  scale_fill_manual(values = c("Cell Ranger consensus" = "grey50",
                               "MACS3 per cell type"   = "#C0392B"), name = NULL) +
  labs(x = "Peak width (bp, log scale)", y = "Number of peaks",
       title = "Peak width: Cell Ranger consensus vs MACS3") +
  theme_classic(base_size = 12) +
  theme(legend.position = "top")

ggsave(file.path(plot_dir, "02_peak_width_comparison.png"), p_width,
       width = 7, height = 4.5, dpi = 300)

Reading this figure: two overlaid histograms on a log scale. The grey one is the Part 3 consensus, which Part 3 showed has a single sharp mode just under 1 kb. The red one is the unified MACS3 set.

The two distributions barely overlap. Cell Ranger peaks form one tall, narrow mode at about 900 bp. MACS3 peaks spread broadly from 200 bp to about 1 kb, with most between roughly 250 and 600 bp (quartiles 277 and 594 bp).

The red spike at the far left is the 200 bp floor. 11,017 peaks are exactly 200 bp wide (Step 8), because Signac’s --extsize 200 makes every insertion a 200 bp window. 81.3 percent of them are new, and Step 5 showed they come mostly from dendritic cells. If you see a spike like this in your own data, check which cell type produces it before treating those peaks as regulatory elements.

The right tail is thin. MACS3 peaks above 2 kb are rare and the longest is 3,804 bp, so merging did not chain neighbouring elements into long regions. A second bump in the kilobase range would have been the warning sign.

How Much Do the Two Peak Sets Share?

#-----------------------------------------------
# STEP 13: Overlap between the two peak sets
#-----------------------------------------------

overlap_df <- rbind(
  data.frame(peak_set = "MACS3 per cell type",
             status   = ifelse(unified$in_cellranger, "Overlaps the other set", "Unique to this set")),
  data.frame(peak_set = "Cell Ranger consensus",
             status   = ifelse(cr_peaks$recalled, "Overlaps the other set", "Unique to this set"))
)

p_overlap <- ggplot(overlap_df, aes(x = peak_set, fill = status)) +
  geom_bar() +
  scale_fill_manual(values = c("Overlaps the other set" = "grey70",
                               "Unique to this set"     = "#C0392B"), name = NULL) +
  labs(x = NULL, y = "Number of peaks", title = "How much do the two peak sets share?") +
  theme_classic(base_size = 12) +
  theme(legend.position = "top")

ggsave(file.path(plot_dir, "03_peak_set_overlap.png"), p_overlap,
       width = 6, height = 5, dpi = 300)

Reading this figure: two stacked bars, one per peak set. Grey is the part of each set that overlaps the other; red is the part unique to it.

Why the two grey segments differ in height. Overlap is counted from each set’s point of view. One broad Cell Ranger peak overlapped by two narrow MACS3 peaks counts once in the Cell Ranger bar and twice in the MACS3 bar. That is why we used overlapsAny() in both directions instead of one Venn diagram: a Venn diagram of intervals with different widths has no single correct count for its middle.

On the MACS3 bar, 43,063 of 155,853 peaks (27.6 percent) are red: new elements the series has not had a feature for since Part 1. On the Cell Ranger bar, 73,626 of 185,581 peaks (39.7 percent) are red: peaks no cell type recalled with MACS3. Each method keeps a substantial set the other does not, so the choice of peak set changes the features every downstream analysis sees.

Where Do the New Peaks Come From?

#-----------------------------------------------
# STEP 14: New cell-type-private peaks versus cell number
#-----------------------------------------------

p_private <- ggplot(type_summary, aes(x = n_cells, y = pct_private_new)) +
  geom_point(size = 3, color = "#C0392B") +
  geom_text_repel(aes(label = cell_type), size = 3, min.segment.length = 0) +
  scale_x_log10(expand = expansion(mult = 0.15)) +
  labs(x = "Cells in the cell type (log scale)",
       y = "Private peaks absent from Cell Ranger set (%)",
       title = "Where do the new peaks come from?") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "04_private_new_peaks.png"), p_private,
       width = 7, height = 5, dpi = 300)

Reading this figure: one point per cell type. The x-axis is cell number on a log scale, from 87 dendritic cells on the left to more than 5,000 CD8 effector T cells and CD14 monocytes on the right. The y-axis is the percentage of that cell type’s private peaks that the Cell Ranger consensus does not contain.

This is the figure that tests the Part 1 argument. If aggregate peak calling dilutes rare-population elements, the points should fall from left to right: small cell types should have the largest share of private peaks that Cell Ranger missed.

The points do fall from left to right (Spearman rho = -0.76, Step 9). Dendritic cells sit alone in the top-left corner at 90.4 percent. CD16 monocytes, CD8 naive T cells and NK cells cluster at 62 to 66 percent. The three largest cell types sit at 49 to 50 percent.

Two points sit off the trend. CD4 naive T cells (62.0 percent) are higher than their size predicts. B cells (47.3 percent) are lower, so more of their private peaks were already in the Cell Ranger consensus. These outputs do not say why either cell type departs from the trend.

FRiP per Cell Type Under Each Peak Set

#-----------------------------------------------
# STEP 15: FRiP under each peak set, by cell type
#-----------------------------------------------

frip_df <- data.frame(
  cell_type = rep(atac$cell_type, 2),
  peak_set  = rep(c("Cell Ranger consensus", "MACS3 per cell type"), each = ncol(atac)),
  FRiP      = c(atac$FRiP_cellranger, atac$FRiP_macs3)
)

p_frip <- ggplot(frip_df, aes(x = cell_type, y = FRiP, fill = peak_set)) +
  geom_boxplot(outlier.size = 0.2) +
  scale_fill_manual(values = c("Cell Ranger consensus" = "grey60",
                               "MACS3 per cell type"   = "#C0392B"), name = NULL) +
  labs(x = NULL, y = "Fraction of fragments in peaks",
       title = "FRiP per cell type under each peak set") +
  theme_classic(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "top")

ggsave(file.path(plot_dir, "05_frip_by_celltype.png"), p_frip,
       width = 9, height = 5, dpi = 300)

Reading this figure: one pair of boxes per cell type, grey for the Part 3 consensus and red for the MACS3 set. Each box summarizes per-cell FRiP for that cell type, so this is the distribution behind Step 10’s medians.

Compare the pairs, not the boxes. Cell types differ in FRiP for biological and technical reasons (Part 3 saw this for cluster 7, now labelled CD16 monocytes), and that difference is the same under both peak sets. What changes between the grey and red box of the same cell type is the peak set alone.

For most cell types the red and grey boxes are almost identical, whole distributions as well as medians. The visible shifts are the ones in Step 10: dendritic cells and CD14 monocytes move up, and CD16 monocytes move down.

The ranking of cell types does not change. NK/T cell unresolved and CD16 monocytes have the lowest FRiP under both sets. Their low values come from the cells, not from the peak set; for CD16 monocytes, this is what Part 3 found when it flagged cluster 7.

A Closer Look: The Strongest New Dendritic-Cell Peak

Numbers say how many peaks are new. A coverage plot shows what one of them looks like. We pick the dendritic-cell peak that is private to dendritic cells, absent from the Cell Ranger consensus, and has the highest MACS3 summit q-value — the strongest example the data offer, chosen by a rule rather than by eye.

#-----------------------------------------------
# STEP 16: Look at the strongest new dendritic-cell peak
#-----------------------------------------------

# Dendritic-cell MACS3 peaks that are private to dendritic cells and absent from Cell Ranger
dc       <- peaks_by_type[["Dendritic_cell"]]
dc_new   <- unified[unified$n_cell_types == 1 & any(unified$called_in == "Dendritic_cell") &
                    !unified$in_cellranger]
dc_hits  <- dc[overlapsAny(dc, dc_new)]
top_peak <- dc_hits[which.max(dc_hits$neg_log10qvalue_summit)]

DefaultAssay(atac) <- "macs3"   # the peak track now shows the MACS3 peaks

p_dc <- CoveragePlot(
  object            = atac,
  region            = top_peak,
  group.by          = "cell_type",
  ranges            = cr_peaks,
  ranges.title      = "Cell Ranger",
  extend.upstream   = 5000,
  extend.downstream = 5000
)

ggsave(file.path(plot_dir, "06_coverage_new_dc_peak.png"), p_dc,
       width = 9, height = 8, dpi = 300)

top_peak

Output:

GRanges object with 1 range and 7 metadata columns:
      seqnames              ranges strand |                   name     score
         <Rle>           <IRanges>  <Rle> |            <character> <integer>
  [1]    chr12 121520617-121521014      * | Dendritic_cell_peak_..       150
      fold_change neg_log10pvalue_summit neg_log10qvalue_summit
        <numeric>              <numeric>              <numeric>
  [1]     10.8974                19.0614                15.0636
      relative_summit_position          ident
                     <integer>    <character>
  [1]                      199 Dendritic_cell
  -------
  seqinfo: 28 sequences from an unspecified genome; no seqlengths

R also prints two warnings from ggplot2 (Removed 156 rows containing missing values or values outside the scale range). They come from drawing elements that fall outside the plotted window, and the figure is complete.

Reading this figure: each coverage track is one cell type, normalized for cell number and depth, across 10 kb around the selected peak. Below them are gene models, the Peaks track (the MACS3 unified peaks, because the default assay is now macs3), and the Cell Ranger track (the Part 3 consensus, added with ranges).

The selected peak is chr12:121,520,617-121,521,014 (398 bp), inside the KDM2B gene as drawn in the gene track. Its fold enrichment is 10.9 and its summit -log10 q-value is 15.06 — about three times the median dendritic-cell peak (4.95, Step 5). Because we chose the strongest qualifying peak, it is the best case, not a typical one.

The figure shows exactly the pattern the dilution argument predicts. The dendritic-cell track has a clear pileup at the peak. The other nine cell types are flat, or nearly flat, at that position. The Peaks track has a bar there, and the Cell Ranger track has none.

The Cell Ranger peak about 2 kb to the right is a different element. It overlaps the second MACS3 peak in the window, where B cells carry weak signal. A dendritic-cell element 2 kb from a Cell Ranger peak was still invisible to the consensus, because Cell Ranger peaks are defined by where pooled signal is enriched, not by proximity.

What this figure does not show is what this element does. Whether it regulates KDM2B in dendritic cells, or anything else, is a question for other data.

🔗 Co-Accessibility With Cicero

With a peak set that includes cell-type-specific elements, we can ask which of them work together. We run Cicero in CD14 monocytes, for continuity: they are the population Parts 5, 6 and 7 studied, and Part 7’s strongest candidates were C/EBP-family motifs. To keep the run short, we analyze chromosome 20, which carries CEBPB, the gene encoding the transcription factor whose motif was the leading candidate in Parts 6 and 7. Signac’s Cicero vignette restricts its example to a single chromosome for the same reason; Best Practice 14 explains how to scale up.

If you restarted R after Step 10, reload what the co-accessibility steps need:

# Only if you restarted R after Step 10: rerun Steps 2 and 3, then reload
atac     <- readRDS(file.path(obj_dir, "atac_macs3_peaks.rds"))
unified  <- readRDS(file.path(obj_dir, "unified_macs3_peaks.rds"))
cr_peaks <- readRDS(file.path(integ_dir, "objects", "consensus_peaks.rds"))

Step 17: Preparing CD14 Monocyte Input for Cicero

Cicero (the Bioconductor release) works on its own data object, a CellDataSet. Its make_atac_cds() function builds one from a simple three-column table — peak, cell, count — so we convert the monocyte counts into that form. We use only the unified peaks that MACS3 called in CD14 monocytes, which is precisely the cell-type-specific input co-accessibility needs.

#-----------------------------------------------
# STEP 17: Prepare CD14 monocyte input for Cicero
#-----------------------------------------------

library(cicero)   # loaded here: it attaches monocle and its dependencies

# CD14 monocytes, and the unified peaks MACS3 called in them, on chromosome 20
mono       <- subset(atac, subset = cell_type == "CD14 monocyte")
mono_peaks <- unified[any(unified$called_in == "CD14_monocyte") & seqnames(unified) == "chr20"]

# Peak IDs in the same "chr-start-end" form as the assay's row names
peak_id <- function(gr) paste(seqnames(gr), start(gr), end(gr), sep = "-")

mono_counts <- LayerData(mono[["macs3"]], layer = "counts")[peak_id(mono_peaks), ]

# Cicero takes a long table: one row per non-zero (peak, cell, count)
nz <- Matrix::summary(mono_counts)
cicero_input <- data.frame(peak  = rownames(mono_counts)[nz$i],
                           cell  = colnames(mono_counts)[nz$j],
                           count = nz$x)

input_cds <- make_atac_cds(cicero_input, binarize = TRUE)

# Group similar cells into overlapping pools of 50, using the Part 3 UMAP
umap_coords <- Embeddings(mono, reduction = "umap.harmony")[colnames(input_cds), ]
cicero_cds  <- make_cicero_cds(input_cds, reduced_coordinates = umap_coords, silent = TRUE)

What each piece does:

• peak_id() writes peak coordinates as chr20-1000-1500, the form Signac uses for row names. Signac’s own GRangesToString() does the same thing but has been deprecated since Signac 1.17.0, so we use a one-line helper instead.
• Matrix::summary() lists only the non-zero entries of the sparse count matrix: row index i, column index j and value x. That is exactly the long table make_atac_cds() expects, and it never builds the full, mostly-zero matrix.
• binarize = TRUE turns every count into 0 or 1 — accessible or not — as in the Cicero vignette. With at most two copies of each locus per nucleus, presence or absence is the information scATAC-seq carries per cell.
• make_cicero_cds() builds the aggregate cells described in the introduction: groups of 50 neighbouring cells in the UMAP, summed and normalized by size factors. We use umap.harmony, the batch-corrected UMAP from Part 3, so cells are grouped by chromatin state rather than by donor batch.
• [colnames(input_cds), ] matters. make_atac_cds() keeps only cells with at least one non-zero count in the chosen peaks, and make_cicero_cds() requires the coordinates to have exactly the same cells. Subsetting the UMAP to colnames(input_cds) guarantees that.

Why one cell type? Co-accessibility comes from covariation across cells. Across all PBMCs, the strongest covariation is between lineages: every B-cell element is open in B cells and closed in everything else, so B-cell elements would look co-accessible with each other simply for being B-cell elements. Within monocytes, what co-varies is the variation between monocyte states — closer to the question of which elements work together in a monocyte. Note that this monocyte population includes Part 3’s severe1-dominated clusters 12 and 13, so some of that variation is a difference between those states and the rest. Cicero’s links pool healthy and severe monocytes; they are not a comparison between conditions.

Step 18: Running Cicero and Building Co-Accessibility Networks

#-----------------------------------------------
# STEP 18: Run Cicero and build co-accessibility networks
#-----------------------------------------------

genome_df <- data.frame(chr    = "chr20",
                        length = seqlengths(BSgenome.Hsapiens.UCSC.hg38)[["chr20"]])

set.seed(1234)
conns <- run_cicero(cicero_cds, genomic_coords = genome_df, sample_num = 100, silent = TRUE)

# Cis-co-accessibility networks (CCANs): groups of peaks linked to each other
ccans <- generate_ccans(conns)

saveRDS(conns, file.path(obj_dir, "cicero_conns_cd14mono_chr20.rds"))
write.csv(ccans, file.path(table_dir, "cicero_ccans_cd14mono_chr20.csv"), row.names = FALSE)

nrow(conns)
summary(conns$coaccess)
length(unique(ccans$CCAN))

# How many peaks belong to a network, and how large are the networks?
nrow(ccans)
summary(as.vector(table(ccans$CCAN)))

Output:

[1] "Coaccessibility cutoff used: 0.09"
[1] 161006
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's
-0.44236  0.00000  0.00000  0.00201  0.00000  0.62456     2876
[1] 59
[1] 770
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   3.00    4.00    7.00   13.05   17.00   49.00

The first line is printed by generate_ccans() as it runs, which is why it appears before the others.

What happens inside run_cicero(), from its source code:

1. Estimate the distance penalty. Cicero samples sample_num windows (here 100) and finds the penalty strength at which few pairs of peaks more than 250 kb apart remain connected. It uses the mean across the sampled windows for the whole run. set.seed() makes that sampling reproducible.
2. Fit the models. It fits a graphical lasso in each overlapping 500 kb window along the chromosome.
3. Assemble the connections. It turns each model into correlations and reconciles pairs that appear in more than one overlapping window.

genome_df tells Cicero which chromosome to tile and how long it is; we take the length of chromosome 20 from the BSgenome package installed in Part 4, so there is no number to type by hand.

Reading this output.

• conns has 161,006 rows: 80,503 distinct pairs of peaks, each listed twice. Every pair appears once as A-B and once as B-A (we checked this on Cicero’s output), and Step 19 removes the duplicates. They come from the 2,928 CD14-monocyte peaks on chromosome 20 (Step 19 counts them).
• Most scores are exactly zero. The first quartile, the median and the third quartile are all 0.00000. That is the graphical lasso at work: its penalty sets most correlations to exactly zero rather than to small values (Friedman et al. 2008), and Cicero makes the penalty larger for distant pairs. What survives is a short tail of positive scores, up to 0.62, and a shorter tail of negative ones.
• 2,876 pairs (1.8 percent) have no score. These pairs fell in two overlapping windows that disagreed in sign, and Cicero’s reconcile() function reports NA for them. They are ambiguous and are dropped in Step 19.
• generate_ccans() chose a cutoff of 0.09. It picks the cutoff automatically, groups peaks connected above it with Louvain community detection, and keeps networks of at least three peaks. Note that 0.09 is much lower than the 0.25 we use for “strong” links in Step 19; a CCAN can contain peaks joined only by weak links.
• 59 CCANs, containing 770 peaks between them. About a quarter of the 2,928 monocyte peaks on chromosome 20 belong to a network. Most networks are small (median 7 peaks), and the largest has 49.

Step 19: Do the Newly Recovered Peaks Take Part in Co-Accessibility?

This step joins the two halves of the tutorial. If the new peaks are real regulatory elements, some of them should be linked to other elements.

#-----------------------------------------------
# STEP 19: Do the newly recovered peaks take part in co-accessibility?
#-----------------------------------------------

# Each pair appears twice (A-B and B-A); keep one copy, and only strong links
strong <- conns[!is.na(conns$coaccess) & conns$coaccess > 0.25 &
                as.character(conns$Peak1) < as.character(conns$Peak2), ]

new_ids <- peak_id(mono_peaks[!mono_peaks$in_cellranger])
strong$link_type <- ifelse(strong$Peak1 %in% new_ids | strong$Peak2 %in% new_ids,
                           "Involves a new MACS3 peak", "Both peaks in Cell Ranger set")

nrow(strong)
table(strong$link_type)
mean(!mono_peaks$in_cellranger)

Output:

[1] 227

Both peaks in Cell Ranger set     Involves a new MACS3 peak
                          221                             6
[1] 0.2400956

Two details of the filter:

• Peak1 < Peak2 keeps one copy of each pair: of A-B and B-A, exactly one has its first peak name sorting before its second.
• The 0.25 threshold is the co-accessibility cutoff the Cicero vignette uses in its plotting example. It is a convention for “clearly positive”, not a significance test, so report the threshold with your results and check that your conclusions hold at nearby values.

Reading this output. 24.0 percent of the monocyte peaks on chromosome 20 are new. But only 6 of the 227 strong links (2.6 percent) involve one of them. The new peaks take part in strong co-accessibility about ten times less often than their numbers would suggest.

Why? One more check compares the two groups of peaks directly: how strong their MACS3 calls were, and how often each group appears in a strong link or a CCAN.

# Strength of each peak: the best MACS3 summit q-value among the CD14 monocyte calls inside it
cd14_calls <- peaks_by_type[["CD14_monocyte"]]
hits   <- findOverlaps(mono_peaks, cd14_calls)
best_q <- tapply(cd14_calls$neg_log10qvalue_summit[subjectHits(hits)], queryHits(hits), max)

peak_check <- data.frame(
  peak_type = ifelse(mono_peaks$in_cellranger, "In Cell Ranger set", "New"),
  width     = width(mono_peaks),
  log10q    = as.numeric(best_q[as.character(seq_along(mono_peaks))]),
  in_link   = peak_id(mono_peaks) %in% c(strong$Peak1, strong$Peak2),
  in_ccan   = peak_id(mono_peaks) %in% ccans$Peak
)

aggregate(cbind(width, log10q) ~ peak_type, data = peak_check, FUN = median)
aggregate(cbind(in_link, in_ccan) ~ peak_type, data = peak_check, FUN = mean)
table(peak_check$peak_type)

Output:

           peak_type width  log10q
1 In Cell Ranger set   513 69.1943
2                New   276 11.4149

           peak_type     in_link    in_ccan
1 In Cell Ranger set 0.082247191 0.32449438
2                New 0.008534851 0.06827881

In Cell Ranger set                New
              2225                703

Reading this output.

• The new peaks are much weaker calls. Their median summit -log10 q-value is 11.4, against 69.2 for peaks the Cell Ranger set also has, and they are about half as wide (276 versus 513 bp). Every one of them still passed MACS3’s q < 0.05, but they are the faint end of the monocyte peak set.
• They join links and networks far less often. 0.9 percent of new peaks (6 of 703) are in a strong link, against 8.2 percent of the others; 6.8 percent are in a CCAN, against 32.4 percent.
• The two facts are probably connected. Cicero estimates co-accessibility from how accessibility varies together across aggregated cells. A weak peak contributes few fragments to each aggregate, which leaves little variation to correlate. So the low link rate fits weak signal, and it does not by itself show that the new peaks are noise or that they are isolated elements. These outputs cannot separate those possibilities.

What this means for the series. Per-cell-type calling added 703 peaks to the monocyte map of chromosome 20, and 6 of the strong links found there could not have been found with the peak set Parts 3 to 7 used. But in this dataset, most of the co-accessibility structure Cicero detects is carried by peaks the Cell Ranger consensus already had.

Step 20: Viewing Co-Accessibility Links at CEBPB

Signac can store co-accessibility links in the assay and draw them under a coverage plot. ConnectionsToLinks() converts Cicero’s table into Signac’s link format, keeping the network labels from ccans.

#-----------------------------------------------
# STEP 20: Add the links to the object and plot them at CEBPB
#-----------------------------------------------

DefaultAssay(mono) <- "macs3"
Links(mono) <- ConnectionsToLinks(conns = conns, ccans = ccans, threshold = 0.25)

p_links <- CoveragePlot(
  object            = mono,
  region            = "CEBPB",
  group.by          = "condition",
  ranges            = cr_peaks,
  ranges.title      = "Cell Ranger",
  extend.upstream   = 150000,
  extend.downstream = 150000
)

ggsave(file.path(plot_dir, "07_cebpb_coaccessibility.png"), p_links,
       width = 10, height = 7, dpi = 300)

Reading this figure: from top to bottom — monocyte accessibility split into Healthy and Severe tracks, the gene models with CEBPB in the middle, the MACS3 peaks (Peaks), the Part 3 consensus (Cell Ranger), and the Links track. Each arc joins the midpoints of two co-accessible peaks; its colour shows the co-accessibility score. Only links above 0.25 are drawn, because we passed that threshold to ConnectionsToLinks().

The window is CEBPB plus 150 kb on each side. On the left are UBE2V1, PEDS1 and the read-through PEDS1-UBE2V1; CEBPB is in the middle.

Five arcs are drawn, with scores between 0.25 and 0.29.

The Healthy and Severe tracks show accessibility only. The links were computed from all monocytes together, so the figure cannot tell you whether a link is stronger in one condition.

The track labels “Peaks” and “Cell Ranger” overlap slightly on the left. This is cosmetic and does not affect the data shown.

💾 Saving Results for Reuse and Supplementary Material

#-----------------------------------------------
# STEP 21: Save tables, peaks and session information
#-----------------------------------------------

unified_df <- as.data.frame(unified)[, c("seqnames", "start", "end", "width",
                                         "peak_called_in", "n_cell_types", "in_cellranger")]
write.csv(unified_df, file.path(table_dir, "unified_macs3_peaks.csv"), row.names = FALSE)
export(unified, file.path(recall_dir, "unified_macs3_peaks.bed"), format = "bed")

write.csv(frip_by_type, file.path(table_dir, "frip_by_celltype.csv"), row.names = FALSE)
write.csv(strong, file.path(table_dir, "cicero_strong_links_cd14mono_chr20.csv"), row.names = FALSE)

writeLines(capture.output(sessionInfo()), file.path(recall_dir, "sessionInfo_part8.txt"))

Together with the objects saved along the way, these files let anyone reproduce every number and figure without rerunning MACS3 or recounting fragments:

Using the new peaks downstream. The macs3 assay is a complete feature set: you can run TF-IDF, LSI, differential accessibility or motif analysis on it exactly as Parts 3 to 7 did on peaks. Because the features are different, recompute everything that depends on them — TF-IDF and LSI, motif matches, gene activity — rather than mixing results from the two assays.

✅ Best Practices for Cell-Type-Specific Peak Calling and Co-Accessibility

1. Recall peaks after annotation, not instead of it. Aggregate peaks are the right features for QC, clustering and a first annotation, because you need cell groups before you can call per-group peaks. Per-cell-type calling is the second pass that the first pass makes possible.

2. Call peaks in cell types, not in quality-control groups. Select the cell types you will analyze before peak calling (Step 4). A pseudobulk of damaged nuclei calls peaks that describe background, and those peaks then enter every downstream analysis.

3. Check the strength of peaks from small cell types before using them. A small pseudobulk can call many peaks, but they are weaker. Here 87 dendritic cells called 30,962 peaks, more than 1,585 CD16 monocytes did, with a median summit -log10 q-value of 4.95 against 16.6 to 69.0 for every other cell type, and 40 percent at the 200 bp minimum width. Summarize width, fold enrichment and q-value per cell type (Step 5). Consider a stricter q-value for small groups, or report their peaks separately, rather than giving them equal weight.

4. Remember that each cell-type pseudobulk pools donors. Pooling gives rare types the signal they need, but a donor with many cells of a type shapes its peak set more than the others. Before calling a peak “cell-type-specific”, check that more than one donor supports it; CoveragePlot(..., split.by = "sample_id") shows this directly.

5. Keep the per-cell-type peak sets. combine.peaks = FALSE costs one line of merging code and preserves each peak’s MACS3 statistics. You can then rebuild the unified set with a different merging rule, or study one cell type’s peaks on their own, without rerunning MACS3.

6. Pass the MACS3 path explicitly. CallPeaks() looks for macs2 unless told otherwise. macs2.path = Sys.which("macs3") makes the choice visible in your code and in your methods section.

7. Think about duplicate handling. In BED mode MACS3 keeps one tag per position and strand by default. Fragments from different nuclei in a pseudobulk are not PCR duplicates of each other; --keep-dup all matches MACS3’s own FRAG-mode behaviour and Part 1.

8. Put temporary files on storage that can hold them. The split fragment files can approach the uncompressed size of all your fragments files. Set fragment.tempdir to project storage rather than a compute node’s local temporary directory.

9. Filter the new peak set exactly as you filtered the old one. Same width bounds, same chromosomes, same blacklist. Otherwise a difference between the two sets may be a difference in filters.

10. Check the merge for chaining. reduce() merges any overlapping peaks, so look at the width distribution after merging. If long merged regions appear, consider ArchR-style iterative overlap merging with fixed-width peaks (Corces et al. 2018; Granja et al. 2021).

11. Requantify from the fragments. New peaks need new counts. Never reuse a count matrix built on one peak set for another; FeatureMatrix() on the fragments files is the only correct source.

12. Compare genome coverage, not just peak counts. A FRiP increase can come entirely from covering more of the genome. Report the total peak width of each set next to any FRiP comparison.

13. Run co-accessibility within a cell type. Across mixed cell types, lineage differences dominate the covariation and almost every lineage-specific element looks linked to every other. Within one cell type, links reflect variation between that cell type’s states.

14. Scale Cicero by chromosome. Signac’s vignette runs one chromosome to save time and notes that the restriction can be dropped to run genome-wide. On a cluster, a practical route is one job per chromosome — each with its own genome_df row — combining the saved conns tables afterwards. Cicero’s windows are 500 kb, so no link crosses a chromosome boundary and nothing is lost by splitting.

15. Report the co-accessibility threshold, and test its sensitivity. 0.25 is a convention from the Cicero vignette’s plotting example, not a statistical cutoff. Check that conclusions such as “new peaks take part in links” hold at 0.2 and 0.3 as well.

16. Treat links as hypotheses. The Cicero authors describe their connections as putative and needing experimental validation. Use links to prioritize regions for follow-up, and look for independent support — chromatin conformation data, eQTLs or perturbation experiments — before claiming that an element regulates a gene. And expect weak peaks to be under-represented: in this run, new monocyte peaks were 24 percent of the peaks on chromosome 20 but appeared in only 2.6 percent of strong links, and their MACS3 calls were far weaker (Step 19).

17. Record versions and parameters. MACS3 version and arguments, the cell types used, the merge rule, the filters and the Cicero settings all shape the result. Save sessionInfo() (Step 21) and keep the CallPeaks call in your methods.

⚠️ Common Pitfalls and How to Avoid Them

🎯 Key Takeaways

• Peak calling is a choice that follows the analysis. Aggregate peaks get you to cell types; per-cell-type peaks are what you analyze cell types with.
• Rare populations pay the highest price for aggregate calling, and the data show it. Nine in ten dendritic-cell private peaks were absent from the Cell Ranger consensus.
• The rarest populations also give the weakest peaks. The same 87 cells that produced the most novel peaks produced the least significant ones. Report both facts together, and check peak strength per cell type before relying on a small group’s peaks.
• Narrower peaks hold the signal more densely. Less than half the genome captured almost the same fraction of fragments, which is the practical case for MACS3 peaks in motif and footprint work.
• New features mean new counts. Every peak set needs its own requantification from the fragments, as Part 3 taught.
• Compare like with like. Same filters, same merge rule, genome coverage alongside peak counts, overlap counted from both sides.
• Co-accessibility links regions, not genes, and predicts rather than measures. It is the strongest RNA-free way to connect distal elements to promoters. In this dataset most strong links were between peaks the Cell Ranger set already had, and weak new peaks rarely reached a strong link.

📚 References

1. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS. Model-based Analysis of ChIP-Seq (MACS). Genome Biology. 2008;9(9):R137. doi:10.1186/gb-2008-9-9-r137
2. MACS3 Project. MACS3 Documentation. https://macs3-project.github.io/MACS/
3. Stuart T, Srivastava A, Madad S, Lareau CA, Satija R. Single-cell chromatin state analysis with Signac. Nature Methods. 2021;18(11):1333-1341. doi:10.1038/s41592-021-01282-5
4. Signac documentation: CallPeaks() reference, peak calling vignette and Cicero vignette. https://stuartlab.org/signac/articles/peak_calling and https://stuartlab.org/signac/articles/cicero
5. Pliner HA, Packer JS, McFaline-Figueroa JL, Cusanovich DA, Daza RM, Aghamirzaie D, Srivatsan S, Qiu X, Jackson D, Minkina A, Adey AC, Steemers FJ, Shendure J, Trapnell C. Cicero predicts cis-regulatory DNA interactions from single-cell chromatin accessibility data. Molecular Cell. 2018;71(5):858-871.e8. doi:10.1016/j.molcel.2018.06.044
6. Cicero Bioconductor package (release 1.28). https://bioconductor.org/packages/release/bioc/html/cicero.html
7. Friedman J, Hastie T, Tibshirani R. Sparse inverse covariance estimation with the graphical lasso. Biostatistics. 2008;9(3):432-441. doi:10.1093/biostatistics/kxm045
8. Granja JM, Corces MR, Pierce SE, Bagdatli ST, Choudhry H, Chang HY, Greenleaf WJ. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nature Genetics. 2021;53(3):403-411. doi:10.1038/s41588-021-00790-6
9. ArchR documentation: The Iterative Overlap Peak Merging Procedure. https://www.archrproject.com/bookdown/the-iterative-overlap-peak-merging-procedure.html
10. Corces MR, Granja JM, Shams S, et al. The chromatin accessibility landscape of primary human cancers. Science. 2018;362(6413):eaav1898. doi:10.1126/science.aav1898
11. Amemiya HM, Kundaje A, Boyle AP. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Scientific Reports. 2019;9(1):9354. doi:10.1038/s41598-019-45839-z
12. Lawrence M, Huber W, Pages H, Aboyoun P, Carlson M, Gentleman R, Morgan MT, Carey VJ. Software for computing and annotating genomic ranges. PLoS Computational Biology. 2013;9(8):e1003118. doi:10.1371/journal.pcbi.1003118
13. deepTools documentation: Effective Genome Size. https://deeptools.readthedocs.io/en/develop/content/feature/effectiveGenomeSize.html
14. 10x Genomics. Cell Ranger ATAC documentation. https://www.10xgenomics.com/support/software/cell-ranger-atac/latest
15. Remodeling of the chromatin landscape in peripheral blood cells in patients with severe Delta COVID-19. Frontiers in Immunology. 2024;15:1415317. doi:10.3389/fimmu.2024.1415317 (scATAC-seq data: NCBI GEO GSE282769)

This tutorial is part of the comprehensive NGS101.com single-cell ATAC-seq analysis series for beginners.