Expanding the DNA Motif Lexicon of the Transcriptional Regulatory Code

Expanding the DNA Motif Lexicon of the Transcriptional Regulatory Code | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Expanding the DNA Motif Lexicon of the Transcriptional Regulatory Code View ORCID Profile Jingyu Fan , Virendra K. Chaudhri , Deepa Bisht , Nicholas A. Pease , David Reed Hall , Peter Gerges , Vivian Fang Yi , Susan Kales , Ching-Huang Ho , Jishnu Das , John P. Ray , Ryan Tewhey , Clifford A. Meyer , Nidhi Sahni , View ORCID Profile Harinder Singh doi: https://doi.org/10.1101/2025.07.09.662874 Jingyu Fan 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jingyu Fan Virendra K. Chaudhri 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Deepa Bisht 2 Department of Epigenetics and Molecular Carcinogenesis, Department of Bioinformatics and Computational Biology, The University of Texas, MD Anderson Cancer Center , Houston, TX, USA 3 Department of Genitourinary Medical Oncology, Division of Cancer Medicine, The University of Texas, MD Anderson Cancer Center , Houston, TX, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Nicholas A. Pease 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site David Reed Hall 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Peter Gerges 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Vivian Fang Yi 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Susan Kales 4 The Jackson Laboratory , Bar Harbor, ME, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ching-Huang Ho 5 Center for Systems Immunology, Benaroya Research Institute , Seattle, WA, USA 6 Department of Immunology, University of Washington , Seattle, WA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Jishnu Das 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site John P. Ray 5 Center for Systems Immunology, Benaroya Research Institute , Seattle, WA, USA 6 Department of Immunology, University of Washington , Seattle, WA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ryan Tewhey 4 The Jackson Laboratory , Bar Harbor, ME, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Clifford A. Meyer 7 Department of Data Sciences, Dana-Farber Cancer Institute and Harvard T.H. Chan School of Public Health , Boston, MA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: cliff_meyer{at}ds.dfci.harvard.edu nidhi.sahni.2025{at}gmail.com harinder{at}pitt.edu Nidhi Sahni 2 Department of Epigenetics and Molecular Carcinogenesis, Department of Bioinformatics and Computational Biology, The University of Texas, MD Anderson Cancer Center , Houston, TX, USA 8 Quantitative and Computational Biosciences Program, Baylor College of Medicine , Houston, TX, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: cliff_meyer{at}ds.dfci.harvard.edu nidhi.sahni.2025{at}gmail.com harinder{at}pitt.edu Harinder Singh 1 Center for Systems Immunology and the Departments of Immunology and Computational and Systems Biology, University of Pittsburgh , PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Harinder Singh For correspondence: cliff_meyer{at}ds.dfci.harvard.edu nidhi.sahni.2025{at}gmail.com harinder{at}pitt.edu Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Transcriptional regulatory sequences in metazoans contain intricate combinations of transcription factor (TF) motifs. Stereospecific arrangements of simple motifs constitute composite elements (CEs) that enhance DNA-protein interaction specificity and enable combinatorial regulatory logic. Despite their importance, CEs remain underexplored. We advance CE discovery and functional characterization by developing an integrated framework that combines computational prediction, experimental testing and deep learning. The extended TF motif catalog comprises both synergistic and counteracting CEs, which are supported by evidence of TF binding in vivo and in vitro . A deep learning model GRACE trained on customized massively parallel reporter assays learns the lexicon of CEs at single-nucleotide resolution. Comparative analysis with a neural network model trained on chromatin accessibility demonstrates striking convergence and distinctions within the expanded regulatory lexicon, enabling joint predictions of motif contributions and the impact of variants on chromatin structure and transcriptional activity in diverse cellular contexts. Introduction DNA‒protein interactions mediated by transcription factors (TFs) are central for decoding regulatory information in metazoan genomes 1 – 5 . Structurally related TFs bind simple DNA motifs (6–10 bp) that have been extensively cataloged using a combination of experimental and computational approaches 3 , 6 . These motifs are the fundamental building blocks of the gene regulatory code and constitute its primary lexicon. However, metazoan promoters and enhancers consist of intricate combinations of such simple TF motifs, thereby generating a contextually rich regulatory syntax that has proven challenging to decipher 2 , 5 . The syntax determines the diverse spatiotemporal as well as signaling-responsive properties of such regulatory DNA sequences. The poor decipherability of the regulatory syntax is compounded by the widespread use of simple TF motif catalogs to generate TF-to-target gene inferences and to assemble and analyze cell type- or cell state-specific gene regulatory networks (GRNs) 7 – 10 . As simple TF motifs display considerable redundancy and have low information content, this impairs the quality of predictions of GRN models based on such motif catalogs. Transcriptional regulatory sequences have been shown to contain striking instances of composite elements (CEs) that are constituted by precise juxta positioning of simple motifs that enable cooperative binding of TFs representing distinct structural families. Exemplars include the ETS::IRF CE (EICE) motif recognized by specific members of the ETS and IRF families of TFs, namely, PU.1/SPIB and IRF4/IRF8 11 ; the AP-1::IRF CE (AICE) motif recognized by AP-1 (BATF/JUN) and IRF4/IRF8 12 ; the NFAT::AP1 CE bound by NFAT and AP-1 paralogs 13 , 14 ; and the SOX::OCT CEs involving OCT4 and various SOX family members 15 . CEs are stereo-specifically constrained by the relative orientation and spacing of their constituent TF motifs, thereby manifesting greater information content and specificity for encoding protein-DNA interactions. Thus, CEs could represent a higher-order lexicon in the grammar of the gene regulatory code underlying combinatorial as well as spatiotemporal logic. Their bipartite structure is suited for the integration of cell-extrinsic or cell-intrinsic inputs into the genome by specific combinations of signaling-regulated TFs. Despite their potential biological importance, the prevalence and diversity of CEs in genomic regulatory sequences remain underexplored. Previous studies examining the prevalence of CEs have focused on either high-throughput in vitro DNA binding assays with recombinant proteins (CAP-SELEX) or computational analyses of accessible or TF-bound chromatin regions 16 , 17 . Deep learning computational methods are proving to be instrumental in analyzing regulatory DNA sequences, predicting TF binding sites, and modeling chromatin accessibility with unprecedented accuracy 18 – 20 . The hierarchical structure of deep neural networks mirrors the complexity of genomic sequences, facilitating the learning of both simple TF motifs and higher-order motif interactions, such as CEs, thereby accelerating the elucidation of the transcriptional regulatory lexicon and syntax. Despite this progress, CEs remain to be systematically analyzed by coupling structural and functional genomic analyses with sequence-based deep learning. Advances in genomic analyses, particularly chromatin accessibility profiles as well as massively parallel reporter assays (MPRAs), now provide an unprecedented opportunity to systematically explore the diversity of CEs and their functional activities across cell types 21 . These DNA sequence datasets, coupled with their chromatin accessibility and transcriptional activity profiles, would in turn enable training of deep learning models that could predict the context-specific molecular activities of CEs in regulating both chromatin structure and transcription at endogenous regulatory sequences. We hypothesized that CEs are diverse, widespread and represent a distinctive higher-order lexicon of the gene regulatory code. To explore this hypothesis systematically, we pursued an integrated analytical and experimental framework ( Figure 1 ). We first developed and utilized a computational pipeline, CEseek, for CE discovery and then experimentally tested the large set of candidate CEs (cCEs) by performing MPRAs 21 with a custom-designed TF-motif library. The transcriptionally active CEs were validated by analyses of TF co-binding in vivo (TF ChIP-seq) and in vitro (CAP-SELEX). Next, we developed and trained a deep learning model, GRACE, using the TF-motif MPRA datasets to accurately predict the transcriptional activity of CEs and simple motifs at single-nucleotide resolution. GRACE enabled predictions of active transcription at endogenous promoters and enhancers, as well as allele-specific regulatory effects. The transcription activity-based CE motif lexicon learned by GRACE was used to systematically assess its contributions to chromatin accessibility by utilizing ChromBPNet 20 , an orthogonal deep learning model. Despite being trained on distinct DNA-associated features, the two models remarkably converged on an extensive set of higher-order TF motif configurations, thereby demonstrating the importance of the expanded transcriptional regulatory lexicon in regulating chromatin structure and transcription. The integrated framework enables predictions of context-specific molecular activities of CEs and distinguishes the contributions of constituent motifs in regulating chromatin structure and/or transcription at endogenous regulatory sequences. The CE regulatory lexicon will greatly facilitate analyses of combinatorial control of spatiotemporal patterns of gene activity as well as the integration of signaling inputs into the genome by distinct signaling-regulated TFs. Download figure Open in new tab Figure 1. Schematic framework for discovering, testing, and modeling composite elements (A) Discover -CEseek identifies candidate composite elements (cCEs) by searching for statistically enriched binary motif configurations by systematically varying the spacing and orientation of simple TF motifs in open chromatin regions (OCRs). (B) Test - Massively Parallel Reporter Assays (MPRAs) measure transcriptional activity of cCEs and the simple motifs comprising them along with their mutated counterparts. (C) Learn and Predict - a deep learning model, GRACE, trained on MPRA datasets that predicts the transcriptional activity of CEs and simple motifs at single-nucleotide resolution. Cellular context-specific models are used to predict transcriptional activity of simple motifs and CEs within endogenous open chromatin regions (OCRs). Results Discovery of Composite Elements using CEseek We developed a computational pipeline named CEseek to identify statistically enriched candidate CEs (cCEs) within cis-regulatory elements (cCREs) delineated by genomic profiling methods such as ATAC-seq, DNase-seq or ChIP-seq 22 , 23 . Heretofore, we will distinguish between cCEs and CEs on the basis that the former display stereo-specific sequence conservation, whereas the latter satisfy one or more additional conditions, such as cooperative or anti-cooperative TF binding and/or synergistic or counteracting activity. CEseek scans a test set of regulatory sequences alongside a background set to identify statistically significant binary configurations using pairs of TF motifs, systematically analyzing occurrences of distinct orientations and user-defined spacings. Known position weight matrices (PWMs) for simple TF motifs are utilized to identify enriched pairs among variable configurations of spacing and orientations of the individual motifs. Specifically, for non-palindromic TF motifs, each spacing yields four potential configurations, and this number is reduced to two for palindromic motifs ( Figure 1 , left panel). Statistically enriched CE configurations are identified by comparing their occurrence in the input set (e.g., open chromatin regions (OCRs)) with that in a background set (e.g., dinucleotide-shuffled sequences, non-OCRs or OCRs from other cell types/states). Sequences of cCEs that are statistically enriched in the input set in relation to the background set are then used to generate an output file of cCEs along with their PWMs. To evaluate CEseek, murine B-cell CREs, previously characterized as active enhancers via STARR-seq 24 , were analyzed for the recovery of EICE and AICE motifs. These CEs have been shown to occur in B-cell CREs and are bound by the transcription factors SPI1 (ETS) and IRF4 or BATF (AP-1) and IRF4, respectively 12 , 25 . For this test, ETS, AP-1, and IRF motifs (with the IRF motif represented as a half-site based on prior knowledge) were used as inputs for CEseek, with analysis of 16 possible motif spacings (default setting, −5 to +10 bp) ( Figure 2A ). A 10-fold dinucleotide-shuffled set of sequences was used as the background set. The outputs of CEseek, represented as a 4×16 matrix of statistically enriched configurations along with the PWMs of the most enriched CE configurations, are displayed ( Figures 2B and 2C ). The pipeline successfully identified the ETS::IRF (EICE) 11 and AP-1::IRF (AICE) 12 CEs. The top two EICE configurations maintained relative orientations and a −2 bp spacing but reversed the order of the ETS and IRF half-sites ( Figure 2B ), whereas the AICE motifs displayed altered spacing (−1 and +3 bp) and reversed the orientation of their IRF motif ( Figure 2C ). These distinctive EICE and AICE configurations have been previously validated by genomic, functional and biochemical analyses 11 , 12 , 26 . Thus, CEseek recovered known CEs along with their multiple stereospecific configurations from genomic regulatory sequences. Download figure Open in new tab Figure 2. CEseek enabled discovery of cCEs in diverse immune and hematopoietic cell types (A) Proof-of-concept workflow using CEseek to recover known CEs composed of ETS, AP-1, and IRF motifs present in murine B-cell cis-regulatory elements (CREs). (B, C) Heatmaps depicting the statistical enrichment of the indicated binary motif configurations (varying in orientation and spacing of the indicated simple motifs) in B-cell CREs. The enrichment was calculated using dinucleotide-shuffled sequences as the background. The top panel (B) shows ETS::IRF (EICE) configurations, and the bottom panel (C) shows AP-1::IRF (AICE) configurations. Examples of the top two statistically enriched configurations of ETS::IRF and AP-1::IRF CEs that include well-characterized EICE and AICE composite motifs are displayed alongside the heatmaps. (D) Workflow for the discovery of cCEs in diverse immune and hematopoietic cell types as well as states. A total of 512,595 cCREs across 86 immune and hematopoietic cell types from the ImmGen ATAC-seq peak database were scanned for 153 non-redundant simple TF motifs using CEseek. All four pairwise orientations allowing for spacings of −5 to +10 bp were analyzed for occurrence in cell type/state-specific cCREs, resulting in 22,698 cell type-enriched cCEs (see Methods ). Further filtering based on stereospecific constraints yielded 5,063 cCEs. (E) Histogram displaying the distribution of TF motif pair counts with their numbers of statistically enriched configurations (n = 1 to 64). Significant motif pairs with ≤6 enriched configurations out of a possible of 64 were deemed to be stereospecifically constrained ( Figure 2D ) and selected as cCEs for further analysis. (F) Histogram showing the false discovery rate (FDR) of the similarity between cCEs and their best-matched de novo TF ChIP-seq motifs recovered from CistromeDB. The fraction of cCEs with matches in the TF ChIP-seq datasets is indicated (FDR < 0.05). A representative example of a cCE discovered by CEseek, along with its matched TF ChIP-seq motif, is displayed. (G) Histogram showing the FDR of the similarity between CAP-SELEX heterodimeric motifs and their best-matched cCEs. The fraction of CAP-SELEX heterodimeric motifs with matches to cCEs is indicated (FDR < 0.05). A representative example of a cCE discovered by CEseek, along with its matched CAP-SELEX motif, is displayed. Download figure Open in new tab Figure S2. CEseek enabled discovery of cCEs in diverse immune and hematopoietic cell types (A ) Histogram displaying the distribution of binary motif pair configuration counts (n=22,698, Figure 2D ) with their numbers of statistically enriched configurations (n=1 to 64). The left end of the distribution highlights motif pairs with ≤6 unique configurations out of a possible of 64 that were designated as cCEs ( Figure 2E ). The right end of the distribution highlights motif pairs with >50 statistically enriched configurations out of a possible of 64. (B, C) Histograms displaying the length distributions of 153 simple TF motifs (B) and 5,063 cCEs (C) . To systematically analyze the identities and occurrence of cCEs in a large set of presumptive transcriptional regulatory sequences, CEseek was used to perform an “all-by-all” CE scan of the ImmGen cis-Regulatory atlas 22 . This database catalogs the OCRs of 86 murine immune and hematopoietic cell types/states along with those of 4 non-hematopoietic cell types. The genomic dataset comprised 512,595 OCRs, each 180 bp in length. A consolidated set of 153 core TF motifs compiled from CIS-BP 6 and HOMER 27 , which are representative of structurally diverse families of mammalian TFs (Table S1), was used to analyze a total of 11,628 binary TF motif combinations. The CEseek scan encompassed 16 possible spacings (−5 to +10 bp) and four orientations for each TF pair, thereby generating 744,192 possible CE configurations ( Figure 2D ). We allowed for partial overlap of the individual motifs based on three considerations: (i) the base pairs flanking core TF motif residues have lower information content 28 ; (ii) cooperative binding can be associated with alterations in the DNA binding specificity of one or both of the TFs 16 ; and (iii) partially overlapping motifs can lead to anti-cooperative binding by the two TFs 29 , 30 . CEseek was used to assess the enrichment of each of these configurations in OCRs that exhibited preferential chromatin accessibility within specific cell types/states (top 5%) compared with all other states, employing stringent statistical criteria (p-value < 1e-20). This analysis identified 22,698 statistically significant binary motif configurations out of a possible 744,192 ( Figure 2D ). Further analysis of the enriched configurations revealed that most of the TF motif pairs displayed significant stereospecific constraints, typically limited to six or fewer permissible configurations out of the total 64 possibilities per pair ( Figure 2E ). In contrast, a smaller group of motif pairs exhibited greater configurational flexibility, spanning numerous possible spacings. Consequently, this yielded a bimodal distribution of TF motif pair counts, distinguishing motif pairs with stereo-specifically constrained configurations (digital) from those with extensive configurational flexibility (fuzzy) ( Figure S2A ). The left end of the distribution (motif pairs with six or fewer configurations; n=2,418 representing 5,063 configurations) was designated as cCEs and was retained for downstream analysis (Table S1). Conversely, the right end of the distribution (motif pairs with 49 or more configurations; n=170 representing 9,399 configurations) was characterized by TF motif pairs lacking strict stereospecific constraints. Notably, known CEs such as EICE, AICE, and NFAT-AP1 were recovered within the stereospecific subset of cCEs. As expected, cCEs exhibited longer sequence lengths compared to their simple constituent motifs ( Figures S2B and S2C ). To analyze whether cCEs predicted by CEseek were enriched in cognate TF-bound genomic regions, the PWMs of cCEs were systematically compared with de novo motifs derived from TF chromatin immunoprecipitation sequencing (ChIP-seq) datasets in CistromeDB 31 – 33 . Given the extensive curated TF ChIP-seq datasets for hematopoietic (K562) and immune (GM12878) cell lines, a comprehensive collection of uniformly processed TF ChIP-seq datasets was assembled, comprising 962 datasets encompassing 337 TFs in K562 cells and 370 datasets encompassing 160 TFs in GM12878 cells (Table S2). We reasoned that, as is the case for simple murine and human TF motifs, cCEs would be conserved across murine and human genomes. The systematic comparison revealed that 41.5% of the cCE PWMs displayed significant similarity (FDR < 0.05) to de novo motifs that were enriched within the TF ChIP-seq datasets ( Figure 2F ). An illustrative example of a CE identified by CEseek and its matched de novo motif enriched in a relevant ChIP-seq dataset is shown. These findings suggest that nearly half of the 5,063 cCEs identified by CEseek likely represent cognate TF binding sites in vivo . These numbers are underestimates given that not all TFs have a corresponding ChIP-seq dataset or have been assayed in the correct cellular context. Nevertheless, the comparative analysis confirmed our expectations that cCEs identified in murine regulatory DNA sequences are conserved within the human genome and likely to be functionally relevant. To assess whether the predicted cCEs enable cooperative DNA binding by their cognate TFs in vitro , we leveraged a comprehensive CAP-SELEX dataset 16 . This dataset comprises in vitro analyses of cooperative binding interactions for multiple pairs of TFs using synthetic DNA sequences and recombinant proteins. Cooperative binding was assessed through consecutive affinity purification coupled with the systematic evolution of ligands by exponential enrichment. The biochemical screen identified 315 distinct TF-TF pairs recognizing 618 CEs, each characterized by a specific orientation and spacing of the constituent motifs. We hypothesized that if the TF-TF cooperative interactions identified by CAP-SELEX were biologically meaningful, a significant portion of these CEs should be contained within the cCEs predicted by CEseek, representing conserved motifs in endogenous genomic regulatory sequences. Therefore, a comparative analysis was conducted between CAP-SELEX-derived CEs and CEseek-predicted cCEs based on their PWMs. Nearly half (44.0%) of the CAP-SELEX CEs exhibited significant similarity (FDR < 0.05) to their corresponding CEseek-predicted cCEs ( Figure 2G ). An illustrative example of a CAP-SELEX CE that is highly similar to a CEseek-predicted CE is shown. Remarkably, these matches, which are based on similarity scores, revealed strong alignments with CEseek-predicted CEs that displayed both gapped and overlapping TF motifs. We note that CAP-SELEX CEs were based on successive affinity selection of DNA sequences through in vitro binding assays, resulting in highly constrained, high-affinity DNA sequences. In contrast, CEseek-derived cCEs identified using the ImmGen dataset represent naturally variable motifs enriched in endogenous regulatory regions of immune and hematopoietic cells. The CAP-SELEX screen did not incorporate any biological constraints or cellular context. Despite these methodological differences, CEseek effectively recovered a substantial fraction of CAP-SELEX-delineated CEs. These CEs represent naturally occurring genomic regulatory sequences that are cooperatively bound by cognate TFs in vitro . Thus, a large fraction of the cCEs discovered by the CEseek scan qualify as CEs based on evidence of TF binding in vivo or in vitro . Functional characterization of cCEs using a TF-motif MPRA library While the integration of CEseek with ChIP-seq, and CAP-SELEX provides substantial support for an expansive set of CEs within regulatory DNA sequences, they do not involve direct measurement of their transcriptional activity in vivo or reveal their binary logic. To comprehensively test the large set of cCEs identified by CEseek, we designed a custom TF-motif Massively Parallel Reporter Assay (MPRA) library 21 that enabled quantitative testing of the transcriptional activities of the cCEs alongside their constituent simple TF motifs ( Figure 3A , Table S3). The TF-motif MPRA library included a total of 60,714 sequences for testing, each 30 base pairs in length, facilitating comprehensive and systematic analysis of simple and CE TF motif activities in vivo . The design of sequences in the TF-motif MPRA library involved several steps. Initially, the set of 5,063 cCEs ( Figures 2D and 2E ) was scanned across ImmGen ATAC-seq peaks, retaining the top 10 motif hits with the highest matching scores for each cCE. These motif hits were then extended to 30 bp. Genomic sequences that showed matches with more than one cCE, primarily represented by cCEs with partially overlapping motif pairs, were excluded, resulting in 3,795 unique cCEs to be tested. For each cCE, the four highest-scoring genomic sequences were selected for inclusion in the library. For each of these 30 bp sequences, three additional mutated variants were designed: one with the left motif mutated, one with the right motif mutated, and one with both motifs mutated at two consecutive positions with the highest information content (nucleotide A changed to C, T changed to G, and vice versa) ( Figure 3A ). This design enabled the quantification of the individual TF motif contributions to the overall transcriptional activity of each cCE. We note that the TF library additionally included the 153 core simple TF motifs that were used to discover the cCEs by CEseek. For such TF motif sequences to be assayed in the MPRA library, we first identified 20 unique genomic hits based on their highest motif scores and then extended these to 30 bp with selected randomized sequences that avoided inadvertent inclusion of alternative simple motifs. For each of these simple TF motif oligonucleotide sequences, we designed corresponding mutant pairs with substitutions at two consecutive base pairs located at positions with the highest information content, as was the case for cCE oligonucleotides. The TF-motif library also contained 1,000 randomly generated 30 bp sequences without recognizable TF motifs, which served as controls for background transcriptional activity. Each 30 bp oligonucleotide sequence was associated with an average of approximately 100 unique barcodes in the construction of the library for robust statistical testing of its transcriptional activity, ensuring reliable measurements of DNA sequence activity. Together, this custom TF-motif MPRA library contained diverse cCE and simple motif (n=3,795+153) sequences along with their mutant counterparts, providing a comprehensive, high-throughput platform for testing cCE activities and quantitatively assessing the contributions of the simple motifs (synergistic, additive or antagonistic) comprising them. Download figure Open in new tab Figure 3. Functional testing of cCEs using a custom TF-motif MPRA library (A) MPRA library design. The MPRA library comprises 30 bp oligonucleotides containing either simple TF motifs (n=153, with 20 instances of each) or CEs (n=3,795, with 4 instances of each). For simple motifs, each wildtype (WT) sequence is paired with a mutated version (M). For cCEs, each wildtype sequence is paired with three mutated variants: left motif mutated (M1), right motif mutated (M2), and both motifs mutated (M1, M2). Mutations were introduced at two consecutive positions with the highest information content in each TF binding site. Specifically, A was changed to C, T to G, and vice versa. The complexity of the library was n=60,714 sequences, and the activity of each sequence was assayable with >100 linked barcodes (see Methods for additional details). (B) MPRA analysis pipeline. The analysis of the TF-motif MPRA datasets is performed at three levels: barcode, sequence, and motif. In the barcode-level analysis, anomalous barcodes with excessively high read counts are removed. The remaining barcodes are then aggregated to their respective TF motif sequences, enabling quantification of sequence-level activity using DESeq2 to compare total barcode counts in RNA vs. DNA. At the motif level, simple motif activities are assessed by comparing the RNA/DNA log2 fold change (log2FC) of the wildtype motif-containing oligonucleotides with the log2FC of matched mutated control oligonucleotides. The activities of composite motifs are similarly estimated by comparing the log2FC values of the RNA/DNA ratios of sequences containing the wild type cCE with those of their double-mutant counterparts (see Methods ). (C) MPRA assays of simple and cCE motif activities in the indicated human cell lines. The MPRA library was transfected into K562 erythro-myeloid cells (6 replicates), GM12878 B cells (5 replicates), and Jurkat T cells (5 replicates) (see Methods ). Box plots showing the distributions of the transcriptional activities of simple motifs and cCEs in the three cell lines. Activating and repressive motifs are distinguished by red and blue color codes, respectively. (D) Counts of activating motifs. Bar plots depict the number of activating simple motifs and cCEs in K562, GM12878, and Jurkat cells. The shaded bars highlight activating cCEs that are composed of simple motifs that lack detectable transcriptional activity on their own in the indicated cellular context. (E) Shared and cell type-specific activating cCEs. Pie chart illustrating the numbers of shared and cell-type specific activating cCEs among the K562, GM12878, and Jurkat cell lines. (F–H) Examples of activating cCEs. Sequence logos (PWMs) represent cCEs, while bar plots show the transcriptional activities of the paired wild-type oligonucleotides and their mutated variants. The dashed lines relate activity of each wildtype oligonucleotide with its mutant counterparts. (I–K) Synergistic and antagonistic interactions within TF motifs comprising cCEs. Volcano plots highlighting cCEs that exhibit significant synergy or antagonism between their two constituent simple TF motifs in the three cell lines. Download figure Open in new tab Figure S3. Functional testing of cCEs using a custom TF motif MPRA library (A) Barcode count distribution. Histograms showing the distribution of unique barcode counts recovered from RNA for TF motif library sequences after transfection in K562 (left), GM12878 (middle), and Jurkat cells (right). Sequences with ≥30 unique barcodes (dashed line) were used for downstream analyses. (B) Reproducibility of the biological replicates. Heatmap displaying pairwise Pearson correlation coefficients for aggregated barcode counts linked with a given TF motif sequence across RNA and DNA biological replicates. (C) Schematic illustrating the determination of simple and composite motif activities from the MPRA datasets. For simple motifs, motif activity was calculated using DEseq2 by comparing the RNA/DNA log2 fold change (log2FC) of motif-containing sequences (sequence activity) as test values and the RNA/DNA log2 fold change (log2FC) of matched mutated control sequences as reference values. For composite motifs, activities were similarly assessed by using the RNA/DNA log2FC of CE motif-containing sequences with those of matched double-mutant controls. To determine individual motif contributions within a CE, the RNA/DNA log2FC values of sequences carrying single-motif mutations were analyzed, with those from matched double-mutant sequences used as controls. (D) Numbers of functional simple and composite TF motifs revealed by MPRA assays. Table summarizing the counts of functional (activating or repressing) simple and composite motifs identified in K562, GM12878, and Jurkat cells, along with their median sequence and motif activities. (E) Comparisons of cCE activities across cell types. Scatter plots of cCE motif activity (log2FC) measured in GM12878 vs. K562 (left), GM12878 vs. Jurkat (middle), and Jurkat vs. K562 (right) cells. The green and blue dots indicate cell type-specific activating or repressing cCEs in the indicated cell types. The red dots are cCEs that are activating or repressing in both cell types. The gray dots are cCEs with no detectable activity in the indicated cell types. Representative simple TF motifs that are constituents of cell-type specific CEs are shown in the left panel. The TF-motif MPRA library was introduced into three human cell lines: K562 (erythro-myeloid), GM12878 (B cell), and Jurkat (T cell). These cell lines were selected for the analyses of the activities of CEs, as K562 and GM12878 have served as ENCODE Tier 1 models and have been extensively profiled using diverse genomic methodologies 23 , 31 – 33 . Their associated datasets enabled downstream analyses of CE containing genomic sequences with chromatin and transcriptional features in the context of endogenous loci. The inclusion of Jurkat T cells enabled an exploration of CEs that were differentially active in human B versus T-cell states. After transfection of the plasmid TF-motif MPRA library into each of the three cell lines, the GFP reporter transcripts were isolated with biotin-labeled probes, converted into cDNA and sequenced to quantify the distribution of the motif-associated barcodes ( Figure S3A ) (see Methods ). Analyses of the MPRA datasets, involving at least 5 biological replicates for each cell line, were performed at three levels: barcode, sequence, and motif ( Figure 3B ). Barcode-level analysis was initially performed via DEseq2 34 to identify and remove barcodes that were statistical outliers characterized by excessively high counts that were not consistent across replicate experiments. These anomalous barcodes were detected by applying Cook’s distance, an established diagnostic metric for assessing the influence of data points. Barcodes identified as statistical outliers at the 99th percentile threshold were excluded from subsequent analysis, ensuring the robustness and accuracy of downstream oligo-level assessments. After such filtering, barcode counts were aggregated to their corresponding TF motif sequences (CE or simple motif), retaining only sequences with unique barcodes equal to or greater than 30 for subsequent analysis ( Figure S3A ). This enabled the activities of the 30 bp oligonucleotide sequences to be quantified by comparing total barcode counts in RNA with corresponding counts in DNA using DESeq2 differential analysis ( Figures 3B and S3B ). The activities of simple and composite TF motifs were then specifically analyzed as follows. For simple motifs, activity was determined by comparing the RNA/DNA log2-fold change (log2FC) of sequences containing the parent motif (sequence activity) with that of matched mutant control sequences (motif activity). Similarly, for composite motifs, motif activities were estimated by comparing the RNA/DNA ratios of parent sequences containing the cCE to those of matched double-mutant control sequences ( Figure S3C ). Simple and composite motifs were classified as activating or repressing only when both sequence-level and motif-level activities reached statistical significance and exhibited consistent directionality (above or below that of the minimal promoter). Overall, the TF-motif MPRA library, with its analytical framework, integrated barcode-level, sequence-level, and motif-level data, enabling comprehensive and precise quantification of TF motif activities and their synergistic, additive, or antagonistic interactions. Systematic analysis of TF motif activities with the TF-motif MPRA library across the three cell lines revealed 893 unique activating CEs: 775 in K562 cells, 321 in GM12878 cells, and 91 in Jurkat cells ( Figure S3D , Table S4). In contrast, substantially fewer activating simple motifs were detected, with 35 in K562 cells, 19 in GM12878 cells, and 5 in Jurkat cells. As anticipated, both the average sequence-level and motif-level activities of CEs were substantially higher than those of their simple motif counterparts ( Figures 3C and S3D ). Further analysis of the constituent simple motifs within these activating CEs revealed that many were composed of simple motifs that lacked detectable transcriptional activity on their own (48.4% in K562, 53.3% in GM12878, and 80.2% in Jurkat; Figure 3D ). This observation suggested that simple motifs were functioning synergistically to activate transcription when configured as CEs. Interestingly, despite the use of a minimal promoter within the MPRA library design that displayed low basal activity, the analysis also uncovered repressing CEs: specifically, 56 in K562 cells, 7 in GM12878 cells, and 1 in Jurkat cells. Next, the CEs were examined for contextual activities that were dictated by the cell types in which they were assayed. A considerable proportion of activating CEs (75.9%) demonstrated cell type-specific activities, most notably in K562 and GM12878 cells, with the remaining manifesting activity in at least two of the three cell types in which they were analyzed ( Figures 3E and S3E , Table S5). The reduced sensitivity of the MPRA assay in Jurkat T cells likely impacted the discovery of cell type-specific CEs in this context. The GM12878 B-cell type-specific cCEs included motifs for EBF, E2A, PAX5, NF-κB, Ikaros and ARID3A, whereas the K562 erythro-myeloid specific CEs encompassed motifs for GATA, RUNX, SPI1, NFE2 and PPAR. These motifs are reflective of TFs that have been shown to regulate the distinctive gene expression programs of these cell types 35 . Although the above analysis demonstrated that many CE sequences discovered by CEseek were more highly active than their counterparts bearing constituent simple TF motifs, it did not involve quantitative testing of their synergistic activities. To demonstrate such testing, we used several previously noted, well-characterized CEs, including their constituent simple TF motifs, along with mutant versions of each parent sequence, all of which were included in the design of the TF-motif MPRA library. For example, in GM12878 cells, the ETS::IRF composite motif (EICE) 11 , which is composed of ETS (SPI1) and IRF binding sites, displayed substantial activity ( Figure 3F ). Assessment of the transcriptional activities of the EICE, with its simple motifs, and their corresponding control sequences revealed that mutation of either the SPI1 or the IRF binding site resulted in activity reductions to levels observed with the isolated simple motifs. Importantly, the additive activity of the individual ETS and IRF motifs was lower than the activity of the intact CE. This was consistent with the expectation that the synergistic action of these motifs was dictated by the cooperative binding of the cognate transcription factors 11 . Similarly, the previously characterized NFAT::AP-1 composite motif 13 , 14 , comprising NFAT and AP-1 motifs separated by two base pairs, also demonstrated significant synergistic activity ( Figure 3G ). This CE showed robust transcriptional activity (motif activity log2FC of 1.02), whereas the isolated NFAT motif exhibited no detectable activity, and the AP-1 partner motif displayed low activity. An example of a novel RUNX::CREB CE that displayed a similar strong synergistic effect is shown in Figure 3H . In K562 cells, neither RUNX nor CREB simple constituent motifs exhibited significant activity. However, their CE demonstrated a remarkable increase in transcriptional activity (motif activity log2FC of 2.49). Based on the above exemplars, we systematically evaluated CEs for synergistic activities in one or more of the cell types. For these analyses, the individual motif contributions within CEs were analyzed by comparing the RNA/DNA ratios of oligos containing singly mutated motifs to those in their matched double-mutant counterparts ( Figures 3B and S3C ). The synergistic activity of a CE was then quantified by contrasting the motif activity of the intact CE with the additive contributions (log scale) of the component motifs. Activating CEs exhibited mean log2-fold changes of 0.33 in K562, 0.16 in GM12878, and 0.39 in Jurkat cells. Using a p-value cutoff of 0.05, we identified 215 non-additively functioning CEs in the three cell lines. ( Figures 3I-K , Table S5). Among these, 190 CEs manifested synergistic transcriptional activation promoted by their constituent TF motifs. The remaining 25 CEs exhibited varying degrees of transcriptional attenuation or antagonism between the two motifs. Thus, these CEs satisfied both structural (stereospecific juxtaposition of TF motifs) and functional (cooperative transcriptional activation or repression) criteria. The quantitative analysis of CE transcriptional activities along with that of their simple motifs enabled illustration of the underlying combinatorial regulatory logic. Specifically, we examined the activity patterns of CEs to determine their correspondence with the following canonical two-input logic gates: AND, AND-NOT, and XOR (Table S6). A CE was classified as an AND gate if it displayed significant positive synergy between its two constituent motifs ( Figures 3F-K ). In all, 190 AND gates were identified (169 in K562, 36 in GM12878, and 21 in Jurkat). A CE was assigned as an AND-NOT gate if one motif was active, but its pairing with an inactive or weakly active motif resulted in loss of CE activity with significant negative synergy. In total, 79 CEs were observed to demonstrate AND-NOT regulatory gate characteristics (72 in K562, 9 in GM12878, and 4 in Jurkat). Examples of CEs demonstrating this regulatory logic included the HSF1[+]:5:ZEB1[+] CE, which exhibited reduced combined activity (CE motif activity log2FC = −0.175; FDR = 0.563) despite strong activity of the HSF motif (log2FC = 0.892; FDR = 1.49E-05) and negligible activity of the ZEB1 motif (log2FC = 0.072; FDR = 0.941). This combination resulted in significant negative synergy (log2FC = −1.14; p-value = 6.28E-06). Another example was the MYC[-]:5:REST[+] CE, which showed modest combined activity (CE motif activity log2FC = 0.163; FDR = 0.754) despite strong activation by the MYC motif (log2FC = 1.112; FDR = 0.00144) and minimal activity contribution from the REST motif (log2FC = −0.067; FDR = 0.969). This combination similarly exhibited significant negative synergy (log2FC = −0.882; p value = 0.0267). Finally, an XOR gate was assigned to CEs in which both simple motifs were individually active, yet the CE displayed diminished activity, also resulting in negative synergy. A single XOR gate-classified CE was observed to function in K562 cells (Table S6). In the XOR-classified CE, PPARD[-]:-3:NRF1[+], the individual simple motifs showed high activity (PPARD motif log2FC = 1.27; NRF1 motif log2FC = 1.53). However, the CE displayed diminished overall activity (log2FC = 0.87), which did not achieve statistical significance (FDR = 0.084), thereby reflecting strong negative synergy (log2FC = −1.93, synergy p-value = 0.000482). Thus, the activities of most CEs corresponded to AND or AND-NOT regulatory logic gates. CEs are enriched in ChIP-seq datasets and associated with co-binding of cognate TFs in vivo To determine whether the MPRA-delineated CEs are specifically associated with co-binding of pairs of cognate TFs in vivo , we utilized the following three-step approach. First, TF motifs within CEs were matched to relevant TF ChIP-seq datasets. A motif was considered matched to a TF ChIP-seq dataset if the TF analyzed by ChIP-seq belonged to the same TF family as the one identified by the simple motif or if the simple motif matched one of the top 10 motifs enriched within a particular TF ChIP-seq dataset. In a second step after generating correspondence between CEs and relevant ChIP-seq datasets, all CEs demonstrating synergistic transcriptional activity (either activating or repressing) were evaluated for their enrichment within peaks from the corresponding TF ChIP-seq datasets. This comprehensive analysis identified 154 CEs significantly enriched within TF-bound regions in K562 and GM12878 cells (Table S7). Finally, those CEs were tested for their statistical enrichment specifically within the co-bound regions of the corresponding pair of TFs by performing CEseek, using the singly bound regions as the background sets. The overall testing framework is illustrated with the SPI1::JUN CE which in K562 cells demonstrated a pronounced synergistic interaction (log2FC of 1.55). ChIP-seq analysis revealed substantial overlap between the SPI1- and JUN-binding regions, with 22.2% of the SPI1-binding regions overlapping with JUN and 30.5% of the JUN-binding regions overlapping with SPI1 ( Figure 4A ). Crucially, CEseek-based enrichment analysis of binary SP1 and JUN configurations within SPI1 and JUN co-bound genomic regions versus their singly bound counterparts revealed that the identified CE configuration was the most significantly enriched among all potential 64 configurations that were evaluated. We note that the JUN[-]:-3:SPI1[-] and JUN[+]:-2:SPI1[-] configurations were effectively identical because of the palindromic nature of the JUN motif. We note that certain CEs displayed significant synergistic interactions despite extended spacing of the simple TF motifs. For example, the NFY::KLF CE in K562 cells exhibited robust synergistic interactions with a 10 bp spacing (synergistic interaction log2FC of 0.65; Figure 4B ). CEseek analysis confirmed this CE as one of the most enriched configurations in the NFYA and KLF13 co-bound regions. Notably, NF-Y and SP/KLF family TFs have been shown to co-bind genomic regions in a significant proportion of human promoters 36 , 37 . Download figure Open in new tab Figure 4. CEs are enriched in TF-bound genomic regions and evidence co-binding of cognate TFs (A–C) Representative CEs demonstrating enrichment in co-bound genomic regions of cognate TFs. The top panels display sequence logos of the CEs (PWMs), whereas the bar plots show the activities of the wild-type oligonucleotides and their mutated variants (connected by dashed lines) in the indicated cells. For each CE, the simple TF motif was used to query CistromeDB for a cognate TF ChIP-seq dataset. Venn diagrams showing the numbers of co-bound versus singly bound regions for the relevant pairs of TFs. Heatmaps depict the enrichment of CE configurations (n=64) in co-bound TF ChIP-seq peak regions relative to singly bound peak regions. The CE configurations assayed in the MPRA library are highlighted with green boxes in the heatmaps and correspond to the PWMs displayed in the top panels. (D, E) Diversity of CEs comprising AP-1 or KLF motifs. Sequence logos of the AP-1 ( D ) or KLF ( E ) families of CEs are shown on the left. The bar plots show the measured sequence, motif, and synergistic activities of each CE in K562, GM12878, and Jurkat cells. These CEs were enriched in corresponding TF ChIP-seq datasets (see panels A-C and Figure S4 ). Download figure Open in new tab Figure S4. CEs are enriched in cognate TF-bound genomic regions (A) Prevalent simple TF motifs contained within CEs. Bar plot showing the counts of simple TF motifs within CEs that are enriched in cognate ChIP-seq datasets (CistromeDB). (B) Diversity of CEs. Grid displaying the pairs of simple TF motifs that comprise experimentally demonstrable CEs. Rows are the top 15 simple TF motifs with the highest CE counts, displayed in panel A . Columns are the more diverse simple TF motifs that the prevalent set (panel A) is paired with. In addition to the large set of synergistically activating CEs that were statistically enriched in co-bound genomic regions of cognate TFs, we also uncovered instances of antagonistic TF motif interactions within novel CEs similarly enriched in co-bound regions. For example, in GM12878 cells, the MYC::REL CE exhibited an antagonistic interaction (log2FC of −0.58; CE activity increased upon mutation of the MYC motif; Figure 4C ). CEseek analysis demonstrated that this CE was the most enriched configuration within the MYC and RELA co-bound regions. We note that co-bound genomic regions in the ChIP-seq datasets could reflect either cooperative or anti-cooperative binding by the cognate TF pairs. We next examined the diversity of the discovered CEs by analyzing the occurrence of simple motifs. The AP-1, KLF and ETS motifs were the most frequently observed simple motifs within the 154 activating CEs ( Figure S4A ). By combining with distinct partner motifs, these prevalent simple motifs greatly diversified the CE repertoire ( Figure S4B ). AP-1 family TFs modulate gene expression in response to diverse stimuli, such as growth factors, hormones, and cytokines, and play crucial roles in regulating key cellular processes, including apoptosis, proliferation, and differentiation 38 . Analysis of CEs that contained the AP-1 motif revealed a wide range of absolute and relative activities that were specific to each cell type ( Figure 4D ). AP-1 CEs involving ZEB, KLF, and NFAT motifs displayed distinct cell type-specific activity patterns. Notably, the SPI1::AP1 CE exhibited the highest activity across all the tested cell lines. In contrast, AP-1 CEs involving TBP and SMAD motifs manifested repressive effects in K562 cells. Importantly, all AP-1 CEs were significantly enriched at genomic regions where AP-1 ChIP-seq peaks intersected with those that were enriched for the partner motif. Collectively, these findings underscore that the AP-1 motif containing CEs can enable contextually specific activating or repressive inputs in regulating transcription. The Krüppel-like transcription factor (KLF) family comprises TFs characterized by three carboxyl-terminal (C-terminal) C2H2-type zinc finger domains 39 . These domains facilitate binding to GC-rich DNA regions. KLFs also regulate diverse cellular processes, including proliferation, differentiation, and apoptosis, as well as tissue development and homeostasis. Analysis of the KLF motif-containing CEs also revealed distinct transcriptional activities that varied significantly depending on the cell type ( Figure 4E ). These CEs exhibited a broad spectrum of both sequence and motif transcriptional activities. Among the identified CEs, NFY::KLF and SPDEF::KLF displayed consistently high activity across all three analyzed cell lines, whereas the KLF::AP-1 CE uniquely demonstrated activity in K562 cells. Interestingly, the KLF:-3:RXRA CE maintained stable activity across cell types; however, altering the spacing and orientation of the same simple motifs in the KLF:-4:RXRA CE led to selective loss of its activity in Jurkat cells. In contrast, the KLF::RBPJ CE displayed repressive activity in all the examined cell types, with the RBPJ motif antagonizing the activity of the KLF motif. All KLF motif-containing CEs were notably enriched at genomic regions where the KLF ChIP-seq peaks overlapped with those from other relevant TF ChIP-seq datasets. Thus, the discovery, testing and validation of a diverse set of AP-1- and KLF-motif-containing CEs suggest the enormous versatility of members of these TF families in controlling diverse patterns of gene activity. Using GRACE to predict TF motif and CE activities across cell types To decode the TF motif lexicon embedded in short DNA sequences and understand how composite elements and their constituent simple motifs control transcription across diverse cellular contexts, we developed a deep learning model, GRACE ( G enome R egulatory A nalysis of simple and C omposite transcriptional E lements). The architecture of the model was designed to be trained on the TF-motif MPRA datasets and to predict context-specific transcriptional activities of short (50 bp) DNA sequences at single-nucleotide resolution. The GRACE model integrates convolutional layers to learn TF binding motifs and employs self-attention mechanisms with positional encoding to capture stereospecific interactions between these motifs 40 , 41 ( Figure 5A ). Download figure Open in new tab Figure 5. Using GRACE to predict CE and TF motif activities (A) Model overview. GRACE is a neural network model that integrates convolutional layers to learn TF binding motifs and employs self-attention with positional encoding to capture interactions among these motifs. See Methods for details. (B) Prediction accuracy of GRACE. Scatterplot comparing GRACE-predicted log2FC transcriptional activities for the withheld set of sequences (n=3,900) from the training set (n=49,042) with their experimentally measured MPRA activities in GM12878 cells. The Pearson correlation coefficient is indicated. (C) Correlations of the GRACE predictions with the MPRA replicates. Heatmap showing pairwise correlations between GRACE predictions and MPRA measurements across individual biological replicates, as well as with the aggregated set of MPRAs in GM12878 cells. (D) Generating transcriptional activity contribution scores with GRACE. Illustration of GRACE assigned activity contribution scores to each base pair using in silico saturation mutagenesis. The indicated 30 bp oligonucleotide, corresponding to the displayed CE position weight matrix (PWM), was systematically mutated to evaluate the impact of individual base pairs on its transcriptional activity using GRACE model trained with the GM12878 MPRA dataset (heatmap). The resulting contribution weight matrix (CWM) is displayed below the heatmap. (E) Histogram showing the false discovery rate (FDR) of the similarity between CEseek PWMs and GRACE CWMs for activating cCEs. All the activating cCE sequences (30 bp) in the MPRA library were used to generate CWMs via the GRACE model of GM12878 cells, which were then compared to their corresponding PWMs. The histogram shows the distribution of median similarity false discovery rate (FDR) values between the PWMs and CWMs, with the fraction of matching cCEs (FDR < 0.05) indicated in the panel. (F) Examples of matched CEseek PWMs with GRACE predicted CWMs. The FDR values for each pair are indicated. Download figure Open in new tab Figure S5. Using GRACE to predict CE and TF motif activities (A, B) Prediction accuracy of GRACE models trained on the K562 (A) and Jurkat (B) MPRA datasets, as shown in Figure 5B . Scatter plots comparing GRACE-predicted transcriptional activities for withheld test sequences to experimentally measured MPRA activities (log2FC). Pearson correlations are indicated. (C, D) Correlations of GRACE predictions with MPRA replicates. Heatmap showing pairwise correlations between GRACE predictions and MPRA measurements across individual biological replicates, as well as with the aggregated set of MPRAs in K562 (C) and Jurkat cells (D) as in Figure 5C . (E, F) Histograms showing the false discovery rate (FDR) of the similarity between CEseek PWMs and GRACE CWMs for activating cCEs in K562 (E) and Jurkat cells (F), as shown in Figure 5E . The fraction of matching cCEs (FDR < 0.05) is indicated in each panel. (G) GRACE activity predictions of different cCE configurations. Heatmap of GRACE-predicted activity (log2FC) for different NFY::KLF CE configurations discovered by CEseek (see Figure 4B ). The median activity for each configuration is shown. Two notable configurations (boxed) represent one tested in the TF motif MPRA library and the other learned by GRACE. (H) Example of the GRACE learned CE configuration. CEseek PWM and GRACE CWM for an NFY::KLF configuration are predicted to exhibit strong transcriptional activity and cooperative binding (see Figures 4B and S5G ). GRACE predictions of the activities of simple motifs comprising the CE along with their PWMs are shown to the left and right of the CE. To enhance the reliability and generalizability of the GRACE model predictions, several considerations were incorporated during model training, informed by the design of the TF-motif MPRA library. Specifically, to mitigate the central positioning bias inherent in the 30 bp oligonucleotide design, we implemented a sequence augmentation strategy. This involved extension of the MPRA assayed oligos with random 20 bp sequences, randomly split into segments with unequal length and appended to the ends of the assayed 30-mers for each GRACE training batch. This strategy effectively minimized positional bias and improved model generalizability. Notably, GRACE’s attention mechanism was restricted to 30 bp blocks, allowing the model to focus on local regulatory features while maintaining computational efficiency. As a result, GRACE can flexibly accommodate input sequences ranging from 30 to 50 bp or longer. This flexibility was enabled by its transformer-based architecture and localized attention mechanism, which allowed the model to learn sequence features in fixed-size windows while maintaining the capacity to process variable-length inputs without loss of resolution or context. Furthermore, to ensure that the GRACE model could accurately predict the activities of variant sequences, we partitioned the dataset into training, validation, and test sets. We assigned variant sequences that were derived from a given parent sequence to the same subset, preventing data leakage, thereby enhancing the model’s ability to generalize predictions to all possible single nucleotide variants. The GRACE models were trained in a context-specific manner with the MPRA datasets from the three cell lines. The predictive performance of the GRACE models was evaluated using held-out test sequences, demonstrating a strong overall Pearson correlation of 0.802 between the predicted and experimentally measured activities in GM12878 cells (0.813 in K562 cells and 0.766 in Jurkat cells; Figures 5B , S5A , and S5B ). This predictive accuracy was comparable to correlations observed among biological replicates of the TF-motif MPRAs. Specifically, Pearson correlation coefficients comparing GRACE predictions to biological replicate measurements for withheld test sequences, versus correlations among biological replicates, were 0.730 compared to 0.707 in K562 cells, 0.642 compared to 0.661 in GM12878 cells, and 0.644 compared to 0.634 in Jurkat cells ( Figures 5C , S5C , and S5D ). To elucidate the underlying motifs contributing to the transcriptional activity of any given DNA sequence, we performed in silico saturation mutagenesis and used GRACE to predict the corresponding transcriptional activities ( Figure 5D ). For each 30-bp oligonucleotide, single-base pair variants were generated, and their transcriptional activities (log2FC) were predicted and displayed as 4 × L matrices (rows = A, C, G, T; columns = positions). To generate contribution weight matrices (CWMs), the mean activity values (log2FC) of all four nucleotides for a given position were subtracted from the activity of a particular base pair at that position, thereby removing position-specific biases and yielding normalized CWMs. We reasoned that if the CEs discovered by CEseek reflected transcriptionally activating sequences at single-nucleotide resolution, then their position weight matrices (PWMs) would closely resemble the corresponding contribution weight matrices (CWMs) learned by GRACE. An example of such a comparison is shown for the SPI1[+]:-2:JUN[-] CE in Figure 5D . In support of the overall hypothesis, very large proportions of active CEs exhibited strong similarity of their PWMs and CWMs, achieving concordance rates of 83.5% in GM12878, 80.3% in K562, and 95.6% in Jurkat cell lines ( Figures 5E , 5F , S5E , and S5F , Table S8). These results demonstrated that the observed transcriptional activities in the CE-containing MPRA sequences are indeed dictated by the motifs identified by CEseek. Furthermore, they highlighted the robustness of motif detection by GRACE despite the limited sequence space that was assayed in the TF-motif MPRA library. They establish that the GRACE-learned CWMs correspond to the CEseek-generated PWMs and illustrate the power of GRACE in making single-nucleotide resolution activity predictions despite being trained on a small set (<50,000) of all possible variants (4 30 ) of a given 30-mer sequence. Next, GRACE was examined for its ability to predict the impact of spacing and orientation on CE transcriptional activities. A focused analysis of NFYA and KLF motifs in K562 cells predicted active NFY::KLF CE configurations beyond those that were designed for and assayed in the MPRA library. Notably, GRACE predicted substantial activities for several untested functional configurations (spacer lengths from 5 to 9 base pairs and different relative orientations of NFY and KLF motifs), including the NFY[-]:8:KLF[+] CE, which had the highest predicted activity ( Figure S5G ). Further analysis of the NFY[-]:8:KLF[+] CE revealed that although the NFYA and KLF motifs displayed detectable independent activities, their combination within the CE manifested a strong synergistic interaction ( Figure S5H ). Importantly, these findings were corroborated by CEseek’s configuration enrichment analysis in open chromatin regions ( Figure 4B ), illustrating that of GRACE can also predict novel stereospecific regulatory interactions within CEs. We note that although variably spaced CE configurations were not included in the MPRA training dataset, GRACE likely captured relevant sequence patterns through various simple and CE motif mutations that were included in the design of the TF-motif MPRA library. These mutations could have generated sequences resembling naturally occurring variably spaced motifs, thereby enabling GRACE to infer and learn patterns beyond those of explicitly provided configurations. These results exemplify the ability of GRACE to learn the CE motif lexicon from measured transcriptional activities of a limited set of short but sufficiently diverse DNA sequences. GRACE learns genome-wide transcription and chromatin features To investigate the relationship between the transcriptional activities of short DNA sequences (30-mers) learned by GRACE and endogenous chromatin features, we applied the GRACE model to predict transcriptional activities across OCRs in GM12878 cells. The transcriptional activity of each OCR was quantified by summing the GRACE-predicted activities (log2FC) of the constituent 30-mer sequences, which exceeded an activity threshold (log2FC=0.2; see Methods ). Notably, OCRs ordered based on their GRACE-predicted transcriptional activities displayed a similar ordering in their levels of the activating histone modification H3K27ac, which is strongly associated with active transcription ( Figures 6A and 6B ). Thus, the GRACE predictions of the transcriptional activities of short DNA sequences generated by training on the MPRA datasets were in accordance with activating histone modifications at endogenous chromatin loci. Download figure Open in new tab Figure 6. Predicting transcriptional activity of open chromatin regions using GRACE (A) Ranking of the OCRs according to their GRACE predicted transcriptional activities. Scatter plot displaying OCRs in GM12878 cells ranked by their GRACE-predicted transcriptional activities (see Methods ). (B) Box plots illustrating the distributions of H3K27ac ChIP-seq signals (reads per million, RPMs) across the GRACE-ranked OCRs in GM12878 cells, as ordered in panel (A) . See Table S9 for the H3K27Ac dataset. (C, D) GRACE predicted transcriptional activities in the promoter and enhancer regions of GM12878 cells. Promoters were defined by having an annotated TSS and an activity measured by GRO-cap in GM12878 cells of >1 TPM (see Table S9 for the dataset). Promoter regions were aligned based on annotated TSSs ( C ). Enhancers lacked an annotated TSS within a 1000 bp and were aligned based on the center of one of the GRACE-predicted transcriptionally active sequences (seqlets) ( D ) (see Methods ). The aligned promoter or enhancer sequences were clustered using K-means, and the heatmaps display the median GRACE-predicted transcriptional activities. The ridge plots above and below each heatmap represent the median GRO-cap signals generated from the promoter sequences in each cluster on the sense/forward and antisense/reverse strands, respectively. The numbers (#) of aligned promoter or enhancer sequences in each K-means cluster are indicated on the right of each panel. (E, F) Analysis of TF motif contributions to chromatin accessibility and transcriptional activity using ChromBPNet and GRACE models. Scatter plots depict the contributions of simple (gray) and CE (red) motifs to chromatin accessibility (x-axis, derived using ChromBPNet) and transcriptional activity (y-axis, derived using GRACE) in K562 ( E ) and GM12878 ( F ) cells. Motif contributions to chromatin accessibility or transcriptional activity are defined as fold changes between contribution scores of motif sequences (ChromBPNet or GRACE) compared with the contribution scores of similarly modeled 50 bp flanking regions (see Methods ). Representative simple TF motifs that are contained either within highly active CEs or within CEs that differentially contribute to chromatin accessibility and transcriptional activity in indicated cell contexts are displayed on plots. (G–I) Examples of activating and repressing CEs with their motif contributions to chromatin accessibility and transcriptional activity. The left and right columns in each panel illustrate the simple TF motifs that comprise the corresponding CEs in the middle column. The first row in each panel displays sequence logos (PWMs) based on CEseek. The second and third rows provide the aggregated contribution weight matrices (CWMs) for chromatin accessibility (second row) and transcriptional activity (third row) using ChromBPNet and GRACE models, respectively. The dashed gray lines indicate baseline contributions estimated for 50 bp flanking regions (see Methods ). Download figure Open in new tab Figure S6. Predicting transcriptional activity of open chromatin regions using GRACE (A, B) GRACE predicted transcriptional activities in the promoter and enhancer regions of GM12878 cells. These panels show the remaining promoter and enhancer K-means clusters from the analysis detailed in Figure 6C , 6D . (C) Box plots illustrating the distribution of ATAC-seq signals (RPMs) across the GRACE-ranked open chromatin regions displayed in Figure 6A . (D) Concordance of GRACE and ChromBPNet predictions at TSSs. Plots show the mean ChromBPNet and GRACE contribution scores ±500 bp of transcription start sites (TSSs) of active promoters in GM12878 cells. Profiles are compared by normalizing to maximum contribution scores for transcriptional activity or chromatin accessibility. (E–G) Examples of activating simple motifs or CEs with their chromatin accessibility and transcriptional activity contributions (See Figures 6G-I for details). To assess whether GRACE predictions align with transcription at endogenous loci, we took advantage of nascent RNA data from Global Run-On sequencing (GRO-cap). This method enables quantitative analysis of transcription initiated at promoters and within active enhancers. GRO-cap analysis has demonstrated bidirectional transcription at both active TSSs and enhancers 42 . Strikingly, OCR clusters identified through K-means clustering based on their GRACE-predicted transcriptional activities demonstrated strong concordance with GRO-cap signals at promoters and enhancers ( Figures 6C and 6D , Table S9). At promoters, GRACE predictions around TSSs revealed distinct spatial patterns between predicted active regions and transcription initiation sites. Specifically, on the sense strand, the intensity of the GRO-cap signals at the TSSs was correlated with the GRACE-predicted transcriptional activity ( Figures 6C and S6A ). Intriguingly, on the antisense strand, GRO-cap signals precisely mapped downstream of the GRACE-predicted active sequences. The gaps between the sense and antisense GRO-cap signals ranged from 50 to 300 bp. Strikingly, at enhancers, GRACE-predicted transcriptionally activating sequences could be aligned with bidirectional GRO-cap signals that manifested stereotypic patterns based on the number and relative positioning of their GRACE active regions. For enhancers containing a single predicted transcriptional activating region, GRO-cap signals on forward and reverse strands precisely flanked the predicted sequence, with GRO-cap signal intensities correlating with GRACE-predicted activities ( Figures 6D and S6B ). For enhancers with two predicted activating sequence regions, GRO-cap signals were positioned in opposite directions relative to the paired elements, likely reflecting their coordinated transcriptional activating functions ( Figure 6C ). Notably, increased spacing between the two activating sequence regions predicted by GRACE correspondingly resulted in greater distancing of the bidirectional transcription initiation sites. Thus, context-specific GRACE models, after training on TF-motif MPRA datasets, can be used to perform genome-wide scans and accurately predict bidirectional nascent transcription at endogenous promoters and enhancers at base pair resolution. Given that GRACE can predict DNA sequences that contribute to nascent transcription at endogenous enhancers, a feature that is recognized as a strong indicator of enhancer activity 42 , we analyzed these regions for the occurrence of CEs. Remarkably, 85.7% of the CEs (281 of 328) delineated by the TF-motif MPRAs in GM12878 cells were enriched in the GRACE predicted transcriptionally active enhancer regions (FDR<0.01, Table S9). Notably, these included known CEs such as EICE and a large set of newly discovered CEs, including stereospecific pairings of the SPI1::VDR, EBF (half-site)::MYC, NFkB::MYC, NFkB::HSF and E2F::NFAT motifs. Strikingly, OCRs ordered based on their GRACE-predicted transcriptional activities also displayed a similar ordering in terms of their levels of chromatin accessibility in ATAC-seq data, suggesting a relationship between GRACE-learned transcriptional features and chromatin accessibility ( Figure S6C ). To explore this relationship at single-nucleotide resolution, we next compared the expanded TF motif lexicon learned by GRACE with that derived from an orthogonal deep learning model, ChromBPNet 20 , which is trained on chromatin accessibility features using ATAC-seq datasets. ChromBPNet is a convolutional neural network model that employs dilated convolutions with residual connections to effectively model and remove the Tn5 enzymatic bias inherent in ATAC-seq profiles and then predicts single-nucleotide resolution contributions to chromatin accessibility. As a first step, we compared the GRACE predictions with those of ChromBPNet at active TSSs. Remarkably, the GRACE model exhibited substantial concordance with the ChromBPNet predictions ( Figure S6D ). The aggregated ChromBPNet chromatin accessibility contribution scores mirrored the GRACE-predicted transcriptional activity profiles, with the latter declining rapidly to minimal levels precisely at the TSS. This finding suggested that both models identify similar underlying DNA sequence motifs despite being trained on orthogonal datasets. Next, a systematic comparative analysis was performed to evaluate the contributions of simple and composite motifs to transcriptional activity (using GRACE) and to chromatin accessibility (using ChromBPNet) ( Figures 6E and 6F , Table S10). Here, the simple motif collection consisted of 879 motifs from the JASPAR 2024 vertebrate core collection 43 , representing a rigorously curated, experimentally validated dataset widely utilized in genomics research to annotate TF binding sites. The CE motif collection included 951 unique functional CEs identified by our MPRA experiments in the K562, GM12878, and Jurkat cell lines ( Figure 3E ). Motif contributions in the GRACE and ChromBPNet models were quantified as the fold difference between aggregated contribution scores on motif hits versus those of flanking regions (±50 bp) within OCRs (see Methods ). Strikingly, many simple motifs and CEs exhibited strong correlations between their chromatin accessibility and transcriptional activity contribution values ( Figures 6E and 6F ; Pearson’s correlation R = 0.722 in K562, R = 0.816 in GM12878). However, notable differences were apparent for some motifs. For instance, the CTCF motif demonstrated a substantially greater contribution to chromatin accessibility than to transcriptional activity in both cell lines (K562: ChromBPNet fold 33.8 [max 43.5], GRACE fold 38.6 [max 231.2]; GM12878: ChromBPNet fold 32.9 [max 44.6], GRACE fold 11.8 [max 144.6]). This finding aligns with the known role of CTCF as an architectural transcription factor essential for maintaining genome organization and 3D chromatin structure 44 . Although the CTCF motif was not included in the original MPRA library design, GRACE successfully learned and identified the CTCF motif through measured activities derived from other sequences, highlighting the model’s robust motif inference capability ( Figure S6E ). Additionally, the GATA motif showed enhanced accessibility contributions in K562 (ChromBPNet fold 28.9 [max 43.5], GRACE fold 22.6 [max 231.2]) but minimal contributions in GM12878 (ChromBPNet fold −0.4 [max 44.6], GRACE fold 1.5 [max 144.6]). This pattern is consistent with the established role of GATA factors as pioneer transcription factors 45 involved in erythroid differentiation 46 . Conversely, some motifs displayed transcriptional activity without significant contributions to chromatin accessibility in a cell type-specific manner. POU-domain factors exhibited strong transcriptional activation but limited chromatin-opening activity in both cell lines (K562: ChromBPNet fold 0.0 [max 43.5], GRACE fold 58.9 [max 231.2]; GM12878: ChromBPNet fold 11.1 [max 44.6], GRACE fold 46.8 [max 144.6]). This observation aligns with prior reports demonstrating that OCT4, a POU-domain factor, lacks intrinsic chromatin-remodeling activity and requires cooperation with BRG1 to activate transcription, with BRG1 rather than OCT4 primarily driving chromatin accessibility 47 . Moreover, the TP53 motif displayed GM12878-specific transcriptional activity (K562: ChromBPNet fold 0.0, GRACE fold −1.4; GM12878: ChromBPNet fold 1.9, GRACE fold 41.5), which was consistent with elevated TP53 expression in GM12878 cells 48 . Notably, repressive motifs were also identified; however, their dynamic ranges were lower than those of their activating counterparts. SNAIL motifs exhibited transcriptional repression with minimal effects on chromatin accessibility (K562: ChromBPNet fold −2.1, GRACE fold −30.0; GM12878: ChromBPNet fold 2.4, GRACE fold −12.8), aligning with known repressive functions of SNAIL family TFs 49 . This finding suggests that SNAIL family members recruit co-repressors to pre-accessible chromatin regions without altering their accessibility and is consistent with the known interaction of RCOR/LSD1 with SNAG repressor domain-containing TFs 50 . In contrast with simple motifs, CEs generally contributed to both transcriptional activity and chromatin accessibility and to a higher extent. Notably, the highest-ranking motifs based on transcriptional activity and chromatin accessibility contributions were CEs. Specifically, AP-1-related CEs showed the greatest contributions in K562 cells, whereas the NFKB- and IRF-associated CEs displayed highest contributions in GM12878 cells ( Figures 6E and 6F ). We note that the JASPAR simple motif collection contained many closely related TF family motifs, whereas our CE collection featured distinct pairings and stereospecific configurations and therefore displayed much greater overall motif diversity. To directly evaluate the relative contributions of CEs compared with their constituent simple motifs, we performed a comparative analysis of their transcriptional activity and chromatin accessibility via CWMs alongside corresponding PWMs. The motif instances for both simple (n=153) and composite motifs (n=893) were identified using the two collections described above and their PWMs were generated from aligned motif sequences within OCRs of GM12878 and K562 cells using CEseek. GRACE CWMs were generated by aggregating contribution score tracks for each motif instance as in Figure 5D . ChromBPNet CWMs were similarly represented as aggregated contribution scores on the same motif hits within the accessible chromatin regions. The functional CEs exhibited distinct contribution patterns when compared with their constituent simple motifs. For example, in K562 cells, we have demonstrated that the SPI1[+]:-2:JUN[-] CE was enriched specifically within regions co-bound by SPI1 and JUN and exhibited synergistic transcriptional activity in MPRA assays ( Figure 4A ). Consistent with these findings, our analysis revealed that the SPI1 simple motif individually contributed minimally to chromatin accessibility and only modestly to transcriptional activity, whereas the JUN simple motif showed higher contributions to both features ( Figure 6G ). Remarkably, these contributions were substantially enhanced within the SPI1[+]:-2:JUN[-] CE context, underscoring the ability of CEs to amplify transcriptional and chromatin accessibility effects. Another illustrative example is the PRDM1[+]:-1:IRF-halfsite[+] CE which resembles an extended interferon response regulatory element (ISRE) with three IRF halfsites each separated by 2 bp. The IRF and PRDM1 motifs displayed minimal accessibility and transcriptional activity contributions but within the context of the CE, both transcriptional and accessibility contributions were markedly increased, highlighting the synergistic action of the CE ( Figure S6F ). Similarly, the NFYA[-]:10:SP3[+] CE demonstrated modest individual contributions from both the NFY and the KLF/SP simple motifs, primarily in transcriptional activity, with minimal effects on chromatin accessibility. Notably, the predicted contributions of both motifs were significantly amplified within the CE context, despite a substantial 10-base-pair spacing between the motifs ( Figure 6H ). A few CEs were predicted to selectively enhance transcriptional activity without affecting chromatin accessibility. For example, the contribution of OCT4[-]:-2:NFE2L1[+] CE to the predicted transcriptional activity was substantially greater than that of its constituent motifs individually; however, neither the simple motif nor the CE itself was predicted to contribute to chromatin accessibility ( Figure S6G ). Repressive CEs were also identified by the comparative analysis of the two deep learning models. The AR-halfsite and ZEB1 simple motifs individually demonstrated modest repressive activity contributions and minimal impacts on accessibility. However, the AR-halfsite[+]:-2:ZEB1[+] CE significantly increased repressive transcriptional activity and had a modest, yet detectable, repressive effect on chromatin accessibility relative to background levels. Together, these findings suggest that CEs can significantly increase transcriptional activity and chromatin accessibility contributions beyond the capacities of individual simple motifs, thereby amplifying their combinatorial regulatory logic. Most strikingly, the strong similarities between GRACE and ChromBPNet learned CWMs of CEs, using orthogonal experimental measurements and distinct deep learning computational methods, testified to the biological validity of the expanded TF motif lexicon. GRACE predicts impact of GWAS variants at TF motifs and CEs We next utilized GRACE to predict and interpret the functional impact of disease-associated variants, with a particular focus on their potential effects on CEs. To do so, we made use of a recent study on autoimmune disease-associated variants involving MPRAs performed with 150 bp reference and variant sequences in an Epstein-Barr virus-transformed B-cell line from a lupus patient 51 . From this MPRA dataset, a total of 1,719 expression-modulating variants (emVARs) were compiled, among which 319 exhibited significant allele-specific effects. We reasoned that functional variants would be predicted by GRACE to have altered transcriptional activity because of the single-nucleotide variant residing within a CE or a simple TF motif. Indeed, GRACE successfully predicted allele-specific effects with a Pearson correlation of 0.512 while focusing on the activities of 50 bp sequences flanking the variants ( Figure 7A ). Furthermore, the precision-recall curve showed that GRACE predictions can generate a precision of 73.6% at 20% recall ( Figure 7B ). Notably, for the 68 functional variants that were successfully predicted and interpreted by GRACE, a large fraction (30.9%, 21/68) appeared to represent loss- or gain-of-function mutations in CEs ( Figure 7C , Table S11). A striking example involves variants (rs3936134 and rs6011530) with a pair of CREB and ETS motifs ( Figures 7D and 7E ). In the case of the rs3936134 variant, the C→G change in the CREB motif was predicted by GRACE to be a loss-of-function mutation but did not result in alteration of the activity of the nearby ETS motif that was separated by 8 bp. In contrast, in the case of the rs9263880 variant, the G→A change in the ETS motif, which was again predicted to be a loss-of-function, also substantially impacted the contribution of the CREB motif that was juxtaposed. Importantly, the juxtaposed CREB and ETS motifs represented a CE, as the stereospecific configuration was enriched in co-bound regions of SPI1 and CREB1 in the GM12878 B-cell ChIP-seq datasets ( Figure S7A ). Additional examples of variants that were predicted to impact the transcriptional activities of either simple TF motifs or CEs are shown (Figures S7B-E). Thus, CEs, which are shown herein to represent an expanded TF motif lexicon, are also predicted by GRACE to account for the functional impact of many disease-associated GWAS variants. Given that CEs have increased contextual information and greater impact on transcription and chromatin structure than the simple TF motifs that comprise them ( Figures 3C , 6E , and 6F ), CEs will enhance the prioritization of variants in non-coding regions of human genomes for phenotypic analyses. Download figure Open in new tab Figure 7. GRACE model predicts and interprets GWAS variants in autoimmune diseases (A) GRACE predictions of GWAS variants. Scatter plot comparing GRACE-predicted and MPRA-measured allele-specific effects (log2 fold change) for autoimmune disease– associated GWAS variants (N=1719) (see Table S11 for the source of the MPRA dataset in the EBV-transformed B-cell line from a lupus patient). The Pearson correlation coefficient (PCC) with its p-value is indicated. (B) Precision-recall curve testing the performance of the GRACE model in identifying GWAS variants with significant allele-specific effects in the MPRA dataset. (C) Sankey diagram linking GRACE-identified simple and CE motifs that were predicted to account for changes in the transcriptional activities (gain or loss of function) of reference and variant allelic sequences (N=86) (from Figure 7A ). CEs and simple TF motifs are linked to their corresponding alleles with red and gray lines, respectively. (D, E) Examples of autoimmune disease-associated GWAS variants. Sequence logos showing nucleotide contributions for the reference (top row) and alternative alleles (bottom row) with their annotated TF motifs. The red dashed line indicates the variant position and nucleotide substitution. The heatmap in the middle row of each panel displays in silico saturation mutagenesis predictions generated by GRACE for the reference sequences. Both alternative alleles are predicted to result in loss-of-function. The ETS::CREB configuration in panel E is a CE (see Figure S7A ). Download figure Open in new tab Figure S7. GRACE model predicts and interprets GWAS variants in autoimmune diseases (A) Identification of a CREB::ETS CE using CEseek and ChIP-seq data for CREB1 and SPI1. The Venn diagram showing the overlap of the CREB1 and SPI1 binding regions (ChIP-seq) in GM12878 cells. Co-binding regions were further analyzed with CEseek, using singly bound regions as controls, to identify statistically enriched CREB::ETS configurations (heatmap). One such configuration corresponding to the motif impacted by the GWAS allele rs6011530 ( Figure 7E ) is highlighted with an orange box. The corresponding CREB::ETS CE PWM is shown below the heatmap. (B-E) Examples of autoimmune disease-associated GWAS variants with their GRACE interpretations. Sequence logos showing GRACE nucleotide contributions for the reference and alternative alleles along with their annotated simple motifs ( B, C ) or CEs ( D, E ). All alternative alleles are predicted to result in gain-of-function. The heatmaps in the middle row depict in silico saturation mutagenesis predictions generated by GRACE for the reference sequences. Discussion The combinatorial complexity of transcriptional regulation in metazoans arises from structurally diverse transcription factors (TFs) that act on densely clustered binding motifs within CREs (enhancers and silencers). While individual TF motifs within CREs have been extensively characterized, higher-order assemblies of simple motifs that are stereo-specifically juxtaposed to constitute composite elements (CEs) and may represent an additional layer in the gene regulatory code have remained underexplored. Here, we advance CE discovery, functional characterization and modeling by developing an integrated framework that combines computational prediction, high-throughput experimental validation, and deep learning. This framework systematically expands the transcriptional regulatory lexicon and demonstrates the prevalence of CEs as higher order fundamental units of the gene regulatory code in mammalian genomes. Prior CE discovery efforts have largely relied on focused molecular biological studies using electrophoretic mobility shift assays, low-throughput reporter assays, and structural analyses, which have delineated a small number of paradigmatic CEs, such as EICE 11 , AICE 12 , NFAT::AP-1 13 , 14 , and OCT::SOX 15 motifs. More recent high-throughput in vitro approaches, notably CAP-SELEX 16 , 52 , have expanded the catalog of TF-TF cooperative binding interactions via synthetic DNA sequences and E. coli -expressed TF domains. Despite its high biochemical throughput, CAP-SELEX is limited by the expression of full length or truncated mammalian proteins in E. coli and the in vitro selection of high-affinity synthetic CEs that may not recapitulate more varied endogenous sequence contexts. Furthermore, the transcriptional activities of such synthetic cCE sequences remain to be determined. Conversely, in vivo TF binding assays such as ChIP-seq and CUT&RUN provide valuable insights into the co-occupancy of TFs at genomic regions and can be used to identify cCEs; however, these methods are restricted by the availability of large-scale pairwise TF binding datasets in diverse cellular contexts 31 – 33 , 23 . Moreover, co-occupancy alone does not establish the presence of functional CEs or resolve their motif syntax or transcriptional activity. Finally, genome-wide discovery of CEs has remained technically challenging due to the absence of combinatorial TF transcriptional activity datasets within diverse mammalian cell types and states. To overcome these limitations, we developed a three-tiered integrative framework. First, the CEseek computational pipeline enabled systematic genome-wide discovery of statistically enriched binary TF motif configurations directly from open chromatin regions without requiring pairwise ChIP-seq or CUT&RUN data. Second, we functionally tested cCEs using a custom-designed TF-motif MPRA library containing short (30 bp) synthetic sequences embedding thousands of cCEs, enabling high-throughput quantification of their transcriptional activities in vivo . Third, we trained a newly developed GRACE deep learning model on the resulting MPRA datasets, enabling the model to learn contextually specific transcriptional regulatory lexicon, comprised of simple motifs and composite elements, at single-nucleotide resolution. Notably, the GRACE-derived contribution weight matrices strongly aligned with the CEseek-predicted position weight matrices, further validating the lexicon’s fidelity. The combined CEseek → MPRA → GRACE pipeline thus provides a scalable, generalizable framework for comprehensive CE discovery. Importantly, this integrated framework functions independently of CAP-SELEX and ChIP-Seq or CUT&RUN analyses to discover and analyze CEs. The latter datasets are utilized in subsequent steps in the analytical framework to provide evidence for cooperative DNA binding of relevant TFs in vitro and/or co-binding of the same TFs in vivo at genomic regions containing the discovered CEs. Our systematic analyses revealed a large and diverse CE repertoire, identifying 3,795 cCEs in immune and hematopoietic cell CREs, of which 893 demonstrated transcriptional activity in MPRA assays in three human cell types. The CE catalog is likely to undergo significant expansion as functional testing is expanded to even more diverse cell contexts. The AP-1 and KLF motifs were shown to be constituents of diverse sets of CEs. On the one hand, these findings demonstrate the varied functional interplay of AP-1 and KLF family members with distinct families of TFs; on the other hand, they suggest the selectivity of action of specific AP-1 and KLF family members with their cognate partners. Structural studies support direct physical interactions between TFs assembled at composite sites, suggesting that CEs encode protein complex assembly rules alongside DNA recognition 13 , 14 . We note that CEs remarkably constrain TF-TF interactions, thereby manifesting marked TF-family specificity. For example, while multiple ETS family members recognize the core ETS motif, only SPI1 and SPIB cooperatively bind to the canonical ETS::IRF composite element (EICE) with IRF4 or IRF8 11 . Similarly, OCT4 recognizes distinct CEs with SOX2 and SOX17, with spacing constraints that likely reflect differential protein-protein interaction interfaces 53 . These findings demonstrate that CEs encode not only sequence-specific information but also stereochemical and structural rules of TF complex assembly. CE function is highly sensitive to motif spacing and orientation. Subtle variations in spacing within NFY::KLF composite elements significantly alter predicted transcriptional activity, likely reflecting steric incompatibility of DNA-bound TFs. Although our TF-motif MPRA library was primarily designed to analyze CEs that synergistically activate transcription in a manner dependent on the non-additive interaction of their stereo-specifically juxtaposed simple motifs reflective of AND regulatory logic gates, it also revealed CEs with alternative activity patterns. In this regard, CEs were detected in which one of the simple motifs was activating whereas the second motif functioned to antagonize the activating function of its partner, thereby representing classical AND-NOT regulatory gates that could underlie alternative or bifurcating cell states. Finally, one CE displayed NOR logic, as it was composed of two activating simple motifs that manifested diminished combinatorial activity. Such regulatory logic may enable mutual attenuation of two gene expression programs that are counterproductive when simultaneously activated. Overall, the expanded lexicon of functionally distinct CE configurations can be used to combinatorially control genome activity via AND, AND-NOT and NOR activity gates. The distinct types of combinatorial regulatory logic enabled by CEs would allow for the integration of distinct signaling inputs into the genome by signaling-induced TF pairings whose assembly in turn is dictated by stereospecific positioning of their binding sites. The existence of widespread CE grammar has far-reaching implications for lineage specification, cellular state transitions and disease susceptibility. A recent study provided in vivo evidence that a single CE, an E-box–homeodomain “Coordinator” motif, modulates chromatin accessibility and gene expression programs that define the embryonic face and limb mesenchyme, underscoring the developmental importance of the CE 54 . As discussed above, CEs not only enable conditional integration of paired regulatory inputs but also may confer both robustness and evolvability to gene regulatory circuits. Small sequence changes that modulate CE configurations may allow rapid evolutionary rewiring of enhancer activity without altering individual TF binding specificities. Notably, immune-cell-specific CEs such as AICE and EICE exemplify how composite logic governs cell type-specific transcriptional programs and immune function 11 , 12 . Extensive dissection of CE configurations across developmental and pathological contexts will further illuminate their role in human disease mechanisms, particularly where mutations disrupt combinatorial rather than individual motif recognition. The synthetic potential of CEs is also substantial. Their defined grammar provides a modular design language for engineering synthetic enhancers with tunable transcriptional outputs, offering new opportunities in gene therapy, synthetic biology, and cell engineering. We note that the CEseek analysis also revealed a smaller set of TF motif pairs, despite their enrichment within cCREs, that did not demonstrate stereospecific configurations. These motif pairs likely utilize molecular mechanisms that do not involve protein or DNA conformation-directed TF-TF interactions to enable their combinatorial activities. They may depend on additional non-DNA-binding transcriptional complexes to enable functional interactions between the DNA-bound TFs that are not constrained by the spacing and orientation of the binding sites. These TF motif pairs appear to represent a complementary second-order lexicon based on a fuzzy or soft logic 18 , 55 . While the theoretical sequence space for CEs is vast, biological constraints appear to limit the realizable lexicon. Vertebrate TFs recognize ∼150-200 unique simple motifs; however, the CEseek → MPRA pipeline experimentally validated 893 candidate CEs, suggesting that functional composite grammar operates within a restricted subset of possible configurations. The strategic design of our MPRA library (∼50,000 sequences) enabled efficient sampling of this constrained space and permitted generalizable model training. Nonetheless, our approach was limited by the analysis of CEs within immune and hematopoietic cCREs, and their activities were assayed in three human cell contexts. Similar TF-motif MPRA libraries that focus on cCEs drawn from CREs in other cellular contexts and assay their activities in such contexts will enable the overall size of the expanded CE lexicon to be determined. Nevertheless, the work herein establishes the prevalence of the CE lexicon in the genome regulatory code and points to its evolutionary adaptability in enabling the control of varied and contextually specific patterns of genome activity. Remarkably, the GRACE model not only demonstrated strong predictive accuracy for the transcriptional activity of synthetic MPRA-tested sequences but also for endogenous promoters and enhancers. Despite training on a limited MPRA dataset, GRACE generalized well to genome-encoded regulatory sequences, which is consistent with the existence of a core regulatory grammar that governs transcriptional output in DNA and chromatin contexts. Importantly, the GRACE outputs converge with an orthogonal chromatin accessibility model, namely, ChromBPNet, revealing overlap but also regulatory divergence between chromatin accessibility and transcriptional activity codes. The coupling of CEseek discovery with deep learning models thus offers a systematic framework to dissect the relative contributions of simple motifs and CEs to transcriptional regulation, chromatin architecture, and higher-order regulatory logic. Approximately 30% of functionally significant nod-coding variants associated with autoimmune diseases map to composite elements, underscoring the importance of CE grammar in interpreting non-coding variation. As many disease-associated variants may disrupt combinatorial rather than individual TF binding, CE-aware variant annotation frameworks may improve the interpretation of genetic risk loci and identify previously unrecognized regulatory mechanisms underlying complex traits. Importantly, our short-sequence-based library design for MPRAs coupled with GRACE modeling offers a reference-independent platform for functional variant annotation, circumventing reliance on population-averaged reference genomes and enabling application to genetically diverse populations or somatic variation in cancer. This capability enhances personalized genomic medicine. Despite these advances, several important limitations remain. First, although our MPRA strategy captures transcriptional activity, it does not directly analyze cooperative or mutually antagonistic TF binding to CEs. Integration with protein binding microarrays 56 and CAP-SELEX 16 , 52 assays is needed to demonstrate cooperative binding of specific TF family members along with structurally guided exploration of molecular surfaces that generate combinatorial specificities on CEs. Second, while we focused on immune and hematopoietic regulatory elements, CE grammar remains to be comprehensively mapped across additional cell types, developmental stages, and many additional disease contexts. Third, perturbation studies targeting specific TF pairs or CE configurations will be essential for dissecting the functional dependencies and mechanistic consequences of CE disruption in vivo . RESOURCE AVAILABILITY Lead Contact Further requests and information concerning this study should be addressed to the lead contact, Harinder Singh ( harinder{at}pitt.edu ). Materials Availability The TF-motif MPRA plasmid library is available upon request. Data and Code Availability All sequencing data have been deposited into the IGVF Data Portal ( https://data.igvf.org/ ): IGVFDS6431ZPLN, IGVFDS7446JCCY, IGVFDS5083RSLD, IGVFDS6230NHQE, IGVFDS2616NTKE, IGVFDS7386LSZJ, IGVFDS4765WKNZ, IGVFDS4645LTYC, IGVFDS2722MXQO, IGVFDS5743OUYJ, IGVFDS1360GZSW, IGVFDS6576ZAJI, IGVFDS9693DMOS, IGVFDS0273TRYL, IGVFDS2337RQZA, IGVFDS5855AVHR, IGVFDS9385FKEJ, IGVFDS2689EODI, IGVFDS5988JIIM, IGVFDS2955VGJT. Supplemental tables have been deposited at GitHub: https://github.com/Jingyu-Fan/TF_motif_lexicon/tables . AUTHOR CONTRIBUTIONS Conceptualization: J.F., V.K.C., H.S. Data Curation: J.F. Formal Analysis: J.F., V.K.C. (CEseek analysis of ImmGen ATAC-seq data) Methodology: J.F., C.A.M., H.S. Project Administration: J.F., N.A.P., N.S., H.S. Software: J.F. (Python implementation of CEseek and development of GRACE), V.K.C. (original R prototype of CEseek) Investigation/ Data Generation: J.F. (CEseek analysis, MPRA design and data analysis, GRACE model development and analysis), D.B. (MPRA assays), N.A.P. (Management and sequencing of MPRA libraries), D.R.H., V.F.Y., C.H.H., J.P.R. (MPRA training), S.K., R.T., N.S. (MPRA library construction). Visualization: J.F. Resources: N.S., H.S. Writing – Original Draft: J.F., H.S. Writing – Review & Editing: J.F., V.K.C., D.B., N.A.P., D.R.H., P.G., V.F.Y., S.K., J.D., J.P.R., R.T., C.A.M., N.S., H.S. Supervision: C.A.M., N.S., H.S. Funding Acquisition: H.S., N.S., J.D. J.F. led the study; he re-implemented CEseek in Python, carried out various CEseek analyses, designed the MPRA library, performed all MPRA data analyses, developed the GRACE framework and performed GRACE analyses, curated and deposited the datasets, generated figures, and drafted the initial manuscript after its logical structuring with H.S.. V.K.C. created the original R prototype of CEseek and produced the initial cCE predictions using the ImmGen dataset and reported in BioRxiv 55 . D.B., N.A.P., D.R.H., S.K., J.P.R., R.T., N.S., and H.S. were involved in the design, construction and characterization of the TF-motif MPRA library and generation of the experimental data involving transient transfections in GM12878, K562 and Jurkat cell lines. N.A.P., V.F.Y generated the ATAC-seq data in GM12878 cells. P.G. trained ChromBPNet models in K562 and GM12878 cells and generated contribution scores on OCRs. C.A.M. guided the development of GRACE, refined the overall methodology together with J.F. and H.S., and provided critical technical guidance and resources. H.S. guided the scientific direction, secured funding, and oversaw all phases of the project. J.F. and H.S. wrote the original draft. V.C., D.B., N.A.P., D.R.H., P.G., V.F.Y., S.K., J.D., J.P.R., R.T., C.A.M., N.S., reviewed the manuscript. J.F. and H.S. made all revisions based on suggested changes. All the authors have read and approved the final version of the manuscript and agree to be accountable for all aspects of the work. DECLARATION OF INTEREST The authors declare no competing interests. SUPPLEMENTAL INFORMATION Table S1. Position probability matrices of the simple TF motifs and cCEs in this study, related to Figure 2 . Table S2. TF ChIP-seq datasets of K562 and GM12878 cells compiled from CistromeDB, related to Figures 2 and 4 . Table S3. Simple TF motif and cCE sequences (30-mers) used in design of the TF-motif MPRA library, related to Figure 3 . Table S4. Transcriptional activities of the designed sequences (30-mers) in the TF-motif MPRA library generated by transfections in K562, GM12878 and Jurkat cells, related to Figure 3 . Table S5. Transcriptional activities of simple motif and cCE containing sequences along with corresponding motif activities in K562, GM12878 and Jurkat cells, related to Figure 3 . Table S6. Synergistic activities of CEs in K562, GM12878 and Jurkat cells, related to Figure 3 . Table S7. Enrichment of CEs in co-bound regions identified in cognate TF ChIP-seq datasets, in K562 and GM12878 cells, related to Figure 4 . Table S8. CWM-PWM similarities of activating cCEs in K562, GM12878 and Jurkat cells, related to Figure 5 . Table S9. Clusters (K-means) of promoters and enhancers based on GRACE activity predictions in GM12878 cells, related to Figure 6 . Table S10. Motif activity and accessibility contribution scores of JASPAR simple motifs and functional CEs in this study, related to Figure 6 . Table S11. GRACE activity predictions and motif interpretations of GWAS variants, related to Figure 7 . Table S12. Primers, probes and sequencing indices used in this study, related to Figure 3 . METHODS CEseek computational framework The CEseek pipeline was developed to systematically identify and analyze composite elements (CEs) within cis-regulatory DNA regions by leveraging position probability matrices (PPMs) to scan genomic sequences for transcription factor (TF) motif hits. Through pairwise motif scanning, CEseek evaluates combinations of all four motif orientations and user-defined spacing intervals, detecting motif co-occurrence patterns suggestive of potential cooperative or antagonistic interactions. Aligned motif hits, including spacer regions, are extracted and used to construct PPMs for the resulting CEs. By quantifying all possible binary configurations for user-defined spacing intervals, CEseek identifies dominant configurations exhibiting orientation and spacing preferences. Statistical enrichment of CEs is assessed via Fisher’s exact test, which compares motif pair frequencies in a test set against those in an appropriate background set (dinucleotide-shuffled sequences or other genomic sequences). CEseek was implemented using Python programming language and incorporates several specialized libraries, including PyTorch for efficient tensor operations and computations, pybedtools 57 , 58 for genomic interval manipulation, and pandas for data management and processing. The integration of GPU acceleration through PyTorch substantially enhances computational performance during motif scanning and sequence analyses. The CEseek tool is publicly available at https://github.com/Jingyu-Fan/CEseek . An archived version implemented in R is accessible at https://github.com/viren-v/CEseek . Compilation of a core simple TF motif set A comprehensive library of TF motifs was compiled for systematic CEseek analysis. Initial collections of mouse and human TF motifs were retrieved from the CIS-BP 6 database and the HOMER 27 motif repository. To reduce redundancy within the library, similarity among all the PPMs was assessed via HOMER’s motif similarity tool. When multiple TFs shared highly similar PPMs (similarity score >0.9), the motif with the highest information content was selected as the representative. This process resulted in a curated, non-redundant final collection consisting of 153 core TF motifs. Delineation of candidate CEs in ImmGen cis-Regulatory atlas Chromatin accessibility data generated by the ImmGen 22 consortium, which includes ATAC-seq profiles for 86 immune and hematopoietic cell types/states along with 4 non-hematopoietic cell types, were utilized for systematic CEseek analysis. The dataset comprises 512,595 open chromatin regions (OCRs), each of which is extended on both sides, yielding standardized 180-base pair (bp) regulatory regions. CEseek was used to scan these OCR sequences systematically, evaluating all possible pairwise motif combinations derived from the library of 153 core TF motifs, with spacer lengths ranging from −5 to +10 base pairs. Each OCR was subsequently assigned to specific cell types or states based on its accessibility, defined as peaks falling within the 95th percentile for that particular cell type or state relative to others in the ImmGen dataset. Fisher’s exact test was employed to identify CE configurations that were statistically enriched in cell type- or state-specific OCRs, and the remaining OCRs were used as a background set (enrichment significance threshold p < 1e-20). Additionally, an independent analysis employing scrambled sequences as a control further validated the prevalence of cCEs, identifying 22,698 significantly enriched cCEs (p < 1e-10). Finally, dominant configurations characterized by preferential spacing and orientation between paired motifs were selected (significant motif pairs with ≤6 enriched configurations out of a possible of 64), resulting in 5,063 high-confidence cCEs. Corresponding PPMs for these selected cCEs were subsequently generated using aligned motif hits identified across the entire ImmGen OCR dataset. Comparison of cCEs with motifs from TF ChIP-seq and CAP-SELEX datasets To provide evidence for TF binding in vivo and in vitro, the cCEs identified by CEseek were compared with motifs obtained from chromatin immunoprecipitation sequencing (ChIP-seq) and CAP-SELEX datasets. Specifically, a comprehensive collection of uniformly processed TF ChIP-seq datasets from CistromeDB 31 – 33 was assembled for hematopoietic (K562) and B cell (GM12878) cell lines. This dataset comprised 962 ChIP-seq experiments corresponding to 337 unique TFs in K562 cells and 370 experiments representing 160 unique TFs in GM12878 cells. From these datasets, a total of 14,897 de novo motifs were identified. Subsequently, each cCE identified by CEseek was systematically compared against this database of de novo motifs using the Tomtom 59 motif comparison tool to assess motif similarity. Additionally, a similar comparative analysis was conducted between cCEs and a set of 618 CAP-SELEX-derived heterodimeric motifs 16 . In this case, each CAP-SELEX motif was compared against the entire database of cCE motifs using the same computational approach (Tomtom). Design of TF-motif MPRA library A Massively Parallel Reporter Assay (MPRA) library comprising simple TF motifs and CEs was designed to enable comprehensive, systematic evaluation of motif activities in vivo . This library consisted of 60,714 sequences, each 30 bp in length. Initially, the set of 5,063 cCEs was scanned across ImmGen ATAC-seq peaks, and the top 10 motif hits with the highest scores for each cCE were selected. These selected hits were subsequently extended to uniform length 30 bp sequences. The sequences assigned to multiple cCEs were removed, resulting in a refined set of 3,795 unique cCEs. For each of these CEs, the four highest-scoring sequences were retained, and three additional mutated variants were systematically generated: one variant with mutations introduced into the left motif, another with mutations in the right motif, and a third variant with mutations in both motifs. Mutations involved substitutions at two consecutive positions of highest information content, converting nucleotide A to C and T to G, and vice versa. This strategy allowed for precise quantification of individual motif contributions within cCEs. Additionally, the MPRA library included 153 core simple motifs recognized by vertebrate TFs. For each simple motif, genomic sequences precisely matching motif length were identified and extended with selected random sequences to prevent unintended creation of alternative motifs. Up to 20 distinct sequences were selected for each simple motif based on their highest motif scores, alongside corresponding mutated sequences with substitutions at two consecutive positions of maximal information content. As controls, the library also contained 1,000 randomly generated 30 bp sequences without any recognizable motifs, serving as baseline measurements of background activity. Construction of TF-motif MPRA library The TF motif library sequences were cloned into the pGL4:23:ΔlucΔxbaI vector as described previously 21 with modifications. To create the mpraΔorf library, barcoded oligos were inserted into Sfi I -digested pGL4:23:ΔlucΔxbaI by Gibson assembly (NEB, E2611) via the use of 1.1 µg of oligos and 1 µg of digested vector in a 50 µL reaction mixture and incubated for 60 min at 50°C. The reaction mixture was purified with 1.2x SPRI and eluted in 20 µL of elution buffer. Test transformation was performed to determine library coverage via the electroporation of 1 µL of ligated vector into 50 µL of electrocompetent 10-beta E. coli (NEB, C3020K) (2kV, 200 ohm, 25 µF). Similarly, for the scaled-up transformation, 1 µL of the library was transformed into 50 µL of electrocompetent cells. The electroporated cells were recovered in 950 µL of SOC media and then split into ten 1 mL aliquots of SOC. These cultures were subsequently recovered for an hour at 37°C, after which each culture was individually expanded in 20 mL of LB supplemented with 100 μg/mL carbenicillin (Sigma-Aldrich, C1389) in a shaker at 37°C for 6.5 hours. After incubation, all the cultures were pooled, and plasmid isolation was performed (Qiagen, 12963). To validate library complexity and connect barcodes to CE sequences, Illumina libraries were prepared for the mpra:Δorf library as described previously 60 and sequenced via 2×150 PE reads on the Illumina NovaSeq platform with 15% PhiX spike-in. To create the final mpra:gfp library, 10 µg of mpra:Δorf library plasmid was linearized with 100 units of Asi SI (NEB, R0630) in a final 400 µL reaction and incubated overnight at 37°C. The linearized product was then column purified (Qiagen, 28104) and eluted in 30 µL. The GFP amplicon was amplified from pMPRAv3:minP-GFP via Q5 NEBNEXT Hot-Start (NEB M0494), 0.5 µM primer 200 and 0.5 µM primer 201 (Table S12), following cycle conditions of 98°C for 30 s, 20 cycles (98°C for 10 s, 60°C for 15 s, 72°C for 45 s), and 72°C for 5 min. The amplified product was incubated with Dpn1 for 30 min at 37°C, followed by 0.5x reverse SPRI and 1.5x forward SPRI purification and elution in 40 µL EB. A second PCR was then performed using the 1:100 diluted purified PCR1 GFP amplicon as PCR1 and then column purified. This amplicon containing a minimal promoter, a GFP open reading frame and a partial 3’ UTR was then inserted by Gibson assembly using 1.6 µg of Asi SI -linearized mpraΔorf plasmid and 5.28 µg of the GFP amplicon in a 400 µL reaction for 90 min at 50°C, followed by 1.5x SPRI purification. The total recovered volume was redigested to remove the remaining uncut vectors by incubation with 50 units of Asi SI , 5 units of RecBCD (NEB, M0345), 10 µg of BSA, 1 mM ATP, and 1x NEB Buffer 4 in a 100-µL reaction mixture incubated overnight at 37°C, followed by 1.5x SPRI purification and elution with 40 µL of EB. To generate the final transfection-ready MPRA library, a test transformation was first performed using 1 µL of the final eluted mpra:gfp plasmid library and 50 µL of 10-beta cells to determine the library coverage. Finally, for scale-up transformation, 8 µL of mpra:gfp plasmid was electroporated (2kV, 200 ohm, 25 µF) into 200 µL of 10-beta cells. The electroporated bacteria were recovered in 12 mL of SOC, split across 6, 2 mL aliquots and incubated for 1 hour at 37°C. Then, each 2 mL culture was expanded to 500 mL of LB with 100 µg/mL carbenicillin and incubated for 16 hours at 37°C, followed by plasmid isolation via the Qiagen Gigaprep Kit (Qiagen, 12991). Transient transfections of TF-motif library in human cell lines GM12878 cells (Coriell Institute for Medical Research) were cultured in RPMI medium (Corning, 10-040-CV) supplemented with 15% FBS (Gibco A52568-01), 1% GlutaMAX (Gibco 35050-061) and 1% Pen-Strep (Corning, 30-002-CI). For each transfection, 1*10 7 cells were mixed with 10 μg of the TF-motif MPRA library in 100 μl of RPMI. These cells were transfected with a Neon transfection system (Thermo Fisher Scientific, MPK5000) using a kit (MPK10096B) with 3 pulses of 1200 V for 20 ms. A total of five biological replicates were performed, each using cells at a density of ∼1 million cells/ml for the transfections. For each replicate, a total of 1.5*10 8 cells were pelleted at 200 × g and resuspended in 1.5 ml of RPMI medium containing 150 μg of the TF motif library. After transfection, each replicate was recovered at a density of 5*10 5 cells/ml in 300 ml of RPMI medium containing 15% FBS, 1% GlutaMAX and 1% Pen-Strep. After 24 hours, the cells were pelleted at 200 × g and washed once with PBS. The cells were resuspended in 15 ml of RLT buffer provided with Qiagen Maxi RNeasy (Qiagen, 75162) and 300 mM DTT (G-Biosciences, 786227). The cells were then homogenized with a tissue homogenizer, TH (OMNI international), for 1 min at maximum speed and the lysates stored at –80°C. Jurkat T cells were cultured in RPMI medium (Gibco 61870-036) supplemented with 10% FBS (Gibco A52568-01) and 1% Pen-Strep (Corning, 30-002-CI). All other steps were the same as for GM12878 cells with the exception that the transfections using the Neon system involved 3 pulses of 1550 V for 10 ms. K562 cells (CIMR) were cultured in Iscove’s modified Dulbecco’s medium (Corning 10-016-CV) supplemented with 10% FBS and 1% Pen-Strep. For each transfection, 1*10 7 cells were mixed with 5 μg of the MPRA library in 100 μl of RPMI and transfected with the Neon transfection system using a kit with three pulses of 1450 V for 10 ms. Five biological replicates were performed as with GM12878 cells. For each replicate, a total of 1.5*10 8 cells were pelleted at 200 × g and resuspended in 1.5 ml of buffer R provided in the Neon Transfection Kit containing 150 μg of the TF motif library. After transfection, each replicate was recovered at a density of 5*10 5 cells/ml in 300 ml of RPMI medium containing 15% FBS and 1% Pen-Strep. After 24 hours, the cells were pelleted at 200 × g and washed once with PBS. The cells were resuspended in 15 ml of RLT buffer and 300 mM DTT. The cells were then homogenized for 1 min at maximum speed and lysates were stored at –80°C. RNA extraction and cDNA synthesis Total RNA was isolated from cells via Qiagen Maxi RNeasy (Qiagen, 75162) following the manufacturer’s protocol, including on-column DNase digestion by RNase-free DNase (Qiagen 79254). The final purified RNA was treated with 5 µl of SUPERase-In (Thermo Fisher Scientific, AM2696). Another DNase treatment was performed on total RNA via the addition of 5 µL of Turbo DNase (Thermo Fisher Scientific, AM2238) in 750 µL of total volume for 1 hour at 37°C. The reaction was terminated by the addition of 7.5 µL of 10% SDS (Life Technologies, 15553--035) and 75 µL of 0.5 M EDTA (Thermo Fisher Scientific, AM9260G), followed by incubation at 70°C for 5 minutes. All the DNase-treated RNA was then subjected to GFP mRNA pulldown. For this reaction, 900 µL of 20X SSC mixture (Life Technologies, 15557-044), 1800 µL of formamide mixture (Sigma Aldrich, 75--12--7) and 2 µL of 100 µM biotin-labeled GFP probe mixture (GFP_BiotinCapture_1--3, Table S12) were added, and the volume was adjusted to 3600 µL. The mixture was then incubated at 65°C for 2.5 hours with intermittent inverting of the tubes every 30 min. The biotin probes were captured via 400 µL of prewashed streptavidin beads (Thermo Fisher Scientific, 65002). The streptavidin beads were washed twice with Buffer-A containing 0.1 M NaOH (Sigma-Aldrich, S2770) and 0.05 M NaCl (Thermo Fisher Scientific, AM9760G) and once with Buffer-B containing 0.1 M NaCl. The beads were then eluted in 500 µL of 20X SSC and added to the RNA probe sample. The hybridized RNA-probe-bead mixture was mixed with a HulaMixer (Thermo Fisher Scientific, V.3A01) at room temperature for 15 minutes. The beads were captured with a DynaMag magnet (Thermo Fisher Scientific, 12321D) and washed once with 1x SSC and twice with 0.1x SSC. The extraction of RNA was performed by the addition of 25 µL of DEPC-treated water (Ambion, AM9906) and heating of the mixture for 2 minutes at 80°C, followed by immediate collection of the eluent on a magnet. A second elution was performed by incubating the beads with an additional 25 µL of water at 80°C. The eluted RNA was then processed for the final DNase treatment with 1 µL of Turbo DNase in a total of 56 µL of reaction. The mixture was incubated for 60 minutes at 37°C, followed by inactivation with 1 µL of 10% SDS. The final DNase-treated GFP mRNA was eluted via RNAClean XP SPRI beads (Beckman, A63987) in 35 µL of DEPC-treated water. First-strand cDNA was synthesized via the use of 30 µL of DNase-treated GFP mRNA with SuperScript III and primer 19 (Table S12) according to the manufacturer’s recommended protocol. Single-stranded cDNA was purified via AMPure SPRI (Beckman, A63881) beads and eluted in 30 µL of EB. Generation and sequencing of Tag-seq libraries To minimize amplification bias during the creation of cDNA tag sequencing libraries, samples were amplified via qRT-PCR to estimate the relative concentrations of GFP cDNA. qRT-PCR was performed via the use of a 1 µL cDNA sample in a 10 µL PCR mixture containing 5 µL of Q5 NEBNext Ultra II master mix, 1.7 µL of SYBR Green I diluted 1:10,000 (Life Technologies, S-7567) and 0.5 µM of TruSeq_Universal_Adapter and MPRA_Illumina_GFP_F primers (Table S12). The samples were amplified with the following conditions: 98°C for 20 s; 40 cycles of 98°C for 10 s, 62°C for 15 s, and 72°C for 30 s; and 72°C for 2 min, followed by melt curve analysis. For sequencing of plasmid library, serial dilutions of the plasmid library were prepared from 1000 pg to 1 fg via 10-fold dilutions. A standard curve was plotted with the plasmid library dilutions and based on the threshold cycles for the cDNA samples, a plasmid dilution with the same cycle threshold value was used for adaptor ligation. To add Illumina sequencing adapters, cDNA samples and 5 plasmid library controls were diluted to normalize the cDNA replicate with the lowest concentration or highest threshold cycle, and 10 µL of normalized sample was amplified using the reaction conditions from the qRT-PCR scaled to 50 µL, without adding SYBR Green to the reaction and by using only n-1 amplification cycles (where n= the cycle number obtained from qRT-PCR). The amplified cDNA was 2x SPRI purified and eluted in 30 µL of EB. A second PCR was performed to add the individual index primers to each sample. The PCR was performed with 20 µL of purified PCR 1 elute in a 50 µL Q5 NEBNext Ultra II reaction with 0.5 µM of TruSeq_Universal_Adapter primer and Illumina_Multiplex primer containing a unique 8 bp index for sample demultiplexing postsequencing (Table S12). The samples were amplified at 98°C for 20 seconds, 6 cycles of 98°C for 10 sec, 62°C for 15 sec, and 72°C for 30 sec, and 72°C for 2 minutes. Indexed libraries were 2x SPRI purified and quantified on a Qubit Flex (Thermo Fisher Scientific, Q33326) using a dsDNA HS Assay Kit (Thermo Fisher Scientific, Q33230) and pooled according to molar estimates from TapeStation (Agilent Technologies, 4200 TapeStation) quantification. The samples were sequenced via 1×20 bp reads on the Illumina NextSeq 2000 platform, with 80 M reads per cDNA library replicate and 60 M reads per plasmid library replicate. MPRA data processing To process the MPRA data, we employed an established pipeline ( https://github.com/tewhey-lab/MPRA_oligo_barcode_pipeline ). The barcode-oligo sequencing FASTQ files were parsed using MPRAmatch function to create the barcode-oligo dictionary, retaining only barcode-oligo links with exact 30 bp matches. Barcode count tables were subsequently generated using MPRAcount function. To quantify the activities of simple motifs and CEs, the analysis was conducted at three levels: barcode, oligo, and motif. First, barcode-level count tables were processed with DESeq2 34 to remove barcodes with excessively large counts inconsistent with those of the biological replicates. This filtering used Cook’s distance, applying a 99% quantile threshold to exclude outliers. After filtering, the barcodes were aggregated into their corresponding oligos, and only those with at least 30 unique barcodes were retained for further analysis. The activities of the 30 bp oligos were assessed by performing DESeq2 differential analysis, which involved comparing the total barcode counts associated with each oligo in the RNA and DNA samples. Normalization and size factor estimation were performed using DESeq2’s internal methods. Additionally, cell type-specific dispersion values were calculated by analyzing only RNA replicates, and these refined dispersions were used to improve the accuracy of the final analysis. To further increase the reliability of the differential activity measurements, log fold change shrinkage was applied via the lfcShrink function from DESeq2, which helps reduce the impact of extreme log fold changes in oligos with low counts and provides more stable effect size estimates. To assess the activity of simple and composite TF motifs accurately, a comprehensive approach was employed. For simple motifs, activities were determined by comparing the RNA/DNA log2-fold change in oligos containing the simple motif (sequence activity) with matched mutated control oligos (motif activity). For composite motifs, sequence and motif activities were estimated similarly by comparing RNA/DNA ratios for oligos containing the CE to those in matched double-mutant control oligos. Additionally, the contributions of individual motifs (motif1 and motif2) within a CE were examined by comparing the RNA/DNA ratios in oligos with one mutated motif to those in matched double-mutant control oligos. Synergistic effects between motifs within a CE were assessed by comparing the motif activity of the CE to the sum of the contributions from motif 1 and motif 2. Simple and composite motifs were classified as activating or repressing based on both sequence and motif activity measures. Specifically, a motif was considered activating or repressing if both activities were significant (false discovery rate, FDR < 0.05) and exhibited the same directionality (activating or repressing). Synergistic effects were considered statistically significant when the p-value was < 0.05. Enrichment analysis of CEs in TF ChIP-seq datasets To determine whether CEs are specifically associated with co-binding of cognate TFs in vivo , a three-step analysis pipeline was implemented. First, simple TF motifs within each CE were matched with molecularly relevant ChIP-seq datasets. A simple motif was considered matched to a TF ChIP-seq dataset if the TF targeted by the ChIP-seq experiment belonged to the same TF family as the simple motif or if the simple motif significantly matched one of the top 10 motifs from the TF ChIP-seq dataset (FDR < 0.05, assessed using Tomtom). Second, all CEs exhibiting significant transcriptional activity (activating or repressing) and synergistic effects were evaluated for enrichment within peaks from matched TF ChIP-seq datasets. This enrichment analysis was conducted using CEseek, with dinucleotide-shuffled sequences (four times the number of test sequences) serving as the background. Statistical significance was assessed using Fisher’s exact test, with a stringent enrichment threshold of p ≤ 1e-10. Finally, a selected set of CEs was assessed for enrichment specifically within regions co-bound by both relevant TFs. This analysis was conducted by comparing the occurrence of the CE within co-bound ChIP-seq peak regions to its occurrence within singly bound regions as the background, employing CEseek and Fisher’s exact tests for statistical evaluation. GRACE model To predict cell type-specific regulatory activities in short DNA sequences, we developed a predictive model named GRACE (Genome Regulatory Analysis of Simple and Composite Transcriptional Elements). The GRACE model architecture starts with a convolutional layer consisting of multiple convolution blocks applied to DNA sequences that are transformed into one-hot encoded matrices. These blocks use different kernel sizes to capture motifs of varying lengths. This is followed by a transformer-based attention mechanism that analyzes motif interactions. The convolutional layer helps identify transcription factor (TF) binding sites of various sizes, whereas the self-attention mechanism, which is enhanced with positional encoding, captures the complex interactions between these motifs. This combination allows the model to effectively represent and analyze intricate regulatory elements within sequences. To mitigate overfitting and enhance the model’s generalizability, we employed dropout layers and data augmentation techniques, such as random sequence additions and reverse complement usage, to diversify the training data. Additionally, early stopping was used to halt training once validation loss failed to improve for a certain number of epochs, ensuring optimal model performance without overfitting. The GRACE model was implemented via PyTorch and PyTorch Lightning, with a custom DataModule that handles data loading, splitting, augmentation, and normalization of activity values, ensuring streamlined and efficient data handling. The model was trained on MPRA experimental datasets containing oligo activities, which were preprocessed and divided into training, validation, and test sets. We assigned variant sequences that were derived from a given parent sequence to the same subset, preventing data leakage, thereby enhancing the model’s ability to generalize predictions to all possible single nucleotide variants. Preprocessing involves filtering sequences by barcode count, augmenting them with random segments to mitigate positional bias, and incorporating reverse complement sequences to improve robustness and generalizability. Sequence padding was applied to ensure uniform input lengths for efficient batch processing. The model’s attention mechanism was restricted to 30 bp blocks, allowing the model to focus on local regulatory features while maintaining computational efficiency and avoiding artificial long-distance attention created by random sequence augmentation. Training was conducted via the AdamW optimizer, which incorporates weight decay to prevent overfitting, along with a learning rate scheduler to dynamically adjust the learning rate throughout training, improving convergence. Early stopping was also employed to prevent overfitting by monitoring the validation loss and halting training once the performance ceased to improve. The evaluation metrics included the mean squared error (MSE), Spearman’s correlation, and Pearson’s correlation, providing a comprehensive assessment of the model’s ability to predict regulatory activity. To validate the model, predictions on withheld test sequences were compared directly to experimentally measured activities, ensuring reliable model performance. The code for the GRACE model can be found at https://github.com/Jingyu-Fan/GRACE . GRACE contribution weight matrix generation and comparison with CEseek PWMs To elucidate the motifs contributing to transcriptional activity, in silico saturation mutagenesis was performed, systematically mutating each nucleotide position within the input DNA sequences to evaluate their individual contributions to transcriptional activity. For each 30-bp oligonucleotide, single-base substitutions were introduced at every position, and the transcriptional activities (log2-fold changes) of these variants were predicted using GRACE. The predictions were organized into 4 × L matrices (rows corresponding to nucleotides A, C, G, and T; columns corresponding to sequence positions). Contribution weight matrices (CWMs) were then generated by subtracting the mean log2 activity of each position from the activity of each nucleotide variant at that position. This normalization step removed position-specific biases. To produce a sparse representation highlighting key nucleotide contributions, only the scores corresponding to nucleotides present in the original (wildtype) sequences were retained, with all other entries set to zero. The resulting matrices were visualized as contribution scores. To compare these GRACE-derived CWMs systematically with CEseek-generated PWMs, the CWMs were first converted into position probability matrices (PPMs) via the Softmax function, ensuring that the probabilities at each position were summed to 1. To assess similarity, each resulting PPM was compared to its corresponding CEseek-derived PWM via Tomtom. For robust statistical comparison, 1,000 randomized background PWMs were generated for each PWM by shuffling positions and nucleotides. The similarity between each CWM-derived PPM and its corresponding CEseek-derived PWM, as well as the background PWMs, was evaluated, and the false discovery rate (FDR) was calculated. Multiple oligonucleotides representing each CE were analyzed, and the median of the negative log10-transformed FDR values was computed. A threshold of FDR < 0.05 was employed to define significant similarity between a CE’s GRACE-derived CWM and its corresponding CEseek-derived PWM. GRACE predictions of transcriptional activity of open chromatin regions The GRACE model was applied to predict the transcriptional activities of the OCRs. Because GRACE requires 50 bp sequences as inputs and OCRs typically exceed this length, a stepwise procedure was implemented. First, each OCR was segmented into overlapping 50 bp subsequences using a 1 bp sliding window, and the activity (log2-fold change, log2FC) of each subsequence was predicted using GRACE. All the predicted activities for the 50-mer sequences within each OCR were subsequently ranked from highest to lowest. The subsequence with the highest predicted activity was selected and included in the overall activity score of the OCR if it exceeded a predefined threshold (log2FC = 0.2). This process was iteratively repeated, excluding subsequences overlapping previously selected segments, until no remaining 50-mer sequences surpassed the activity threshold. The sum of the selected subsequences’ activities represented the final GRACE-predicted transcriptional activity of the OCR. The H3K27ac ChIP-seq signal of an OCR is calculated as the average read per million (RPM) for the 1 kb size window centered at the OCR and the ATAC-seq signal. GRACE predictions of transcriptional activity at promoters and enhancers The GRACE model was used to predict transcriptional activities at promoters and enhancers in GM12878 cells. Promoters included in the analysis were those that were associated with genes that manifested at least 1 transcript per million (TPM) nascent RNA signal from Global Run-On sequencing (GRO-cap). Enhancers used in the analysis were OCRs that did not overlap with any annotated transcriptional start sites (+/- 500bp). Contribution score peaks (seqlets) were called from the genome-wide GRACE contribution score track on the OCRs via macs3’s bdgpeakcall function 61 . The GRACE activity predictions for enhancers centered on the GRACE-defined seqlets were used for k-means clustering. Finally, the median GRO-cap read count at each nucleotide position relative to the center of the TSS or enhancer seqlet was calculated for each cluster. Systematic analysis of contributions of simple motifs and CEs to transcriptional activity and chromatin accessibility The contributions of simple motifs and CEs to transcriptional activity and chromatin accessibility were systematically evaluated using genome-wide contribution scores derived from the GRACE and ChromBPNet models, respectively. For this analysis, the collection of simple motifs included 879 motifs sourced from the JASPAR 2024 vertebrate core collection 43 . The CE collection consisted of 951 unique CEs identified as being active in the MPRA experiments conducted in the K562, GM12878, and Jurkat cell lines. To identify instances of simple and composite motif hits within OCRs, CEseek was employed. The contribution scores obtained from GRACE and ChromBPNet were organized into CWMs, which were formatted as 4 × L matrices (rows representing nucleotides A, C, G, and T; columns representing sequence positions). For each motif occurrence, the contribution scores were summed and averaged by dividing by the number of motif hits. Each averaged CWM was subsequently normalized through elementwise multiplication with the corresponding PPM derived from CEseek-detected motif hits. The resulting sum was divided by the motif length to yield the motif’s normalized contribution score. The baseline contribution scores were computed similarly using ±50 bp flanking regions around motif hits. The final motif contribution was calculated as the fold difference between the motif hit region and the baseline flanking region. GRACE analysis of GWAS variants associated with autoimmune diseases The GRACE model was employed to predict and interpret the functional impact of disease-associated genetic variants identified through genome-wide association studies (GWASs). A recently published MPRA dataset comprising functional assays conducted with 150 bp reference and variant sequences in an Epstein-Barr virus-transformed B-cell line derived from a lupus patient was used 51 . From this MPRA dataset, expression-modulating variants demonstrating ambiguous allele-specific effects without statistical significance (|log2FC| > 0.1 and FDR > 0.05) were excluded. A total of 1,719 expression-modulating variants were subsequently compiled, 319 of which exhibited significant allele-specific effects (|log2FC| > 0.6 and FDR < 0.05). To predict the variant effects, GRACE activity predictions were performed on 50 bp sequences centered on each variant site, comparing reference and alternative alleles. The difference in predicted activity between the two alleles was utilized as a measure of the variant effect. An absolute value of 0.2 of a GRACE-predicted variant effect was used as the threshold for determining significant variant effects. Applying this criterion resulted in the successful prediction of 68 functional variants. GRACE-generated contribution scores for these 50 bp sequences were then analyzed to interpret the TF motifs (simple or CE) underlying the observed variant effects. ACKNOWLEDGEMENTS This work was primarily supported by the National Institutes of Health (NIH) grant 5U01HG012041 (to N.S., J.D., and H.S.). This work was also funded by NIH grant RC2DK122376 (H.S.). Next-generation sequencing was performed by either the High-Throughput Genomics Core (UPMC) (formally the UPMC Genome Center) or the Health Sciences Sequencing Core (University of Pittsburgh). This research was supported in part by the University of Pittsburgh Center for Research Computing and Data, RRID:SCR_022735, through the resources provided. Specifically, this work used the HTC cluster, which is supported by the NIH award number S10OD028483. This work used bridge-2 (EM, RM, GPU) and ocean systems at the Pittsburgh Super Computing Center through the allocation of CIS240372 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by the U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. We thank Anshul Kundaje (Stanford), Lee Grimes and Nathan Salomonis (Cincinnati Children’s Hospital Medical Center) and Song (Stephen) Yi (Baylor) for valuable discussions and feedback. Funder Information Declared National Institutes of Health, https://ror.org/01cwqze88 , 5U01HG012041 , RC2DK122376 REFERENCES 1. ↵ Johnson , P.F. , and McKnight , S.L . ( 1989 ). Eukaryotic transcriptional regulatory proteins . Annu Rev Biochem 58 , 799 – 839 . doi: 10.1146/annurev.bi.58.070189.004055 . OpenUrl CrossRef PubMed Web of Science 2. ↵ Levine , M. , and Tjian , R . ( 2003 ). Transcription regulation and animal diversity . Nature 424 , 147 – 151 . doi: 10.1038/nature01763 . OpenUrl CrossRef PubMed Web of Science 3. ↵ Kadonaga , J.T . ( 2004 ). Regulation of RNA Polymerase II Transcription by Sequence-Specific DNA Binding Factors . Cell 116 , 247 – 257 . doi: 10.1016/S0092-8674(03)01078-X . OpenUrl CrossRef PubMed Web of Science 4. Young , R.A . ( 2011 ). Control of the embryonic stem cell state . Cell 144 , 940 – 954 . OpenUrl CrossRef PubMed Web of Science 5. ↵ Farley , E.K. , Olson , K.M. , and Levine , M.S . ( 2015 ). Regulatory Principles Governing Tissue Specificity of Developmental Enhancers . Cold Spring Harb Symp Quant Biol 80 , 27 – 32 . doi: 10.1101/sqb.2015.80.027227 . OpenUrl Abstract / FREE Full Text 6. ↵ Weirauch , M.T. , Yang , A. , Albu , M. , Cote , A.G. , Montenegro-Montero , A. , Drewe , P. , Najafabadi , H.S. , Lambert , S.A. , Mann , I. , Cook , K. , et al. ( 2014 ). Determination and inference of eukaryotic transcription factor sequence specificity . Cell 158 , 1431 – 1443 . OpenUrl CrossRef PubMed Web of Science 7. ↵ Qin , Q. , Fan , J. , Zheng , R. , Wan , C. , Mei , S. , Wu , Q. , Sun , H. , Brown , M. , Zhang , J. , Meyer , C.A. , et al. ( 2020 ). Lisa: inferring transcriptional regulators through integrative modeling of public chromatin accessibility and ChIP-seq data . Genome Biology 21 , 32 . doi: 10.1186/s13059-020-1934-6 . OpenUrl CrossRef PubMed 8. Lynch , A.W. , Theodoris , C.V. , Long , H.W. , Brown , M. , Liu , X.S. , and Meyer , C.A . ( 2022 ). MIRA: joint regulatory modeling of multimodal expression and chromatin accessibility in single cells . Nat Methods 19 , 1097 – 1108 . doi: 10.1038/s41592-022-01595-z . OpenUrl CrossRef PubMed 9. Wang , L. , Trasanidis , N. , Wu , T. , Dong , G. , Hu , M. , Bauer , D.E. , and Pinello , L . ( 2023 ). Dictys: dynamic gene regulatory network dissects developmental continuum with single-cell multiomics . Nat Methods 20 , 1368 – 1378 . doi: 10.1038/s41592-023-01971-3 . OpenUrl CrossRef 10. ↵ Kamimoto , K. , Stringa , B. , Hoffmann , C.M. , Jindal , K. , Solnica-Krezel , L. , and Morris , S.A . ( 2023 ). Dissecting cell identity via network inference and in silico gene perturbation . Nature 614 , 742 – 751 . doi: 10.1038/s41586-022-05688-9 . OpenUrl CrossRef PubMed 11. ↵ Brass , A.L. , Zhu , A.Q. , and Singh , H . ( 1999 ). Assembly requirements of PU.1-Pip (IRF-4) activator complexes: inhibiting function in vivo using fused dimers . EMBO J 18 , 977 – 991 . OpenUrl Abstract / FREE Full Text 12. ↵ Glasmacher , E. , Agrawal , S. , Chang , A.B. , Murphy , T.L. , Zeng , W. , Vander Lugt , B. , Khan , A.A. , Ciofani , M. , Spooner , C.J. , Rutz , S. , et al. ( 2012 ). A genomic regulatory element that directs assembly and function of immune-specific AP-1-IRF complexes . Science 338 , 975 – 980 . doi: 10.1126/science.1228309 . OpenUrl Abstract / FREE Full Text 13. ↵ Jain , J. , McCaffrey , P.G. , Valge-Archer , V.E. , and Rao , A . ( 1992 ). Nuclear factor of activated T cells contains Fos and Jun . Nature 356 , 801 – 804 . doi: 10.1038/356801a0 . OpenUrl CrossRef PubMed 14. ↵ Chen , L. , Glover , J.N.M. , Hogan , P.G. , Rao , A. , and Harrison , S.C . ( 1998 ). Structure of the DNA-binding domains from NFAT, Fos and Jun bound specifically to DNA . Nature 392 , 42 – 48 . doi: 10.1038/32100 . OpenUrl CrossRef PubMed Web of Science 15. ↵ Ng , C.K.L. , Li , N.X. , Chee , S. , Prabhakar , S. , Kolatkar , P.R. , and Jauch , R . ( 2012 ). Deciphering the Sox-Oct partner code by quantitative cooperativity measurements . Nucleic Acids Research 40 , 4933 – 4941 . doi: 10.1093/nar/gks153 . OpenUrl CrossRef PubMed Web of Science 16. ↵ Jolma , A. , Yin , Y. , Nitta , K.R. , Dave , K. , Popov , A. , Taipale , M. , Enge , M. , Kivioja , T. , Morgunova , E. , and Taipale , J . ( 2015 ). DNA-dependent formation of transcription factor pairs alters their binding specificity . Nature 527 , 384 – 388 . doi: 10.1038/nature15518 . OpenUrl CrossRef PubMed 17. ↵ Datta , V. , Siddharthan , R. , and Krishna , S . ( 2018 ). Detection of cooperatively bound transcription factor pairs using ChIP-seq peak intensities and expectation maximization . PLoS One 13 , e0199771 . OpenUrl CrossRef PubMed 18. ↵ Avsec , Z. , Weilert , M. , Shrikumar , A. , Krueger , S. , Alexandari , A. , Dalal , K. , Fropf , R. , McAnany , C. , Gagneur , J. , Kundaje , A. , et al. ( 2021 ). Base-resolution models of transcription-factor binding reveal soft motif syntax . Nat Genet 53 , 354 – 366 . OpenUrl CrossRef PubMed 19. Avsec , Ž. , Agarwal , V. , Visentin , D. , Ledsam , J.R. , Grabska-Barwinska , A. , Taylor , K.R. , Assael , Y. , Jumper , J. , Kohli , P. , and Kelley , D.R. ( 2021 ). Effective gene expression prediction from sequence by integrating long-range interactions . Nat Methods 18 , 1196 – 1203 . doi: 10.1038/s41592-021-01252-x . OpenUrl CrossRef PubMed 20. ↵ Pampari , A. , Shcherbina , A. , Kvon , E.Z. , Kosicki , M. , Nair , S. , Kundu , S. , Kathiria , A.S. , Risca , V.I. , Kuningas , K. , Alasoo , K. , et al. ( 2025 ). ChromBPNet: bias factorized, base-resolution deep learning models of chromatin accessibility reveal cis-regulatory sequence syntax, transcription factor footprints and regulatory variants . Preprint at bioRxiv , doi: 10.1101/2024.12.25.630221 10.1101/2024.12.25.630221. OpenUrl Abstract / FREE Full Text 21. ↵ Tewhey , R. , Kotliar , D. , Park , D.S. , Liu , B. , Winnicki , S. , Reilly , S.K. , Andersen , K.G. , Mikkelsen , T.S. , Lander , E.S. , Schaffner , S.F. , et al. ( 2016 ). Direct Identification of Hundreds of Expression-Modulating Variants using a Multiplexed Reporter Assay . Cell 165 , 1519 – 1529 . doi: 10.1016/j.cell.2016.04.027 . OpenUrl CrossRef PubMed 22. ↵ Yoshida , H. , Lareau , C.A. , Ramirez , R.N. , Rose , S.A. , Maier , B. , Wroblewska , A. , Desland , F. , Chudnovskiy , A. , Mortha , A. , Dominguez , C. , et al. ( 2019 ). The cis-Regulatory Atlas of the Mouse Immune System . Cell 176 , 897 – 912 e20. doi: 10.1016/j.cell.2018.12.036 . OpenUrl CrossRef PubMed 23. ↵ Dunham , I. , Kundaje , A. , Aldred , S.F. , Collins , P.J. , Davis , C.A. , Doyle , F. , Epstein , C.B. , Frietze , S. , Harrow , J. , Kaul , R. , et al. ( 2012 ). An integrated encyclopedia of DNA elements in the human genome . Nature 489 , 57 – 74 . doi: 10.1038/nature11247 . OpenUrl CrossRef PubMed Web of Science 24. ↵ Chaudhri , V.K. , Dienger-Stambaugh , K. , Wu , Z. , Shrestha , M. , and Singh , H . ( 2020 ). Charting the cis-regulome of activated B cells by coupling structural and functional genomics . Nat Immunol 21 , 210 – 220 . OpenUrl CrossRef PubMed 25. ↵ Ochiai , K. , Maienschein-Cline , M. , Simonetti , G. , Chen , J. , Rosenthal , R. , Brink , R. , Chong , A.S. , Klein , U. , Dinner , A.R. , Singh , H. , et al. ( 2013 ). Transcriptional Regulation of Germinal Center B and Plasma Cell Fates by Dynamical Control of IRF4 . Immunity 38 , 918 – 929 . doi: 10.1016/j.immuni.2013.04.009 . OpenUrl CrossRef PubMed Web of Science 26. ↵ Tamura , T. , Thotakura , P. , Tanaka , T.S. , Ko , M.S.H. , and Ozato , K . ( 2005 ). Identification of target genes and a unique cis element regulated by IRF-8 in developing macrophages . Blood 106 , 1938 – 1947 . doi: 10.1182/blood-2005-01-0080 . OpenUrl Abstract / FREE Full Text 27. ↵ Heinz , S. , Benner , C. , Spann , N. , Bertolino , E. , Lin , Y.C. , Laslo , P. , Cheng , J.X. , Murre , C. , Singh , H. , and Glass , C.K . ( 2010 ). Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities . Mol Cell 38 , 576 – 589 . doi: 10.1016/j.molcel.2010.05.004 . OpenUrl CrossRef PubMed Web of Science 28. ↵ Lambert , S.A. , Jolma , A. , Campitelli , L.F. , Das , P.K. , Yin , Y. , Albu , M. , Chen , X. , Taipale , J. , Hughes , T.R. , and Weirauch , M.T . ( 2018 ). The Human Transcription Factors . Cell 175 , 598 – 599 . doi: 10.1016/j.cell.2018.09.045 . OpenUrl CrossRef PubMed 29. ↵ Spooner , C.J. , Cheng , J.X. , Pujadas , E. , Laslo , P. , and Singh , H . ( 2009 ). A recurrent network involving the transcription factors PU.1 and Gfi1 orchestrates innate and adaptive immune cell fates . Immunity 31 , 576 – 586 . doi: 10.1016/j.immuni.2009.07.011 . OpenUrl CrossRef PubMed Web of Science 30. ↵ Madison , B.J. , Clark , K.A. , Bhachech , N. , Hollenhorst , P.C. , Graves , B.J. , and Currie , S.L . ( 2018 ). Electrostatic repulsion causes anticooperative DNA binding between tumor suppressor ETS transcription factors and JUN-FOS at composite DNA sites . J Biol Chem 293 , 18624 – 18635 . doi: 10.1074/jbc.RA118.003352 . OpenUrl Abstract / FREE Full Text 31. ↵ Mei , S. , Qin , Q. , Wu , Q. , Sun , H. , Zheng , R. , Zang , C. , Zhu , M. , Wu , J. , Shi , X. , Taing , L. , et al. ( 2017 ). Cistrome Data Browser: a data portal for ChIP-Seq and chromatin accessibility data in human and mouse . Nucleic Acids Research 45 , D658 – D662 . doi: 10.1093/nar/gkw983 . OpenUrl CrossRef PubMed 32. Zheng , R. , Wan , C. , Mei , S. , Qin , Q. , Wu , Q. , Sun , H. , Chen , C.-H. , Brown , M. , Zhang , X. , Meyer , C.A. , et al. ( 2019 ). Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis . Nucleic Acids Research 47 , D729 – D735 . doi: 10.1093/nar/gky1094 . OpenUrl CrossRef PubMed 33. ↵ Taing , L. , Dandawate , A. , L’Yi , S. , Gehlenborg , N. , Brown , M. , and Meyer , C.A . ( 2024 ). Cistrome Data Browser: integrated search, analysis and visualization of chromatin data . Nucleic Acids Research 52 , D61 – D66 . doi: 10.1093/nar/gkad1069 . OpenUrl CrossRef 34. ↵ Love , M.I. , Huber , W. , and Anders , S . ( 2014 ). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biology 15 , 550 . doi: 10.1186/s13059-014-0550-8 . OpenUrl CrossRef PubMed 35. ↵ Orkin , S.H. , and Zon , L.I . ( 2008 ). Hematopoiesis: An Evolving Paradigm for Stem Cell Biology . Cell 132 , 631 – 644 . doi: 10.1016/j.cell.2008.01.025 . OpenUrl CrossRef PubMed Web of Science 36. ↵ Benner , C. , Konovalov , S. , Mackintosh , C. , Hutt , K.R. , Stunnenberg , R. , and Garcia-Bassets , I . ( 2013 ). Decoding a Signature-Based Model of Transcription Cofactor Recruitment Dictated by Cardinal Cis-Regulatory Elements in Proximal Promoter Regions . PLOS Genetics 9 , e1003906 . doi: 10.1371/journal.pgen.1003906 . OpenUrl CrossRef PubMed 37. ↵ Suske , G . ( 2017 ). NF-Y and SP transcription factors — New insights in a long-standing liaison . Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms 1860 , 590 – 597 . doi: 10.1016/j.bbagrm.2016.08.011 . OpenUrl CrossRef 38. ↵ Hess , J. , Angel , P. , and Schorpp-Kistner , M . ( 2004 ). AP-1 subunits: quarrel and harmony among siblings . Journal of Cell Science 117 , 5965 – 5973 . doi: 10.1242/jcs.01589 . OpenUrl Abstract / FREE Full Text 39. ↵ Presnell , J.S. , Schnitzler , C.E. , and Browne , W.E. KLF/SP Transcription Factor Family Evolution: Expansion, Diversification, and Innovation in Eukaryotes . 40. ↵ Vaswani , A. , Shazeer , N. , Parmar , N. , Uszkoreit , J. , Jones , L. , Gomez , A.N. , Kaiser , Ł. ukasz , and Polosukhin , I. ( 2017 ). Attention is All you Need . In Advances in Neural Information Processing Systems ( Curran Associates, Inc .). 41. ↵ Su , J. , Lu , Y. , Pan , S. , Murtadha , A. , Wen , B. , and Liu , Y . ( 2023 ). RoFormer: Enhanced Transformer with Rotary Position Embedding . Preprint at arXiv , doi: 10.48550/arXiv.2104.09864 10.48550/arXiv.2104.09864. OpenUrl CrossRef 42. ↵ Core , L.J. , Martins , A.L. , Danko , C.G. , Waters , C.T. , Siepel , A. , and Lis , J.T . ( 2014 ). Analysis of nascent RNA identifies a unified architecture of initiation regions at mammalian promoters and enhancers . Nat Genet 46 , 1311 – 1320 . doi: 10.1038/ng.3142 . OpenUrl CrossRef PubMed 43. ↵ Rauluseviciute , I. , Riudavets-Puig , R. , Blanc-Mathieu , R. , Castro-Mondragon , J.A. , Ferenc , K. , Kumar , V. , Lemma , R.B. , Lucas , J. , Chèneby , J. , Baranasic , D. , et al. ( 2024 ). JASPAR 2024: 20th anniversary of the open-access database of transcription factor binding profiles . Nucleic Acids Research 52 , D174 – D182 . doi: 10.1093/nar/gkad1059 . OpenUrl CrossRef PubMed 44. ↵ Ong , C.-T. , and Corces , V.G . ( 2014 ). CTCF: an architectural protein bridging genome topology and function . Nat Rev Genet 15 , 234 – 246 . doi: 10.1038/nrg3663 . OpenUrl CrossRef PubMed 45. ↵ Tanaka , H. , Takizawa , Y. , Takaku , M. , Kato , D. , Kumagawa , Y. , Grimm , S.A. , Wade , P.A. , and Kurumizaka , H . ( 2020 ). Interaction of the pioneer transcription factor GATA3 with nucleosomes . Nat Commun 11 , 4136 . doi: 10.1038/s41467-020-17959-y . OpenUrl CrossRef PubMed 46. ↵ Yu , C.-H. , Yang , S.-Q. , Zhang , Y.-J. , Rong , L. , and Yi , Z.-C . ( 2023 ). The role of GATA switch in benzene metabolite hydroquinone inhibiting erythroid differentiation in K562 cells . Arch Toxicol 97 , 2169 – 2181 . doi: 10.1007/s00204-023-03541-0 . OpenUrl CrossRef PubMed 47. ↵ King , H.W. , and Klose , R.J . ( 2017 ). The pioneer factor OCT4 requires the chromatin remodeller BRG1 to support gene regulatory element function in mouse embryonic stem cells . eLife 6 , e22631 . doi: 10.7554/eLife.22631 . OpenUrl CrossRef PubMed 48. ↵ An integrated encyclopedia of DNA elements in the human genome | Nature https://www-nature-com.pitt.idm.oclc.org/articles/nature11247 . 49. ↵ Cano , A. , Pérez-Moreno , M.A. , Rodrigo , I. , Locascio , A. , Blanco , M.J. , del Barrio , M.G. , Portillo , F. , and Nieto , M.A . ( 2000 ). The transcription factor Snail controls epithelial– mesenchymal transitions by repressing E-cadherin expression . Nat Cell Biol 2 , 76 – 83 . doi: 10.1038/35000025 . OpenUrl CrossRef PubMed Web of Science 50. ↵ Lin , Y. , Wu , Y. , Li , J. , Dong , C. , Ye , X. , Chi , Y. , Evers , B.M. , and Zhou , B.P . ( 2010 ). The SNAG domain of Snail1 functions as a molecular hook for recruiting lysine-specific demethylase 1 . The EMBO Journal 29 , 1803 – 1816 . doi: 10.1038/emboj.2010.63 . OpenUrl Abstract / FREE Full Text 51. ↵ Fu , Y. , Kelly , J.A. , Gopalakrishnan , J. , Pelikan , R.C. , Tessneer , K.L. , Pasula , S. , Grundahl , K. , Murphy , D.A. , and Gaffney , P.M . ( 2024 ). Massively parallel reporter assay confirms regulatory potential of hQTLs and reveals important variants in lupus and other autoimmune diseases . HGGADVANCE 5 . doi: 10.1016/j.xhgg.2024.100279 . OpenUrl CrossRef 52. ↵ Xie , Z. , Sokolov , I. , Osmala , M. , Yue , X. , Bower , G. , Pett , J.P. , Chen , Y. , Wang , K. , Cavga , A.D. , Popov , A. , et al. ( 2025 ). DNA-guided transcription factor interactions extend human gene regulatory code . Nature 641 , 1329 – 1338 . doi: 10.1038/s41586-025-08844-z . OpenUrl CrossRef 53. ↵ Merino , F. , Ng , C.K.L. , Veerapandian , V. , Schöler , H.R. , Jauch , R. , and Cojocaru , V . ( 2014 ). Structural Basis for the SOX-Dependent Genomic Redistribution of OCT4 in Stem Cell Differentiation . Structure 22 , 1274 – 1286 . doi: 10.1016/j.str.2014.06.014 . OpenUrl CrossRef 54. ↵ Kim , S. , Morgunova , E. , Naqvi , S. , Goovaerts , S. , Bader , M. , Koska , M. , Popov , A. , Luong , C. , Pogson , A. , Swigut , T. , et al. ( 2024 ). DNA-guided transcription factor cooperativity shapes face and limb mesenchyme . Cell 187 , 692 – 711 .e26. doi: 10.1016/j.cell.2023.12.032 . OpenUrl CrossRef PubMed 55. ↵ Chaudhri , V.K. , and Singh , H . ( 2021 ). Diverse digital and fuzzy composite transcriptional elements are prevalent features of mammalian cis-regulomes . Preprint at bioRxiv , doi: 10.1101/2021.11.26.470154 10.1101/2021.11.26.470154. OpenUrl Abstract / FREE Full Text 56. ↵ Inge , M.M. , Miller , R. , Hook , H. , Bray , D. , Keenan , J.L. , Zhao , R. , Gilmore , T.D. , and Siggers , T . ( 2024 ). Rapid profiling of transcription factor–cofactor interaction networks reveals principles of epigenetic regulation . Nucleic Acids Research 52 , 10276 – 10296 . doi: 10.1093/nar/gkae706 . OpenUrl CrossRef PubMed 57. ↵ Quinlan , A.R. , and Hall , I.M . ( 2010 ). BEDTools: a flexible suite of utilities for comparing genomic features . Bioinformatics 26 , 841 – 842 . doi: 10.1093/bioinformatics/btq033 . OpenUrl CrossRef PubMed Web of Science 58. ↵ Dale , R.K. , Pedersen , B.S. , and Quinlan , A.R . ( 2011 ). Pybedtools: a flexible Python library for manipulating genomic datasets and annotations . Bioinformatics 27 , 3423 – 3424 . doi: 10.1093/bioinformatics/btr539 . OpenUrl CrossRef PubMed Web of Science 59. ↵ Gupta , S. , Stamatoyannopoulos , J.A. , Bailey , T.L. , and Noble , W.S . ( 2007 ). Quantifying similarity between motifs . Genome Biology 8 , R24 . doi: 10.1186/gb-2007-8-2-r24 . OpenUrl CrossRef PubMed 60. ↵ Mouri , K. , Guo , M.H. , de Boer , C.G. , Lissner , M.M. , Harten , I.A. , Newby , G.A. , DeBerg , H.A. , Platt , W.F. , Gentili , M. , Liu , D.R. , et al. ( 2022 ). Prioritization of autoimmune disease-associated genetic variants that perturb regulatory element activity in T cells . Nat Genet 54 , 603 – 612 . doi: 10.1038/s41588-022-01056-5 . OpenUrl CrossRef PubMed 61. ↵ Zhang , Y. , Liu , T. , Meyer , C.A. , Eeckhoute , J. , Johnson , D.S. , Bernstein , B.E. , Nusbaum , C. , Myers , R.M. , Brown , M. , Li , W. , et al. ( 2008 ). Model-based Analysis of ChIP-Seq (MACS) . Genome Biology 9 , R137 . doi: 10.1186/gb-2008-9-9-r137 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted July 15, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Expanding the DNA Motif Lexicon of the Transcriptional Regulatory Code Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Expanding the DNA Motif Lexicon of the Transcriptional Regulatory Code Jingyu Fan , Virendra K. Chaudhri , Deepa Bisht , Nicholas A. Pease , David Reed Hall , Peter Gerges , Vivian Fang Yi , Susan Kales , Ching-Huang Ho , Jishnu Das , John P. Ray , Ryan Tewhey , Clifford A. Meyer , Nidhi Sahni , Harinder Singh bioRxiv 2025.07.09.662874; doi: https://doi.org/10.1101/2025.07.09.662874 Share This Article: Copy Citation Tools Expanding the DNA Motif Lexicon of the Transcriptional Regulatory Code Jingyu Fan , Virendra K. Chaudhri , Deepa Bisht , Nicholas A. Pease , David Reed Hall , Peter Gerges , Vivian Fang Yi , Susan Kales , Ching-Huang Ho , Jishnu Das , John P. Ray , Ryan Tewhey , Clifford A. Meyer , Nidhi Sahni , Harinder Singh bioRxiv 2025.07.09.662874; doi: https://doi.org/10.1101/2025.07.09.662874 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7619) Biochemistry (17639) Bioengineering (13865) Bioinformatics (41854) Biophysics (21405) Cancer Biology (18544) Cell Biology (25432) Clinical Trials (138) Developmental Biology (13356) Ecology (19863) Epidemiology (2067) Evolutionary Biology (24287) Genetics (15587) Genomics (22465) Immunology (17701) Microbiology (40300) Molecular Biology (17142) Neuroscience (88440) Paleontology (666) Pathology (2825) Pharmacology and Toxicology (4814) Physiology (7633) Plant Biology (15107) Scientific Communication and Education (2042) Synthetic Biology (4285) Systems Biology (9809) Zoology (2268)