Motif Enrichment Names a Family, Not a Factor

Zetobit · Bioinformatics Insight Series

Regulatory Genomics

Motif Enrichment Names a Family, Not a Factor

The last figure in most accessibility analyses is a list of transcription factors. It is a list of sequence patterns with protein names attached by a database — and nothing upstream checks whether those proteins are in the cells.

Kanna Nandakumar, PhD · Zetobit, LLC

An ATAC-seq analysis usually ends with transcription factors. A ranked table, or a heatmap of deviation scores, populated with names: GATA1, CEBPA, SPI1, an ETS factor or two. The output reads as a statement about proteins doing things in cells. What was computed is a statement about short degenerate sequence patterns and their distribution across a peak set. The protein names arrived from a database at the annotation step, and they are the part of the result that travels.

The arithmetic of redundancy

There are on the order of sixteen hundred human transcription factors, of which about 1,107 have a known motif.² But motifs are short and tolerant — typically eight to fourteen bases with several positions accepting multiple nucleotides⁹ — and they are nowhere near as distinguishable as the proteins are distinct. Computing pairwise similarity across 2,179 motif models from human and mouse and clustering them yields 286 distinct motif clusters, each summarised by an archetypal consensus.¹ Roughly eight models per separable pattern.

That number reframes the familiar output. When five related factors occupy the top ten rows of an enrichment table, that is not five converging lines of evidence. It is one pattern reported five times, and the ordering among them is determined by small differences between position weight matrices rather than by anything measured in the experiment. A reader who takes the top hit as "the" factor is reading a ranking that the data did not produce.

Paralogs are not a database artefact

It would be convenient if the redundancy were purely a cataloguing problem. It is not. Paralogous transcription factors are frequently reported to have identical binding motifs while performing distinct, sometimes opposing, regulatory functions, and they are often co-expressed in the same cells.⁴ The androgen and glucocorticoid receptors are the textbook case: nearly identical DNA-binding specificities, entirely different transcriptional programmes in vivo.⁵

Careful in vitro work shows that closely related paralogs do bind genomic sites differently, contrary to the assumption that their specificities are indistinguishable. But the way they differ is exactly wrong for our purposes: across eleven factors from four structural families, the divergence was concentrated at medium- and low-affinity sites, which current motif models capture poorly.⁴

Paralogs diverge most at the sites the motif model represents worst. The discrimination you need lives in the part of the signal a PWM discards.

What a deviation score actually is

chromVAR is a well-built method and worth describing precisely rather than dismissively. It aggregates chromatin accessibility across peaks sharing a common feature — a motif match — and applies bias correction using background peak sets matched for GC content and average accessibility, producing a per-cell deviation from the expectation under equal accessibility profiles.³ The bias correction is the substantive contribution; without it the scores track GC content and depth.

Note the shape of the resulting quantity. It is a property of a set of peaks defined by a sequence match, expressed relative to matched background. It is not an estimate of protein abundance, nuclear occupancy, or regulatory activity. The row label is the motif's name. Two documented caveats deserve to be surfaced more often than they are. ArchR's own documentation notes that the absolute value of the deviation z-score is correlated with per-cell read depth⁸ — so in sparse data the magnitude carries a technical component. And chromVAR reports fractionMatches, the fraction of peaks containing a match for each motif; a pattern present in 60% of your peaks and one present in 2% produce scores in the same column, with very different meanings, and that diagnostic is almost never shown.

The demonstration is in the method's own paper

The strongest argument here is not a criticism of the tool. It is something the tool's authors did deliberately: chromVAR computes deviation scores for seven-mers, displaying the closest matching known motif underneath for interpretation.³ The method works on raw k-mers. It produces interpretable, cell-type-discriminating scores for sequence features that have no protein attached to them at all.

Which settles what the factor name is. If the identical pipeline yields a useful result for a seven-base string with no name, then in a normal run the transcription factor's name is a label applied to a sequence feature after the computation, for the convenience of interpretation. The designers understood this perfectly well; the k-mer mode exists precisely because the annotation is separable from the measurement. Downstream, in a figure legend and then in a slide, the convenience quietly becomes the claim.

Motif presence is not binding, by two orders of magnitude

A typical PWM scan with default parameters yields over a million potential binding sites across three billion bases, often with several high-scoring matches per gene, while the number actually occupied is far lower.⁹ One worked example: the SP1 motif occurs more than 3.6 million times in the human genome, and fewer than 1% of those occurrences are bound in a lymphoblastoid cell line.⁶

Restricting the scan to accessible peaks is exactly the right response, and it is why doing this in ATAC data is more informative than scanning the genome. But it narrows the gap rather than closing it: occupancy also depends on nucleosome positioning, cooperative and competitive binding with other factors, local sequence context and cell state.⁶ And the reference assay has its own version of the problem. Across hundreds of TF ChIP-seq datasets, a small set of profiles — CTCF-like, ETS-like, JUN-like and THAP11 — turn up significantly enriched in experiments that did not target them, and peaks carrying these patterns while lacking the immunoprecipitated factor's own motif can make up as much as 45% of a dataset.⁷ If the direct measurement of occupancy carries that much off-target structure, an inference from sequence alone deserves correspondingly wide error bars.

The check nobody runs, and why it is better than it looks

Nothing in the standard workflow verifies that the named factor is expressed in the cells being analysed. This is the cheapest available filter, it is almost never applied, and in a multiome or paired experiment the expression matrix is sitting in the same object.

It is also better than a sanity check, for a reason that is easy to miss. Analysis of human paralogous TF pairs found that when the binding motifs of two paralogs are similar, the expression of the genes encoding them has typically diverged — the two modes of functional diversification substitute for one another.¹⁰ That is an exploitable asymmetry rather than a coincidence. The redundancy that defeats sequence-based discrimination is anticorrelated with redundancy in expression. Precisely where the motif cannot separate family members, the expression pattern most often can. Filtering an archetype's membership to the factors actually expressed will not identify the responsible protein, but it frequently removes most of the candidates that cannot be responsible — and it does so using the axis along which those candidates differ.

The claim as written	What was computed	Standing
“Motif X is enriched in these peaks”	A set-level comparison against a background, matched for GC and accessibility	Supported, if the background is stated
“TF X is active in these cells”	The same set-level quantity, with a database name attached	Family-level at best; expression unchecked
“TF X is the top regulator”	A ranking among near-identical models within one archetype	The ordering is not a finding
“TF X activity increased 2-fold”	A change in aggregate accessibility of a motif-defined peak set	Magnitude partly tracks read depth
“These peaks are X binding sites”	Sequence matches within accessible regions	Presence is not occupancy, by orders of magnitude

Only the first row describes what the pipeline measured. Every row below it adds a noun the arithmetic never supplied.

Schematic. Members of a motif archetype are not separable by sequence, so they appear as consecutive rows in an enrichment table with near-identical statistics — four p-values reporting one observation. The ordering among them is set by small differences between position weight matrices. Expression, which tends to diverge in paralogs precisely when the motifs have not, is the axis on which the family can often be narrowed.

What to do about it

1. Use clustered motif archetypes rather than raw databases. ArchR supports the Vierstra archetype sets natively; the redundancy is a solved problem with tool support.
2. Report the archetype first and its members second. “An ETS-family archetype, whose members include…” is honest, still useful, and survives review.
3. Filter to expressed factors. In multiome or paired data this is one join, and it is the single highest-value addition to a motif result.
4. Show fractionMatches next to every score. A motif in most of your peaks and a motif in a few percent are not comparable quantities.
5. Do not rank paralogs against each other on motif score. That ordering is an artefact of matrix construction.
6. Treat deviation magnitude cautiously in sparse data — check whether it correlates with per-cell depth in your own dataset before reading it as activity.
7. State the background explicitly. GC- and accessibility-matched background is what makes the number interpretable; an unstated background makes it uninterpretable.
8. Where a factor name drives a decision — a target, a mechanism claim, a client deliverable — plan orthogonal evidence: knockdown, ChIP, or at minimum expression plus a literature argument.
9. Write the claim at the level the data supports. Most of the scientific value of a motif result survives the downgrade from factor to family; almost none of the risk does.

The name is the last step and the only one nobody checks

Three failures in this series' recent chromatin entries share a structure worth naming. A zero acquired the meaning “closed” from a data structure inherited from transcriptomics. A regulatory link acquired a mechanism from being reported as a genomic coordinate. Here, a set-level accessibility statistic acquires a protein from a row label. In each case the arithmetic is sound and the noun is doing work the arithmetic never did — and in each case nobody makes the decision, because it arrives with the output format.

This is not an argument against motif enrichment, which is among the most useful things you can do with an accessibility experiment and often the only route from peaks to mechanism. The methods are careful, the bias correction is real, the archetype resources exist and are supported in the standard tools, and the expression filter costs a single line. The gap between what the analysis computes and what the figure says is narrower than in most problems this series covers — which is exactly why leaving it open is hard to defend.

The motif was enriched. The factor was named. A database did the naming, and the database has never seen your cells.

References

1. Vierstra J, Lazar J, Sandstrom R, et al. Global reference mapping of human transcription factor footprints. Nature. 2020;583:729–736. doi:10.1038/s41586-020-2528-x (Motif archetype clustering resource: vierstra.org/resources/motif_clustering.)
2. Lambert SA, Jolma A, Campitelli LF, et al. The human transcription factors. Cell. 2018;172(4):650–665. doi:10.1016/j.cell.2018.01.029
3. 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
4. Shen N, Zhao J, Schipper JL, et al. Divergence in DNA specificity among paralogous transcription factors contributes to their differential in vivo binding. Cell Systems. 2018;6(4):470–483.e8. doi:10.1016/j.cels.2018.02.009
5. Androgen and glucocorticoid receptor direct distinct transcriptional programs by receptor-specific and shared DNA binding sites. Nucleic Acids Research. 2021. PMC8053126 (Author list and volume not confirmed at time of writing.)
6. Zheng A, Lamkin M, Zhao H, et al. Deep neural networks identify context-specific determinants of transcription factor binding affinity. bioRxiv. 2020:2020.02.26.965343. doi:10.1101/2020.02.26.965343 (Check for the peer-reviewed version.)
7. Worsley Hunt R, Wasserman WW. Non-targeted transcription factors motifs are a systemic component of ChIP-seq datasets. Genome Biology. 2014;15:412. doi:10.1186/s13059-014-0412-4
8. Granja JM, Corces MR, Pierce SE, et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nature Genetics. 2021;53(3):403–411. doi:10.1038/s41588-021-00790-6
9. Jolma A, Yan J, Whitington T, et al. DNA-binding specificities of human transcription factors. Cell. 2013;152(1–2):327–339. doi:10.1016/j.cell.2012.12.009
10. Singh LN, Hannenhalli S. Functional diversification of paralogous transcription factors via divergence in DNA binding site motif and in expression. PLoS ONE. 2008;3(6):e2345. doi:10.1371/journal.pone.0002345

Zetobit, LLC · Bioinformatics consulting · Lexington, KY · zetobit.com
Related in this series: Sparsity in Single-Cell ATAC-seq, Peak-to-Gene Assignment, Cell Composition in Bulk ATAC-seq.