How to Analyze Single-Cell ATAC-seq Data — A Complete Beginner’s Guide Part 6: Footprinting Analysis

Go beyond “which motifs are enriched” and ask “is the transcription factor actually sitting there?” — run scATAC-seq footprinting two ways, with Signac Footprint() on fragments and TOBIAS on cell-type pseudobulk BAMs, and learn how far each answer can be trusted.

Part 5 ended with a list of transcription factor families. In CD14 monocytes, the 46 peaks that opened in severe COVID-19 were enriched for the RFX X-box motif, and the 50 peaks that closed were enriched for AP-1 motifs (FOS, JUN, FOSL2, BATF and relatives). That is a statement about DNA sequence: those motifs are over-represented in a set of regions. It says nothing about whether any protein is bound to them in these cells.

scATAC-seq footprinting analysis asks the next question. Instead of counting motifs in peaks, it goes back to the individual Tn5 insertion events around every instance of a motif and looks at their shape. A bound transcription factor changes that shape: the chromatin around its site is kept open, and the site itself can be partially shielded from Tn5. Footprinting reads that pattern, aggregated over thousands of motif sites, and turns it into evidence about transcription factor occupancy.

This tutorial runs footprinting on the same CD14 monocytes, severe versus healthy, with two tools that take very different routes to the answer. Signac’s Footprint() works directly on the fragments file inside your Seurat object. TOBIAS, the tool from our bulk ATAC-seq footprinting tutorial, needs a BAM file, so we build one per condition from the monocyte barcodes. Then we put the two answers side by side, next to the Part 5 motif enrichment results, and ask whether three independent lines of evidence tell the same story.

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, Part 4: Cell Type Identification and Part 5: Differential Accessibility Analysis. You need the atac_annotated.rds object from Part 4, the Part 5 motif tables, the Cell Ranger ATAC outs/ folders from Part 1 (the TOBIAS half uses possorted_bam.bam), the same Pixi workspace, and an HPC account for the TOBIAS steps. No new raw data is downloaded.

🧭 Introduction: What Footprinting Adds to scATAC-seq Analysis

What Is Footprinting Analysis and When Do You Need It?

Tn5 transposase cuts and tags DNA where chromatin is open. When a transcription factor binds DNA, two things can happen to the Tn5 cuts around its site:

1. The flanks open up. Many transcription factors displace or reposition nearby nucleosomes, so the DNA on either side of the motif becomes more accessible. Tn5 cuts there more often.
2. The motif itself can be protected. A protein sitting on the motif physically blocks Tn5 from those few bases, producing a dip in cuts exactly over the motif.

Averaged over thousands of motif sites, those two effects produce the classic “peak-dip-peak” footprint profile. Buenrostro et al. (2013) showed it for CTCF in the original ATAC-seq paper: a clear notch in ATAC-seq signal at the CTCF motif, aggregated over the genome-wide CTCF motif sites.

Beginner’s Tip: Think of motif enrichment as reading the address book — which addresses exist in a neighborhood. Footprinting is looking at the doorsteps — which houses show signs that someone is home. The address book is the same in every cell; the doorsteps are not.

When do you need footprinting? Three situations come up again and again:

• You have a motif family and want evidence for binding. Part 5 gave you “AP-1 motifs are enriched in regions that close in severe disease.” Footprinting tests whether AP-1 sites genome-wide in monocytes show less binding-associated signal in severe patients.
• You want to compare transcription factor activity between conditions or cell types using all motif sites, not only the handful inside differentially accessible peaks. Part 5’s motif test used 46 and 50 regions; footprinting uses every motif instance in accessible chromatin.
• You want candidate binding sites for follow-up. TOBIAS classifies individual motif sites as bound or unbound, which gives you site-level hypotheses to test with ChIP-seq, CUT&RUN or reporter assays.

Motif Enrichment and Footprinting: Two Different Questions

These two analyses are easily confused, because they start from the same motif database and produce lists of the same transcription factor names.

The difference that matters: motif enrichment is blind to your cells. Give FindMotifs() the same peak set from a different experiment and it returns the same answer, because it only ever looks at DNA sequence. The cells enter the analysis solely through which peaks you chose to hand it — in Part 5, the 46 peaks that opened and the 50 that closed.

Footprinting is computed from the Tn5 insertions of the cells you choose. At a single identical motif site, severe and healthy monocytes can give different answers. That is what makes it an independent test of a motif hypothesis rather than a restatement of it.

Why Footprinting Is Harder Than It Looks

Before running anything, you need three facts that shape every interpretation later in this tutorial.

Tn5 has its own sequence preference. Tn5 does not cut DNA uniformly; its insertion rate depends on the local sequence. A motif’s own sequence therefore creates a characteristic cut pattern even on naked DNA. Sung et al. (2014) showed, for DNase-seq, that the cleavage profile inside a footprint originates largely from the DNA sequence of the binding site rather than from the protein. Every footprinting method therefore estimates an expected cut pattern from sequence alone and compares the observed pattern against it.

Many transcription factors never leave a measurable footprint. Sung et al. (2014) found that factors with short DNA residence times show no footprint at bound motif elements. Baek et al. (2017) reported that 80 percent of transcription factor motifs showed no measurable footprint and proposed measuring motif-flanking accessibility alongside footprint depth for exactly this reason. A missing dip does not mean a missing factor.

In scATAC-seq the central dip is the least reliable part of the plot. The Signac maintainer, asked how to read Signac footprint plots, answered that the height of the motif-flanking region is what to look at, and wrote: “I would not really interpret the dip in the center to mean anything” (Signac GitHub Discussion #968). His reasoning: scATAC-seq data are sparse, and Tn5 insertion bias is strong, so the simplest explanation for a central dip is bias. This tutorial follows that advice. For Signac we quantify the flank; for TOBIAS we use its footprint score, which combines depletion with flank accessibility.

How scATAC-seq Footprinting Compares to Bulk ATAC-seq Footprinting

If you have worked through bulk ATAC-seq Part 3: Footprinting Analysis, the TOBIAS commands in this tutorial will look familiar. What changes is how you get to a BAM file, and how much signal that BAM contains.

The row that should worry you is signal per condition. In Part 5, Step 8, the CD14 monocyte pseudobulk held 16.8 million peak counts for the two healthy donors combined (4,838,159 + 11,979,419) and 9.4 million for the two severe donors (7,019,769 + 2,363,335). That is a fraction of one cell type from a fraction of each library. Single-cell footprinting always works with less signal than bulk footprinting of a sorted population sequenced to the same depth, and that is why this tutorial footprints the largest well-annotated population, aggregates many motif sites, and checks every result donor by donor.

Signac Footprint() Versus TOBIAS: Two Different Approaches

Both tools correct for Tn5 bias and aggregate over motif sites, but almost everything else differs.

They are complementary rather than competing. TOBIAS screens; Signac lets you check a small number of motifs donor by donor without rebuilding anything.

The Part 6 Workflow

atac_annotated.rds (Part 4)
        |
        v
[1] Subset CD14 monocytes, confirm motif positions and fragment files
        |
        +------------------------------+-------------------------------+
        v                                                              v
[2] Signac Footprint()                                   [3] TOBIAS on pseudobulk BAMs
    5 chosen motifs, fragments-based                         export barcodes, peaks, motifs (R)
    PlotFootprint by condition and by donor                  samtools: select cells by CB tag
    quantify motif-flank enrichment per donor                samtools: merge per condition, filter
        |                                                    ATACorrect -> ScoreBigwig -> BINDetect
        |                                                              |
        +------------------------------+-------------------------------+
                                       v
[4] Compare: Signac flank change vs TOBIAS differential score vs Part 5 motif enrichment
                                       |
                                       v
[5] Visualization: footprint profiles, donor-level flank plot,
    BINDetect volcano, method concordance

Which Steps Run Where

This tutorial moves between R and the Linux command line, and the two halves can be days apart because the TOBIAS stage usually runs as a batch job. The work is therefore arranged as three uninterrupted stretches, with everything the next stretch needs written to disk at each boundary.

You can close R after Step 11 and close the terminal after Step 18. Step 19 reloads the saved objects in two lines, so nothing before it has to be repeated.

⚙️ Setting Up the Analysis Environment

What Is New in Part 6

Nothing new is needed in R. Footprint(), PlotFootprint() and GetFootprintData() are part of Signac, and BSgenome.Hsapiens.UCSC.hg38 was installed and repaired in Part 4 (if you skipped that repair, see the Bioconda stub fix in Part 4 — library(BSgenome.Hsapiens.UCSC.hg38) failing with “there is no package called” is that problem).

Two command-line tools are new, and only for the TOBIAS half:

Adding a Dedicated TOBIAS Environment to the Same Pixi Workspace

TOBIAS is a Python command-line program and your analysis environment is a large R stack: Signac, Seurat and several hundred Bioconductor packages. Mixing them is possible, but it makes every future pixi add in that project solve a bigger, more constrained problem, and a failed solve then appears in the middle of the environment your R work depends on.

Pixi handles this cleanly with one workspace and several environments. We add a footprint feature and build an environment from that feature alone (--no-default-feature), so it neither inherits nor disturbs the R stack. Both environments share one pixi.toml and one pixi.lock, so the project stays reproducible as a whole.

#-----------------------------------------------
# STEP 1: Add a separate footprinting environment to the Pixi workspace
#-----------------------------------------------

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

# Put the tools in a feature called "footprint"
pixi add --feature footprint tobias samtools

# Create an environment named "footprint" that contains ONLY this feature
pixi workspace environment add footprint --feature footprint --no-default-feature

# Install it and confirm the tools resolve
pixi install -e footprint
pixi run -e footprint TOBIAS --version
pixi run -e footprint samtools --version

Your default environment (R, Signac, Seurat) is untouched: pixi run R works exactly as before. For the TOBIAS steps, enter the new environment with pixi shell -e footprint.

A warning you can ignore on a cluster. If your home directory is on shared storage, Pixi prints WARN cache for Repodata ... is on a network/parallel filesystem ... redirected to /tmp/.... That is not an error: Pixi moves its package cache to local disk for that run and continues. Set PIXI_CACHE_DIR to a local path in your shell profile if you would rather keep one cache across runs.

Loading Libraries and Creating the Part 6 Directory Structure

Launch R in the default environment:

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

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

library(Signac)
library(Seurat)
library(BSgenome.Hsapiens.UCSC.hg38)  # genome sequence for Tn5 bias and expected insertions
library(GenomicRanges)                # seqnames(), start(), end() for the peak BED export
library(ggplot2)
library(ggrepel)
library(patchwork)

# Same seed as every other part of the series
set.seed(1234)

# Footprinting builds large per-cell insertion matrices
options(future.globals.maxSize = 60 * 1024^3)

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

base_dir   <- "/projects/mylab/shared/scatac_tutorial"
annot_dir  <- file.path(base_dir, "annotation")    # Part 4 output
da_dir     <- file.path(base_dir, "differential")  # Part 5 output
fp_dir     <- file.path(base_dir, "footprinting")  # this tutorial

obj_dir    <- file.path(fp_dir, "objects")
plot_dir   <- file.path(fp_dir, "plots")
table_dir  <- file.path(fp_dir, "tables")
tobias_dir <- file.path(fp_dir, "tobias")

# TOBIAS gets one subfolder per stage so intermediate files never mix
tobias_sub <- file.path(tobias_dir,
                        c("barcodes", "bam", "atacorrect",
                          "footprints", "bindetect", "aggregate"))

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

Your project tree now looks like this:

/projects/mylab/shared/scatac_tutorial/
├── results/               # Cell Ranger ATAC output (Part 1) -- BAMs and fragments
├── downstream_qc/         # per-sample QC objects (Part 2)
├── integrated/            # consensus peaks, integrated object (Part 3)
├── annotation/            # annotated object (Part 4)
├── differential/          # DA results and motif tables (Part 5)
└── footprinting/          # this tutorial
    ├── objects/
    ├── plots/
    ├── tables/
    └── tobias/
        ├── barcodes/      # one barcode list per donor
        ├── bam/           # CD14 monocyte pseudobulk BAMs
        ├── atacorrect/    # bias-corrected cut-site tracks
        ├── footprints/    # footprint score tracks
        ├── bindetect/     # per-motif binding results
        └── aggregate/     # TOBIAS aggregate plots

📦 Example Data: The Annotated Object from Part 4

Loading the Object and Selecting CD14 Monocytes

We footprint the same cell type as Part 5, for the same reasons: it is the largest well-annotated population, it has useful numbers in both conditions, and Part 5 already gave us motif-family hypotheses for it.

#-----------------------------------------------
# STEP 4: Load the Part 4 object and subset CD14 monocytes
#-----------------------------------------------

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

# subset() keeps the peaks assay, the motif object and the fragment file links
mono <- subset(atac, subset = cell_type == "CD14 monocyte")
Idents(mono) <- "condition"

# Free memory: the full object is no longer needed
rm(atac)
invisible(gc())

table(mono$sample_id, mono$condition)

Output:

         Healthy Severe
healthy1     807      0
healthy2    2575      0
severe1        0   1323
severe2        0    508

Reading this table. The same design as Part 5, and the same warning. There are 5,213 cells but only four donors, and the donors are unbalanced: healthy2 contributes 2,575 of the 3,382 healthy cells and severe2 only 508. Any footprint computed on “Healthy” is dominated by one donor. That is why every result below is also computed per donor.

Confirming the Motif Positions and Fragment Files Are in Place

Footprint() needs two things that you did not create in this tutorial: the exact genomic positions of each motif match, and readable fragments files. Both should already be there. Check before spending compute.

#-----------------------------------------------
# STEP 5: Confirm motif positions and fragment files
#-----------------------------------------------

# AddMotifs() in Part 4 stored one GRanges of motif matches per motif.
# These are matches inside the consensus peaks, found by motifmatchr.
motif_pos <- GetMotifData(mono, slot = "positions")
length(motif_pos)

# Every fragments file the object points to must still exist on disk
sapply(Fragments(mono), function(f) file.exists(GetFragmentData(f, slot = "path")))

Output:

[1] 633

[1] TRUE TRUE TRUE TRUE

Reading this output. 633 motifs is the same JASPAR2020 CORE vertebrate set that FindMotifs() tested in Parts 4 and 5, so the footprinting here and the enrichment there are drawn from an identical motif universe. Four TRUE values mean all four fragments files are still where Part 2 left them.

If length(motif_pos) fails or returns 0, the motif object was built without positions and Footprint() will stop with “Motif positions not present in Motif object.” Re-run AddMotifs() from Part 4 on this object. If any fragments check is FALSE, the Cell Ranger output has moved since Part 2. Point the object at the new location with Signac’s UpdatePath() before continuing.

Choosing Which Motifs to Footprint

Signac footprints only the motifs you name, and each one costs real compute. Choose them for a reason. We use five, each with a job:

Two reminders from Parts 4 and 5 apply fully here. First, motifs identify families, not proteins: FOSL2 stands for the AP-1 site that FOS, JUN, FOSL1, FOSL2, BATF and others all bind, and RFX2 stands for the X-box bound by the whole RFX family. Footprinting cannot tell family members apart either, because they bind the same DNA. Second, the Part 5 RFX result rested on only 46 regions; footprinting is a genuinely independent test because it uses every RFX2 motif site in accessible chromatin, not just those 46.

#-----------------------------------------------
# STEP 6: Record the chosen motifs, their IDs and motif widths
#-----------------------------------------------

fp_motifs <- c("CTCF", "SPI1", "CEBPB", "FOSL2", "RFX2")

# Motif positions are stored under JASPAR IDs; footprints are stored under names
fp_ids <- ConvertMotifID(mono, name = fp_motifs)

# Motif width is needed later to define the flanking window
motif_width <- setNames(sapply(motif_pos[fp_ids], function(g) width(g)[1]), fp_motifs)

data.frame(motif = fp_motifs, id = fp_ids,
           width = motif_width, n_sites = lengths(motif_pos[fp_ids]))

Output:

      motif       id width n_sites
CTCF   CTCF MA0139.1    19   55393
SPI1   SPI1 MA0080.5    20  117781
CEBPB CEBPB MA0466.2    10    7127
FOSL2 FOSL2 MA0478.1    11   39118
RFX2   RFX2 MA0600.2    16   15832

Reading this table, because it predicts how noisy each footprint will be. n_sites is how many motif matches the profile will be averaged over, and it spans a 17-fold range here. SPI1 has 117,781 sites and will give the smoothest curve; CEBPB has 7,127 and will be the jumpiest. When you look at the profiles in Step 22, expect the visual noise to track this column, and do not mistake a jumpy CEBPB line for a biological difference.

Width matters for the flank window. Step 8 defines the flank as the region from the motif edge out to 50 bp beyond it, so a 10 bp motif (CEBPB) and a 20 bp motif (SPI1) measure slightly different windows. That is intentional — it is how PlotFootprint() defines the flank internally — but it means flank values are comparable between conditions for the same motif, not between motifs.

🔬 Footprinting Analysis

Method 1: Signac Footprint() on the Fragments File

Footprint() does four things for each motif:

1. Extends every motif match by 250 bp on each side (upstream and downstream defaults).
2. Counts Tn5 insertions at each base of that window, per cell, from the fragments file.
3. Estimates Tn5 insertion bias once per object — the frequency of each hexamer at insertion sites on chr1 divided by its frequency in the chr1 sequence — and uses it to compute the expected insertion profile from the sequence of the motif windows.
4. Stores the per-cell matrix and the expected profile inside the peaks assay.

Because the per-cell matrix is stored, grouping happens later: you run Footprint() once and can plot it by condition, by donor, or by anything else in the metadata.

#-----------------------------------------------
# STEP 7: Compute footprints for the five motifs
#-----------------------------------------------

# in.peaks = TRUE keeps only motif matches inside peaks of this assay.
# AddMotifs() already scanned only the peaks, so this changes nothing here --
# it is set explicitly so the code stays correct if the motif object was
# ever built on a wider set of regions.
mono <- Footprint(
  object     = mono,
  motif.name = fp_motifs,
  genome     = BSgenome.Hsapiens.UCSC.hg38,
  in.peaks   = TRUE
)

Expect this to take a while. The first call reads every chr1 fragment for all 5,213 cells to estimate the insertion bias, and each motif then reads the fragments around all of its sites. Run it on a compute node, and save the object in Step 11 so you never repeat it.

Turning the Profile Into a Number: Motif-Flank Enrichment

PlotFootprint() draws the profiles, but a picture is hard to compare across donors. Following the Signac maintainer’s advice, we extract the height of the motif flank, computed exactly the way PlotFootprint() computes it internally:

• GetFootprintData() returns, for each group, the mean insertion count at each position, divided by the mean of the outermost 50 bp on each side. The same scaling is applied to the expected profile.
• PlotFootprint() then subtracts expected from observed (its default normalization = "subtract").
• It defines the flank as the region from the edge of the motif out to 50 bp beyond it, on both sides.

So a flank value of 0 means “no more Tn5 insertion next to the motif than the sequence alone predicts, relative to 200 bp away.” Positive values mean the region around the motif is more open than its surroundings.

#-----------------------------------------------
# STEP 8: Quantify motif-flank enrichment per group
#-----------------------------------------------

# Reproduces the flank calculation inside PlotFootprint():
#   observed minus expected, averaged from the motif edge out to +50 bp
flank_enrichment <- function(object, group.by) {
  fp   <- GetFootprintData(object, features = fp_motifs, group.by = group.by)
  obs  <- fp[fp$class == "Observed", ]
  expd <- fp[fp$class == "Expected", ]

  # Subtract the sequence-expected Tn5 profile at the same position
  obs$corrected <- obs$norm.value -
    expd$norm.value[match(paste(obs$feature, obs$position),
                          paste(expd$feature, expd$position))]

  # Flank window: beyond half the motif width, within 50 bp of the motif edge
  half     <- ceiling(motif_width[obs$feature] / 2)
  in_flank <- abs(obs$position) > half & abs(obs$position) < half + 50

  aggregate(corrected ~ feature + group, data = obs[in_flank, ], FUN = mean)
}

# Pooled by condition, and separately for each donor
flank_cond  <- flank_enrichment(mono, group.by = "condition")
flank_donor <- flank_enrichment(mono, group.by = "sample_id")
flank_donor$condition <- ifelse(grepl("^severe", flank_donor$group), "Severe", "Healthy")

# One row per motif: Healthy, Severe and the difference
flank_wide <- reshape(flank_cond, idvar = "feature", timevar = "group",
                      direction = "wide")
flank_wide$signac_delta <- flank_wide$corrected.Severe - flank_wide$corrected.Healthy
flank_wide

Output:

  feature corrected.Healthy corrected.Severe signac_delta
1   CEBPB         1.0926130        1.3896869   0.29707392
2    CTCF         1.5205405        1.3572774  -0.16326301
3   FOSL2         2.4527200        2.3792277  -0.07349232
4    RFX2         0.9129154        0.9608078   0.04789249
5    SPI1         0.8613500        0.9003179   0.03896793

Reading this table, starting with the column you did not come for. Look at the two condition columns before the difference. Every value is well above zero, which is the reassuring part: at all five motifs, the 50 bp flanking the motif carries more Tn5 insertion than the sequence alone predicts. The motifs are sitting in genuinely open, actively transposed chromatin.

The ranking between motifs is not a ranking of importance. FOSL2 has the highest flank enrichment (2.45) and SPI1 the lowest (0.86), yet SPI1 is the myeloid master regulator and FOSL2 is one family among many. Flank enrichment measures how open the neighbourhood of a motif is relative to 200 bp away, so it rewards motifs that sit in small, sharply defined accessible islands and penalizes motifs spread through broad open domains. Compare a motif to itself across conditions; do not read the column top to bottom.

Now the difference. CEBPB moves the most by a wide margin (+0.297), and it moves up in severe disease. CTCF is second (-0.163) — and CTCF is the control, which should give you pause rather than satisfaction. FOSL2 (-0.073) and RFX2 (+0.048) move in the directions Part 5 predicted, but by a third to a fifth as much as CEBPB. SPI1 barely moves (+0.039), which is what you want from a lineage factor that both groups still are.

Why the normalization matters here. Each group is scaled to its own outer flanks before anything is compared, which removes the effect of sequencing depth and cell number on overall signal height — necessary when severe2 contributes 508 cells and healthy2 contributes 2,575. What it does not remove is noise: a group built from fewer insertions still gives a jumpier profile, and a jumpier profile can shift a windowed average.

Testing Across Donors

The pooled comparison in Step 8 has no error bars, because it collapses each condition into a single profile. The four donor values are the data that can carry a test, and the donor is the correct unit: cells within a donor share that person’s genetics, batch and clinical course, exactly as in Part 5.

So we do the ordinary thing and run a t-test on them.

#-----------------------------------------------
# STEP 9: Test the flank enrichment across donors
#-----------------------------------------------

donor_test <- do.call(rbind, lapply(split(flank_donor, flank_donor$feature), function(d) {
  h <- d$corrected[d$condition == "Healthy"]
  s <- d$corrected[d$condition == "Severe"]

  # Two donors per group, so Student's t with equal variance: df = 2.
  tt <- t.test(s, h, var.equal = TRUE)

  data.frame(feature = d$feature[1],
             mean_H  = mean(h),
             mean_S  = mean(s),
             diff    = mean(s) - mean(h),
             pct     = 100 * (mean(s) - mean(h)) / mean(h),
             t       = unname(tt$statistic),
             p       = tt$p.value)
}))

# Five motifs tested, so correct across them
donor_test$p_BH <- p.adjust(donor_test$p, method = "BH")

donor_test[order(donor_test$p), ]

Output:

      feature   mean_H   mean_S     diff    pct       t       p    p_BH
CEBPB   CEBPB   1.0861   1.3744   0.2883  26.55   8.734 0.01286 0.06429
SPI1     SPI1   0.8657   0.8961   0.0304   3.51   2.276 0.15061 0.37652
CTCF     CTCF   1.5071   1.3899  -0.1172  -7.78  -1.601 0.25048 0.41747
FOSL2   FOSL2   2.4321   2.3816  -0.0505  -2.08  -1.047 0.40484 0.50604
RFX2     RFX2   0.9065   0.9204   0.0140   1.54   0.176 0.87629 0.87629

Reading this table. One motif clears p < 0.05 before correction: CEBPB, at 0.013, with a t of 8.7. After BH correction across five motifs it sits at 0.064 — above the line, and you should say so rather than quoting the raw value. Nothing else comes close; the next smallest raw p is 0.15.

The controls behave. CTCF at p = 0.25 and SPI1 at 0.15 are where control motifs belong. That matters because CTCF has the second largest effect in the table (-7.8 percent): the t-test correctly declines to call it, because its two severe donors differ from each other by more than the groups differ.

The Part 5 hypotheses do not survive. FOSL2 (p = 0.40) and RFX2 (p = 0.88) move in the predicted directions and fail the test. RFX2 is the clearest failure in the set — a 1.5 percent difference against donors that interleave completely.

Use the pct column when comparing motifs to each other. The baselines span 0.87 (SPI1) to 2.43 (FOSL2), nearly threefold, so a raw difference of 0.03 means something different at each end. CEBPB’s +26.5 percent is the number to put in a sentence; its +0.288 is not comparable to FOSL2’s -0.051.

Why diff here is not the signac_delta from Step 8. CEBPB reads +0.2883 above and +0.2971 in Step 8; RFX2 reads +0.0140 here against +0.0479 there, more than three times smaller. Step 8 pools all cells into one profile per condition, so it is cell-weighted: healthy2 contributes 76 percent of the healthy cells and severe1 72 percent of the severe ones, and those two donors effectively set the pooled values. The donor means here weight each person equally. For RFX2 that matters a great deal, because severe1 is both the largest severe donor and its highest RFX2 value.

This is the same pseudoreplication that Part 5 spent a whole section on, wearing a different hat. Report the donor-weighted difference, and treat the pooled profile as a picture rather than a measurement.

What a t-test can and cannot do at n = 2 per group. You need |t| > 4.303 for p < 0.05 at df = 2, so the test is close to all-or-nothing: CEBPB clears it at 8.7 and nothing else exceeds 2.3. Each standard deviation rests on a single degree of freedom, which makes them unstable — look at FOSL2, where sd is 0.068 in healthy and 0.0069 in severe, a tenfold difference that is itself noise. Had the healthy pair happened to land as tightly as the severe pair, a trivial difference would have come out significant.

A rank-based test cannot help here. There are only six ways to label four donors as two-and-two, so the smallest two-sided p a permutation or Wilcoxon test can return is 2/6 = 0.33. At this sample size a distribution-free test has no power at all, which is why the t-test’s normality assumption is worth making. It is the same trade Part 5 made when it handed DESeq2 a 2-versus-2 design.

If you footprint many motifs, moderate the variance. With one degree of freedom per motif the variance estimates are barely estimates. limma::eBayes shrinks each motif’s variance toward a trend fitted across all motifs — the same borrowing-strength logic DESeq2 uses for dispersion — and is the standard treatment for this shape of data. With five motifs there is almost nothing to borrow from, so a plain t-test is the honest choice here.

Where the Test Came From, and Why Signac Does Not Provide One

It is worth being explicit about something you may have noticed: Signac has no differential footprinting function. Footprint() computes the profiles, PlotFootprint() draws them and GetFootprintData() returns the numbers. There is no test, no p-value and no recommended statistic anywhere in the footprinting API.

That is not an oversight, and two design decisions explain it. Footprint() sums insertions across every instance of a motif as it builds the pileup, so the stored matrix is cells by positions — there is no per-site dimension left to test. And grouping happens at plot time through group.by, so the object never knows what your experimental design is. The package gives you a profile and leaves the statistics to you.

So the flank statistic in Step 8 and the t-test in Step 9 are a choice this tutorial makes, not something Signac prescribes. The flank measure follows the Signac maintainer’s advice in GitHub Discussion #968 that the height of the motif flank is the informative part; the t-test is the ordinary way to compare a per-replicate summary across two groups. Anyone reproducing this is entitled to make a different choice, and should say which one they made.

The Two Methods Measure Significance of Completely Different Things

TOBIAS does report a p-value, and it is tempting to line it up beside the one above. Do not. They are not the same kind of number.

The row that matters is the last one. Signac’s smallest p is 0.013; TOBIAS’s largest p among our five motifs is 2.7e-72. That gap is not a disagreement about biology. TOBIAS’s degrees of freedom come from an iteration count set in software, and its effect size is driven by how many binding sites a motif has — SPI1’s 125,604 sites buy it a p of 6.9e-104 for a change of 0.027. Nothing in that column responds to how many people you sequenced.

TOBIAS’s own documentation says as much. The wiki notes that the p-value “can be very small due to the large number of transcription factor binding sites found, so this should always be considered in combination with the change column,” and the BINDetect output description states that when bound-site counts and the change score disagree, “the _change score is the more correct metric to use.” Read TOBIAS through change, and treat its p-value as a statement about consistency across sites within these two libraries.

How to write this up. Quote the Signac p-value for the replicate-level claim and the TOBIAS change for the site-level effect, and never put the two p-values in the same sentence. A reader who sees 1e-138 next to 0.013 will assume the first is the stronger evidence, when it is the one that knows nothing about biological variation.

This is also why TOBIAS cannot be made replicate-aware by arranging its inputs differently. It takes one bigWig per condition; there is no replicate dimension in its data model at all. Its FAQ is direct about this: “TOBIAS does not deal with individual biological replicates, and it is therefore recommended to merge replicate .bam-files prior to correction and footprinting.” Merging is what Step 14 does, and the consequence is that the donor-level question can only be answered on the Signac side.

Method 2: TOBIAS on Cell-Type Pseudobulk BAMs

TOBIAS was built for bulk ATAC-seq and reads BAM files. The job here is to build a “bulk” CD14 monocyte BAM for each condition from the Cell Ranger BAMs, which contain every cell type. R exports three inputs; everything after that happens on the command line.

Exporting Barcodes, Peaks and Motifs From R

#-----------------------------------------------
# STEP 10: Export the three TOBIAS inputs from the Seurat object
#-----------------------------------------------

# (a) One barcode list per donor, for samtools.
# Object barcodes carry a sample prefix added at merging ("severe1_AAAC...-1");
# the CB tag in the Cell Ranger BAM does not ("AAAC...-1"), so strip it.
# One barcode per line, no header: that is what "samtools view -D" expects.
bc_by_donor <- split(sub("^[^_]+_", "", colnames(mono)), mono$sample_id)

invisible(Map(function(barcodes, donor) {
  writeLines(barcodes,
             file.path(tobias_dir, "barcodes", paste0(donor, "_CD14mono.txt")))
}, bc_by_donor, names(bc_by_donor)))

# (b) Peaks that are actually open in CD14 monocytes, as a BED file.
# Footprinting all 185,581 consensus peaks would include regions these cells
# never open. AccessiblePeaks() keeps peaks detected in >= 10 monocytes.
open_peaks <- AccessiblePeaks(mono, idents = c("Healthy", "Severe"))
peaks_gr   <- sort(StringToGRanges(open_peaks))
write.table(
  data.frame(as.character(seqnames(peaks_gr)), start(peaks_gr) - 1, end(peaks_gr)),
  file.path(tobias_dir, "cd14mono_peaks.bed"),   # BED starts are 0-based
  sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE
)
length(peaks_gr)

# (c) The same 633 JASPAR2020 matrices that AddMotifs() stored in Part 4,
# written in JASPAR format so TOBIAS and Signac use identical motifs.
pfm_list    <- GetMotifData(mono, slot = "pwm")          # raw count matrices, rows A/C/G/T
motif_names <- GetMotifData(mono, slot = "motif.names")
jaspar_txt  <- unlist(lapply(names(pfm_list), function(id) {
  m <- pfm_list[[id]]
  c(paste0(">", id, "\t", motif_names[[id]]),
    paste0(rownames(m), " [ ", apply(m, 1, paste, collapse = " "), " ]"))
}))
writeLines(jaspar_txt, file.path(tobias_dir, "JASPAR2020_CORE_vertebrates.jaspar"))

Output of length(peaks_gr):

[1] 115324

Reading this number. 115,324 of the 185,581 consensus peaks — 62 percent — are detected in at least 10 CD14 monocytes and go into the TOBIAS peak set. The remaining 38 percent are peaks other cell types contributed to the consensus, and footprinting inside them would be scanning motifs in chromatin these cells keep closed. Every TOBIAS number below is computed within these 115,324 regions, which is also why the motif-site counts BINDetect reports differ from the n_sites column in Step 6: Signac counted matches in all consensus peaks, BINDetect rescans this monocyte-restricted set.

Why export the motifs instead of downloading them? The bulk tutorial downloaded a JASPAR file from the web. Here that would be a mistake: if TOBIAS used a different JASPAR release from Signac, any disagreement between the two methods could simply be a motif-database difference. Writing the matrices out of the object guarantees that both methods scan for exactly the same 633 motifs.

Now switch to a terminal. The rest of Method 2 runs in the footprint Pixi environment. These steps read and write whole BAM files, so on a shared cluster, run them inside a batch job rather than on a login node; our HPC and SLURM beginner’s guide shows how. Keep --cores equal to the CPUs you request.

Saving the R Session Before You Switch to the Terminal

This is the first boundary between R and the command line, and everything below runs in a terminal. The footprinting you just did is expensive, so write it to disk now rather than relying on the R session still being open when you come back.

#-----------------------------------------------
# STEP 11: Save everything the second R session will need
#-----------------------------------------------

# The object carries the per-cell footprint matrices and the Tn5 bias model,
# so reloading it is far cheaper than re-running Footprint().
saveRDS(mono, file.path(obj_dir, "cd14mono_footprint.rds"))

# The small results from Steps 6, 8 and 9, in one file.
saveRDS(list(fp_motifs    = fp_motifs,
             fp_ids       = fp_ids,
             motif_width  = motif_width,
             flank_cond   = flank_cond,
             flank_donor  = flank_donor,
             flank_wide   = flank_wide,
             donor_test   = donor_test),
        file.path(obj_dir, "signac_footprint_results.rds"))

You can now close R safely. The TOBIAS half takes hours and is usually submitted as a batch job, so assume this session will not survive it. Step 19 reloads both files in one line each. If you do keep R open, Step 19 tells you which part to skip.

Setting the Paths Every Command Below Uses

Steps 13 to 18 all refer to the same five paths. Define them once here.

#-----------------------------------------------
# STEP 12: Paths for every command-line step in this tutorial
#-----------------------------------------------

cd /projects/mylab/shared/scatac-analysis
pixi shell -e footprint

RESULTS=/projects/mylab/shared/scatac_tutorial/results
TOBIAS_DIR=/projects/mylab/shared/scatac_tutorial/footprinting/tobias
GENOME=/projects/mylab/shared/reference/refdata-cellranger-arc-GRCh38-2024-A/fasta/genome.fa
BLACKLIST=/projects/mylab/shared/reference/blacklist/ENCFF356LFX.bed
PEAKS=${TOBIAS_DIR}/cd14mono_peaks.bed

cd ${TOBIAS_DIR}

# Confirm all five resolve before running anything expensive.
# A "No such file or directory" for an empty name means that variable is unset.
ls -d "${RESULTS}" "${TOBIAS_DIR}" "${GENOME}" "${BLACKLIST}" "${PEAKS}"

The genome FASTA must be the one the reads were aligned to — the Cell Ranger reference from Part 1. ATACorrect checks that the BAM and FASTA chromosome names match and stops if they do not. The blacklist is the ENCODE file downloaded in Part 2.

Selecting One Cell Type’s Reads by Barcode

Cell Ranger writes each read’s corrected cell barcode into the CB tag of the BAM. Pulling out one cell type is therefore a tag filter: keep the reads whose CB value appears in our barcode list. samtools view -D CB:FILE does exactly that — it keeps alignments whose CB tag matches one of the values listed in a file.

The same command can apply the quality filters at the same time, so each donor’s BAM is built in a single pass. Those filters matter for method comparison: Signac reads the fragments file, and TOBIAS reads the BAM, so we filter the BAM the way Cell Ranger ATAC filters fragments — properly paired, not a PCR duplicate, mapping quality above 30. Cell Ranger keeps pairs with MAPQ > 30, and samtools view -q keeps reads with MAPQ greater than or equal to the value given, so -q 31 is the match.

#-----------------------------------------------
# STEP 13: Pull each donor's CD14 monocyte reads out of its Cell Ranger BAM
#-----------------------------------------------

# Check the barcode lists before touching any BAM
wc -l barcodes/*.txt

# -D CB:FILE  keep reads whose CB tag is listed in FILE (one barcode per line)
# -f 2        properly paired
# -F 1024     not a PCR duplicate
# -q 31       MAPQ > 30, as in the Cell Ranger fragments file
for s in healthy1 healthy2 severe1 severe2; do
  samtools view -@ 8 -b \
    -D CB:barcodes/${s}_CD14mono.txt \
    -f 2 -F 1024 -q 31 \
    -o bam/${s}_CD14mono.bam \
    ${RESULTS}/${s}/outs/possorted_bam.bam
done

Output:

   807 barcodes/healthy1_CD14mono.txt
  2575 barcodes/healthy2_CD14mono.txt
  1323 barcodes/severe1_CD14mono.txt
   508 barcodes/severe2_CD14mono.txt
  5213 total

Check this against the table in Step 4. The line counts must match the cells per donor exactly. If a file is empty or the numbers differ, the export in Step 10 went wrong, and every downstream TOBIAS result would be built on the wrong cells.

Merging Donors Into One BAM Per Condition

#-----------------------------------------------
# STEP 14: Merge the two donors of each condition and index
#-----------------------------------------------

for cond in Healthy Severe; do
  lower=$(echo ${cond} | tr 'A-Z' 'a-z')   # Healthy -> healthy

  # Both inputs are coordinate-sorted, so the merge is too
  samtools merge -@ 8 -f -o bam/${cond}_CD14mono.bam \
    bam/${lower}1_CD14mono.bam bam/${lower}2_CD14mono.bam

  # ATACorrect needs the index
  samtools index bam/${cond}_CD14mono.bam
done

# Reads available to TOBIAS in each condition
samtools view -c bam/Healthy_CD14mono.bam
samtools view -c bam/Severe_CD14mono.bam

Output:

51419481
33158974

Write these two numbers down. They are the most important pair for judging everything TOBIAS reports. The healthy pseudobulk holds 51.4 million reads from 3,382 cells, the severe one 33.2 million from 1,831 cells — a 1.55-fold difference in reads against a 1.85-fold difference in cells.

Per cell, the severe monocytes are actually slightly deeper: roughly 18,100 reads per cell versus 15,200 in healthy. That matches Part 5, Step 5, where the severe cells had marginally higher median fragment counts. So the gap here is driven by how many monocytes each group contributed, not by libraries of different quality.

It still leaves the severe track noisier in absolute terms. ATACorrect normalizes for read number and BINDetect normalizes footprint scores across conditions, but normalization cannot manufacture observations: an aggregate profile built from 33 million reads is grainier than one built from 51 million. Keep that in mind whenever a difference looks like “severe is lower” — and note that on this dataset the largest change, CEBPB, runs the other way, which makes a pure depth artifact an unlikely explanation for it.

ATACorrect removes duplicates and chrM reads itself, so the filtering in Step 13 is about matching the fragments file and making the pseudobulk definition explicit rather than about protecting TOBIAS.

Correcting Tn5 Sequence Bias With ATACorrect

#-----------------------------------------------
# STEP 15: Correct Tn5 bias
#-----------------------------------------------

# Uses GENOME, PEAKS and BLACKLIST from Step 12
for cond in Healthy Severe; do
  TOBIAS ATACorrect \
    --bam bam/${cond}_CD14mono.bam \
    --genome ${GENOME} \
    --peaks ${PEAKS} \
    --blacklist ${BLACKLIST} \
    --outdir atacorrect \
    --cores 16
done

ATACorrect shifts every read by +4/-5 bp to the Tn5 insertion point, learns the Tn5 sequence preference from reads outside the peaks, and writes four tracks per condition, prefixed with the BAM name:

Open the _atacorrect.pdf for both conditions. It shows the learned Tn5 insertion preference and how much of it remains after correction. After correction the pattern should be largely flattened. If the severe PDF looks much noisier than the healthy one, that is the lower read count showing itself.

Scoring Footprints With ScoreBigwig

#-----------------------------------------------
# STEP 16: Calculate footprint scores
#-----------------------------------------------

# ScoreBigwig replaces the older FootprintScores command.
# Uses PEAKS from Step 12 -- re-run that block if this is a new shell.
for cond in Healthy Severe; do
  TOBIAS ScoreBigwig \
    --signal atacorrect/${cond}_CD14mono_corrected.bw \
    --regions ${PEAKS} \
    --output footprints/${cond}_CD14mono_footprints.bw \
    --cores 16
done

The TOBIAS footprint score looks at each position for a depletion of corrected signal 20-50 bp wide (--fp-min, --fp-max) flanked by 10-30 bp of accessible signal (--flank-min, --flank-max). It rewards both depletion and flank accessibility, so, as the TOBIAS documentation notes, when there is no depletion the score converges to the accessibility of the flanks. A high TOBIAS score does not require a visible dip. That makes it closer to Signac’s flank measure than you might expect — and more robust to the weak-footprint problem described in the introduction.

Detecting Differential Binding With BINDetect

#-----------------------------------------------
# STEP 17: Predict bound sites and differential binding for all 633 motifs
#-----------------------------------------------

# Uses GENOME and PEAKS from Step 12.
# --cond-names   labels in the same order as --signals
# --skip-excel   skip the .xlsx copies; much faster with hundreds of motifs
TOBIAS BINDetect \
  --motifs JASPAR2020_CORE_vertebrates.jaspar \
  --signals footprints/Healthy_CD14mono_footprints.bw \
            footprints/Severe_CD14mono_footprints.bw \
  --genome ${GENOME} \
  --peaks ${PEAKS} \
  --cond-names Healthy Severe \
  --outdir bindetect \
  --skip-excel \
  --cores 16

Note the hyphens. Current TOBIAS options are --cond-names, --signal-labels, --share-y and --plot-boundaries. Some older examples, including our bulk ATAC-seq tutorial, write them with underscores; the current command-line parser rejects those with “unrecognized arguments.”

For each motif, BINDetect finds every match in the peaks, takes its footprint score in each condition, splits sites into bound and unbound (default bound p-value 0.001), and summarizes the change between conditions. The main outputs are:

• bindetect/bindetect_results.txt — one row per motif, with total_tfbs, Healthy_mean_score, Severe_mean_score, Healthy_bound, Severe_bound, Healthy_Severe_change, Healthy_Severe_pvalue and Healthy_Severe_highlighted.
• bindetect/<TF>/beds/ — BED files of all sites and of bound/unbound sites per condition, for example bindetect/FOSL2_MA0478.1/beds/FOSL2_MA0478.1_all.bed.
• bindetect/<TF>/<TF>_overview.txt — every site with its score in each condition and a per-site log2 fold change.
• bindetect/bindetect_figures.pdf and an interactive .html volcano plot.

Two details of the _change and _pvalue columns decide how you read everything below:

Direction. The TOBIAS documentation defines <condition1>_<condition2>_change so that negative values mean more binding in condition 2. With --cond-names Healthy Severe, a negative Healthy_Severe_change means more binding in Severe.

What the p-value tests. In the TOBIAS preprint, the p-value is obtained by repeatedly subsampling the background and testing the observed change against it. The units being resampled are binding sites, not donors. Like the single-cell test in Part 5, it can become very small because there are thousands of sites, and the TOBIAS documentation itself advises reading it together with the change score. It is not a test across biological replicates. With one merged BAM per condition, TOBIAS cannot know that the healthy BAM is 76 percent healthy2.

Aggregate Footprint Plots From TOBIAS

To look at the bias-corrected TOBIAS signal itself, aggregate it over all binding sites of a motif. The _all.bed file from BINDetect contains every motif site in the peaks, so both conditions are compared at identical positions.

#-----------------------------------------------
# STEP 18: Aggregate corrected signal at CTCF and FOSL2 sites
#-----------------------------------------------

# All paths are relative to the tobias directory (Step 12)
cd ${TOBIAS_DIR}

for tf in CTCF_MA0139.1 FOSL2_MA0478.1; do
  TOBIAS PlotAggregate \
    --TFBS bindetect/${tf}/beds/${tf}_all.bed \
    --signals atacorrect/Healthy_CD14mono_corrected.bw \
              atacorrect/Severe_CD14mono_corrected.bw \
    --signal-labels Healthy Severe \
    --share-y both \
    --plot-boundaries \
    --flank 100 \
    --title ${tf} \
    --output aggregate/${tf}_aggregate.png
done

Reading these figures: three panels each — healthy alone, severe alone, then both overlaid. The y-axis is the ATACorrect bias-corrected signal averaged over every site of the motif (59,954 for CTCF, 45,873 for FOSL2), with the motif boundaries as dashed lines. Because the expected signal has already been subtracted, zero means “cut exactly as much as the sequence predicts,” and negative means protected.

FOSL2 is the textbook footprint this tutorial warned you might never see. Two sharp shoulders at about +/- 20 bp, and between them a deep central depression that goes negative, down to -0.05, precisely over the motif. Fewer Tn5 insertions than sequence alone predicts, at the exact bases a protein would occupy. This is what the Signac panels could not show, because Signac’s hexamer bias model and 250 bp window smear it, while ATACorrect’s dinucleotide model and 100 bp window resolve it.

CTCF has no such dip, despite being the canonical footprinting example. Its middle is a forest of spikes reaching 0.22 with no sustained protection. Footprint shape is motif-specific and depends on the bias model, the depth and the site set, so “CTCF always footprints cleanly” is a statement about deeply sequenced bulk data, not about a cell type extracted from four scATAC libraries.

Now compare how the two conditions differ in each, because the difference is not the same kind of difference. In the CTCF overlay, the healthy track sits above the severe track everywhere — at the shoulders, in the middle, and out at +/- 100 bp where no CTCF motif has any influence. A uniform vertical offset across the whole window is the signature of a global difference in signal level between the two tracks, not of transcription factor binding.

In the FOSL2 overlay, the two tracks converge at +/- 100 bp and inside the central dip, and separate only at the shoulders — healthy about 0.096 against severe 0.075 at the -20 bp peak. A difference localized to the part of the profile that reflects binding, with agreement everywhere else, is the shape a real change in occupancy makes.

This is the cleanest single explanation of the CTCF problem in this run. Signac scored CTCF at -0.163, the second largest change in the five-motif panel, because its flank metric averages a window that a uniform offset shifts just as effectively as a real footprint change. TOBIAS scored it at -0.027 because its footprint score is built from local depletion relative to nearby flanks, which a uniform offset largely cancels out of. When the two methods disagree about a motif, plot the aggregate and look at the far flanks. If the tracks differ 100 bp from the motif, you are looking at normalization, not binding.

And note what this does to the FOSL2 conclusion. The TOBIAS aggregate shows healthy shoulders above severe — AP-1 signal down in severe disease, agreeing with the negative tobias_delta, the negative signac_delta, and Part 5’s enrichment of AP-1 motifs in peaks that closed. Four consistent observations. The donors still do not separate, and that still governs what you may claim.

Back in R: Reloading and Reading the BINDetect Results

The command-line work is finished. Everything from here runs in R again.

#-----------------------------------------------
# STEP 19: Reopen R and reload what Step 11 saved
#-----------------------------------------------
# SKIP THIS BLOCK if your R session from Step 11 is still open.

library(Signac)
library(Seurat)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomicRanges)
library(ggplot2)
library(ggrepel)
library(patchwork)

set.seed(1234)
options(future.globals.maxSize = 60 * 1024^3)

# Same paths as Step 3
base_dir   <- "/projects/mylab/shared/scatac_tutorial"
annot_dir  <- file.path(base_dir, "annotation")
da_dir     <- file.path(base_dir, "differential")
fp_dir     <- file.path(base_dir, "footprinting")
obj_dir    <- file.path(fp_dir, "objects")
plot_dir   <- file.path(fp_dir, "plots")
table_dir  <- file.path(fp_dir, "tables")
tobias_dir <- file.path(fp_dir, "tobias")

mono <- readRDS(file.path(obj_dir, "cd14mono_footprint.rds"))

# list2env() unpacks the saved list into fp_motifs, motif_width, flank_cond,
# flank_donor, flank_wide and donor_test as separate objects, exactly as
# they were named in the first session.
invisible(list2env(readRDS(file.path(obj_dir, "signac_footprint_results.rds")),
                   envir = .GlobalEnv))

ls()

What you should see: ls() lists mono, donor_test, flank_cond, flank_donor, flank_wide, fp_ids, fp_motifs and motif_width, plus the path variables. If any is missing, Step 11 did not complete.

Reading the BINDetect Results

#-----------------------------------------------
# STEP 20: Load BINDetect results
#-----------------------------------------------

bd <- read.delim(file.path(tobias_dir, "bindetect", "bindetect_results.txt"))

# Find the comparison columns by suffix rather than hard-coding their names
chg_col <- grep("_change$",      colnames(bd), value = TRUE)
p_col   <- grep("_pvalue$",      colnames(bd), value = TRUE)
hl_col  <- grep("_highlighted$", colnames(bd), value = TRUE)

# BINDetect reports Healthy relative to Severe (negative = more bound in Severe).
# Flip the sign so that, as in Parts 5 and 6, positive means "higher in severe".
bd$tobias_delta <- -bd[[chg_col]]
bd$tobias_p     <- bd[[p_col]]
bd$highlighted  <- as.logical(bd[[hl_col]])

nrow(bd)
table(bd$highlighted)
head(bd[order(bd$tobias_p), c("name", "motif_id", "total_tfbs", "tobias_delta", "tobias_p")], 10)

Output:

[1] 633

FALSE  TRUE
  565    68

      name motif_id total_tfbs tobias_delta     tobias_p
401  KLF15 MA1513.1      84511     -0.14993 1.18481e-163
557  CEBPA MA0102.4      30057      0.16903 2.27831e-148
545 ZBTB14 MA1650.1      34763     -0.13726 1.57216e-144
252  CEBPE MA0837.1      10957      0.23083 3.86461e-144
558  CEBPD MA0836.2      26938      0.16097 4.73100e-144
40    NRF1 MA0506.1      32652     -0.15757 7.96636e-143
410    MAZ MA1522.1      98275     -0.07284 1.42448e-141
548 ZNF148 MA1653.1     134261     -0.07614 2.95927e-141
554   ATF4 MA0833.2      27649      0.13475 3.93947e-139
251  CEBPB MA0466.2      11167      0.23037 9.82676e-138

Start with the p-values, so they stop distracting you. The smallest is 1e-163. These come from resampling thousands of binding sites, not from comparing donors, so with 84,511 KLF15 sites even a change of -0.15 is enormously “significant.” Nothing in this column speaks to whether four people differ. Sort by it to see what TOBIAS is most confident about within these two libraries; never quote it as evidence of a disease effect.

Now read the table as families rather than rows. Four of the top ten are C/EBP proteins — CEBPA, CEBPE, CEBPD, CEBPB — and all four move the same way, up in severe disease. That is one finding reported four times, because these factors bind nearly identical sites. ATF4, also positive, is a bZIP factor binding a related composite element.

The negative side is a family too, and a less obvious one. KLF15, ZBTB14, NRF1, MAZ and ZNF148 are all GC-rich, CpG-island, promoter-proximal binders. They are not a coordinated group of transcription factors responding to disease; they are what you get when signal shifts away from GC-rich promoters. Before writing “NRF1 activity falls in severe COVID-19,” consider the simpler description: the accessible landscape tilts from promoters toward C/EBP-bearing elements.

highlighted is 68 of 633, and it is not a significance call. BINDetect marks motifs above the 95th percentile of -log10(p) and/or in the top 5 percent of change in either direction. Those are percentiles, so a run with no real differences would also highlight roughly this many. Treat the highlighted set as “the extremes of this comparison.”

Looking at the Families Directly

The top-ten view is misleading in one direction and honest in another, so it is worth pulling out one family in full.

#-----------------------------------------------
# STEP 20b: The C/EBP family, and what the top of each direction looks like
#-----------------------------------------------

# BINDetect's own clustering of motifs by binding-site overlap
bd[grep("^CEBP", bd$name), c("name", "cluster", "total_tfbs", "tobias_delta", "tobias_p")]

# The strongest movers in each direction, by effect size rather than p-value
head(bd[order(-bd$tobias_delta), c("name", "cluster", "tobias_delta")], 12)
head(bd[order( bd$tobias_delta), c("name", "cluster", "tobias_delta")], 12)

Output:

             name       cluster total_tfbs tobias_delta     tobias_p
252         CEBPE       C_CEBPA      10957      0.23083 3.86461e-144
251         CEBPB       C_CEBPA      11167      0.23037 9.82676e-138
253         CEBPG       C_CEBPA       9691      0.21948 2.08792e-134
557         CEBPA       C_CEBPA      30057      0.16903 2.27831e-148
558         CEBPD       C_CEBPA      26938      0.16097 4.73100e-144
556  CEBPG(var.2) C_CEBPG(var.2)     31747      0.15064 2.03159e-137

# strongest up in severe
             name             cluster tobias_delta
           CEBPE             C_CEBPA      0.23083
           CEBPB             C_CEBPA      0.23037
           CEBPG             C_CEBPA      0.21948
           CEBPA             C_CEBPA      0.16903
             DBP             C_NFIL3      0.16788
           CEBPD             C_CEBPA      0.16097
    CEBPG(var.2)      C_CEBPG(var.2)     0.15064
             HLF             C_NFIL3      0.14791
             TEF             C_NFIL3      0.14418
            ATF4      C_CEBPG(var.2)     0.13475
           NFIL3             C_NFIL3      0.12283
             JUN C_FOSL1::JUN(var.2)     0.09777

# strongest up in healthy
             name             cluster tobias_delta
            NRF1              C_NRF1     -0.15757
           KLF15             C_KLF15    -0.14993
           TFDP1              C_E2F6    -0.13806
          ZBTB14            C_ZBTB14    -0.13726
          ZBTB33            C_ZBTB33    -0.11863
           TCFL5             C_TCFL5    -0.11619
            HES1               C_MNT   -0.09716
           HINFP             C_HINFP   -0.09452
             YY2               C_YY2   -0.08408
          TFAP2E    C_TFAP2C(var.2)   -0.08292
            ZIC5              C_ZIC5   -0.07958
            EGR3              C_EGR1   -0.07928

The cluster column is BINDetect doing the family grouping for you. It clusters motifs by how much their predicted binding sites overlap in the genome, and names each cluster after one member. All six C/EBP entries fall into two clusters, C_CEBPA and C_CEBPG(var.2). Report them as “C/EBP-family motifs,” and cite the cluster rather than listing six proteins.

The severe side is one structural class, not twelve findings. Eleven of the top twelve are bZIP factors: the C/EBP cluster, the PAR-bZIP cluster (C_NFIL3: DBP, HLF, TEF, NFIL3), and ATF4. These all bind variations on the same palindromic element, so the assay cannot separate them, and the honest statement is that bZIP/C/EBP-type elements gain footprint signal in severe disease.

JUN is the interesting twelfth. It sits at +0.098, up in severe — while FOSL2, which Part 5 flagged in peaks that closed in severe, sits at -0.022. Both are AP-1. That is not a contradiction the data can resolve for you: AP-1 is a family of heterodimers with overlapping but non-identical sites, and different members can move in different directions. It does mean that “AP-1 closes in severe monocytes” is too coarse a statement for this result, and Part 5’s conclusion needs the qualifier that it was based on 50 peaks and on the FOS/FOSL-type matrices specifically.

The healthy side is a promoter signature. NRF1, KLF15, TFDP1/E2F6, ZBTB33, HINFP, YY2, ZNF148 and EGR3 are GC-rich binders concentrated at CpG-island promoters. A redistribution of accessible signal from promoters toward distal bZIP-bearing elements would produce exactly this list, without any of these factors changing their behaviour.

How this sits against the literature. Giroux et al. (2022) profiled PBMC chromatin accessibility in COVID-19 and reported that CD14+ monocytes from patients with moderate symptoms were characterized by up-regulation of the CEBPB regulatory region. That is a different severity group, a different cohort and a different analysis — accessibility at the gene’s own locus rather than footprints at its motif — so it is not confirmation. It is a reason to take the C/EBP result seriously enough to test properly, with more donors and with expression data.

Read the top of the list by families, not rows. Expect groups of near-identical motifs to move together — several AP-1 dimers, several CEBP members, several RFX members. BINDetect’s cluster column groups motifs whose binding sites overlap. As in Part 5, ten rows from one cluster are one finding, not ten.

Comparing the Signac and TOBIAS Results

We now have three independent kinds of evidence about the same five motifs in the same cells:

• Part 5 motif enrichment — sequence only, on 46 opened and 50 closed differential peaks.
• Signac flank enrichment — Tn5 insertions from the fragments file, every motif site in the peaks, bias from a hexamer model, per donor.
• TOBIAS differential binding — Tn5 insertions from filtered BAMs, every motif site in monocyte-accessible peaks, bias from a dinucleotide model, per condition.

#-----------------------------------------------
# STEP 21: Put the three lines of evidence side by side
#-----------------------------------------------

cmp <- merge(flank_wide[, c("feature", "signac_delta")],
             bd[, c("name", "tobias_delta", "tobias_p", "total_tfbs")],
             by.x = "feature", by.y = "name")

# Do the two footprinting methods move in the same direction?
cmp$same_direction <- sign(cmp$signac_delta) == sign(cmp$tobias_delta)

# The replicate-level result from Step 9
cmp$signac_p    <- donor_test$p[match(cmp$feature, donor_test$feature)]
cmp$signac_pct  <- donor_test$pct[match(cmp$feature, donor_test$feature)]

# Part 5 motif enrichment in peaks opened / closed in severe disease
m_up <- read.csv(file.path(da_dir, "tables", "motifs_opened_severe.csv"))
m_dn <- read.csv(file.path(da_dir, "tables", "motifs_closed_severe.csv"))
cmp$part5_padj_opened <- m_up$p.adjust[match(cmp$feature, m_up$motif.name)]
cmp$part5_padj_closed <- m_dn$p.adjust[match(cmp$feature, m_dn$motif.name)]

cmp

Output:

  feature signac_delta tobias_delta     tobias_p total_tfbs same_direction
1   CEBPB   0.29707392      0.23037 9.82676e-138      11167           TRUE
2    CTCF  -0.16326301     -0.02669  9.76539e-82      59954           TRUE
3   FOSL2  -0.07349232     -0.02214  2.74649e-72      45873           TRUE
4    RFX2   0.04789249      0.09382  3.78557e-99       9610           TRUE
5    SPI1   0.03896793      0.02664 6.92667e-104     125604           TRUE
  part5_padj_opened part5_padj_closed signac_p signac_pct
1      1.0000000000      1.000000e+00  0.01286      26.55
2      1.0000000000      1.000000e+00  0.25048      -7.78
3      1.0000000000      8.803895e-06  0.40484      -2.08
4      0.0006599133      1.000000e+00  0.87629       1.54
5      1.0000000000      1.000000e+00  0.15061       3.51

The first thing to notice is that same_direction is TRUE five times out of five. Two tools that read different files (fragments versus a filtered BAM), model Tn5 bias differently (a hexamer table versus a dinucleotide weight matrix), and score entirely different quantities (windowed insertion enrichment versus a depletion-plus-flank footprint score) agree on the sign of every motif. That is the single most reassuring result in this tutorial. It says the pipeline is measuring something real rather than each tool’s own artifacts.

It also means direction alone is cheap. Five for five is what you would expect if both methods are sensitive to the same underlying signal, including the same underlying technical signal. Agreement in direction is a necessary check, not evidence of biology.

CEBPB is the only row that survives every column. Largest change in both methods, and the only one the donor t-test calls at all (p = 0.013, +26.5 percent). Part 5 gave it nothing (p.adjust = 1 in both directions), which is not a contradiction: Part 5 tested only the 46 opened and 50 closed peaks, while footprinting used all 11,167 CEBPB sites in monocyte-accessible chromatin. This is precisely the case footprinting exists to catch — a change spread thinly across thousands of sites, invisible to a test restricted to the handful of peaks that reached significance on four donors.

FOSL2 and RFX2 are the honest disappointments. Both are the motifs Part 5 pointed at, both agree in direction between the two footprinting methods, and both have Part 5 enrichment backing that direction — FOSL2 in closed peaks (p.adjust 8.8e-06), RFX2 in opened peaks (p.adjust 6.6e-04). Three lines of evidence, aligned in direction. And the donor test returns p = 0.40 and p = 0.88. The correct write-up is “directionally consistent across three analyses, not supported at the level of biological replicates,” and then a request for more donors.

CTCF and SPI1 are the controls doing their job. CTCF has the second-largest Signac effect in the table (-7.8 percent) and a p of 0.25, because its two severe donors differ from each other more than the two groups do. TOBIAS independently puts it at -0.027, the second smallest change of the five, and the aggregate plot in Step 18 shows why. Three separate pieces of the analysis agree that the largest-looking control effect is not one.

Note the two p-value columns side by side, and do not compare them. signac_p ranges from 0.013 to 0.88 and counts four donors; tobias_p ranges from 1e-72 to 1e-138 and counts binding sites and Monte Carlo draws. They are in the table together for convenience, not for comparison — see the statistics table at the end of Method 1.

Scales are not comparable between the two methods. Signac’s CTCF change (-0.163) is six times its TOBIAS change (-0.027), while for RFX2 the TOBIAS change is twice the Signac one. The units are different and the ranks only roughly agree. Compare signs and rough ordering; never plot one against the other and fit a line.

How to read each row. The table is designed to be read motif by motif, with a small set of explicit rules:

What this run actually supports. One motif family, C/EBP, gained footprint signal in severe-disease CD14 monocytes, by both methods, with a between-condition gap four times the donor-to-donor spread, and with the whole C/EBP cluster plus the related PAR-bZIP cluster moving together in the genome-wide BINDetect screen. Everything else in the five-motif panel is either directionally suggestive without replicate support (FOSL2, RFX2) or a control that moved for reasons the controls exist to expose (CTCF, SPI1).

What it does not support. It does not establish that C/EBP factors drive severe COVID-19 monocyte biology. Four donors, two per group, cannot establish that. It does not identify which C/EBP protein is responsible, because CEBPA, CEBPB, CEBPD, CEBPE and CEBPG bind the same element and BINDetect places them in one cluster. And it does not resolve the AP-1 question from Part 5: FOSL2 moves down while JUN moves up, which is a real complication rather than a tidy confirmation.

How this connects to the literature. Part 5 noted Brauns et al. (2022), who found increased AP-1 and MAF accessibility in monocytes of patients recovering from severe COVID-19, and argued that reduced AP-1 accessibility in acute disease would be directionally consistent. The footprinting result is messier than that: the FOS/FOSL-type matrices move down and JUN moves up. The C/EBP result has a separate and closer anchor in Giroux et al. (2022), who reported up-regulation of the CEBPB regulatory region in CD14+ monocytes from COVID-19 patients — though in a moderate-symptom outpatient cohort, and measured as accessibility at the gene locus rather than as footprints at its motif. Neither paper confirms anything here. Both make the C/EBP direction worth a properly powered experiment.

📊 Visualization

All of these run in the R session you reopened at Step 19. The TOBIAS aggregate plots were already produced on the command line in Step 18, alongside the tool that made them.

Signac Footprint Profiles by Condition

#-----------------------------------------------
# STEP 22: Footprint profiles, severe vs healthy
#-----------------------------------------------

# Default normalization = "subtract": observed minus sequence-expected.
# show.expected = TRUE draws the expected Tn5 profile under each panel.
p_fp_condition <- PlotFootprint(mono, features = fp_motifs, group.by = "condition") +
  plot_layout(ncol = 2)

ggsave(file.path(plot_dir, "01_signac_footprint_by_condition.png"), p_fp_condition,
       width = 11, height = 13, dpi = 300)

Reading this figure: each panel is one motif. The x-axis is distance from the motif centre in base pairs; the upper track is Tn5 insertion enrichment (observed minus sequence-expected) for each condition, and the lower track is the expected profile from sequence alone.

The first thing to see is that none of the five shows a peak-dip-peak. All five are a single broad hump, 300 to 400 bp wide, centred on the motif, with a violent spike-and-crash region in the middle 20 bp. If you came looking for the textbook footprint, it is not here — and this is the normal appearance of sparse single-cell pseudobulk data, not a failed run.

Look at the expected track underneath to understand why. It is flat at 1.0 across the entire window except for a sharp excursion in exactly that middle 20 bp: up to 2.5 for CTCF and CEBPB, down to 0.4 for SPI1. Almost everything happening at the motif itself is sequence composition, which is precisely why this tutorial measures the flank. The centre of a Signac footprint plot on scATAC-seq data is a bias readout with a little biology buried in it.

Noise tracks n_sites from Step 6, exactly. SPI1 (117,781 sites) is the smoothest curve; CEBPB (7,127) is visibly the jaggiest, with excursions of +/- 0.3 that are pure sampling. Before reading any difference between two lines, check the motif’s site count. A CEBPB wiggle and an SPI1 wiggle do not mean the same thing.

CEBPB is the one panel where the two lines genuinely part company. The severe line sits above the healthy line across roughly 200 bp of the rise and fall, not at one or two positions. That sustained offset is what produced the +0.297 in Step 8, and it is the only separation in this figure you can see without being told where to look.

CTCF, the control, is where the eye deceives. The healthy line runs slightly above severe on the left shoulder — visible, but a fraction of the CEBPB gap, and the two lines are indistinguishable beyond +/- 150 bp. Step 8 scored this as -0.163, the second largest change in the panel, and Step 9 explains why that does not survive: the gap is real but smaller than the spread between the two severe donors, so the t-test returns p = 0.25.

SPI1, FOSL2 and RFX2 are effectively overlapping. Their Step 8 deltas (+0.039, -0.073, +0.048) are the residue of lines that trace each other.

Ignore the label ordering. PlotFootprint() labels groups by flank height, so with two groups both get labelled and the order carries no statistical meaning.

Donor-Level Footprints and Flank Enrichment

#-----------------------------------------------
# STEP 23: The same footprints, one line per donor
#-----------------------------------------------

# One control and one hypothesis motif, split by donor
p_fp_donor <- PlotFootprint(mono, features = c("CTCF", "FOSL2"),
                            group.by = "sample_id", label.top = 4) +
  plot_layout(ncol = 1)

ggsave(file.path(plot_dir, "02_signac_footprint_by_donor.png"), p_fp_donor,
       width = 8, height = 10, dpi = 300)

# Flank enrichment per donor for all five motifs
p_flank_donor <- ggplot(flank_donor, aes(x = condition, y = corrected, color = condition)) +
  geom_point(size = 3) +
  geom_text_repel(aes(label = group), size = 3, show.legend = FALSE) +
  facet_wrap(~ feature, scales = "free_y", nrow = 1) +
  scale_color_manual(values = c(Healthy = "#2E86AB", Severe = "#C0392B")) +
  labs(x = NULL,
       y = "Motif-flank Tn5 enrichment\n(observed - expected)",
       title = "Signac motif-flank enrichment per donor -- CD14 monocytes") +
  theme_classic(base_size = 11) +
  theme(legend.position = "none")

ggsave(file.path(plot_dir, "03_signac_flank_by_donor.png"), p_flank_donor,
       width = 12, height = 4, dpi = 300)

Reading the donor footprints (first figure): the four CTCF lines lie almost exactly on top of each other across the whole window. That is the noise floor of the method, and it is impressively low — which is worth knowing, because it means the small CTCF difference in Step 8 is not random jitter. It is a small, consistent offset. Small and consistent is still small.

FOSL2 tells a different story, and not the one the numbers told. Here severe1 rides visibly above the other three donors across the left rise, from about -150 to -50 bp. Yet in the flank table severe1 has the lowest FOSL2 value of all four donors. Both are true: the flank window only spans the motif edge out to +50 bp, so the elevation that dominates the picture sits almost entirely outside the region being measured. The eye is drawn to the 300 bp hump; the metric reads a 50 bp shoulder. When a figure and a number disagree, check which piece of the x-axis each one is using.

Reading the per-donor dot plot (second figure), starting with a warning about the axes. Each panel has its own y-scale, because scales = "free_y" was needed to see anything at all. SPI1’s axis spans 0.05 units and CEBPB’s spans 0.35 — seven times wider. The SPI1 panel therefore looks as cleanly separated as CEBPB while representing a gap one twentieth the size. Always read the axis numbers on a free-scale facet before reading the pattern.

CEBPB is the only panel that separates on any scale. Both severe donors sit far above both healthy donors, and the gap between the groups is visibly larger than the distance between the two donors within either group. This is what p = 0.013 and t = 8.7 look like.

CTCF and SPI1 show why the eye is not a test. In CTCF the four points do fall in condition order, but severe2 (1.456) sits just below healthy1 (1.475) while severe1 sits far below severe2 — the within-severe spread swamps the between-group gap, and the t-test returns 0.25. In SPI1 the four points are almost evenly spaced across a 0.05 axis and happen to alternate correctly; p = 0.15. Both look convincing at a glance and neither is an effect.

FOSL2 and RFX2 show honest overlap. In FOSL2, healthy2 is an outlier at 2.48 while healthy1 sits between the two severe donors. In RFX2 the four donors interleave completely, with severe1 highest and severe2 lowest — the two severe donors bracket the entire healthy pair. A pooled RFX2 difference computed from that is not describing a condition.

TOBIAS BINDetect Volcano Plot

BINDetect writes its own volcano plots, but a ggplot version lets us mark the five motifs and keep the direction convention used throughout this series.

#-----------------------------------------------
# STEP 24: BINDetect volcano plot with the five motifs labelled
#-----------------------------------------------

bd$label <- ifelse(bd$name %in% fp_motifs, bd$name, NA)

p_volcano_tobias <- ggplot(bd, aes(x = tobias_delta, y = -log10(tobias_p))) +
  geom_point(aes(color = highlighted), size = 1.2, alpha = 0.7) +
  geom_text_repel(aes(label = label), size = 3.5, na.rm = TRUE,
                  min.segment.length = 0, max.overlaps = Inf) +
  scale_color_manual(values = c(`FALSE` = "grey75", `TRUE` = "#8E44AD"),
                     labels = c(`FALSE` = "Not highlighted", `TRUE` = "Highlighted by BINDetect"),
                     name = NULL) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  labs(x = "Differential binding score (Severe - Healthy)",
       y = "-log10 BINDetect p-value",
       title = "TOBIAS BINDetect -- CD14 monocytes, 633 JASPAR2020 motifs") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "04_tobias_bindetect_volcano.png"), p_volcano_tobias,
       width = 8, height = 6, dpi = 300)

Reading this figure: each point is one of the 633 motifs. Right of the dashed line means a higher footprint score in severe disease, left means higher in healthy. The y-axis reaches 165, because p-values come from resampling tens of thousands of binding sites — it measures consistency across sites, never across donors.

The V shape is the point. Effect size and significance rise together because both are driven by the same site counts, so the plot is really a single axis wearing two. Read the x-axis and treat the y-axis as a rough confidence ordering within this one comparison.

CEBPB sits alone at the right edge, at +0.23, further from zero than all but a handful of the 633. The cluster of purple points around it at +0.14 to +0.23 is the rest of the C/EBP and PAR-bZIP family. A single labelled motif out on a limb would deserve suspicion; a labelled motif with its whole family beside it is the pattern you want.

RFX2 at +0.094 sits inside the highlighted group on the right shoulder, respectable but far short of CEBPB.

SPI1, CTCF and FOSL2 are all crowded against the dashed line, grey rather than purple, indistinguishable from the bulk of the 633 motifs. For CTCF, that is the control behaving correctly — which directly contradicts what the Signac flank metric said about it, and is the reason Step 21 trusts TOBIAS on this point.

The two wings are not symmetric. The left wing is a dense purple cluster around -0.05 to -0.16 — the GC-rich promoter factors. The right wing is sparser but reaches further. That asymmetry is the redistribution described in Step 20: many promoter-type motifs drifting modestly down, a smaller family moving sharply up.

Do the Two Methods Agree?

#-----------------------------------------------
# STEP 25: Method concordance for the five motifs
#-----------------------------------------------

p_concordance <- ggplot(cmp, aes(x = signac_delta, y = tobias_delta)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(aes(shape = signac_p < 0.05), size = 3.5, color = "#34495E") +
  geom_text_repel(aes(label = feature), size = 4) +
  scale_shape_manual(values = c(`TRUE` = 16, `FALSE` = 1),
                     name = "Signac donor\nt-test p < 0.05") +
  labs(x = "Signac: change in motif-flank enrichment (Severe - Healthy)",
       y = "TOBIAS: differential binding score (Severe - Healthy)",
       title = "Signac vs TOBIAS -- CD14 monocytes, severe vs healthy") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "05_method_concordance.png"), p_concordance,
       width = 7, height = 6, dpi = 300)

Reading this figure: all five points land in the upper-right or lower-left quadrant, so the two methods agree on direction for every motif. Filled circles passed the donor-separation check in Step 9.

CEBPB is far out on its own in the upper right, roughly three times further from the origin on both axes than anything else, and filled. That is the shape of the one reportable result.

The two open circles are the Part 5 hypotheses. RFX2 sits up and to the right of the origin, FOSL2 down and to the left — both in the direction Part 5 predicted, both hollow. The figure shows agreement and withholds support in the same glance, which is exactly what it is for.

CTCF in the lower left is the figure’s most useful point. It has the most negative Signac value of the five, yet its TOBIAS value is barely below zero, and it is hollow because the donor t-test returns p = 0.25. A point that is extreme on one axis, unremarkable on the other, and non-significant across donors is a measurement disagreement rather than a biological finding.

Five points is not a correlation. Do not fit a line or quote a coefficient. The axes are in different units and the ranks only roughly agree — CTCF is second by Signac and fourth by TOBIAS. For a genome-wide method comparison, footprint the top BINDetect motifs in Signac and repeat, at the compute cost warned about above.

💾 Saving Results for Reuse and Supplementary Material

#-----------------------------------------------
# STEP 26: Save the object, tables and session information
#-----------------------------------------------

# The object was already saved in Step 11 and has not changed since,
# so only the tables are written here.
write.csv(flank_cond,   file.path(table_dir, "signac_flank_by_condition.csv"), row.names = FALSE)
write.csv(flank_donor,  file.path(table_dir, "signac_flank_by_donor.csv"),     row.names = FALSE)
write.csv(donor_test,   file.path(table_dir, "signac_donor_ttest.csv"),        row.names = FALSE)
write.csv(cmp,         file.path(table_dir, "method_comparison_5_motifs.csv"), row.names = FALSE)
write.csv(bd,          file.path(table_dir, "tobias_bindetect_results.csv"),  row.names = FALSE)

# Full Signac footprint profiles, for anyone who wants to replot them
write.csv(GetFootprintData(mono, features = fp_motifs, group.by = "sample_id"),
          file.path(table_dir, "signac_footprint_profiles_by_donor.csv"), row.names = FALSE)

writeLines(capture.output(sessionInfo()), file.path(fp_dir, "sessionInfo_part6.txt"))

Record the TOBIAS and samtools versions too — they live in the other environment, so sessionInfo() cannot see them:

cd /projects/mylab/shared/scatac-analysis
pixi list -e footprint > /projects/mylab/shared/scatac_tutorial/footprinting/pixi_footprint_env.txt

✅ Best Practices for scATAC-seq Footprinting Analysis

1. Footprint one cell type at a time. Footprints averaged across cell types mostly reflect which cell types are present. Subset first, as in Step 4, exactly as you did for differential accessibility.

2. Screen broadly, then footprint selectively. Signac footprints only the motifs you name, and each one costs real compute, so it is the wrong tool for discovery. Let BINDetect rank all 633 motifs first, then footprint the top hits plus your controls in Signac, where you can split them by donor. Motif enrichment (Part 5) and prior biology are the other legitimate ways to pick that shortlist.

3. Always include a control motif, and expect to learn from it. CTCF is the standard choice. Its job is not only to confirm the pipeline works — in this run it produced the second-largest Signac change in the panel, and tracing why (a uniform offset between the two aggregate tracks, visible 100 bp from any motif) is what calibrated every other number. A control that moves is information, not a failure.

4. Read the flank, not the dip. In scATAC-seq, the central dip is strongly shaped by Tn5 sequence bias. Motif-flank accessibility is the more robust measure, and it is the one Baek et al. (2017) showed carries information even for factors without a measurable footprint.

5. Test across donors, not across cells or sites. In Signac this costs one argument (group.by = "sample_id") and one t.test(). Resist eyeballing whether the donor values separate: with two per group they separate by chance one time in three, and in this dataset they separated for both control motifs and for neither hypothesis. Report the p-value and the percent change together, and correct across the motifs you tested.

6. Use identical motifs in every method. Export the motif matrices from the object, as in Step 10, rather than downloading a different JASPAR release. Otherwise method disagreements may just be database differences.

7. Match the input reads between methods. Filter the BAM to properly paired, non-duplicate reads with MAPQ > 30 so that TOBIAS sees what the fragments file contains.

8. Use the genome FASTA the reads were aligned to. For TOBIAS that is the Cell Ranger reference fasta/genome.fa. ATACorrect checks that the BAM and FASTA chromosomes match.

9. Report depth for every group. State the cells and reads behind each pseudobulk. The severe CD14 monocyte BAM is built from 1,831 cells and the healthy from 3,382; a reader needs those numbers to judge how noisy each footprint is. If your conditions are very unequal, consider downsampling the larger BAM (for example samtools view -s) and checking that the conclusions hold.

10. Interpret TOBIAS p-values as site-level evidence. They measure how consistently binding sites shift, not whether donors differ. For donor-level evidence, build one BAM per donor and compare, or use Signac by donor.

11. Talk about families, not proteins. Motifs from the same family are bound by the same DNA sequence, and footprinting cannot separate them. Write “AP-1 family footprints,” not “FOSL2 binding.”

12. Define your paths in one block and re-run it in every new shell. Shell variables do not survive a new SSH session or a batch script. An unset variable does not stop TOBIAS with a clear message — it produces expected one argument, which reads like a syntax error in your command. The ls -d check in Step 12 turns that into an obvious answer in one second.

13. Diagnose method disagreements with the aggregate plot. When Signac and TOBIAS disagree about a motif, run PlotAggregate and look at the edges of the window. Tracks separated at +/- 100 bp differ globally; tracks that converge there and separate only at the shoulders differ in binding.

14. Record every version. Signac, TOBIAS and the motif database all shape the result. Save sessionInfo() and the Pixi environment listing with every analysis. This tutorial’s results were produced with R 4.5.3, Signac 1.17.1, Seurat 5.5.1, BSgenome.Hsapiens.UCSC.hg38 1.4.5 and TOBIAS 0.17.5.

⚠️ Common Pitfalls and How to Avoid Them

🎯 Conclusion

You have taken the motif families from Part 5 and tested them with an entirely different kind of evidence. Along the way you have:

1. Separated two questions that are often confused: motif enrichment, which reads DNA sequence in a chosen set of peaks, and footprinting, which reads the base-pair pattern of Tn5 insertions from your own cells.
2. Learned the limits of footprinting: Tn5 sequence bias, factors that leave no footprint, and why the flank height is more trustworthy than the central dip in sparse single-cell data.
3. Run Signac Footprint() directly on the fragments file, and quantified motif-flank enrichment by condition and by donor.
4. Built cell-type pseudobulk BAMs from Cell Ranger output with samtools, selecting cells by CB tag and filtering to match the fragments file.
5. Run the full TOBIAS pipeline — ATACorrect, ScoreBigwig and BINDetect — on those BAMs, with exactly the same 633 motifs Signac uses.
6. Compared the two methods and the Part 5 enrichment with explicit rules for what agreement and disagreement mean.

What the analysis found. Both methods agreed on the direction of all five motifs. One result survived every check: C/EBP-family motifs gain footprint signal in severe-disease CD14 monocytes — the largest change in both methods, the only motif the donor-level t-test calls (p = 0.013, +26.5 percent), and the entire C/EBP cluster plus the related PAR-bZIP cluster moving together at the top of the genome-wide screen, against a counter-moving set of GC-rich promoter motifs. The Part 5 hypotheses, AP-1 down and X-box up, were directionally reproduced by both methods but returned p = 0.40 and p = 0.88 across donors. With two donors per group, that is where the evidence stops.

Key Takeaways:

• Footprinting uses your cells’ reads; motif enrichment does not. That is what makes it an independent test of a motif hypothesis — and why CEBPB, invisible to Part 5’s 96 differential peaks, surfaced here across 11,167 sites.
• Signac and TOBIAS are complementary. Signac is fast, flexible and replicate-aware for a few motifs; TOBIAS screens every motif and predicts site-level binding.
• Controls earn their place by failing. CTCF produced the second-largest Signac change in this run. TOBIAS, the donor margin and the aggregate plot all said it was not real. Without a control motif, it would have been a finding.
• Direction is cheap; a replicate-level test is not. Five out of five motifs agreed in direction between two independent tools. Only one of the five survived a t-test across four donors.
• Test at the donor level, and say which test you used. Signac provides no differential footprinting function, so the summary statistic and the test are the analyst’s choice, and must be stated. Eyeballing whether donor values separate is not a test: with two per group they separate by chance one time in three, and here they separated for both controls and neither hypothesis.
• Families, not factors. Neither motif enrichment nor footprinting can tell apart proteins that bind the same DNA. BINDetect’s cluster column does that grouping for you — use it.
• When two methods disagree, look at the far flanks of the aggregate plot. Tracks that differ 100 bp from any motif differ for reasons that have nothing to do with binding.

Moving Forward

Everything in this tutorial is a function of the chosen cell type. The same workflow applies to any population with enough cells in both conditions — B cells and CD4 T cells are the obvious next candidates from Part 5 — by changing the subset in Step 4 and the barcode groups in Step 10. For a replicate-aware TOBIAS analysis, build one BAM per donor instead of per condition and pass all four to BINDetect, which compares every pair.

The biggest remaining limitation is that chromatin can only tell you what a cell is prepared to do. Whether AP-1 or RFX target genes are actually expressed differently in severe-disease monocytes is a question for the matched scRNA-seq from this study, and linking accessible regions to the genes they control is where the series goes next.

📚 References

1. Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nature Methods. 2013;10:1213-1218. doi:10.1038/nmeth.2688
2. Stuart T, Srivastava A, Madad S, Lareau CA, Satija R. Single-cell chromatin state analysis with Signac. Nature Methods. 2021;18:1333-1341. doi:10.1038/s41592-021-01282-5
3. Bentsen M, Goymann P, Schultheis H, et al. ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation. Nature Communications. 2020;11:4267. doi:10.1038/s41467-020-18035-1
4. Bentsen M, Goymann P, Schultheis H, et al. Beyond accessibility: ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation. bioRxiv. 2019. doi:10.1101/869560
5. Sung MH, Guertin MJ, Baek S, Hager GL. DNase footprint signatures are dictated by factor dynamics and DNA sequence. Molecular Cell. 2014;56(2):275-285. doi:10.1016/j.molcel.2014.08.016
6. Baek S, Goldstein I, Hager GL. Bivariate genomic footprinting detects changes in transcription factor activity. Cell Reports. 2017;19(8):1710-1722. doi:10.1016/j.celrep.2017.05.003
7. 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
8. Li Z, Schulz MH, Look T, Begemann M, Zenke M, Costa IG. Identification of transcription factor binding sites using ATAC-seq. Genome Biology. 2019;20:45. doi:10.1186/s13059-019-1642-2
9. Fornes O, Castro-Mondragon JA, Khan A, et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Research. 2020;48(D1):D87-D92. doi:10.1093/nar/gkz1001
10. Teo AYY, Squair JW, Courtine G, Skinnider MA. Best practices for differential accessibility analysis in single-cell epigenomics. Nature Communications. 2024;15:8805. doi:10.1038/s41467-024-53089-5
11. Squair JW, Gautier M, Kathe C, et al. Confronting false discoveries in single-cell differential expression. Nature Communications. 2021;12:5692. doi:10.1038/s41467-021-25960-2
12. Brauns E, Azouz A, Grimaldi D, et al. Functional reprogramming of monocytes in patients with acute and convalescent severe COVID-19. JCI Insight. 2022;7(9):e154183. doi:10.1172/jci.insight.154183
13. Giroux NS, Ding S, McClain MT, et al. Differential chromatin accessibility in peripheral blood mononuclear cells underlies COVID-19 disease severity prior to seroconversion. Scientific Reports. 2022;12:11714. doi:10.1038/s41598-022-15668-8
14. Signac documentation: Transcription factor footprinting vignette. https://stuartlab.org/signac/articles/footprint
15. Signac GitHub Discussion #968: TF footprinting analysis and interpretation (answer by the Signac maintainer). https://github.com/stuart-lab/signac/discussions/968
16. TOBIAS documentation wiki (ATACorrect, ScoreBigwig, BINDetect). https://github.com/loosolab/TOBIAS/wiki
17. samtools view manual, -D STR:FILE tag filtering. http://www.htslib.org/doc/samtools-view.html
18. 10x Genomics. Cell Ranger ATAC algorithms overview (duplicate marking and fragment filtering). https://software.10xgenomics.com/single-cell-atac/software/pipelines/2.1/algorithms/overview
19. Pixi documentation: Multi environment tutorial. https://pixi.prefix.dev/latest/tutorials/multi_environment/
20. Gene Expression Omnibus accession GSE282769. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE282769

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