Paper deep dive
Causal dictionary learning reveals and validates transcription-factor binding features in genomic language models
Sarwan Ali
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/23/2026, 2:26:28 AM
Summary
The paper introduces a framework combining sparse dictionary learning with causal intervention to extract and validate interpretable features in genomic language models (Nucleotide Transformer and DNABERT-2). It demonstrates that naive validation of features against transcription factor (TF) motifs is confounded by sequence composition (GC content, repeats). By using a composition-matched binding test and causal ablation (measuring KL divergence shifts upon feature removal), the authors identify features that causally represent cell-type-specific TF binding (CTCF, GATA1, REST) rather than just motif presence.
Entities (10)
Relation Signals (7)
Sparse Autoencoder → extractsfeaturesfrom → Nucleotide Transformer
confidence 95% · Training top-k sparse autoencoders on the hidden activations of two architecturally distinct models, Nucleotide Transformer...
Sparse Autoencoder → extractsfeaturesfrom → DNABERT-2
confidence 95% · Training top-k sparse autoencoders on the hidden activations of two architecturally distinct models... DNABERT-2...
Causal Ablation → validates → CTCF
confidence 92% · Across three transcription factors (CTCF, GATA1, REST)... causally validated binding features emerge reproducibly
Causal Ablation → validates → GATA1
confidence 92% · Across three transcription factors (CTCF, GATA1, REST)... causally validated binding features emerge reproducibly
Causal Ablation → validates → REST
confidence 92% · Across three transcription factors (CTCF, GATA1, REST)... causally validated binding features emerge reproducibly
GC Composition → confounds → Naive Validation
confidence 90% · naive validation of such features against position weight matrices is severely confounded by GC composition
ChIP-seq → providesgroundtruthfor → TranscriptionFactor
confidence 88% · experimental ChIP-seq binding, is available... ENCODE ChIP-seq IDR peak sets for CTCF... GATA1... REST
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Genomic language models achieve strong performance across regulatory-genomics tasks, yet what these models internally represent remains opaque, and the field lacks a principled procedure for verifying that an apparent ``concept'' inside a model is real rather than an artifact of sequence composition. We introduce a framework that combines sparse dictionary learning with causal intervention to extract, validate, and causally test interpretable features in genomic foundation models. Training top-$k$ sparse autoencoders on the hidden activations of two architecturally distinct models, Nucleotide Transformer ($6$-mer tokenization) and DNABERT-2 (byte-pair encoding), we recover thousands of monosemantic features that map to transcription-factor (TF) sequence motifs. We show that the naive validation of such features against position weight matrices is severely confounded by GC composition and repetitive elements, producing hundreds of spurious ``TF features'', and we develop a composition-matched, binding-resolved protocol that removes these confounds. Critically, we move beyond correlation: by ablating individual dictionary directions during the model's forward pass and measuring the induced shift in the model's own predictive distribution, we establish that specific features are \emph{causally} used to represent cell-type-specific TF binding, not merely motif presence. Across three transcription factors (CTCF, GATA1, REST) and both architectures, causally validated binding features emerge reproducibly ($7$--$14$ of $15$ tested features per condition), while two classes of negative control, scrambled binding labels and randomly selected features, yield no detectable signal. The framework is purely computational, uses only public data, and provides a reusable standard for interpretability claims in genomic deep learning.
Tags
Links
- Source: https://arxiv.org/abs/2607.19618v1
- Canonical: https://arxiv.org/abs/2607.19618v1
Trouble viewing inline? Open PDF directly →
Full Text
41,842 characters extracted from source content.
Expand or collapse full text
Causal dictionary learning reveals and validates transcription-factor binding features in genomic language models Sarwan Ali Columbia University Irving Medical Center, New York, USA sa4559@cumc.columbia.edu Abstract Genomic language models achieve strong performance across regulatory-genomics tasks, yet what these models internally represent remains opaque, and the field lacks a principled procedure for verifying that an apparent “concept” inside a model is real rather than an artifact of sequence composition. We introduce a framework that combines sparse dictionary learning with causal intervention to extract, validate, and causally test interpretable features in genomic foundation models. Training top-k sparse autoencoders on the hidden activations of two architecturally distinct models, Nucleotide Transformer (66-mer tokenization) and DNABERT-2 (byte-pair encoding), we recover thousands of monosemantic features that map to transcription-factor (TF) sequence motifs. We show that the naive validation of such features against position weight matrices is severely confounded by GC composition and repetitive elements, producing hundreds of spurious “TF features”, and we develop a composition-matched, binding-resolved protocol that removes these confounds. Critically, we move beyond correlation: by ablating individual dictionary directions during the model’s forward pass and measuring the induced shift in the model’s own predictive distribution, we establish that specific features are causally used to represent cell-type-specific TF binding, not merely motif presence. Across three transcription factors (CTCF, GATA1, REST) and both architectures, causally validated binding features emerge reproducibly (77–1414 of 1515 tested features per condition), while two classes of negative control, scrambled binding labels and randomly selected features, yield no detectable signal. The framework is purely computational, uses only public data, and provides a reusable standard for interpretability claims in genomic deep learning. 1 Introduction DNA language models trained by self-supervision on reference genomes have rapidly become general-purpose tools in regulatory genomics, supporting variant-effect prediction, chromatin-state annotation and regulatory-element discovery [6, 31, 23, 19, 18, 2]. Their predictive utility is now well documented, but their interpretability is not. Unlike supervised models built around designed features [30, 12, 1], foundation models distribute information across high-dimensional activations in which individual neurons are polysemantic, responding to many unrelated sequence properties at once [8]. Consequently, it is rarely possible to state what a genomic language model has actually learned, which limits both scientific trust and the use of these models for hypothesis generation about genome regulation [20]. In natural-language processing, sparse autoencoders (SAEs) have emerged as a powerful remedy for polysemanticity: by reconstructing model activations through an over-complete, sparsely active dictionary, they decompose distributed representations into approximately monosemantic features [5, 3, 27, 9]. Whether this paradigm transfers to genomic models, and whether the resulting features correspond to bona fide biological entities, has not been established. Two obstacles stand in the way. First, genomic models differ fundamentally in how they tokenize DNA, fixed k-mers, learned byte-pair vocabularies, or single nucleotides, and it is unknown whether dictionary learning behaves consistently across these schemes. Second, and more seriously, the genome’s sequence statistics make naive validation treacherous: many transcription-factor motifs are GC-rich, as are CpG islands and the most abundant repetitive elements in the human genome [7, 24], so a feature that merely detects local GC content or an Alu repeat will spuriously appear “enriched” for any GC-rich motif. An interpretability method that does not control for this will systematically over-report success. We address both obstacles and, in doing so, propose a general standard. Our contributions are as follows. (i) We train top-k SAEs on the activations of two genomic language models with different tokenizations and show that monosemantic, motif-aligned features arise robustly in both, with middle layers yielding the most interpretable dictionaries. (i) We demonstrate quantitatively that position-weight-matrix enrichment, the obvious validation, is dominated by compositional confounds, and we introduce a GC-matched, binding-resolved test that isolates genuine signal. (i) Most importantly, we develop a causal validation: ablating a single dictionary direction during the forward pass and measuring the Kullback–Leibler shift in the model’s masked-token predictions at bound versus unbound motif sites. This distinguishes features the model uses to represent binding from features that merely correlate with it. (iv) We validate the whole framework with stringent negative controls and show cross-architecture, cross-TF reproducibility. The result is a purely computational, fully reproducible procedure, built entirely on public models and data, for making, and disciplining, interpretability claims about genomic deep learning. 2 Related Work Genomic language models. Self-supervised models of DNA have progressed rapidly from early k-mer transformers [11] to large multi-species encoders such as the Nucleotide Transformer [6] and the byte-pair-encoded DNABERT-2 [31], and to long-context architectures including HyenaDNA [19], the state-space model Caduceus [23] and the genome-scale Evo model [18, 10]. These models support competitive variant-effect prediction [2] and regulatory annotation, but recent evaluations question how much regulatory signal their representations actually contain [26], underscoring the need for tools that probe what they encode rather than only how well they predict. Supervised sequence models such as DeepSEA [30], Basenji [12] and Enformer [1] remain the performance reference for binding and expression prediction; our aim is not to compete with them but to interpret the self-supervised models. Interpretability of deep models in genomics. Attribution methods such as DeepLIFT [25] and SHAP [13], together with in-silico mutagenesis, are the standard tools for explaining genomic neural networks [20]. These produce per-input importance maps but do not yield a global vocabulary of reusable features, and they inherit the polysemanticity of individual neurons [8]. Our work instead learns an explicit, model-wide dictionary and asks which of its elements the model causally uses. Sparse autoencoders and mechanistic interpretability. Sparse dictionary learning on model activations has become a leading approach to extracting monosemantic features from language models [5, 3, 27], with the top-k formulation offering a clean sparsity–fidelity trade-off [9]. Validating such features is an open problem: activation patching and causal mediation establish that components are used by a model [28, 16], but apparent feature interpretations can be illusory if not causally grounded [14]. We bring this causal-intervention discipline to genomics, where, uniquely, an external ground truth, experimental ChIP-seq binding, is available, and where sequence-composition confounds make naive validation especially hazardous. Transcription-factor binding resources. We ground features against curated motif models from JASPAR [22] and HOCOMOCO [29] and against experimental occupancy from ENCODE ChIP-seq [4, 17]. The biology of our test factors, CTCF as a GC-rich insulator [21] whose motifs are frequently embedded in Alu repeats [7, 24], is precisely what makes composition control essential. 3 Methodology We formalize the three components of the framework, the sparse dictionary, the composition-matched binding test, and the causal ablation, and then state two propositions that explain why the design isolates binding-specific causal use. The full procedure is summarized in Algorithm 1. 3.1 Top-k sparse autoencoders Let ∈ℝdx ^d denote a (mean-centered, scaled) model activation. A top-k sparse autoencoder [9] encodes x into a sparse code ∈ℝmf ^m with m=16dm=16d and decodes it back: =enc(−pre), =W_enc(x-b_pre), (1) =TopKk(ReLU()), =TopK_k\! (ReLU(z) ), (2) x =dec+pre, =W_decf+b_pre, (3) where TopKkTopK_k retains the k largest entries and zeroes the rest, the columns jd_j of dec∈ℝd×mW_dec ^d× m are constrained to unit norm, and pre∈ℝdb_pre ^d is a learned pre-bias. The training objective is reconstruction error with an auxiliary term that revives inactive features, ℒ=‖−^‖22⏟reconstruction+α‖(−^)−^aux‖22⏟dead-feature revival,L= \|x- x \|_2^2_reconstruction+α\, \|(x- x)- x_aux \|_2^2_dead-feature revival, (4) where ^aux x_aux reconstructs the residual using only the top dead features and α=1/16α=1/16. The fixed sparsity ‖0=k\|f\|_0=k is enforced architecturally by TopKkTopK_k, avoiding the shrinkage bias of ℓ1 _1 penalties. We initialized dictionary directions from random data activations, which proved essential for keeping the dictionary alive (raising the live fraction from 45%45\% to 93%93\% at fixed reconstruction error). 3.2 Composition-matched binding test For each transcription factor we build three classes of length-matched genomic window: bound (ChIP-seq peak centers containing a strong motif), unbound-motif (strong motif occurrences outside all peaks), and background (random windows). A position weight matrix M∈ℝL×4M ^L× 4 from JASPAR [22] is converted to a log-odds matrix against a uniform background, and a window is deemed to contain a strong motif if its maximum log-odds score over all offsets exceeds a fixed fraction of the matrix maximum. To remove the compositional confound, the unbound-motif set is GC-matched to the bound set by rejection sampling against the bound GC histogram. For a feature j we max-pool its activation fjf_j over the tokens of each window and test, with a one-sided Mann–Whitney U [15], whether fjbound>fjunboundf_j^bound>f_j^unbound. We summarize effect size by AUCj=U/(nbnu)AUC_j=U/(n_bn_u), the probability that a random bound window activates the feature more strongly than a random unbound-motif window. A feature is motif-selective if fjunbound>fjbackgroundf_j^unbound>f_j^background and binding-sensitive if AUCj>0.55AUC_j>0.55 at Bonferroni-corrected significance; we additionally require |corr(fj,GC)|<0.2|corr(f_j,GC)|<0.2 for “GC-robust” status. 3.3 Causal feature ablation Let hℓ()h_ (x) be the hidden state at the SAE’s layer and gℓ→L(⋅)g_ → L(·) the remainder of the network mapping that hidden state to output logits over the token vocabulary. For dictionary direction jd_j with encoded activation fjf_j, the ablated hidden state removes the feature’s contribution in the (normalized) SAE space, hℓ′()=hℓ()−s−1(fj()j),h_ (x)=h_ (x)-s^-1\, (f_j(x)\,d_j ), (5) where s is the activation scale used during SAE training. We then complete the forward pass with the perturbed state and compare the model’s predictive distributions at each token position t via the Kullback–Leibler divergence, Δj()=1T∑t=1TKL(pt(gℓ→L(hℓ))∥pt(gℓ→L(hℓ′))), _j(x)= 1T _t=1^TKL\! (p_t (g_ → L(h_ ) )\, \|\,p_t (g_ → L(h_ ) ) ), (6) where pt(⋅)=softmaxp_t(·)=softmax of the logits at position t and T is the number of tokens. The binding-specific causal effect of feature j is the one-sided Mann–Whitney comparison of Δj()\ _j(x)\ between bound and unbound-motif windows, summarized again by an AUCAUC. A feature passes if this AUC>0.55AUC>0.55 at Bonferroni-corrected significance. Algorithm 1 Causal dictionary learning and validation 1:Input: genomic model with layer ℓ ; windows; ChIP-seq peaks; motif PWMs 2:Extract per-token activations \x\ at layer ℓ via forward hooks 3:Train top-k SAE on \x\ by minimizing Eq. (4) ⊳ Eqs. (1)–(3) 4:Build bound / unbound-motif (GC-matched) / background window sets 5:for each live feature j do 6: compute AUCjAUC_j (bound vs. unbound) and corr(fj,GC)corr(f_j,GC) 7: mark j binding-sensitive if AUCj>0.55AUC_j>0.55 (Bonferroni) and |corr|<0.2|corr|<0.2 8:end for 9:for each binding-sensitive feature j (and random controls) do 10: ablate jd_j (Eq. (5)); measure Δj _j (Eq. (6)) at bound vs. unbound 11: mark j causally validated if binding-specific AUC>0.55AUC>0.55 (Bonferroni) 12:end for 13:Output: causally validated binding features per factor and model 3.4 Theoretical analysis We formalize why the causal ablation test isolates binding-specific use, and why the composition-matched design is necessary. Throughout, let B∈0,1B∈\0,1\ indicate ChIP-seq binding, let S denote the observable sequence of a window, and let fj(S)f_j(S) be feature j’s (max-pooled) activation. Proposition 1 (Spurious enrichment from composition). Let ρ(S)∈[0,1]ρ(S)∈[0,1] be a compositional summary of S (e.g. GC fraction), and suppose a motif log-odds score σ(S)σ(S) satisfies [σ(S)∣ρ(S)=r]E[σ(S) ρ(S)=r] strictly increasing in r. If a feature is a pure composition detector, fj(S)=ϕ(ρ(S))f_j(S)=φ(ρ(S)) for some increasing ϕφ, then in an unmatched comparison where the high-activation set has higher mean composition than background, the feature’s top windows are enriched for high σ, i.e. [σ(S)|fj(S)large]>[σ(S)],E [σ(S)\, |\,f_j(S)\ large ]>E [σ(S) ], even though fj⟂B∣ρf_j B ρ, i.e. the feature carries no information about binding beyond composition. Proof. Since fj=ϕ(ρ)f_j=φ(ρ) with ϕφ increasing, the event fj large\f_j large\ equals ρ>r0\ρ>r_0\ for a threshold r0r_0 above the marginal median of ρ. By the monotonicity of [σ∣ρ=r]E[σ ρ=r] in r and the law of total expectation, [σ∣ρ>r0]=∫r01[σ∣ρ=r]Fρ∣ρ>r0(r)>∫01[σ∣ρ=r]Fρ(r)=[σ],E[σ ρ>r_0]=\! _r_0^1\!E[σ ρ=r]\,dF_ρ ρ>r_0(r)>\! _0^1\!E[σ ρ=r]\,dF_ρ(r)=E[σ], because conditioning on ρ>r0ρ>r_0 stochastically increases ρ and [σ∣ρ=r]E[σ ρ=r] is increasing. The independence fj⟂B∣ρf_j B ρ holds by construction since fjf_j is a deterministic function of ρ alone. Hence the feature shows motif enrichment while being conditionally uninformative about binding. ∎ Proposition 1 formalizes the empirical failure of naive enrichment and motivates comparing bound against GC-matched unbound windows: matching equalizes the marginal law of ρ across groups, so a pure composition detector satisfies [fj∣B=1]=[fj∣B=0]E[f_j B=1]=E[f_j B=0] in the matched population and its binding-sensitivity AUC→1/2AUC→ 1/2. The composition-matched test is therefore calibrated against this confound. Definition 1 (Causal binding use). Feature j is causally used for binding if its ablation changes the model’s output distribution more on bound than on unbound-motif inputs, i.e. [Δj∣B=1]>[Δj∣B=0]E[ _j B=1]>E[ _j B=0] on motif-matched inputs. Proposition 2 (Separation of importance from binding-specificity). Let the contribution removed by ablation be cj(S)=s−1fj(S)jc_j(S)=s^-1f_j(S)d_j and write the post-layer map’s local sensitivity as J(S)=∂gℓ→L/∂hℓJ(S)=∂ g_ → L/∂ h_ . To first order, Δj(S)≈12cj(S)⊤G(S)cj(S), _j(S)\ ≈\ 12\,c_j(S) G(S)\,c_j(S), (7) where G(S)=J(S)⊤HKLJ(S)⪰0G(S)=J(S) H_KL\,J(S) 0 and HKLH_KL is the Fisher information of the softmax output. Then a feature with large average effect [Δj]E[ _j] but fj⟂Bf_j B (motif-only) has [Δj∣B=1]=[Δj∣B=0]E[ _j B=1]=E[ _j B=0] on motif-matched inputs, hence binding-specificity AUC=1/2AUC=1/2; whereas a feature whose activation is elevated on bound inputs, [fj∣B=1]>[fj∣B=0]E[f_j B=1]>E[f_j B=0], yields [Δj∣B=1]>[Δj∣B=0]E[ _j B=1]>E[ _j B=0] whenever G(S)G(S) does not anti-correlate with binding. Proof. Equation (7) is the second-order Taylor expansion of the KL divergence between pt(g(hℓ))p_t(g(h_ )) and pt(g(hℓ−cj))p_t(g(h_ -c_j)) about cj=0c_j=0; the first-order term vanishes because the KL divergence and its gradient are zero at cj=0c_j=0, and the quadratic form uses the Fisher information HKLH_KL as the Hessian of the KL divergence at the unperturbed distribution, propagated through J(S)J(S) by the chain rule. Taking cj(S)=s−1fj(S)jc_j(S)=s^-1f_j(S)d_j gives Δj(S)≈12s−2fj(S)2j⊤G(S)j _j(S)≈ 12s^-2f_j(S)^2\,d_j G(S)d_j. Write γ(S)=12s−2j⊤G(S)j≥0γ(S)= 12s^-2d_j G(S)d_j≥ 0. For a motif-only feature, fj⟂Bf_j B on motif-matched inputs, and if the geometric factor γ does not anti-correlate with binding (in particular if it is independent of B, true to first order when ablation does not systematically interact with binding), then [fj2γ∣B=1]=[fj2γ∣B=0]E[f_j^2γ B=1]=E[f_j^2γ B=0], giving equal conditional effects and AUC=1/2AUC=1/2. For a binding-sensitive feature, [fj2∣B=1]>[fj2∣B=0]E[f_j^2 B=1]>E[f_j^2 B=0] since fj≥0f_j≥ 0 and its mean is larger on bound inputs; multiplying by the nonnegative, binding-independent γ preserves the inequality, so [Δj∣B=1]>[Δj∣B=0]E[ _j B=1]>E[ _j B=0] and the binding-specificity AUC>1/2AUC>1/2. ∎ Proposition 2 predicts a counterintuitive phenomenon we observe empirically (Section 5): a motif-only feature can have a far larger absolute ablation effect than a binding feature, yet a binding-specificity AUCAUC at chance, because absolute effect is governed by [fj2γ]E[f_j^2γ] while specificity is governed by the difference of this quantity across binding states. The causal test reads the latter. 4 Experimental Setup Models. We used the publicly released checkpoints of Nucleotide Transformer (NT) v2 (500500M parameters, 66-mer tokenization, 2424 transformer blocks) [6] and DNABERT-2 (117117M parameters, byte-pair encoding, 1212 blocks) [31]. Both expose a masked-language-model head, enabling the logit-level causal readout of Eq. (6). Data. The human reference genome (hg38), ENCODE SCREEN candidate cis-regulatory elements [17], JASPAR 2024 CORE vertebrate motifs [22], and ENCODE ChIP-seq IDR peak sets for CTCF (GM12878), GATA1 (K562) and REST (K562) were obtained from their public repositories. For activation extraction we sampled 5×1045× 10^4 cis-regulatory and 5×1045× 10^4 random 200200 bp windows, retaining only windows with unambiguous nucleotides, yielding 3.63.6–4.2×1064.2× 10^6 token activations per layer. Because the two models tokenize DNA differently, NT into non-overlapping 66-mers, DNABERT-2 into variable-length byte-pair tokens, we recorded, for every token, the exact base-pair span it covers, mapping all features to genomic coordinates regardless of tokenization. Special tokens were excluded from all analyses. Training and layer selection. For each model we sampled four evenly spaced layers and trained one top-k SAE per layer (dictionary size 16×16× the hidden width; k=32k=32) for one epoch with Adam (learning rate 4×10−44× 10^-4, batch size 40964096). Layer selection used reconstruction fidelity and live-feature fraction (Section 5, Fig. 1). Binding sets and controls. Per-factor window sets were built as in Methodology, with the unbound-motif set GC-matched to bound. Two negative controls were processed by the identical pipeline: scramble (random windows with random bound/unbound labels) and gata1-scram (real GATA1 motif windows with randomized labels, so that sequence and motif content are unchanged and only the binding label is destroyed). Population causal analyses ablate the top 1515 binding-sensitive, GC-robust features and 1515 randomly selected live features per condition. Statistics. All P-values use the one-sided Mann–Whitney U test [15] and are Bonferroni-corrected for the number of features tested; significance thresholds are stated with each analysis. 5 Results 5.1 Sparse dictionaries are monosemantic and layer-dependent The autoencoders reconstructed activations faithfully while remaining sparse: the best dictionaries explained 7474–76%76\% of activation variance with only 3232 of ∼ 12,000–16,000 features active per token. Reconstruction quality and the fraction of “live” (ever-active) features both peaked at intermediate depth (Fig. 1), mirroring observations in language models [3]: NT layer 1414 (50%50\% depth) and DNABERT-2 layer 66 (50%50\% depth) provided the best balance of fidelity and dictionary utilization and were used in all subsequent analyses. Both tokenization schemes therefore admit high-quality sparse dictionaries, despite their very different vocabularies. Figure 1: Sparse-autoencoder reconstruction across model depth. Variance explained (solid) and fraction of live dictionary features (dashed) for Nucleotide Transformer (blue) and DNABERT-2 (red) as a function of relative layer depth. Both architectures attain their best trade-off of reconstruction fidelity and feature utilization at intermediate depth (shaded), motivating the use of NT layer 1414 and DNABERT-2 layer 66. 5.2 Naive motif enrichment is confounded by sequence composition To ask whether dictionary features correspond to regulatory grammar, we first applied the obvious test: for each feature, we collected the genomic windows on which it fired most strongly and scored them against transcription-factor position weight matrices [22]. Using CTCF, a ubiquitous insulator with a long, information-rich motif [21], this naive procedure flagged 440440 of 6,5766,576 tested features as significantly CTCF-enriched. Inspection revealed this to be largely an artifact. The single most “CTCF-enriched” feature fired almost exclusively on near-identical copies of the Alu consensus, the most abundant repetitive element in the human genome [7], which scores above background on the GC-rich CTCF matrix without representing CTCF binding; other top features detected generic GC-rich or CpG-island sequence. Because CTCF, CpG islands and Alu elements are all GC-rich, a feature that merely tracks base composition is spuriously “enriched” for CTCF. Comparing each feature’s motif scores against a GC-matched background rather than an unmatched one collapsed the apparent enrichment of these features toward zero, identifying them as compositional, as predicted by Proposition 1. The count of motif-selective features for each factor after composition control is summarized in Fig. 2. Figure 2: Motif-selective features per transcription factor after composition control, for NT (blue) and DNABERT-2 (red). CTCF and REST, long, information-rich motifs, yield many motif-selective features; GATA1, with a short WGATAR motif, yields almost none (33 and 0, respectively), foreshadowing that its binding is encoded through context rather than the canonical motif. 5.3 A composition-matched, binding-resolved test isolates genuine signal The deeper question is whether a feature encodes binding, cell-type-specific occupancy measured by ChIP-seq, rather than mere motif presence. Testing whether each feature’s activation is higher on bound than on GC-matched unbound-motif windows, NT layer 1414 contained 265265 motif-selective CTCF features but only a small number that remained binding-sensitive after GC control; the cleanest, feature 80878087, exhibited a monotonic activation gradient from background to unbound-motif to bound windows (Fig. 3a) with negligible GC correlation. A stringent, composition-aware test therefore recovers a specific, interpretable, binding-associated feature where the naive test had returned hundreds of confounded ones. 5.4 Feature ablation establishes causal use of binding representations Correlation between activation and binding does not establish that the model uses a feature to represent binding. We therefore ablated single dictionary directions during the forward pass (Eq. (5)) and measured the Kullback–Leibler shift in masked-token predictions (Eq. (6)). For the exemplar CTCF feature, ablation shifted NT’s predictions significantly more at bound than at unbound sites (AUC=0.63AUC=0.63; P=1.7×10−18P=1.7× 10^-18). Two controls confirmed specificity (Fig. 3b). Ablating a random feature produced no binding-specific effect (AUC=0.51AUC=0.51). More tellingly, ablating a motif-selective but not binding-sensitive feature produced a large overall prediction shift, an order of magnitude larger in absolute terms than the binding feature, yet no bound-versus-unbound differential (AUC=0.51AUC=0.51), exactly the dissociation predicted by Proposition 2. The causal test thus distinguishes features that matter to the model in general from features that matter specifically for binding, which a magnitude-based analysis cannot. Figure 3: A causally validated CTCF binding feature. (a) Mean activation of two exemplar binding features (DNABERT-2 feature 857857, red; NT feature 80878087, blue) across background, unbound-motif and ChIP-seq-bound windows; activation increases monotonically with binding, not merely with motif presence. (b) Causal effect of ablating each of three NT features, measured as the prediction-distribution shift (Kullback–Leibler divergence) at bound (blue) versus unbound-motif (grey) sites. The binding feature shows a binding-specific effect (AUC=0.63AUC=0.63); a random feature shows none; a motif-only feature shows a large overall effect but no binding specificity (AUC=0.51AUC=0.51), demonstrating that the test isolates binding-specific causal use rather than general feature importance. 5.5 Causal binding features are reproducible across factors and architectures We next asked whether causal binding features are a population phenomenon and whether the framework generalizes. For each transcription factor and model we ablated the top 1515 binding-sensitive, GC-robust features and 1515 random features. Across three factors of distinct structural class, CTCF (insulator), GATA1 (lineage-specific activator) and REST (repressor), and both architectures, a substantial fraction of binding features showed significant binding-specific causal effects, whereas random-feature controls showed essentially none (Fig. 4, Fig. 5, Table 1). The effect was, if anything, stronger in DNABERT-2 than in NT, indicating that byte-pair tokenization does not impede, and may sharpen, the emergence of binding representations. A biologically informative subtlety emerged for GATA1. Its binding features were only weakly motif-selective in the correlational test (Fig. 2), which we had initially read as weak representation. The causal test overturned this: GATA1 binding features were robustly causal in both models (9/159/15 and 13/1513/15; Table 1), indicating that the models encode GATA1 occupancy through contextual sequence features beyond the canonical WGATAR motif, signal that the causal test detects but motif-based correlation misses. This illustrates the added value of intervention over association. Figure 4: Causal validation matrix. Number of features (of 1515 tested) with a significant binding-specific causal effect (Bonferroni-corrected P<1.7×10−3P<1.7× 10^-3 and AUC>0.55AUC>0.55), for each transcription factor (rows) and model (columns). Random-feature controls yielded 0/150/15 in every cell. The two bottom rows are negative controls in which binding labels were scrambled; these produced no binding-sensitive features and hence nothing to test (“0 feat.”), confirming that the framework reports a true null when no real binding signal is present. Figure 5: Causal effect sizes separate binding features from controls. Distribution of causal-effect AUCAUC (bound versus unbound motif) for the top 1515 binding features of each transcription factor and for pooled random-feature controls, in NT (left) and DNABERT-2 (right). Target features lie above the significance threshold (AUC=0.55AUC=0.55, red dotted) while random controls cluster at chance (AUC=0.5AUC=0.5, grey dashed). Points are individual features; boxes show median and interquartile range. 5.6 Negative controls confirm that the framework reports true nulls A validation method is only trustworthy if it can return a negative. In both negative controls, scramble (random windows, random labels) and the sharper gata1-scram (real GATA1 motif windows, randomized labels), and in both architectures, the binding-sensitivity test identified zero qualifying features, leaving nothing for the causal stage to test (Table 1, Fig. 4, bottom rows). The framework thus manufactures no signal from structured-but-unlabeled data: the positive results for real transcription factors reflect genuine binding information, not artifacts of window construction or feature ablation. Table 1: Causally validated binding features across transcription factors, models and controls. Entries give the number of features (of 1515 tested) with a significant binding-specific causal effect (Bonferroni P<1.7×10−3P<1.7× 10^-3, AUC>0.55AUC>0.55); random-feature controls were 0/150/15 throughout. Negative controls (scrambled labels) yielded no binding-sensitive features to test. “Motif-sel.” gives the number of motif-selective features identified by the composition-matched test. Dashes ‘—’ denote not-applicable cells: the causal test is undefined for negative controls (no binding-sensitive features exist to test), and motif-selectivity is undefined for randomly chosen features. NT (6-mer, layer 14) DNABERT-2 (BPE, layer 6) Condition Causal / 15 Motif-sel. Causal / 15 Motif-sel. CTCF (insulator) 77 265265 1010 8888 GATA1 (activator) 99 33 1313 0 REST (repressor) 77 367367 1414 135135 scramble (neg.) — 0 — 0 gata1-scram (neg.) — 0 — 0 Random-feature control 0 — 0 — 6 Discussion Three findings have implications beyond our specific models. First, the dominant obstacle to interpreting genomic models is not extracting features but validating them: position-weight-matrix enrichment, the field’s default, is confounded by GC content and repetitive elements to the point of returning hundreds of spurious “TF features”, a failure we both observed empirically and proved is guaranteed under mild conditions (Proposition 1). Composition-matched, binding-resolved testing is therefore not optional but necessary. Second, association and causation diverge in an informative way. The GATA1 case exhibits a factor whose binding is causally encoded yet poorly captured by motif correlation, because the model represents occupancy through contextual sequence beyond the canonical motif; only intervention reveals this, and Proposition 2 explains why magnitude-based attribution would have missed it. Third, the emergence of binding features in a model trained purely on reference sequence, with no cell-type labels, indicates that self-supervised genomic models internalize determinants of cell-type-specific occupancy latent in sequence context, an observation that invites systematic study. Our analysis has limitations that also define natural extensions. We focused on two masked-language-model encoders; autoregressive and state-space genomic models [18, 23, 10] expose different readouts and merit dedicated treatment. We examined three well-characterized factors; the framework scales directly to the hundreds of factors with ENCODE ChIP-seq, enabling a systematic atlas of which regulatory programs genomic models encode and use. Finally, the same causal machinery could be turned from validation to steering, editing binding features to test sequence-design hypotheses in silico. 7 Conclusion We presented a purely computational framework that extracts interpretable features from genomic language models, validates them against experimental binding data while controlling for the compositional confounds that plague naive analyses, and tests by direct intervention whether the model actually uses each feature to represent binding. The framework recovers reproducible, causally validated transcription-factor binding features across two architectures with different tokenizations and across three transcription factors of distinct regulatory function, while returning clean nulls under two classes of negative control. By coupling sparse dictionary learning to causal intervention and disciplined controls, it offers the regulatory-genomics community a reusable standard for interpretability claims: a feature is real not when it correlates with a motif, but when the model demonstrably uses it. References [1] Ž. Avsec, V. Agarwal, D. Visentin, J. R. Ledsam, A. Grabska-Barwinska, K. R. Taylor, Y. Assael, J. Jumper, P. Kohli, and D. R. Kelley (2021) Effective gene expression prediction from sequence by integrating long-range interactions. Nature methods 18 (10), p. 1196–1203. Cited by: §1, §2. [2] G. Benegas, S. S. Batra, and Y. S. Song (2023) DNA language models are powerful predictors of genome-wide variant effects. Proceedings of the National Academy of Sciences 120 (44), p. e2311219120. Cited by: §1, §2. [3] T. Bricken, A. Templeton, J. Batson, B. Chen, A. Jermyn, T. Conerly, N. Turner, C. Anil, C. Denison, A. Askell, R. Lasenby, Y. Wu, S. Kravec, N. Schiefer, T. Maxwell, N. Joseph, Z. Hatfield-Dodds, A. Tamkin, K. Nguyen, B. McLean, J. E. Burke, T. Hume, S. Carter, T. Henighan, and C. Olah (2023) Towards monosemanticity: decomposing language models with dictionary learning. Transformer Circuits Thread. External Links: Link Cited by: §1, §2, §5.1. [4] E. P. Consortium et al. (2012) An integrated encyclopedia of dna elements in the human genome. Nature 489 (7414), p. 57. Cited by: §2. [5] H. Cunningham, A. Ewart, L. Riggs, R. Huben, and L. Sharkey (2023) Sparse autoencoders find highly interpretable features in language models. arXiv preprint arXiv:2309.08600. Cited by: §1, §2. [6] H. Dalla-Torre, L. Gonzalez, J. Mendoza-Revilla, N. Lopez Carranza, A. H. Grzywaczewski, F. Oteri, C. Dallago, E. Trop, B. P. De Almeida, H. Sirelkhatim, et al. (2025) Nucleotide transformer: building and evaluating robust foundation models for human genomics. Nature Methods 22 (2), p. 287–297. Cited by: §1, §2, §4. [7] P. Deininger (2011) Alu elements: know the sines. Genome biology 12 (12), p. 236. Cited by: §1, §2, §5.2. [8] N. Elhage, T. Hume, C. Olsson, N. Schiefer, T. Henighan, S. Kravec, Z. Hatfield-Dodds, R. Lasenby, D. Drain, C. Chen, et al. (2022) Toy models of superposition. arXiv preprint arXiv:2209.10652. Cited by: §1, §2. [9] L. Gao, T. Dupre la Tour, H. Tillman, G. Goh, R. Troll, A. Radford, I. Sutskever, J. Leike, and J. Wu (2025) Scaling and evaluating sparse autoencoders. In International Conference on Learning Representations, Vol. 2025, p. 26721–26754. Cited by: §1, §2, §3.1. [10] A. Gu and T. Dao (2023) Mamba: linear-time sequence modeling with selective state spaces. arXiv preprint arXiv:2312.00752. Cited by: §2, §6. [11] Y. Ji, Z. Zhou, H. Liu, and R. V. Davuluri (2021) DNABERT: pre-trained bidirectional encoder representations from transformers model for dna-language in genome. Bioinformatics 37 (15), p. 2112–2120. Cited by: §2. [12] D. R. Kelley, Y. A. Reshef, M. Bileschi, D. Belanger, C. Y. McLean, and J. Snoek (2018) Sequential regulatory activity prediction across chromosomes with convolutional neural networks. Genome research 28 (5), p. 739. Cited by: §1, §2. [13] S. M. Lundberg and S. Lee (2017) A unified approach to interpreting model predictions. Advances in neural information processing systems 30. Cited by: §2. [14] A. Makelov, G. Lange, and N. Nanda (2023) Is this the subspace you are looking for? an interpretability illusion for subspace activation patching. arXiv preprint arXiv:2311.17030. Cited by: §2. [15] H. B. Mann and D. R. Whitney (1947) On a test of whether one of two random variables is stochastically larger than the other. The annals of mathematical statistics, p. 50–60. Cited by: §3.2, §4. [16] K. Meng, D. Bau, A. Andonian, and Y. Belinkov (2022) Locating and editing factual associations in gpt. Advances in neural information processing systems 35, p. 17359–17372. Cited by: §2. [17] J. E. Moore, M. J. Purcaro, H. E. Pratt, C. B. Epstein, N. Shoresh, J. Adrian, T. Kawli, C. A. Davis, A. Dobin, et al. (2020) Expanded encyclopaedias of dna elements in the human and mouse genomes. Nature 583 (7818), p. 699–710. Cited by: §2, §4. [18] E. Nguyen, M. Poli, M. G. Durrant, B. Kang, D. Katrekar, D. B. Li, L. J. Bartie, A. W. Thomas, S. H. King, G. Brixi, et al. (2024) Sequence modeling and design from molecular to genome scale with evo. Science 386 (6723), p. eado9336. Cited by: §1, §2, §6. [19] E. Nguyen, M. Poli, M. Faizi, A. Thomas, M. Wornow, C. Birch-Sykes, S. Massaroli, A. Patel, C. Rabideau, Y. Bengio, et al. (2023) Hyenadna: long-range genomic sequence modeling at single nucleotide resolution. Advances in neural information processing systems 36, p. 43177–43201. Cited by: §1, §2. [20] G. Novakovsky, N. Dexter, M. W. Libbrecht, W. W. Wasserman, and S. Mostafavi (2023) Obtaining genetics insights from deep learning via explainable artificial intelligence. Nature Reviews Genetics 24 (2), p. 125–137. Cited by: §1, §2. [21] C. Ong and V. G. Corces (2014) CTCF: an architectural protein bridging genome topology and function. Nature Reviews Genetics 15 (4), p. 234–246. Cited by: §2, §5.2. [22] I. Rauluseviciute, R. Riudavets-Puig, R. Blanc-Mathieu, J. A. Castro-Mondragon, K. Ferenc, V. Kumar, R. B. Lemma, J. Lucas, J. Chèneby, D. Baranasic, et al. (2024) JASPAR 2024: 20th anniversary of the open-access database of transcription factor binding profiles. Nucleic acids research 52 (D1), p. D174–D182. Cited by: §2, §3.2, §4, §5.2. [23] Y. Schiff, C. Kao, A. Gokaslan, T. Dao, A. Gu, and V. Kuleshov (2024) Caduceus: bi-directional equivariant long-range dna sequence modeling. Proceedings of machine learning research 235, p. 43632. Cited by: §1, §2, §6. [24] F. W. Schmitges, E. Radovani, H. S. Najafabadi, M. Barazandeh, L. F. Campitelli, Y. Yin, A. Jolma, G. Zhong, H. Guo, T. Kanagalingam, et al. (2016) Multiparameter functional diversity of human c2h2 zinc finger proteins. Genome research 26 (12), p. 1742. Cited by: §1, §2. [25] A. Shrikumar, P. Greenside, and A. Kundaje (2017) Learning important features through propagating activation differences. In International conference on machine learning, p. 3145–3153. Cited by: §2. [26] Z. Tang, N. Somia, Y. Yu, and P. K. Koo (2025) Evaluating the representational power of pre-trained dna language models for regulatory genomics. Genome Biology 26 (1), p. 203. Cited by: §2. [27] A. Templeton, T. Conerly, J. Marcus, J. Lindsey, T. Bricken, B. Chen, A. Pearce, C. Citro, E. Ameisen, A. Jones, et al. (2024) Scaling monosemanticity: extracting interpretable features from Claude 3 Sonnet. Transformer Circuits Thread. External Links: Link Cited by: §1, §2. [28] J. Vig, S. Gehrmann, Y. Belinkov, S. Qian, D. Nevo, Y. Singer, and S. Shieber (2020) Investigating gender bias in language models using causal mediation analysis. Advances in neural information processing systems 33, p. 12388–12401. Cited by: §2. [29] I. E. Vorontsov, I. A. Eliseeva, A. Zinkevich, M. Nikonov, S. Abramov, A. Boytsov, V. Kamenets, A. Kasianova, S. Kolmykov, I. S. Yevshin, et al. (2024) HOCOMOCO in 2024: a rebuild of the curated collection of binding models for human and mouse transcription factors. Nucleic Acids Research 52 (D1), p. D154–D163. Cited by: §2. [30] J. Zhou and O. G. Troyanskaya (2015) Predicting effects of noncoding variants with deep learning–based sequence model. Nature methods 12 (10), p. 931–934. Cited by: §1, §2. [31] Z. Zhou, Y. Ji, W. Li, P. Dutta, R. Davuluri, and H. Liu (2024) DNABERT-2: efficient foundation model and benchmark for multi-species genomes. In International Conference on Learning Representations, Vol. 2024, p. 41642–41665. Cited by: §1, §2, §4.