How to Analyze Single-Cell ATAC-seq Data — A Complete Beginner’s Guide Part 7: Genome-Wide Differential TF Activity with chromVAR and limma

Meta description: Learn genome-wide differential transcription factor activity analysis for scATAC-seq: donor pseudobulk, chromVAR motif deviations and limma moderated t-tests, with R code.

Stop hand-picking motifs. Score all 633 JASPAR2020 transcription factor motifs in CD14 monocytes with chromVAR, test severe versus healthy at the level of donors with limma, and learn how to read the result next to the differential accessibility of Part 5 and the footprints of Part 6.

Part 6 ended with one motif family standing. C/EBP-family motifs gained footprint signal in severe-disease CD14 monocytes: CEBPB showed the largest change in both Signac and TOBIAS, and it was the only one of five hand-picked motifs that a donor-level t-test supported (+26.5 percent, p = 0.013, which became 0.064 after correction across the five motifs). It also ended with a limitation that this tutorial exists to remove. We chose five motifs because Signac footprinting is expensive, and TOBIAS, which did screen all 633 motifs, cannot see biological replicates at all.

So two questions were left open. Of all 633 motifs, which ones change between severe and healthy monocytes when donors, not cells or binding sites, are the unit of replication? And does C/EBP still stand out when 628 other motifs are tested beside it under the same rules?

That is genome-wide differential transcription factor (TF) activity analysis, and it has a standard, benchmarked recipe: score every motif in every sample with chromVAR, then test the scores across replicates with limma. This tutorial walks through that recipe step by step on the same CD14 monocytes, the same four donors and the same motif database you used in Parts 4 to 6, and then lines the answer up against everything the series has found so far.

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, Part 5: Differential Accessibility Analysis and Part 6: Footprinting Analysis. You need the atac_annotated.rds object from Part 4, the fragments files it points to (Part 1 and Part 2), the motif tables from Part 5, the result tables from Part 6, and the same Pixi workspace. No new raw data is downloaded, and everything in this tutorial runs in R.

🧭 Introduction: What Genome-Wide TF Activity Analysis Adds to scATAC-seq

What Is Differential TF Activity Analysis?

A transcription factor is a protein, and scATAC-seq never measures proteins. What it measures is where Tn5 could reach the DNA. The link between the two is indirect but powerful: when a transcription factor is active in a cell, the regions carrying its binding motif tend, as a group, to be more accessible. Its binding helps keep those regions open, recruits remodelers, and displaces nucleosomes. When the factor is absent or inactive, the same regions tend, as a group, to close.

TF activity analysis turns that idea into a number. For every motif, it pools the accessibility of all peaks that contain the motif and asks whether that pool is more or less accessible than expected in each sample. Differential TF activity analysis then asks whether that number differs between conditions — here, severe COVID-19 versus healthy donors.

Two features make the approach valuable:

1. It aggregates. A single peak in a 2-versus-2 design rarely carries enough signal to pass a statistical test. A motif present in 10,000 peaks pools 10,000 small, coordinated shifts into one score. A factor whose effect is spread thinly across the genome can be invisible peak by peak and obvious motif by motif.
2. It covers every motif. Instead of asking about a short list chosen in advance, you score the entire motif database and let the data rank the factors.

Beginner’s Tip: Think of each motif as a political party and each peak as a voter. Differential accessibility (Part 5) asks which individual voters changed their minds. TF activity analysis asks whether, overall, more supporters of a party turned up to vote this time than last time — even if no single voter’s change was dramatic.

A word on vocabulary before going further. Papers, packages and this tutorial all say “TF activity”, but what chromVAR actually measures is what its authors called transcription-factor-associated accessibility: the accessibility of chromatin that contains a factor’s motif (Schep et al. 2017). It is evidence about a TF, not a direct measurement of it, and every interpretation in this tutorial keeps that distinction.

Why Do We Need It? What Parts 5 and 6 Could Not Tell Us

Each of the last two parts answered a real question and left a gap that only a genome-wide, replicate-aware motif test can fill.

Part 5 found too few peaks to say much about motifs. The donor-level pseudobulk test with DESeq2 found 96 differentially accessible peaks in CD14 monocytes (46 opened, 50 closed in severe disease) out of 77,570 tested. Part 5 then ran motif enrichment on those 46 and 50 peaks and found RFX X-box motifs among the opened peaks and AP-1 motifs among the closed ones. Motif enrichment on such small sets is fragile, and it can only find motifs inside peaks that already reached significance one at a time.

Part 6 measured binding, but only for a handful of motifs, or without replicates. Signac footprinting was replicate-aware but limited to five motifs chosen by hand. TOBIAS screened all 633 motifs, but its p-values resample binding sites, not donors — its documentation recommends merging replicates before footprinting, so the donor dimension disappears.

The gap in one sentence: nothing so far has tested every motif with the donor as the unit of replication. That is exactly what chromVAR plus limma does.

Part 6 also gave us a concrete reason to want it. CEBPB had p.adjust = 1 in both directions in Part 5’s motif enrichment, yet footprinting picked it up across 11,167 motif sites. That is the signature of a change spread thinly across many sites — the situation TF activity analysis is designed for.

How chromVAR Turns Accessibility Into a TF Activity Score

chromVAR (Schep et al. 2017) computes its score in four steps. Each is worth understanding, because each one shapes how you read the results.

1. Expected accessibility. For each peak, chromVAR computes the fraction of all fragments, across all samples, that fall in that peak. The expected count for a sample is that fraction multiplied by the sample’s total fragments in peaks. This corrects for sequencing depth: a deeper sample is expected to have more fragments everywhere.
2. Raw deviation. For a motif, chromVAR adds up the observed fragments in all peaks containing the motif, subtracts the expected count, and divides by the expected count. A positive value means “the peaks with this motif are more accessible in this sample than its depth predicts.”
3. Background peak sets. Some peaks are accessible in every sample simply because of their GC content or because they are strongly open everywhere. To remove that effect, chromVAR builds many background peak sets, each with the same number of peaks as the motif set and matched on GC content and average accessibility, and computes the same deviation for each.
4. Bias-corrected deviation and z-score. The motif’s deviation minus the mean background deviation is the bias-corrected deviation. Dividing that by the standard deviation of the background deviations gives the deviation z-score — how many background standard deviations the motif sits away from what sequence composition alone would produce.

The z-score is what we test. Two consequences follow directly from the definitions above:

• Scores are relative to the samples you give chromVAR. The expectation in step 1 is built from those samples, so a z-score says how one sample compares with the average of the set, not how accessible a motif is in absolute terms. Never compare z-scores between separate chromVAR runs.
• The number of background sets controls how precise the z-score is. The default in getBackgroundPeaks() is 50 background sets. That choice becomes important in a moment.

Differential Accessibility, Footprinting and TF Activity: When to Use Which

All three analyses start from the same fragments and the same motif database, which makes them easy to confuse. They answer different questions.

When to use which:

• Start with differential accessibility (Part 5) when your question is about loci: which enhancers or promoters change, and near which genes.
• Use TF activity analysis (this tutorial) when your question is about regulators: which transcription factor families are associated with the change, across the whole genome, with a replicate-level statistic. It is the right screen when you have no strong prior about which factors matter.
• Use footprinting (Part 6) to follow up a short list — typically the top families from this tutorial — and ask whether the change looks like a change in binding rather than a change in accessibility around the motif.

In practice the three run in this order of increasing specificity: peaks, then motifs genome-wide, then footprints at selected motifs. This series ran footprinting before TF activity because Part 5’s motif hypotheses were already on the table. Part 7 closes the loop by screening everything.

Why Test Donors, Not Cells: The Benchmark Behind This Workflow

chromVAR was designed for single cells, and per-cell z-scores are useful: coloring a UMAP by motif activity shows which clusters use which factors. (Signac versions before 1.17.0 offered RunChromVAR() for exactly that; it was removed in 1.17.0, which is why this tutorial calls chromVAR directly.) But testing per-cell scores between conditions counts every cell as an independent replicate — the pseudoreplication problem from Part 5. Instead, we sum one cell type’s fragments per donor, which turns four donors’ CD14 monocytes into four small “libraries,” and test across those. Two benchmarks support that choice.

Teo et al. (2024) compared differential accessibility methods for single-cell epigenomic data and found that pseudobulk methods, which aggregate cells per sample before testing, consistently ranked at the top, while single-cell methods produced large numbers of false discoveries in null comparisons. Part 5 is built on that result. Their benchmark is about peaks, not motifs.

Gerbaldo et al. (2024) benchmarked methods for finding differentially active transcription factors from ATAC-seq, using datasets in which known factors had been perturbed plus semi-simulated data, with 2 to 3 samples per group — the same scale as our two donors per condition. Their findings set every parameter in this tutorial:

• chromVAR followed by limma was one of the three most promising approaches, alongside monaLisa and a GC-aware normalization combined with multivariate modeling.
• Testing chromVAR z-scores with limma clearly outperformed chromVAR’s own differentialDeviations() test. The authors attribute this to limma’s moderation of the variance, which matters most when each group has only two or three samples.
• The default number of background sets is too small for reproducible results. Differential results changed with the random seed until the number of background iterations was raised substantially; they ran all comparisons with 2000 and recommend 1000 to 2000. They also note that this is not easy with single-cell data — one more reason to work on pseudobulks, where 2000 iterations over four columns is cheap.
• Quantile normalization of the z-scores before testing performed best among the downstream normalizations they tried.
• Peaks were resized to 300 bp, which they describe as recommended by chromVAR and close to their median peak size.
• P-values were not well calibrated: the observed false discovery rate was nearly always higher than the nominal threshold. Treat adjusted p-values as a ranking with a rough threshold, not as an exact error rate.
• They recommend running more than one of the top methods to increase confidence.

An honest gap. No benchmark has yet evaluated this exact combination — cell-type pseudobulk by donor, then chromVAR and limma — on scATAC-seq data. Gerbaldo et al. is mostly built on bulk data, and Teo et al. tests peaks rather than motifs. The workflow below is the best-supported option assembled from both, and the seed check in Step 11 and the cross-method comparison in Step 13 exist precisely because no benchmark has done that checking for us.

The Part 7 Workflow

atac_annotated.rds (Part 4)
        |
        v
[1] Subset CD14 monocytes                                    Step 4
        |
        v
[2] Monocyte-accessible peaks, resized to 300 bp             Step 5
        |
        v
[3] Count fragments per peak per cell, sum per donor         Step 6
    -> peaks x 4 donors
        |
        v
[4] GC content, overlap filter                               Step 7
[5] Match 633 JASPAR2020 motifs to the resized peaks         Step 8
        |
        v
[6] chromVAR: 2000 GC- and accessibility-matched             Step 9
    background sets -> deviation z-scores, 633 motifs x 4 donors
        |
        v
[7] Quantile-normalize z-scores, limma moderated t-test      Step 10
    (Severe vs Healthy, donors as replicates)
        |
        +--> [8]  Seed reproducibility check                 Step 11
        +--> [9]  Collapse motifs into families              Step 12
        +--> [10] Compare with Part 5, Part 6 and motif GC   Step 13
        |
        v
[11] Visualization: volcano, donor heatmap, donor dots,      Steps 14-18
     chromVAR vs TOBIAS, seed reproducibility

⚙️ Setting Up the Analysis Environment

What Is New in Part 7

Almost everything this tutorial needs is already in the Pixi workspace. chromVAR, motifmatchr, TFBSTools, JASPAR2020, BSgenome.Hsapiens.UCSC.hg38 and SummarizedExperiment were installed for Part 4’s motif analysis (the Part 4 session information records chromVAR 1.32.0). Only one package is new: limma – Linear models with empirical Bayes variance moderation — the replicate-level test.

#-----------------------------------------------
# STEP 1: Add limma to the Pixi workspace
#-----------------------------------------------

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

pixi add bioconductor-limma

# Confirm every package this tutorial loads is visible to R
pixi run Rscript -e 'sapply(c("chromVAR", "motifmatchr", "limma", "TFBSTools", "JASPAR2020", "BSgenome.Hsapiens.UCSC.hg38"), requireNamespace, quietly = TRUE)'

Expected output: all six names followed by TRUE.

If JASPAR2020 or BSgenome.Hsapiens.UCSC.hg38 returns FALSE, you have met the Bioconda stub problem described in Part 4: large Bioconductor data packages are shipped as small stubs whose post-link scripts download the real data, and Pixi does not run those scripts. Apply the Part 4 fix before continuing. limma is ordinary software and is not affected.

Loading Libraries and Creating the Part 7 Directory Structure

This tutorial runs entirely in R, in one session. Launch R in the default environment, preferably on a compute node with several CPUs and generous memory (see the note after Step 2):

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

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

library(Signac)
library(Seurat)
library(chromVAR)                     # motif deviation z-scores
library(motifmatchr)                  # scan the resized peaks for motifs
library(JASPAR2020)                   # the same motif database as Parts 4-6
library(TFBSTools)                    # getMatrixSet() to pull motifs from JASPAR2020
library(BSgenome.Hsapiens.UCSC.hg38)  # genome sequence for GC content and motif scanning
library(SummarizedExperiment)
library(GenomicRanges)
library(Matrix)
library(limma)                        # moderated t-test across donors
library(BiocParallel)
library(ggplot2)
library(ggrepel)
library(patchwork)

set.seed(1234)

# computeDeviations() runs motifs in parallel; match workers to the CPUs you requested
register(MulticoreParam(workers = 4))

Why the resource request matters. For every motif, chromVAR builds a sparse matrix with one row per background set and one entry for every peak containing the motif. With 2000 background sets and motifs present in tens of thousands of peaks, that matrix is large, and each parallel worker builds its own. We suggest requesting at least 64 GB of memory for four workers. If the session is killed during Step 9 or Step 11, lower workers, or use register(SerialParam()) to run one motif at a time — slower, but with the smallest memory footprint.

#-----------------------------------------------
# STEP 3: Define and create the Part 7 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")  # Part 6 output
tfa_dir   <- file.path(base_dir, "tf_activity")   # this tutorial

obj_dir   <- file.path(tfa_dir, "objects")
plot_dir  <- file.path(tfa_dir, "plots")
table_dir <- file.path(tfa_dir, "tables")

invisible(lapply(c(obj_dir, plot_dir, table_dir),
                 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)
├── downstream_qc/         # per-sample QC objects (Part 2)
├── integrated/            # consensus peaks, integrated object (Part 3)
├── annotation/            # annotated object (Part 4)        <- input
├── differential/          # DA results and motif tables (Part 5)   <- comparison
├── footprinting/          # Signac and TOBIAS results (Part 6)     <- comparison
└── tf_activity/           # this tutorial
    ├── objects/
    ├── plots/
    └── tables/

📦 Example Data: Choosing the Right Annotated Object

Which Output From Parts 4 to 6 Should You Start From?

Three earlier tutorials left objects on disk, and it is worth being deliberate about which one feeds this analysis.

The Part 4 object is the cleanest starting point: it is the common ancestor of Parts 5 and 6, so all three analyses describe exactly the same cells. We need the fragments files it points to, because Step 6 recounts fragments in resized peaks. The Part 5 and Part 6 result tables are read only at the end, to see whether the three approaches agree.

Loading the Object and Selecting CD14 Monocytes

We analyze CD14 monocytes for the third time running. They are the largest well-annotated population, both conditions contribute meaningful numbers of cells, and Parts 5 and 6 have already produced hypotheses for them.

#-----------------------------------------------
# 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 fragments 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. 5,213 cells, four donors, two per condition — the same design as Parts 5 and 6. The cell numbers are unbalanced (healthy2 contributes 2,575 cells, severe2 only 508), and this tutorial handles that the right way: each donor becomes one column in the analysis, whatever its cell count. A donor with more cells gives a less noisy pseudobulk, but it does not get more votes.

🔬 Differential TF Activity Analysis: Severe vs Healthy CD14 Monocytes

The analysis has three stages: build a peak-by-donor count matrix (Steps 5 to 7), turn it into motif z-scores (Steps 8 and 9), and test them (Step 10). Steps 11 to 13 then check how far the result can be trusted.

Step 5: Choosing Monocyte Peaks and Resizing Them to 300 bp

chromVAR works best with fixed-width peaks. The consensus peaks from Part 3 vary in width, and most are about three times wider than 300 bp (you will see the exact numbers in a moment). A 900 bp peak carries three times as many chances of containing any given motif as a 300 bp peak, and much of that extra sequence lies away from the accessible centre. Resizing every peak to 300 bp around its midpoint follows Gerbaldo et al. (2024) and puts all peaks on an equal footing.

We start from the peaks actually open in CD14 monocytes rather than all 185,581 consensus peaks, for the same reason as Part 6: peaks that only other cell types open would add noise and nothing else.

#-----------------------------------------------
# STEP 5: Choose the peaks and resize them to 300 bp
#-----------------------------------------------

# Peaks detected in at least 10 CD14 monocytes -- the same set Part 6 gave TOBIAS
open_peaks <- AccessiblePeaks(mono, idents = c("Healthy", "Severe"))
peaks_gr   <- sort(granges(mono)[rownames(mono) %in% open_peaks])

# Fixed-width peaks centred on the original peak midpoint (Gerbaldo et al. 2024)
peaks_300 <- resize(peaks_gr, width = 300, fix = "center")

length(peaks_300)
summary(width(peaks_gr))

Output:

[1] 115324

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  200.0   866.0   913.0   913.6   952.0  2598.0

Reading this output. 115,324 peaks — 62 percent of the consensus set, and exactly the peak set Part 6 handed to TOBIAS, because the same function with the same arguments was applied to the same cells. Keeping the two analyses on an identical set of regions makes the comparison in Step 13 fair.

The width summary shows why resizing matters. Half of the monocyte peaks are between 866 and 952 bp wide, with a median of 913 bp — about three times the 300 bp window we resize to. A few are much longer (up to 2,598 bp) and a few shorter (down to 200 bp). Resizing trims the typical peak to its central third and stretches the shortest ones slightly, so every peak contributes the same amount of sequence to the motif scan in Step 8.

Why the midpoint? The ideal centre is the peak summit — the single base with the most Tn5 insertions. The consensus peaks from Part 3 do not carry a summit position, so the midpoint is the best available choice, and the one Gerbaldo et al. (2024) used when resizing.

Step 6: Counting Fragments in the New Peaks and Building a Donor Pseudobulk

Resized peaks need new counts. Signac’s FeatureMatrix() reads the fragments files and counts, for every cell, the fragments overlapping each region. We then sum the cells of each donor, exactly like the pseudobulk in Part 5, but on the 300 bp peaks.

#-----------------------------------------------
# STEP 6: Count fragments in the resized peaks and sum them per donor
#-----------------------------------------------

# Re-reads the fragments files: one row per resized peak, one column per cell
cell_counts <- FeatureMatrix(
  fragments = Fragments(mono),
  features  = peaks_300,
  cells     = colnames(mono)
)

# Pseudobulk: add up the cells of each donor -> peaks x 4 donors
donor_of_cell <- mono$sample_id[colnames(cell_counts)]
pb_counts <- sapply(split(colnames(cell_counts), donor_of_cell),
                    function(cells) Matrix::rowSums(cell_counts[, cells, drop = FALSE]))

dim(pb_counts)
colSums(pb_counts)

# Keep the pseudobulk; the per-cell matrix is no longer needed
saveRDS(pb_counts, file.path(obj_dir, "cd14mono_pseudobulk_counts_300bp.rds"))
rm(cell_counts)
invisible(gc())

Output:

[1] 115324      4

healthy1 healthy2  severe1  severe2
 3896051  9619830  5524531  1801008

Reading this output. One row per resized peak, one column per donor, with columns in alphabetical order. colSums() gives each donor’s total fragments in the 300 bp peaks: 3.9, 9.6, 5.5 and 1.8 million. The ranking of donors is the same as in Part 5’s pseudobulk (4.84, 11.98, 7.02 and 2.36 million peak counts), and the totals are lower because the windows are narrower.

Even after trimming every peak to a third of its width, the pseudobulks keep most of their fragments. That is because Tn5 cuts concentrate near the centre of an accessible region, which is exactly the part the resizing keeps.

The donors differ more than fivefold in depth. severe2 has 1.8 million fragments, healthy2 9.6 million. That follows directly from the cell numbers in Step 4 (508 versus 2,575 monocytes).

Unequal depth is fine here. chromVAR’s expected counts are each donor’s total multiplied by each peak’s overall share, so a donor with three times the fragments is expected to have three times the counts in every peak. What chromVAR tests is the distribution of fragments across peaks, not the total.

Step 7: Building the chromVAR Input and Adding GC Content

chromVAR works with a SummarizedExperiment: a count matrix packaged together with the genomic coordinates of its rows and a table describing its columns. We add each peak’s GC content, which chromVAR uses to choose background peaks.

#-----------------------------------------------
# STEP 7: Build the chromVAR input and add GC content
#-----------------------------------------------

donor_info <- data.frame(
  sample_id = colnames(pb_counts),
  condition = factor(ifelse(grepl("^severe", colnames(pb_counts)), "Severe", "Healthy"),
                     levels = c("Healthy", "Severe")),
  row.names = colnames(pb_counts)
)

se <- SummarizedExperiment(
  assays    = list(counts = pb_counts),
  rowRanges = GRanges(sub("-", ":", rownames(pb_counts))),  # "chr1-100-400" -> "chr1:100-400"
  colData   = donor_info
)

# GC content of every peak -- used to pick matched background peaks
se <- addGCBias(se, genome = BSgenome.Hsapiens.UCSC.hg38)

# Drop peaks with no fragments, and keep one peak where resized peaks overlap
se <- filterPeaks(se, non_overlapping = TRUE)

dim(se)
summary(rowData(se)$bias)

Output:

[1] 115305      4

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 0.2267  0.4367  0.4967  0.5149  0.5767  0.9233

Reading this output. filterPeaks() removes peaks with no fragments in any donor and, where two resized peaks overlap, keeps the one with more fragments. Only 19 of 115,324 peaks were dropped, leaving 115,305. That is expected: every peak was detected in at least 10 monocytes, so the losses are a handful of neighbouring peaks whose 300 bp windows now overlap. The 115,305 rows are the peak set for everything that follows.

The bias column is GC content, the fraction of G and C bases in each 300 bp peak. Here it spans 23 to 92 percent, with half the peaks between 44 and 58 percent. That range is why GC matching matters: a motif found mostly in GC-rich peaks (CpG-island promoters, for example) is compared in Step 9 with background peaks of similar GC content, not with accessible chromatin at large.

Setting factor levels matters. levels = c("Healthy", "Severe") makes Healthy the reference group, so every effect in this tutorial reads as Severe minus Healthy — the same direction convention as Parts 5 and 6.

Step 8: Scanning the Resized Peaks for the 633 JASPAR2020 Motifs

Motif matches stored in the Part 4 object refer to the original, variable-width peaks, so they have to be redone on the 300 bp peaks. We retrieve the motifs with the exact query used in Part 4, which guarantees the same 633 matrices as Parts 4, 5 and 6.

#-----------------------------------------------
# STEP 8: Scan the resized peaks for the 633 JASPAR2020 motifs
#-----------------------------------------------

# Identical query to Part 4, so the motif universe matches Parts 4-6
pfm <- getMatrixSet(
  x    = JASPAR2020,
  opts = list(species = 9606, collection = "CORE", all_versions = FALSE)
)

motif_ix <- matchMotifs(pfm, se, genome = BSgenome.Hsapiens.UCSC.hg38)

# JASPAR ID -> transcription factor name, for readable tables
motif_names <- setNames(TFBSTools::name(pfm), TFBSTools::ID(pfm))

dim(motif_ix)
summary(colSums(motifMatches(motif_ix)))

Output:

[1] 115305    633

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    312    2137    3934    5051    6087   31100

Reading this output. motif_ix is a peaks-by-motifs table of TRUE/FALSE: does this 300 bp peak contain at least one match to this motif? matchMotifs() uses the same default match threshold (p < 5e-05) that Signac’s AddMotifs() used in Part 4. The second line summarizes how many peaks each motif hits.

Every motif has hundreds to tens of thousands of peaks to pool. The rarest hits 312 peaks, the median motif 3,934, and the most frequent 31,100. Pooling is what gives TF activity analysis its power, and even the rarest motif here pools more peaks than the 96 differential peaks of Part 5 combined.

Why the counts differ from Part 6. Part 6 reported, for example, 117,781 SPI1 matches in the original peaks. Matches in 300 bp peaks will be fewer, because a third of the sequence is being scanned, and matchMotifs() also computes its background nucleotide frequencies from the sequences it scans. Both are expected and both affect all motifs alike.

Step 9: Computing chromVAR Deviation Z-Scores With 2000 Background Sets

This is the step Gerbaldo et al. (2024) singled out. The default of 50 background sets gives z-scores that change noticeably from one random seed to the next, which made their differential results irreproducible. We use 2000.

#-----------------------------------------------
# STEP 9: Background peaks and motif deviation z-scores
#-----------------------------------------------

# 2000 background sets matched on GC content and average accessibility
set.seed(1234)
bg_peaks <- getBackgroundPeaks(se, niterations = 2000)

dev <- computeDeviations(se, annotations = motif_ix, background_peaks = bg_peaks)

# One z-score per motif per donor: rows are JASPAR IDs, columns are donors
z <- deviationScores(dev)

dim(z)
round(head(z), 2)

saveRDS(dev, file.path(obj_dir, "cd14mono_chromvar_deviations.rds"))

Output:

[1] 633   4

         healthy1 healthy2 severe1 severe2
MA0030.1     1.80    -2.26    1.90   -3.26
MA0031.1     1.64    -2.14    1.34   -1.41
MA0051.1    -5.50    -3.33    3.42    4.52
MA0057.1     1.15    -0.04   -0.39   -0.14
MA0059.1    -1.04    -0.21    0.64   -0.17
MA0066.1    -1.92    -2.40    3.06   -1.83

Reading this output. A 633-by-4 matrix: one deviation z-score per motif per donor. A z of +5 for a donor means the peaks containing that motif are five background standard deviations more accessible in that donor than the average of the four donors would predict, after matching for GC content and overall accessibility.

Look at the signs across a row. Because each donor is compared with the average of all four, the z-scores of a motif tend to balance out across the row. The six rows printed here are simply the first six motifs by JASPAR ID, and they show the three patterns you will meet:

• A condition pattern. MA0051.1 (IRF2) is negative in both healthy donors (-5.50, -3.33) and positive in both severe donors (3.42, 4.52). This is what a motif that is truly higher in severe disease looks like.
• A donor pattern. MA0030.1 (FOXF2) alternates: positive in healthy1 and severe1, negative in healthy2 and severe2. That split follows the donors, not the conditions, and should not be called.
• A single-donor pattern. MA0066.1 (PPARG) is driven by severe1 alone (3.06, with the other severe donor at -1.83).

Values span roughly -21 to +20 across the whole matrix. The z-score divides by the spread of the background deviations, and that spread shrinks as fragment counts grow, so pseudobulks with millions of fragments give larger z-scores than a single cell would. The size of a single z-score is not what we test; the difference between the two groups, relative to donor-to-donor variation, is. No z-score is NA: every motif had enough expected fragments in every donor.

Step 10: Testing Severe vs Healthy With limma Moderated t-Statistics

With two donors per group, each motif’s ordinary t-test would rest on two degrees of freedom — exactly the instability Part 6 ran into, where one motif’s healthy standard deviation was ten times its severe one by chance. limma solves this by borrowing strength across motifs. Smyth (2004) showed how to shrink each feature’s variance toward a value estimated from all features, which gives a moderated t-statistic whose degrees of freedom are augmented by the information from all the others. How much it borrows is estimated from the data; Step 10’s output lets us work it out for this run.

We first quantile-normalize the z-scores, as recommended by Gerbaldo et al. (2024). Quantile normalization forces the four donors’ z-score distributions to be identical, so a donor whose scores are globally stretched or shifted cannot drive the result.

#-----------------------------------------------
# STEP 10: Quantile-normalize the z-scores and test with limma
#-----------------------------------------------

# Wrapped in a function so Step 11 can rerun it on a second background draw
test_tf_activity <- function(z, condition) {
  z_qn   <- normalizeQuantiles(z)                 # same distribution in every donor
  design <- model.matrix(~ condition)             # intercept = Healthy, slope = Severe - Healthy
  fit    <- eBayes(lmFit(z_qn, design))           # moderated t: variances shrunk across motifs
  res    <- topTable(fit, coef = "conditionSevere", number = Inf, sort.by = "P")
  res$motif_id <- rownames(res)
  res$motif    <- motif_names[res$motif_id]       # readable TF name next to the ID
  res
}

tf_res <- test_tf_activity(z, donor_info$condition)

head(tf_res[, c("motif", "logFC", "AveExpr", "t", "P.Value", "adj.P.Val")], 10)
sum(tf_res$adj.P.Val < 0.05)
table(direction = sign(tf_res$logFC[tf_res$adj.P.Val < 0.05]))

Output:

                      motif      logFC     AveExpr         t      P.Value  adj.P.Val
MA0488.1                JUN  16.395259 -0.74153333 12.680703 0.0001549246 0.09806724
MA0639.1                DBP  11.068870 -0.28810159  8.578428 0.0007725259 0.10064404
MA0605.2               ATF3  12.688800  0.08571070  7.989358 0.0010297165 0.10064404
MA0492.1        JUND(var.2)  15.750800  0.49930077  7.508436 0.0013209379 0.10064404
MA0838.1              CEBPG  14.938510 -1.27896513  7.158655 0.0015976028 0.10064404
MA0837.1              CEBPE  15.523027 -1.60812890  7.099915 0.0016507556 0.10064404
MA0466.2              CEBPB  15.908684 -1.60624290  6.993309 0.0017528427 0.10064404
MA0507.1             POU2F2 -12.888611 -0.52604027 -6.981138 0.0017649817 0.10064404
MA0018.4              CREB1  11.827106  0.13164685  6.979340 0.0017667835 0.10064404
MA1143.1 FOSL1::JUND(var.2)   8.807854  0.05560247  6.875712 0.0018745293 0.10064404

[1] 0

limma was written for gene expression, so some column names need translating:

Start with the headline: no motif passes adj.P.Val < 0.05. The [1] 0 says so, and the empty table on the last line is what R prints when there is nothing to count. The smallest adjusted p-value is 0.098, for JUN, and 20 more motifs sit just above it at 0.1006 to 0.1010. Benjamini-Hochberg adjustment often gives a run of motifs at the top of the list the same value, which is why so many read 0.1006. Say this plainly in a write-up: at a 5 percent false discovery rate, this 2-versus-2 design finds no differentially active motif.

Why so close and yet not there? JUN has a moderated t of 12.7 — the four donors separate cleanly — yet its p-value is only 1.5e-4. With two donors per group, each motif has 2 residual degrees of freedom. limma’s moderation added the equivalent of about 2.2 more, borrowed from the other motifs (back-calculated from JUN’s t and p). Even with about 4.2 in total, a t of 12.7 cannot get much below p = 1e-4. Multiplying by 633 tests, as Benjamini-Hochberg roughly does at the top of the list, lands at 0.1. With two donors per group, a motif needs an almost perfect separation just to reach the edge of significance.

The unadjusted p-values show there is signal. 95 motifs have p < 0.05 and 40 have p < 0.01. If no motif changed at all, you would expect about 32 and 6 by chance. And the excess is lopsided: 35 of those 40 motifs are higher in severe disease. Something real is shifting, even though no single motif clears the multiple-testing bar.

Nine of the top ten are bZIP motifs from three groups, not nine separate factors.

• C/EBP: CEBPG, CEBPE and CEBPB, ranks 5 to 7.
• bZIP motifs with “AP-1” or CREB/ATF names: JUN, ATF3, JUND(var.2), CREB1 and FOSL1::JUND(var.2).
• PAR bZIP: DBP.

The only motif higher in healthy donors is the octamer motif POU2F2. Step 12 sorts these into families properly.

CEBPB, Part 6’s one surviving motif, ranks 7th of 633. It is higher in severe disease (logFC +15.9, t = 7.0, p = 0.0018), and so are its relatives. That is the direction Part 6 found by footprinting, now from a different input (peak counts, not insertion profiles) and a different test (moderated t across 633 motifs, not a t-test across five). Like Part 6’s p = 0.013 (BH 0.064 across five motifs), it does not survive correction (adj.P = 0.10).

Keep the calibration warning in mind. Gerbaldo et al. (2024) found the real false discovery rate of these tests was nearly always higher than the nominal one. Adjusted p-values around 0.1 are therefore a ranking of candidates, not a finding with a known error rate.

Step 11: Checking That the Result Does Not Depend on the Random Seed

getBackgroundPeaks() samples background peaks at random. Gerbaldo et al. (2024) found that, at default settings, which motifs came out significant depended on that randomness. That is why Step 9 used 2000 background sets instead of the default 50. Step 11 keeps the same 2000 and changes only the random seed: it draws a completely new background and reruns the test, to check that 2000 was enough for our data.

#-----------------------------------------------
# STEP 11: Check that the result does not depend on the random seed
#-----------------------------------------------

set.seed(2026)
bg_peaks_2 <- getBackgroundPeaks(se, niterations = 2000)
z_2 <- deviationScores(computeDeviations(se, annotations = motif_ix,
                                         background_peaks = bg_peaks_2))
tf_res_2 <- test_tf_activity(z_2, donor_info$condition)

# Same motifs, same order, then compare
seed_cmp <- data.frame(motif_id = tf_res$motif_id,
                       motif    = tf_res$motif,
                       t_seed1  = tf_res$t,
                       t_seed2  = tf_res_2$t[match(tf_res$motif_id, tf_res_2$motif_id)])

cor(seed_cmp$t_seed1, seed_cmp$t_seed2, method = "spearman")
length(intersect(head(tf_res$motif_id, 20), head(tf_res_2$motif_id, 20)))

Output:

[1] 0.99966
[1] 19

Reading this output. The first number is the rank correlation of all 633 moderated t-statistics between the two background draws: 0.99966, effectively identical. The second says 19 of the top 20 motifs are shared.

The one difference is a swap at the very edge of the list. IRF7 is 20th with seed 1234, and CEBPD takes its place with seed 2026. The largest change in any motif’s t-statistic between the two runs is 0.68 (JUND(var.2), 7.51 versus 8.19). Nothing in the top ten changes rank enough to matter. With 2000 background sets, the ranking is a property of the data, not of the random number generator.

If your correlation is clearly below 1, or the top 20 overlap poorly, the background is still too small for your data. Increase niterations (memory permitting) and repeat, and report only motifs that are stable across draws. This one check is the most direct protection against the irreproducibility Gerbaldo et al. (2024) documented, and it costs one extra run.

Step 12: Grouping Motifs Into Families

JASPAR2020 contains many near-identical motifs, and a factor family that changes will drag every member up the list together. Gerbaldo et al. (2024) note that many motif models in the collection are highly similar to each other. Part 6’s BINDetect run already solved the grouping problem for us: its cluster column groups all 633 motifs by how much their predicted binding sites overlap in the genome. We reuse it here.

#-----------------------------------------------
# STEP 12: Group significant motifs into families
#-----------------------------------------------

# Part 6 BINDetect clustered all 633 motifs by binding-site overlap
bd <- read.csv(file.path(fp_dir, "tables", "tobias_bindetect_results.csv"))
tf_res$cluster <- bd$cluster[match(tf_res$motif_id, bd$motif_id)]

# One row per family: its best motif, and whether the other members agree
family_tab <- do.call(rbind, lapply(split(tf_res, tf_res$cluster), function(d) {
  best <- d[which.min(d$P.Value), ]
  data.frame(cluster     = best$cluster,
             n_motifs    = nrow(d),
             n_sig       = sum(d$adj.P.Val < 0.05),
             same_dir    = mean(sign(d$logFC) == sign(best$logFC)),
             top_motif   = best$motif,
             top_logFC   = best$logFC,
             top_adjP    = best$adj.P.Val)
}))
family_tab <- family_tab[order(family_tab$top_adjP), ]
rownames(family_tab) <- NULL

nrow(family_tab)
head(family_tab, 15)

Output:

[1] 285

                 cluster n_motifs n_sig same_dir          top_motif  top_logFC   top_adjP
1    C_FOSL1::JUN(var.2)       18     0        1                JUN  16.395259 0.09806724
2                C_CEBPA        5     0        1              CEBPG  14.938510 0.10064404
3         C_CEBPG(var.2)        2     0        1       CEBPG(var.2)  20.119877 0.10064404
4                  C_FEV       11     0        1                ERF   8.798639 0.10064404
5   C_FOSL1::JUND(var.2)        1     0        1 FOSL1::JUND(var.2)   8.807854 0.10064404
6                 C_IRF1        5     0        1               IRF9   8.474150 0.10064404
7                 C_IRF3        1     0        1               IRF3   7.878578 0.10064404
8                C_NFIL3        4     0        1                DBP  11.068870 0.10064404
9               C_POU2F3        3     0        1             POU2F2 -12.888611 0.10064404
10                C_ETV2        2     0        1               ETV2   8.291050 0.12429816
11      C_BHLHA15(var.2)        6     0        1     BHLHE22(var.2)  -6.689925 0.12938399
12                C_IRF2        1     0        1               IRF2   6.673528 0.13249943
13                C_EBF1        2     0        1               EBF1  -8.710232 0.13344205
14               C_GMEB2        1     0        1              GMEB2   5.253666 0.18291352
15               C_NR4A1        1     0        1              NR4A1 -12.263105 0.19082608

Reading this table. Count findings in this table, not in Step 10’s. The 633 motifs collapse into 285 families. With n_sig at 0 everywhere, same_dir carries the information: a value of 1 means every member moves the way the best member does. All fifteen families shown here have same_dir = 1.

On the severe side:

• The C/EBP cluster, C_CEBPA. All five members point up: CEBPA, CEBPB, CEBPD, CEBPE and CEBPG, with logFC +14.9 to +20.1 and p 0.0016 to 0.0046. Its sister cluster C_CEBPG(var.2) (CEBPG(var.2) and ATF4) is up too.
• The PAR-bZIP cluster, C_NFIL3. DBP, TEF, HLF and NFIL3 are all up. These are the same two clusters that topped the severe side of Part 6’s BINDetect screen.
• Two interferon regulatory factor clusters, C_IRF1 and C_IRF3. IRF9, IRF7, IRF8 and IRF4 in C_IRF1 are all up with p < 0.012; IRF1 itself is up but weaker (p = 0.16).
• An ETS cluster, C_FEV. All 11 members are up, led by ERF.

On the healthy side: octamer (POU) motifs (C_POU2F3, and C_POU1F1 just below the top 15), EBF motifs (EBF1, EBF3), and single-motif clusters such as NR4A1. EBF1 was one of the B-cell lineage motifs in Part 4’s annotation panel. PAX5, the other B-cell motif in that panel, does not move here (logFC +1.8, adj.P = 0.60), so the EBF and POU signal is not simply a B-cell signature. Still, a lineage-type motif falling in one condition is worth checking against annotation purity before you interpret it.

Name families, not proteins. Write “C/EBP-family motifs” or cite the BINDetect cluster, not “CEBPB activity increases.” Motif analysis cannot tell family members apart; deciding which member is responsible needs expression data or protein-level evidence.

If you skipped TOBIAS in Part 6, the cluster column will be empty. You can still group by eye using the TF names and the JASPAR family annotation, but the BINDetect clusters are data-driven and, can catch distinctions that names hide.

Step 13: Comparing Part 7 With Parts 5 and 6, and Checking for GC Effects

This is the payoff of running three analyses on the same cells. For the five motifs footprinted in Part 6 we can now line up four independent answers, and for all 633 motifs we can compare chromVAR with TOBIAS. Finally, we ask one diagnostic question that Part 6 raised: is the ranking secretly a GC-content ranking?

#-----------------------------------------------
# STEP 13: Line up Part 7 against Parts 5 and 6
#-----------------------------------------------

# Part 6: per-donor Signac flank enrichment for the five footprinted motifs.
# Recompute Part 6's donor t-test (Student's t, 2 vs 2) from those values.
flank_donor <- read.csv(file.path(fp_dir, "tables", "signac_flank_by_donor.csv"))
fp_ttest <- do.call(rbind, lapply(split(flank_donor, flank_donor$feature), function(d) {
  h <- d$corrected[d$condition == "Healthy"]
  s <- d$corrected[d$condition == "Severe"]
  data.frame(feature = d$feature[1],
             pct     = 100 * (mean(s) - mean(h)) / mean(h),
             p       = t.test(s, h, var.equal = TRUE)$p.value)
}))

# Part 6: TOBIAS differential binding score for all 633 motifs
tf_res$tobias_delta <- bd$tobias_delta[match(tf_res$motif_id, bd$motif_id)]

# Part 5: motif enrichment in peaks that opened or 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"))

five <- tf_res[match(fp_ttest$feature, tf_res$motif),
               c("motif", "logFC", "t", "P.Value", "adj.P.Val", "tobias_delta")]
five$signac_pct        <- fp_ttest$pct
five$signac_p          <- fp_ttest$p
five$part5_padj_opened <- m_up$p.adjust[match(five$motif, m_up$motif.name)]
five$part5_padj_closed <- m_dn$p.adjust[match(five$motif, m_dn$motif.name)]
five

# Genome-wide: does chromVAR agree with TOBIAS across all 633 motifs?
cor.test(tf_res$logFC, tf_res$tobias_delta, method = "spearman")

# Is the ranking driven by motif GC content? (fraction of C + G in each matrix)
motif_gc <- sapply(as.list(pfm), function(m) {
  x <- TFBSTools::Matrix(m)
  sum(x[c("C", "G"), ]) / sum(x)
})
tf_res$motif_gc <- motif_gc[tf_res$motif_id]
cor.test(tf_res$logFC, tf_res$motif_gc, method = "spearman")

Output:

         motif      logFC          t     P.Value adj.P.Val tobias_delta
MA0466.2 CEBPB  15.908684  6.9933093 0.001752843 0.1006440      0.23037
MA0139.1  CTCF -14.510009 -2.1788685 0.090876332 0.3900205     -0.02669
MA0478.1 FOSL2  -1.110631 -0.3378268 0.751532341 0.8961142     -0.02214
MA0600.2  RFX2   8.519315  1.8783367 0.129349850 0.4192195      0.09382
MA0080.5  SPI1  -2.387368 -0.5002203 0.641728345 0.8458417      0.02664
         signac_pct   signac_p part5_padj_opened part5_padj_closed
MA0466.2  26.547059 0.01285707      1.0000000000      1.000000e+00
MA0139.1  -7.779012 0.25048344      1.0000000000      1.000000e+00
MA0478.1  -2.076521 0.40483559      1.0000000000      8.803895e-06
MA0600.2   1.540406 0.87628776      0.0006599133      1.000000e+00
MA0080.5   3.506595 0.15060916      1.0000000000      1.000000e+00

	Spearman's rank correlation rho

data:  tf_res$logFC and tf_res$tobias_delta
S = 23902835, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho
0.4345547

	Spearman's rank correlation rho

data:  tf_res$logFC and tf_res$motif_gc
S = 41405236, p-value = 0.6064
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho
0.02051799

R also prints Cannot compute exact p-value with ties after each cor.test(). That warning is harmless: some motifs share identical values, so R uses an approximation for the p-value.

First, check the recomputed Part 6 columns. signac_pct and signac_p match the published Part 6 table (CEBPB +26.5 percent, p = 0.013; CTCF -7.8 percent, p = 0.25). That confirms the per-donor flank file reproduces Part 6’s donor test exactly.

Now read the five motifs one at a time.

• CEBPB: three methods agree in direction, and none clears multiple testing.
  • Part 5’s motif enrichment gave it nothing, because it only looked at the 96 differential peaks.
  • Signac footprinting gave +26.5 percent (p = 0.013).
  • TOBIAS gave the largest change of the five (+0.230).
  • chromVAR now gives +15.9 (p = 0.0018, rank 7 of 633), with the whole C/EBP cluster moving the same way.
  The two replicate-aware tests both have nominal p < 0.05 and neither survives correction (BH 0.064 across five motifs in Part 6; 0.10 across 633 here). That is the strongest candidate this four-donor dataset produces, and it is still a candidate.
• CTCF, the control, moves a lot and is driven by one donor. Its logFC of -14.5 is the second-largest effect of the five, yet p = 0.09. The per-donor plot in Step 16 shows why: severe1 sits far below everyone (-13.8), while severe2 (+4.8) is in the healthy range. Part 6 saw the same thing in Signac, where severe1 was also the low CTCF donor. Two independent measurements agree that severe1 has unusually low CTCF-motif accessibility. The test correctly refuses to call that a disease effect.
• FOSL2 shows nothing (logFC -1.1, p = 0.75). This is the third method in a row that does not support Part 5’s “AP-1 closes in severe disease” at the donor level. Step 12 showed why the label was too coarse: FOSL2’s 7-bp-site cluster barely moves, while the 8-bp CREB/ATF-type cluster goes up.
• RFX2 leans up but rests on one donor (logFC +8.5, p = 0.13). The healthy donors agree closely (about -7). severe1 is high (+8.0) and severe2 is low (-4.8). The whole RFX family leans the same way (RFX1 to RFX5, logFC +7.7 to +8.8), all with p > 0.1. Part 6’s Signac flank result for RFX2 also had severe1 highest. That is consistent, but it is still a single-donor effect, and Part 5’s X-box enrichment remains unconfirmed.
• SPI1 behaves as a lineage control should (logFC -2.4, p = 0.64). Its two healthy donors differ by more than the conditions do, which is exactly why there is no call.

chromVAR and TOBIAS agree moderately across all 633 motifs (Spearman rho = 0.43, p < 2.2e-16). The two methods share no input file, score or statistic, so that much agreement means they are partly tracking the same biology.

The GC correlation is essentially zero (rho = 0.02, p = 0.61). After chromVAR’s GC-matched background, a motif’s GC content does not predict whether it goes up or down. So the chromVAR ranking is not a GC ranking. Part 6 asked whether the promoter-motif signal on the healthy side of TOBIAS was a shift of signal away from GC-rich promoters. The replicate-aware test here gives that signal little support, which fits that explanation better than a coordinated change in many promoter factors.

📊 Visualization

The same R session continues. Five figures cover the result from the whole-genome overview down to individual donors.

Volcano Plot of Differential TF Activity

#-----------------------------------------------
# STEP 14: Volcano plot of differential TF activity
#-----------------------------------------------

fp_motifs <- fp_ttest$feature
top_up    <- head(tf_res$motif[tf_res$adj.P.Val < 0.05 & tf_res$logFC > 0], 8)
top_dn    <- head(tf_res$motif[tf_res$adj.P.Val < 0.05 & tf_res$logFC < 0], 8)

tf_res$status <- ifelse(tf_res$adj.P.Val >= 0.05, "Not significant",
                        ifelse(tf_res$logFC > 0, "Higher in Severe", "Higher in Healthy"))
tf_res$label  <- ifelse(tf_res$motif %in% c(top_up, top_dn, fp_motifs), tf_res$motif, NA)

p_volcano <- ggplot(tf_res, aes(x = logFC, y = -log10(P.Value), color = status)) +
  geom_point(size = 1.2, alpha = 0.7) +
  geom_text_repel(aes(label = label), size = 3, color = "black", na.rm = TRUE,
                  min.segment.length = 0, max.overlaps = Inf) +
  scale_color_manual(values = c("Higher in Severe" = "#C0392B",
                                "Higher in Healthy" = "#2E86AB",
                                "Not significant"  = "grey75"), name = NULL) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  labs(x = "Difference in chromVAR z-score (Severe - Healthy)",
       y = "-log10 limma p-value",
       title = "Differential TF activity -- CD14 monocytes, 633 JASPAR2020 motifs") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "01_tf_activity_volcano.png"), p_volcano,
       width = 8, height = 6, dpi = 300)

Reading this figure: each point is one of the 633 motifs. Right of the dashed line means more motif-associated accessibility in severe disease; left means more in healthy.

Every point is grey, and only the five Part 6 motifs are labelled. No motif passed adj.P.Val < 0.05, so nothing was colored. The code picks extra labels only among significant motifs, so top_up and top_dn are empty. That is the honest picture of this dataset. If you want the strongest motifs named anyway, drop the adj.P.Val < 0.05 condition from the top_up and top_dn lines, and say in the legend that the labelled motifs are not significant.

A large effect is not the same as a small p-value. CTCF sits far out on the left (logFC -14.5) but at -log10 p of only about 1, because one donor carries the whole effect. RFX2 is on the right at +8.5 for the same reason. With two donors per group, how far a point sits from zero says little about how well the donors agree; the height does.

Donor-Level Heatmap of the Top Motifs

#-----------------------------------------------
# STEP 15: Heatmap of the top motifs, one column per donor
#-----------------------------------------------

# The 15 most significant motifs in each direction (tf_res is sorted by p-value)
z_qn    <- normalizeQuantiles(z)
top_ids <- c(head(tf_res$motif_id[tf_res$logFC > 0], 15),
             head(tf_res$motif_id[tf_res$logFC < 0], 15))

# Centre each motif across donors so the color shows severe-vs-healthy contrast
z_top <- z_qn[top_ids, ] - rowMeans(z_qn[top_ids, ])
heat_df <- data.frame(motif = rep(motif_names[top_ids], times = ncol(z_top)),
                      donor = rep(colnames(z_top), each = nrow(z_top)),
                      z     = as.vector(z_top))
heat_df$motif <- factor(heat_df$motif, levels = rev(unique(motif_names[top_ids])))
heat_df$donor <- factor(heat_df$donor, levels = c("healthy1", "healthy2", "severe1", "severe2"))

p_heatmap <- ggplot(heat_df, aes(x = donor, y = motif, fill = z)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "#2E86AB", mid = "white", high = "#C0392B",
                       name = "Centred\nz-score") +
  labs(x = NULL, y = NULL,
       title = "Top 15 motifs in each direction") +
  theme_minimal(base_size = 10) +
  theme(panel.grid = element_blank())

ggsave(file.path(plot_dir, "02_tf_activity_heatmap.png"), p_heatmap,
       width = 5, height = 8, dpi = 300)

Reading this figure: rows are motifs (the 15 most significant with positive logFC on top, the 15 with negative logFC below), columns are donors, and each row is centred so its four values average zero. These are the same quantile-normalized values limma tested.

The upper half is a clean block, blue in both healthy columns and red in both severe columns. Thirteen of its 15 rows are bZIP motifs; IRF3 and ERF are the exceptions. Look closer and the severe columns are not identical: two programs with different donor emphasis are hiding in this block.

• The C/EBP rows are deeper red in severe1.
• The CREB/ATF-type rows are as deep or deeper in severe2.

Both severe donors are above both healthy donors in each program, which is why both reach the top of the list. But with two donors, “C/EBP is strongest in one patient and CREB/ATF-type in the other” cannot be told apart from noise. Keep it as an observation, not a finding.

The lower half is less tidy. Seven of its 15 rows are octamer (POU) motifs, plus EBF1 and EBF3. Here severe2 is the deepest blue and severe1 paler.

• Some rows are carried by one donor.
• This is the pattern to worry about: a row where a single donor carries the color. limma’s variance moderation reduces the weight of such rows but cannot remove it, and a motif driven by one person is not a disease effect.

Remember the donor sizes. severe2 is the smallest pseudobulk (508 cells, 1.8 million fragments) and healthy2 the largest (2,575 cells, 9.6 million). In a pseudobulk analysis each counts once, but the smaller donor’s z-scores are noisier, so an extreme severe2 value alone deserves less weight than one shared by both severe donors.

Per-Donor TF Activity for Selected Motifs

#-----------------------------------------------
# STEP 16: Donor-level z-scores for selected motifs
#-----------------------------------------------

# Four strongest motifs each way, plus the five motifs footprinted in Part 6
show_names <- unique(c(head(top_up, 4), head(top_dn, 4), fp_motifs))
show_ids   <- tf_res$motif_id[match(show_names, tf_res$motif)]

dot_df <- data.frame(motif = rep(show_names, times = ncol(z_qn)),
                     donor = rep(colnames(z_qn), each = length(show_ids)),
                     z     = as.vector(z_qn[show_ids, ]))
dot_df$condition <- donor_info[dot_df$donor, "condition"]
dot_df$motif     <- factor(dot_df$motif, levels = show_names)

p_donor <- ggplot(dot_df, aes(x = condition, y = z, color = condition)) +
  geom_point(size = 3) +
  geom_text_repel(aes(label = donor), size = 2.5, show.legend = FALSE) +
  facet_wrap(~ motif, scales = "free_y", ncol = 4) +
  scale_color_manual(values = c(Healthy = "#2E86AB", Severe = "#C0392B")) +
  labs(x = NULL, y = "Quantile-normalized chromVAR z-score",
       title = "Per-donor TF activity -- CD14 monocytes") +
  theme_classic(base_size = 10) +
  theme(legend.position = "none")

ggsave(file.path(plot_dir, "03_tf_activity_by_donor.png"), p_donor,
       width = 10, height = 2.4 * ceiling(length(show_names) / 4), dpi = 300)

Reading this figure: one panel per motif, four points per panel — the data limma actually saw.

Only five panels, for the same reason as the volcano. top_up and top_dn are empty because no motif passed 0.05, so the figure shows just the five Part 6 motifs. That makes it a direct Part 7 counterpart of Part 6’s per-donor flank plot, panel for panel.

Read the axis numbers before the pattern. Each panel has its own y-scale (scales = "free_y"). Here the ranges are similar (roughly -15 to +12), but a free scale can make a tiny difference look as clean as a large one.

• CEBPB is the only panel where the groups separate.

Compare with Part 6’s Signac donor plot. There too, CEBPB was the only motif whose two severe donors both sat above both healthy donors. That agreement is reassuring about the methods, and it is a reminder that individual donors, not only conditions, shape every result in this series.

chromVAR Versus TOBIAS Across All 633 Motifs

Part 6 could only compare its two methods on five motifs, and warned that five points is not a correlation. With chromVAR scores for every motif, the comparison can now be made genome-wide.

#-----------------------------------------------
# STEP 17: chromVAR versus TOBIAS across all motifs
#-----------------------------------------------

p_concord <- ggplot(tf_res, aes(x = tobias_delta, y = logFC)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(aes(color = status), size = 1.2, alpha = 0.7) +
  geom_text_repel(aes(label = label), size = 3, color = "black", na.rm = TRUE,
                  min.segment.length = 0, max.overlaps = Inf) +
  scale_color_manual(values = c("Higher in Severe" = "#C0392B",
                                "Higher in Healthy" = "#2E86AB",
                                "Not significant"  = "grey75"),
                     name = "chromVAR + limma") +
  labs(x = "TOBIAS differential binding score, Part 6 (Severe - Healthy)",
       y = "chromVAR z-score difference, Part 7 (Severe - Healthy)",
       title = "Accessibility-based versus footprint-based TF activity") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "04_chromvar_vs_tobias.png"), p_concord,
       width = 8, height = 6.5, dpi = 300)

Reading this figure: each point is a motif. The x-axis is TOBIAS’s change from Part 6 (sign flipped so that positive means higher in severe); the y-axis is chromVAR’s logFC. Points in the upper-right and lower-left quadrants are motifs on which the two methods agree in direction — 61 percent of all 633. All points are grey because no motif passed correction in chromVAR.

The upper-right arm is where the methods agree most.

The lower-left quadrant has no matching arm. The motifs most negative in TOBIAS are only mildly negative in chromVAR.

The tall column at TOBIAS around zero is where chromVAR sees something TOBIAS does not.

Disagreement on one axis has two documented explanations:

• Many factors leave no measurable footprint. Part 6 covered this: their motif regions can open or close without any change in footprint shape.
• TOBIAS’s score largely cancels a change in how open the whole region is. It compares depletion at the motif with its immediate flanks, so a uniform shift drops out. chromVAR counts every fragment in the peak and sees the full shift. Part 6 showed exactly this for CTCF: the two condition tracks were offset uniformly, and TOBIAS’s change stayed small. Step 16 adds that the offset comes from one donor, severe1.

The axes are in different units. Compare signs and ranks, as the Spearman correlation in Step 13 does (rho = 0.43). Do not fit a regression line between them.

Seed Reproducibility

#-----------------------------------------------
# STEP 18: Seed reproducibility plot
#-----------------------------------------------

p_seed <- ggplot(seed_cmp, aes(x = t_seed1, y = t_seed2)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(size = 1, alpha = 0.6, color = "#34495E") +
  labs(x = "Moderated t, background seed 1234",
       y = "Moderated t, background seed 2026",
       title = "Same data, two background draws") +
  theme_classic(base_size = 12)

ggsave(file.path(plot_dir, "05_seed_reproducibility.png"), p_seed,
       width = 6, height = 6, dpi = 300)

Reading this figure: every point is a motif, plotted by its moderated t-statistic under two independent background draws. With 2000 background sets the points lie on the dashed diagonal. Points scattered off the line at the extremes — where calls are made — would be the reproducibility problem Gerbaldo et al. (2024) described. Here there are none. This figure belongs in a supplement: it shows reviewers that your top motifs are a property of the data, not of a random number generator.

💾 Saving Results for Reuse and Supplementary Material

#-----------------------------------------------
# STEP 19: Save tables and session information
#-----------------------------------------------

write.csv(tf_res[, c("motif_id", "motif", "cluster", "logFC", "AveExpr", "t", "P.Value",
                     "adj.P.Val", "B", "tobias_delta", "motif_gc")],
          file.path(table_dir, "chromvar_limma_cd14mono.csv"), row.names = FALSE)
write.csv(family_tab, file.path(table_dir, "chromvar_limma_families.csv"), row.names = FALSE)
write.csv(five,       file.path(table_dir, "five_motif_comparison_parts5to7.csv"), row.names = FALSE)
write.csv(data.frame(motif_id = rownames(z), motif = motif_names[rownames(z)], z,
                     check.names = FALSE),
          file.path(table_dir, "chromvar_zscores_by_donor.csv"), row.names = FALSE)
write.csv(seed_cmp,   file.path(table_dir, "seed_reproducibility.csv"), row.names = FALSE)

writeLines(capture.output(sessionInfo()), file.path(tfa_dir, "sessionInfo_part7.txt"))

Together with the two objects saved in Steps 6 and 9, these files let anyone reproduce every number and figure without recounting fragments:

✅ Best Practices for Differential TF Activity Analysis

1. Make the donor the unit of replication. Pseudobulk each cell type per donor before running chromVAR, then test across donors. Per-cell chromVAR scores are excellent for UMAPs and cluster descriptions; testing them between conditions repeats the pseudoreplication error that Teo et al. (2024) documented for peaks.

2. Analyze one cell type at a time. A pseudobulk of all cells measures cell-type composition as much as regulation. If severe donors have proportionally more monocytes, every monocyte motif will look “more active” in severe disease without any monocyte changing.

3. Use fixed-width peaks. Resize to 300 bp around the summit if you have one, the midpoint if you do not, and recount. Variable-width peaks give long peaks more motif matches and more fragments at once.

4. Raise the background iterations to 1000-2000. The default of 50 is fine for clustering cells and not fine for differential testing (Gerbaldo et al. 2024). On pseudobulks the extra cost is modest.

5. Test z-scores with limma, not with differentialDeviations(). Variance moderation is what makes 2- and 3-replicate designs workable, and limma on z-scores clearly outperformed chromVAR’s built-in test in the benchmark.

6. Quantile-normalize before testing. It was the best-performing downstream normalization in the benchmark, and it removes donor-wide shifts that could otherwise masquerade as many coordinated motif changes.

7. Prove seed robustness. Rerun the background with a different seed and compare (Step 11). It is cheap and is the most direct defence against irreproducible hits.

8. Report families, not motifs. Collapse redundant matrices before counting results, using BINDetect clusters, motif clustering or the JASPAR family annotation, and check that family members agree in direction.

9. Treat adjusted p-values as a ranking with a threshold. Gerbaldo et al. (2024) found the observed false discovery rate was nearly always above nominal. Prefer families with several concordant members, strong effects and support from a second method over a single motif just under 0.05.

10. Expect few or no FDR-significant motifs with two donors per group, and report that honestly. In this run the best motif had a moderated t of 12.7 and still only reached adj.P = 0.098, because each motif had about 4 degrees of freedom even after moderation. When nothing passes, report the ranked families, the excess of small unadjusted p-values over chance (95 motifs below 0.05 against about 32 expected), and the agreement with other methods. Then say clearly that more donors are needed.

11. Cross-check with a second, different method. The benchmark recommends combining top methods. In this series that is footprinting from Part 6; within the benchmarked set, monaLisa was another of the top performers.

12. Check for GC effects explicitly. chromVAR matches background on GC, but sample-specific GC biases exist in ATAC-seq data. A correlation between motif GC content and logFC (Step 13) is a quick warning sign.

13. Keep motif universes and peak sets identical across analyses. This tutorial uses the same JASPAR2020 query as Part 4 and the same monocyte peaks as Part 6. Without that, disagreements between methods could simply be database or peak-set differences.

14. Never compare z-scores across separate chromVAR runs. Each run’s expectation is built from the samples in that run. To compare cell types, either run them together and model cell type in limma, or compare test results, not raw z-scores.

15. Say what the score is. chromVAR measures transcription-factor-associated accessibility. Write “C/EBP motif-associated accessibility increases” rather than “C/EBP binding increases,” unless footprinting or ChIP-type data back the stronger claim.

16. Check the matrix, not just the name. A TF name tells you which protein a matrix was measured for, not which DNA sequence it describes. In this run, motifs labelled JUN or FOS::JUN(var.2) were built on the 8-bp TGA(C/T)GTCA core shared with CREB1, while FOS and FOSL2 were built on the 7-bp TGACTCA core, and the two groups moved in opposite directions. Before writing “AP-1 goes up,” look at the matrices or at the motif clusters.

17. Look at the donors behind every top motif. In this dataset severe1 alone drove CTCF, NRF1, STAT3 and THAP11 down and RFX2 up, and the same donor-specific pattern appeared in Part 6. A per-donor plot (Step 16) takes seconds and catches what a p-value summarizes away.

18. Record versions. chromVAR, motifmatchr, JASPAR2020, limma, Signac and the random seeds all shape the result. Save sessionInfo() with every run (Step 19).

⚠️ Common Pitfalls and How to Avoid Them

🎯 Conclusion

You have taken the series from hand-picked motifs to a genome-wide, replicate-aware screen of every transcription factor motif in the database. Along the way you have:

1. Learned what TF activity analysis measures: the pooled accessibility of all peaks containing a motif, relative to GC- and accessibility-matched background peaks — transcription-factor-associated accessibility, not binding.
2. Placed it beside Parts 5 and 6: differential accessibility asks about loci, footprinting asks about binding at chosen motifs, and TF activity asks about regulators across the genome.
3. Built a donor pseudobulk on fixed-width peaks: 115,324 monocyte-accessible peaks resized to 300 bp, recounted from the fragments files, and filtered to 115,305.
4. Run chromVAR with benchmark-based settings: 2000 background sets, z-scores for all 633 JASPAR2020 motifs in each of four donors.
5. Tested severe versus healthy with limma: quantile-normalized z-scores, moderated t-statistics that borrow variance information across motifs.
6. Checked the result four ways: across random seeds, across motif families, against Parts 5 and 6, and against motif GC content.

📚 References

1. Schep AN, Wu B, Buenrostro JD, Greenleaf WJ. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nature Methods. 2017;14(10):975-978. doi:10.1038/nmeth.4401
2. Gerbaldo FE, Sonder E, Fischer V, Frei S, Wang J, Gapp K, Robinson MD, Germain PL. On the identification of differentially-active transcription factors from ATAC-seq data. PLoS Computational Biology. 2024;20(10):e1011971. doi:10.1371/journal.pcbi.1011971
3. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 2015;43(7):e47. doi:10.1093/nar/gkv007
4. Smyth GK. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology. 2004;3(1):Article 3. doi:10.2202/1544-6115.1027
5. 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
6. 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
7. 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
8. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15:550. doi:10.1186/s13059-014-0550-8
9. 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

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