Sequence determinants of human gene regulatory elements

1The temporal and spatial pattern of gene expression is encoded in the DNA sequence; this information is read and interpreted by TFs, which recognize and bind specific short DNA sequence motifs. Major efforts have been undertaken to determine the DNA-binding specificities of TFs in vitro– and map their binding positions in vivo,. TFs regulate gene expression by binding to distal enhancer elements and to promoters located close to the TSS,. Both enhancers and promoters are characterized by RNA transcription, the presence of open chromatin and histone H3 lysine 27 acetylation (H3K27ac). In addition, promoters and enhancers are preferentially marked by histone H3 lysine 4 trimethylation and monomethylation, respectively. Although these features can be mapped genome-wide in a high-throughput manner, they are correlative in nature and do not establish that an element can act as an enhancer, increasing expression from a promoter irrespective of position and orientation. To more directly measure enhancer activity, MPRAs have been developed to study the activity of yeast, Drosophila and human gene regulatory elements on a genome-wide scale. However, unbiased discovery of sequence determinants of human gene expression using only genomic sequences is made difficult by the fact that the genome is repetitive and has evolved to perform multiple functions. Furthermore, the human genome is too short to even encode all combinations, orientations and spacings of approximately 1,639 human TFs in multiple independent sequence contexts. Thus, despite the vast amount of information generated by genome-scale experiments, most sequence determinants that drive the activity of human enhancers and promoters, and the interactions between them, remain unknown.

2To systematically characterize the sequence determinants of human gene regulatory element activity, we developed a set of four MPRA libraries that cover more than 100 times the sequence space of the human genome (Fig. 1a and Methods). The libraries are based on the self-transcribing active regulatory region sequencing (STARR-seq) design, in which putative enhancers are cloned downstream of a minimal promoter to the 3′ untranslated region (UTR) of a reporter gene. The constructs are transfected to cells, and the enhancer activity of the UTRs are then determined by RNA sequencing (RNA-seq) (Fig. 1b). Three libraries were designed to measure enhancer activities of combinations of known TF binding motifs embedded within two different 49-bp sequence contexts, ~500-bp fragments of genomic DNA and synthetic random 170-bp sequences; a fourth library was designed to measure both enhancer and promoter activities of synthetic random 150-bp sequences. Sequencing of the input libraries revealed their ultrahigh complexity, reaching billions of unique fragments (Supplementary Fig. 1a,b and Methods).

3To measure the enhancer activity of the known TF consensus sequences, we transfected GP5d colon carcinoma cells with the motif libraries (Fig. 1a, i) and purified total poly(A)+ RNA from the transfected cells. The synthetic motif sequences that were transcribed to RNA were recovered using reverse-transcription PCR, and the abundance of each sequence was then quantified by massively parallel sequencing (Methods). Comparison of the median activities of the individual TF consensus sequences revealed that several TFs had enhancer activity in GP5d cells (Fig. 1c, Extended Data Fig. 1a–c, Supplementary Note and Supplementary Table 5). The consensus sequence corresponding to the p53 protein family (p53, p63 and p73) displayed the strongest enhancer activity in this assay, suggesting that there is constitutive p53 activity in GP5d cells (Fig. 1c and Supplementary Fig. 1c–f). As the library contained each single-base substitution to the consensus sequences, we were able to generate activity position weight matrices (PWMs) for the motifs. For 11 motifs, the activity PWMs were highly similar to that of the motifs derived from an in vitro binding-specificity assay (high-throughput systematic evolution of ligands by exponential enrichment, HT-SELEX; Fig. 1d; Extended Data Fig. 1d), indicating that the measured enhancer activity originated from the TFs that bound to the motifs, demonstrating that the assay can be used to faithfully measure TF activities in cells.

4Comparison of enhancer activities of motifs with the DNA-binding activities of respective TFs measured from the nuclear extract of GP5d cells by an active TF identification (ATI) assay (Methods) revealed that the transcriptional and DNA-binding activities were only weakly correlated (log2 fold change, Pearson R = 0.032; Fig. 2a and Supplementary Note). These results suggest that largely distinct sets of TFs display strong enhancer activity and strong DNA-binding activity in a cell.

5Apart from simple cellular alarm signals, most transcription is thought to require combinatorial action of many TFs–. Consistent with this, we observed that the average activity of all consensus sequences was very low, and for the majority of the TFs, the enhancer activity increased as a function of the number of consensus sequences (Extended Data Fig. 1e, red horizontal lines). Conversely, for the TFs that can activate transcription alone (e.g., p53 and IRF), two consensus sequences had lower activity than that predicted from an additive model (Fig. 2b, red dotted line), presumably due to saturation of both the occupancy and the downstream transcriptional activation. For TFs with intermediate activity levels (e.g., NFAT and YY), activity increased linearly rather than synergistically as a function of the number of binding sites (Extended Data Fig. 1e). The simplest model consistent with these observations is that human enhancer activation requires overcoming a repressive activity, after which activation is linear (additive) until it starts to saturate as it approaches a maximum level.

6To discover sequence features that contribute to human enhancer activity in an unbiased manner, we used extremely complex random enhancer library (Fig. 1a, iii) in GP5d cells. Motif mapping across replicate experiments indicated that motif activities were highly reproducible (Pearson R = 0.963; Extended Data Fig. 2a), displaying additivity and saturation similar to that observed with the motif library (Extended Data Fig. 2b,c). Enrichment of motifs corresponding to known TFs specific to colon cancer and intestinal lineage, such as TCF/LEF, GRHL and HNF4, was clearly observed (Extended Data Fig. 2b). De novo motif mining identified 22 TF motifs; most of these were for individual TFs or conventional heterodimers (Extended Data Fig. 2d and Supplementary Fig. 2a). One strong de novo ETS-bZIP composite motif was also identified, revealing a potential role for ETS-bZIP combinatorial control in colon cancer cells (Extended Data Fig. 2d). Analysis of spacing between motif matches identified few significantly overrepresented spacing preferences for motif pairs such as p53 family–p53 family, GRHL–ETS class I and GRHL–ATF6 (Fig. 2c); weak overall preference for motifs that were relatively close (<50 bp) was also observed (Extended Data Fig. 2e,f; Supplementary Note). These results suggest that TF grammar is strong at the level of heterodimers (analogous to ‘compound words’) but relatively weak at the level of specific combinations and spacing and orientation preferences between TFs (‘sentences’).

7To determine sequence features present in the de novo enhancers, we used machine learning classifiers. First, we determined the importance of known motif features using a logistic regression model (Methods); we found that only a handful of known TF binding motifs are needed for optimal classification, as only 26 out of 19,150 features had regression coefficient absolute values within 10% of the largest regression coefficient (Fig. 2d, Extended Data Fig. 3a and Methods). The most predictive features were motifs for known TFs important for tumorigenesis and colon development (Fig. 2d). These motifs were enriched, suggesting that the corresponding TFs act as transcriptional activators. The interactions between the motifs were largely additive, as specific pairwise combinations did not add substantially to the predictive power.

8Next, to identify possible novel sequence features that would allow more optimal classification, we trained a convolutional neural network (CNN)–based classifier similar to DeepBind on the sequence data. This method is capable of learning the sequence motifs, their combinations and their relative weights de novo. The CNN classifier performed substantially better than logistic regression using the same training, validation and test sets (11% increase in the area under the precision-recall curve (AUprc); Extended Data Fig. 3a,b and Methods). Analysis of the CNN classifier revealed that it had learned motif features similar to those identified by logistic regression (Fig. 2e, Extended Data Figs. 3c and 4, Supplementary Fig. 2 and Supplementary Note). In conclusion, these results indicate that individual TFs contribute to de novo enhancer function mostly without specific interactions between them.

9To determine whether enhancers are cell-type specific, we used the random enhancer library (Fig. 1a, iii) to identify sequence features important for enhancer activity in HepG2 hepatocellular carcinoma cells. Comparison of enhancer motifs between the GP5d and HepG2 cells revealed that most motifs had similar enhancer activity across the cell lines (Pearson R = 0.78; Fig. 2f). The motifs with differential activity corresponded to lineage-determining TFs (GRHL in GP5d) and TFs important for tissue function (TEAD and ATF4:CEBPB in HepG2 cells). Importantly, the lineage-determining factors showed clear differential expression between the two cell types (Fig. 3a), indicating that activities of individual TFs are commonly affected by their expression level, although the overall correlation between motif activity and expression of corresponding TF family was weak (Extended Data Fig. 1a–c). These results show that the transcriptional landscape of a cell is dominated by cell-biological or ‘housekeeping’ TFs that show comparable activity across cell types and that the largest differences of motif activity between cell types are driven by TFs important for lineage specification.

10To determine how sequence features combine to generate functional genomic enhancers, we assayed genomic enhancer activity in GP5d and HepG2 cells at ~1.5-bp resolution (Fig. 1a, ii, Supplementary Fig. 1a,b, Extended Data Fig. 5a,b and Methods) before and after methylation of the library (Fig. 3b). To determine the role of TP53 in enhancer activity, we performed similar experiments in TP53−/− GP5d cells. The signal was highly specific, as indicated by the fact that loss of TP53 resulted in loss of most enhancer peaks containing its motif (Fig. 3c,d). However, despite being the strongest activator in both cell types (see Fig. 2f), TP53 contributed to a relatively small proportion of the overall enhancer activity in both GP5d and HepG2 cells; only 16% and 4.9% of the genomic STARR-seq peaks overlapped with TP53 chromatin immunoprecipitation sequencing (ChIP-seq) peaks (Extended Data Fig. 5c,d and Supplementary Note). Analysis of the methylated libraries revealed that activities of methylated genomic elements were consistent with the known effect of methylation on TF DNA binding (Fig. 3d and Yin et al.). Consistent with the known association between accessible chromatin, TF binding and enhancer activity, the STARR-seq peaks overlapped significantly with chromatin accessibility; specifically, 30% of the STARR-seq peaks in GP5d and 27% in HepG2 cells overlap with assay for transposase-accessible chromatin using sequencing (ATAC-seq) peaks in the same cell types (Figs. 3e and 4a–c and Extended Data Figs. 5c and 6a,b). Furthermore, ATAC-seq peaks could be predicted by a CNN trained using genomic or random synthetic STARR-seq sequences (AUprc 0.80 and 0.71, respectively; Extended Data Fig. 6c), indicating that the sequence features discovered using STARR-seq correspond partially to the features that are associated with open chromatin in vivo.

11We next used the differential signals for chromatin accessibility (ATAC+ or ATAC−) and classical enhancer activity (STARR+ or STARR−) for defining different classes of gene regulatory elements along with ChIP-seq data for individual TFs; the histone marks H3K27ac, H3K9me3 and H3K27me3; and the structural chromatin protein CTCF (Fig. 4a,c). This analysis revealed six classes of elements: (1) closed-chromatin enhancers (STARR+ and ATAC−), (2) cryptic enhancers (silenced STARR+ and ATAC− regions), (3) promoters (ATAC+ and STARR+/−), (4) chromatin-dependent enhancers (STARR−/low and ATAC+ with active histone mark H3K27ac), (5) structural chromatin elements (STARR−, ATAC+ and CTCF+) and (6) classical enhancers (STARR+ and ATAC+).

12Analysis of the methylated genomic elements revealed that the cryptic enhancers were not silenced by methylation (Fig. 3b). Instead, they were inactive due to the presence of either H3K27me3 (polycomb; ‘poised’ enhancer) or H3K9me3/HP1 repressive chromatin marks. The three other types of enhancers (closed chromatin, chromatin dependent and classical) appeared active based on the fact that inclusion of the corresponding features improved prediction of differential gene expression between GP5d and HepG2 cells (Supplementary Table 6 and Methods). Analysis of ChIP-seq peaks and motifs present in the different classes of elements revealed that classical and closed-chromatin enhancers bound to TFs and contained motifs that were similar to those that were found in active elements selected from random sequences (Extended Data Fig. 7a and Fig. 2f). Classical enhancers were preferentially bound by TFs with strong activator domains (e.g., FOS and JUN), whereas chromatin-dependent enhancers displayed relatively weak preference for HLF and FOXA motifs, and both types of enhancers were bound by HNF4A (Fig. 3e and Extended Data Fig. 7b). These results indicate that cells contain three distinct classes of enhancers (Supplementary Note): (1) classical enhancers that overlap with open chromatin and transactivate a heterologous promoter regardless of position or orientation; (2) chromatin-dependent enhancers that cannot be effectively detected using STARR-seq (see also Inoue et al. ) and have strong signal for open chromatin and the activating histone mark H3K27ac; and (3) closed-chromatin enhancers whose detection requires STARR-seq, as these elements are not strongly enriched for chromatin marks associated with enhancer activity.

13Consistent with few TFs determining the overall transcriptional landscape of a cell, the genomic STARR-seq peaks were enriched for relatively few motifs (Extended Data Fig. 2d). The motifs themselves were similar to known monomeric, dimeric and composite TF motifs determined using HT-SELEX and consecutive affinity-purification systematic evolution of ligands by exponential enrichment (CAP-SELEX) (Extended Data Fig. 2d). The motifs discovered from genomic and random enhancers were also largely similar (Extended Data Fig. 2d). The main difference was the enrichment of pioneer factor motifs such as GATA and SOX in genomic fragments (Extended Data Fig. 7c); these motifs may be specifically associated with classical genomic enhancers because of the ability of the corresponding TFs to displace nucleosomes and/or open higher-order chromatin. Many discovered genomic STARR-seq motifs also displayed strong DNA-binding activity in an ATI assay (Extended Data Fig. 2d), indicating that strong DNA binders are important for in vivo enhancer activity, potentially because they are capable of opening chromatin. In summary, the sequence features of classical genomic enhancers are highly similar to those enriched from random sequence; these motifs define the classical enhancer activity of a cell. In addition to this activity, additional chromatin-dependent enhancers confer tissue specificity to genes; these elements are characterized by motifs for TFs that have lower transactivation activity, suggesting that these TFs act via chromatin to facilitate the activity of promoters and associated classical enhancers. Consistent with this view, the strongest cellular enhancers, superenhancers, typically consist of arrays of chromatin-dependent elements associated with a classical enhancer (Fig. 4b and Extended Data Fig. 7d).

14To identify sequence determinants of human promoter activity, we assayed the activity of the binary STARR-seq library consisting of random sequences placed in the position of both the promoter and the enhancer (Fig. 1a, iv). For this analysis, we used two tumor cell lines (GP5d and HepG2; endodermal origin) and an untransformed cell line derived from retinal pigment epithelium (RPE1; ectodermal origin). Robust promoter activity was observed in all three cell lines from a subset of the random sequences, and motif mapping across replicate experiments in GP5d cells showed that motif activities were highly reproducible (Pearson R = 0.997; Extended Data Fig. 2a). As observed for the motifs at active enhancers, most motifs enriched at promoters were similar in all cell types (Fig. 5a). The motifs that displayed differential activity were linked to lineage determination (e.g., HNF1A) and specialized cell functions (ATF4:CEBP in HepG2 cells; Fig. 5a). Comparison of the active sequences in GP5d cells revealed that many sequence motifs were enriched in both the promoter and enhancer positions (Fig. 5b). However, some specificity in the enrichment was also observed. For example, although p53 and YY motifs were similarly enriched at promoters and enhancers, ETS (promoters vs. enhancers 9.5 versus 2.1 in linear scale) and recently discovered BANP motifs (17.7 versus 2.0) were preferentially, and NRF1 (8.4 versus 1.3) as well as HNF1 motifs (4.3 versus 0.9) almost exclusively enriched at promoters (Fig. 5b,c). No motif enriched only at enhancers, indicating that all motifs with enhancing activity can also act from a proximal position at the promoter (Fig. 5b). Of note, some negative effects were also observed (Fig. 5b), consistent with previously known repressive functions of the corresponding TFs (e.g., OVO-like transcriptional repressor 1 and cut-like homeobox). In summary, these results indicate that human promoters can be enriched from random sequences and that active promoter elements are highly similar among different cell types.

15To evaluate the positioning of the different features relative to the TSS, we first determined the TSS position within the promoters derived from random sequences by recovering the 5′ end of the transcript using a template switch (Fig. 6a), yielding 85,217 unique TSS positions. Alignment of the recovered sequences with respect to the TSS positions (Methods) revealed a relatively high information content feature located at the TSS that corresponded to the classic initiator motif (Fig. 6b). In addition, a clear AT-rich region was observed at the canonical −30 position of the TATA box. However, we did not detect other TSS-proximal motifs that have previously been described (BRE, DPE, MTE, DCE, X-core promoter element and TCT,). The transcript side was characterized by a modest increase in G across a relatively wide region (+10 to +35); this feature is also observed in genomic promoters (Supplementary Fig. 3). To identify interactions between the features, we performed mutual information analysis (Methods). The strongest signal was for short-range interactions located 5ʹ to the TSS, excluding a region just upstream of the TATA box; this signal represents enrichment of individual TF motifs. Two mutually exclusive longer-range interactions were detected: one between the TATA box and the TSS and the other between the TSS and the G-rich downstream sequence (Fig. 6c). This pattern is consistent with the loading of the RNA polymerase II either ‘heel first’ (TFIID) or ‘toe first’ (TFIIH) with respect to the TSS.

16Motif mapping revealed that many TF motifs were also specifically positioned and oriented relative to the TSS (Fig. 6d,e). The strongest positional signals were observed for the TATA box, initiator and YY (YY1). YY1 motifs were mainly enriched on the transcript side (the first C of the CCAT sequence occurring on the minus strand at position +12), oriented in a manner that the YY1 protein can position and orient the RNA polymerase II to direct transcription toward the YY motif (Fig. 6d and Houbaviy et al.). In addition, many TF motifs preferentially enriched close to the TSS (Fig. 6d). On the 5ʹ side, the strongest enrichment occurs close to the TSS, slowly decreasing toward the TATA box; preferential enrichment upstream of the TATA box was also observed for some TFs (e.g., BANP). On the 3ʹ side, the enrichment declines more sharply with very little motif enrichment observed beyond the +20 position from the TSS (Fig. 6d,e). In summary, these results highlight that some, but not all, TFs have positional dependency related to the TSS.

17To determine how well transcription can be predicted based on the de novo promoter sequences, we trained a CNN model (Methods) to predict the TSS positions genome-wide. To test the CNN, we first used it to score wild-type and mutant forms of the TERT promoter, (Methods); the model correctly predicted that known cancer-associated mutations increase the activity of this promoter (Fig. 6f and Extended Data Fig. 8a–d). We next used the CNN to predict the positions of active TSSs in GP5d cells using TSS annotation derived from the Eukaryotic Promoter Database (EPD), and the activity of the TSSs was determined using cap analysis of gene expression (CAGE; Methods). This analysis revealed that promoters enriched from random sequences were more predictive than the genomic sequences themselves; 33% of the positions of unseen genomic TSSs were accurately predicted by the CNN trained on the promoters enriched from random sequences, as opposed to 7% predicted by the CNN trained on the EPD promoters (Fig. 6g and Extended Data Fig. 8e). A mutual information–based analysis of interactions learned by the CNN classifiers (Methods) revealed that the classifiers trained on STARR-seq data learned a stronger position-specific signal than the classifiers trained on the EPD data, which relied more on information present at a relatively short region around the TSS (Extended Data Fig. 8f,g and Methods). These results highlight the power of unbiased interrogation of sequence space that is 100 times larger than that of the human genome.

18The binary STARR-seq approach allows identification of interactions between promoters and enhancers. For this analysis, we counted single motif matches at the promoter and enhancer positions and all pairs of motif matches. When promoters and enhancers were analyzed separately, almost all pairs of TF motifs enriched independently of each other. Strikingly, even across promoters and enhancers, all motifs were independently enriched (Fig. 7a), suggesting that TFs bound to enhancers activate promoters, but in a very nonspecific manner. Some highly enriched TF–TF motif pairs, however, displayed weaker activity than that expected from a model that assumes additive action of the enhancer and promoter (Fig. 7b). In addition, three TF–TF motif pairs displayed stronger transcriptional activity than that expected from independent action of the individual TFs (Fig. 7a and Extended Data Fig. 9a); all three pairs combined a p53 family motif at the promoter with a repressive motif at the enhancer (Fig. 7a). These results are consistent with a model in which enhancer and promoter activities are integrated into total transcriptional activity; the observed saturation is consistent with a strong promoter not needing an enhancer and with a strong enhancer rendering weak and strong promoters equally active.

19Unbiased machine learning analysis also supported a general mechanism of integration of promoter and enhancer activities (Fig. 7c). A CNN classifier using only promoter sequences outperformed a classifier using only enhancer sequences. As expected, combining the promoters with the correct enhancer sequences increased performance substantially. However, permutating the pairings between the promoters and enhancers resulted in similar performance (Fig. 7c), indicating that there was no predictive power in the specific pairing of individual promoters and enhancers. Taken together, our results indicate that the mechanisms that control transcription are very general and that the activities of almost all TFs can independently contribute to transcriptional activity.

20Learning the rules by which DNA sequence determines where and when genes are expressed has proven surprisingly hard, despite the availability of full genome sequences of several mammals, extensive maps of genomic features,, and genome-scale data about TF protein expression levels and TF DNA binding in vitro,,. Direct comparison of activities between TFs has remained difficult, and therefore, we generally lack parameters describing the relative strength of the different sequence features and their interactions—features that are critically important for prediction of transcriptional activity. To address this, we have here defined sequence determinants of human regulatory element activity in an unbiased manner, using an approach in which genomic, designed and random sequences are identified that display promoter or enhancer activity.

21We found that the cellular gene regulatory system is relatively complex, consisting of several distinct kinds of elements. Motif grammar is relatively strong at the level of heterodimers but weaker at the level of spacing and orientation of specific TF motif combinations. In transcriptionally active sequences, precise TF arrangements such as those found in the interferon enhanceosome are rare, with most elements consisting of TFs acting together in a largely additive manner,,–. Our results are consistent with a recent report showing that independent actions of TFs can explain over 92% of the transcriptional activity measured from random yeast promoters. The presence of a weak motif syntax is also consistent with known existence of spacing and orientation preferences of TFs in vitro and at human genomic enhancers and synthetic yeast promoters.

22Our results also show that different cell types have very similar TF activities and that the topology of the gene regulatory network is hierarchical, with few TFs displaying very strong transactivation activity. This is consistent with our previous work showing that relatively few TFs show strong DNA-binding activity in a cell and that many of the strong binders are common to various cell types. Our findings contrast with the known tissue specificity of many putative enhancer elements in vivo. Interestingly, the level of conservation of many endogenous promoters and enhancers appears to be higher than the elements selected in our assay (Extended Data Fig. 9b). The simplest explanation for these two facts is that enhancers in vivo evolve to be specific and that due to the similarity between cells, specificity is more difficult to achieve than activity. Specificity will naturally require specific TF combinations and also fine-tuning using motif number, spacing, orientation and affinity (e.g., Panne et al. and Crocker et al.). Specificity is also required to silence strongly active elements in cell types in which their target protein is not needed due to the substantial fitness cost of protein expression. Further analysis using main cell types representing all three germ layers is needed to determine whether and to what extent differentiated human cell types have retained the regulatory mechanisms that existed in their common unicellular ancestor. Moreover, the contribution of specific TFs to the transcriptional activity in the cell could be further dissected, for example, by testing mutated genomic fragments in MPRAs.

23The original functional definition described enhancers as genetic elements that can activate a promoter from a distance, irrespective of their orientation relative to the TSS. We find here that in addition to these elements, two other types of enhancing gene regulatory elements exist: chromatin-dependent enhancers and closed-chromatin enhancers (Fig. 7d). Chromatin-dependent enhancers are characterized by forkhead motifs, binding of Mediator and p300 protein and a strong signal for H3K27 acetylation. Unlike classical enhancers, chromatin-dependent enhancers do not transactivate a heterologous promoter strongly, most likely due to lack of binding of TFs with strong transactivator domains. Their presence is, however, strongly predictive of tissue-specific gene expression, suggesting that they act to increase gene expression via chromatin modification or structural changes in higher-order chromatin. Several chromatin-dependent enhancers also combine with a single classical enhancer to form superenhancers (see Fig. 4b), indicating that these elements may be required for driving high levels of gene expression from distal promoters. The third element type, closed-chromatin enhancers, are located in regions that show little or no signal for DNase I hypersensitivity or ATAC-seq. They are not silenced by CpG methylation. These elements appear to consist of only a single TF (e.g., p53; see also Peng et al.) or a set of closely bound TFs that fit between or associate directly with well-ordered nucleosomes. The prevalence of both the closed-chromatin enhancers and chromatin-dependent enhancers suggests that they may contribute substantially to regional control of gene expression.

24By using machine learning approaches, we show here that transcriptional activity in human cells can be predicted from sequence features (see also Avsec et al. and Agarwal and Shendure). Interestingly, we found that the promoters enriched from completely random synthetic sequences in a single experimental step are even more predictive of transcriptional activity than the genomic sequences themselves. By analysis of de novo promoters enriched from random sequences, we discovered a G-rich element that interacts with the TSS, potentially positioning RNA polymerase II to the TSS independently of the TATA box. Overall, TF activities could be classified into three groups: TSS position–determining activity (e.g., TATA box and YY), short-range promoting activity (e.g., NRF1) and enhancing activity (many TFs). We did not detect a separate class of distal enhancing activity, suggesting that activities that would allow an enhancer to selectively act at a very long range are likely to be associated with chromatin-dependent enhancers and not classical enhancers,. The three classes of activities detected are not mutually exclusive, suggesting that TFs act at multiple levels and/or scales to regulate transcription (Fig. 7d). For example, YY1 acts as both an enhancing TF and a TSS-determining one, and FOXA motifs are present at both chromatin-dependent and classical enhancer elements. Our results thus indicate that TF motifs are the atomic units of gene expression and should be the ultimate basis of analysis and prediction of genomic elements controlling gene regulatory activity.

25Our random promoter–enhancer design allowed unbiased discovery of features that facilitate interactions between classical enhancers and promoters at a relatively short range. No specific pair of motifs controlling such interactions was found. This, together with the fact that no specific TF that only acts from an enhancer was found, is consistent with a generic and indirect mechanism of action, where the activities of individual TFs bound to an enhancer are aggregated and their total activity then activates the promoter. Molecularly, these results are consistent with mediation of the effect by the least specific type of biochemical interaction, steric hindrance. The simplest mechanism for enhancer action would involve opening of higher-order and local chromatin in such a way that the steric hindrance that prevents large macromolecular complexes such as Mediator or RNA polymerase II from loading to DNA is decreased (Fig. 7d and Maeshima et al.). However, in the highly evolved genomic context, more specific interactions can exist between chromatin-dependent enhancers and particular promoters, as reported in a few cases, such as multichromosome structures that control the expression of the repertoire of olfactory receptor genes or the complex regulatory landscape of HOX genes.

26In summary, we show here that direct experimentation to interrogate transcriptional activities of sequences that represent on aggregate >100 times larger sequence space than that of the human genome can be used to determine mechanisms of action of, and interaction between, gene regulatory elements. The experiments revealed unexpected simplicity of gene regulatory logic. The discovery of the relative simplicity of the interactions, together with the ability to measure transcriptional activities of all TFs in a cell, represents a major step toward achieving the ultimate aim of regulatory genomics: predicting gene expression from a sequence.

27We designed a modified STARR-seq reporter construct pGL4.10-Sasaki-SS (a) based on an earlier published design in the pGL4.10 backbone (Promega, E6651). The sequence between SacI and AfeI was replaced with a sequence containing CG-depleted chicken lens δ1-crystallin gene (Sasaki) promoter, a synthetic intron (pIRESpuro3; Clontech, 631619), an ORF (fusion of Nanoluc-EmGFP), homology arms for library cloning with AgeI and SalI restriction enzyme (RE) sites flanking the ccdB gene, a small 52-bp DNA stuffer (a part of the neomycin resistance cassette) and a 20-bp sequence from the 3ʹ-Illumina adapter for optimally sized final library for Illumina sequencing and the SV40 late poly(A) signal from the pGL3 backbone (Promega, E1751).

28To enable the analysis of CpG methylation on enhancer activity, we designed modified STARR-seq vectors in a CpG-free backbone with Lucia reporter gene (Invivogen, pcpgf-promlc) driven either by the EF1α promoter (b. pCpG-free-EF1α-SS) or the Sasaki promoter (c. pCpG-free-Sasaki-SS-v1) as above. To facilitate the cloning of the synthetic DNA library to the 3ʹ UTR of the reporter gene, the cloning cassette from the pGL4.10-Sasaki-SS vector (a) containing the homology arms with AgeI and SalI RE sites, the 52-bp DNA stuffer and the 20-bp sequence from the 3ʹ-Illumina adapter as above was introduced to the CpG-free vectors using the NheI site.

29Standard Illumina adapters harbor CG dinucleotides, and to make our modified STARR-seq design completely CpG-free, we designed custom adapters for Illumina sequencing (oligos 3 and 4 in Supplementary Table 2). To accommodate the cloning of genomic DNA and random sequence inserts with flanking CpG-free custom adapters, the cloning cassette in CpG-free-Sasaki-SS-v1 was modified by removing the 3ʹ-Illumina adapter and the 52-bp stuffer. In addition, this vector was further improved by replacing the AgeI and SalI RE sites with the AflII and PvuII sites devoid of CG dinucleotides and introducing a DNA stuffer of 1.2 kb between the RE sites to the resulting pCpG-free-Sasaki-SS-v2 vector (d) to unambiguously detect and purify the linearized reporter backbone for downstream cloning.

30For the binary STARR-seq approach in which random sequences were cloned as both promoters and enhancers, the pCpG-free-Sasaki-SS-v2 vector (d) was modified by replacing the Sasaki promoter with a custom CpG-free 5ʹ-adapter sequence and an AgeI RE site and introducing a SalI RE site and a custom CpG-free 3ʹ adapter immediately downstream of the ORF. Moreover, to optimize the random promoter and random enhancer library size for Illumina sequencing, the Lucia reporter gene was replaced by a small 11-amino-acid ORF from Drosophila melanogaster (Dm tal-1A) in the pCpG-free-promoter–enhancer-SS vector (e). The cloned random promoter–random enhancer input library is paired-end sequenced to map the promoter–enhancer pairs, and thus, the random enhancer sequences obtained after sequencing the reporter-specific RNA library can be used to identify the corresponding promoter sequence from the input library. In total, the constant sequence between promoter and enhancer elements is 872 bp in the pCpG-free-Sasaki-SS-v2 construct and 215 bp in the pCpG-free-promoter–enhancer-SS construct.

31The new reporter vectors (a–e) were used in different experiments as follows (their complete sequences are provided in Supplementary Table 1): a, pGL4.10-Sasaki-SS (5,754 bp) was used for experiments with the synthetic motif library shown in Extended Data Fig. 1b; b, pCpG-free-EF1α-SS (3,497 bp) was used for all experiments with the synthetic motif library; c, pCpG-free-Sasaki-SS-v1 (3,388 bp) intermediate plasmid was not used in the experiments; d, pCpG-free-Sasaki-SS-v2 (4,458 bp) was used for all experiments with genomic fragments and random enhancer (N170) sequences; and e, pCpG-free-promoter–enhancer-SS (2,551 bp) was used for all experiments with random promoter (N150)–random enhancer (N150) sequences.

32STARR-seq reporter libraries were generated from rationally designed oligonucleotides harboring TF binding motifs, from fragmented human genomic DNA and from synthetic oligonucleotide with completely random DNA sequences as detailed in the Supplementary Methods. All the oligonucleotides that were used for cloning of the libraries were purchased from Integrated DNA Technologies, and their sequences are provided in Supplementary Table 2.

33The genomic DNA library was methylated using M.SssI (New England Biolabs) for 4 h at 37 °C with the reaction volumes scaled for 62.5 µg plasmid DNA per reaction and inactivated for 20 min at 65 °C, followed by purification and ethanol precipitation of the methylated library.

34The cell lines used in this study were the colon cancer cell line GP5d (Sigma, 95090715), the liver cancer cell line HepG2 (ATCC, HB-8065) and the retinal pigment epithelial cell line hTERT-RPE1 (ATCC, CRL-4000). The cells were maintained in their respective media (GP5d in DMEM, HepG2 in MEM and RPE1 in DMEM/F12) supplemented with 10% fetal bovine serum, 2 nM l-glutamine and 1% penicillin–streptomycin.

35The TP53-null GP5d cell line was generated by CRISPR-Cas9 targeting of exon 4 of the TP53 gene using Alt-R CRISPR-Cas9 from Integrated DNA Technologies. Briefly, annealed sgRNA duplex from crRNA (oligo 12; Supplementary Table 2) and tracrRNA with atto550 were used for ribonucleoprotein complex formation with Cas9-HiFi protein, and the ribonucleoprotein complex was transfected to GP5d cells using CRISPRMAX (Invitrogen). The next day, atto550+ cells were FACS sorted, and single-cell colonies were cultured to produce a clonal TP53-null cell line. The clonal cells lines were screened for TP53 depletion by western blotting, and clones were verified by Sanger sequencing using oligos 13 and 14 (Supplementary Table 2).

36In STARR-seq experiments, 1 µg of each input library DNA was transfected per million cells. For TF motif DNA libraries, a total of 50 and 35 million GP5d cells were transfected for the libraries in the pGL4.10-Sasaki-SS (a) and pCpG-free-EF1α-SS (b) vectors, respectively. Experiments were performed in two replicates with random enhancer libraries in GP5d and HepG2 cells and random promoter–enhancer libraries in GP5d cells (250 million cells per each replicate). Genomic STARR-seq experiments were performed in two replicates in HepG2 cells (170 million cells per replicate) and four different conditions in GP5d cells (wild-type and TP53-null GP5d cells using both methylated and nonmethylated input DNA libraries; 500 million cells per condition). For random promoter–enhancer libraries in HepG2 and RPE1 cells, a total of 400 and 480 million cells were transfected, respectively. Briefly, a day before transfection, 6.7–10 million cells were plated per 15-cm dish in their respective media without antibiotics. The next morning, plasmid DNA was mixed with transfection reagent optimized for each cell line (Transfex (ATCC) for GP5d, Transfectin (Bio-Rad) for HepG2 and FuGENE HD (Promega) for RPE1) at a 1:3 ratio in Opti-MEM medium (Gibco), incubated for 15 min at room temperature and added dropwise to the cells. The cells were incubated for 24 h in a 37 °C incubator with 5% CO2.

37Cells were harvested and total RNA isolated 24 h after transfection using the RNeasy Maxi kit (Qiagen) with on-column DNase I digestion. The poly(A)+ RNA fraction was purified using the Dynabeads mRNA DIRECT Purification kit (Invitrogen) followed by DNase treatment using TurboDNase (Ambion) and purification using RNeasy Minelute kit (Qiagen) as previously described.

38The library preparation protocol was adapted from Arnold et al. essentially in all steps but with primers matching our modified STARR-seq vectors, and the exact protocol is described in Supplementary Methods.

39To generate a sequencing library using a template-switch strategy, a 40-µg aliquot of total RNA from the random promoter–random enhancer STARR-seq experiment in the GP5d cell line was used. See Supplementary Methods for the detailed protocol and section ‘Mapping TSS positions based on template switching’ for respective data analysis.

40ChIP-seq was performed as previously described using the following antibodies (2 μg per reaction): H3K27ac, H3K9me3 and H3K27me3 (Diagenode, C15410196, C15410193 and C15410195, respectively); FOXA1 (Abcam, ab23738); p53, HNF4a, IRF3 and CTCF (Santa Cruz Biotechnology, sc-126x, sc-8987x, sc-33641x and sc-15914x, respectively); SMC1 (Bethyl Laboratories, 300-055A); and normal rabbit, mouse and goat IgG (Santa Cruz Biotechnology, sc‐2027, sc‐2025 and sc‐2028, respectively). To analyze the genomic occupancy of TP53, GP5d cells were treated with 350 µM 5-fluorouracil (Sigma) 24 h before harvesting the cells. To analyze the effect of STARR-seq plasmid transfection on cellular alarm signals by ChIP-seq, RNA-seq and ATAC-seq, HepG2 cells were collected 24 h after the following treatments: mock (DMSO), 350 µM 5-fluorouracil treatment and transfection of the genomic STARR-seq library using similar conditions as in the STARR-seq experiments described above. The details of conditions, protocols and analysis parameters are described in Supplementary Methods.

41For RNA-seq, total RNA was isolated using RNeasy Mini kit (Qiagen) and RNA-seq libraries were generated using KAPA stranded mRNA-seq kit for Illumina (Roche). CAGE library was prepared from total RNA isolated from GP5d cells as previously described from 1 µg total RNA. The bisulfite sequencing data for DNA methylation in GP5d cells were obtained from Yin et al.. The ATAC-seq libraries were prepared from 50,000 cells as previously described for GP5d and HepG2 cells. All samples were paired-end sequenced using Illumina platforms.

42The ATI assay in GP5d cells and processing of the data were done as previously described. Transcribed enhancer regions defined using TT-seq data are based on Lidschreiber et al.. The experiments were performed from GP5d cells in two biological replicates as previously described. The details of the protocols and analyses are described in Supplementary Methods.

43For testing activities of known TF motifs, a set of 3,226 HT-SELEX motifs were collected (refs.,, and unpublished draft motifs). The motif collection and library design, measurement of enhancer activity of TF motif consensus sequences, and generation of activity PWM are described in detail in the Supplementary Methods. Moreover, previously described promoter motifs were used in TSS analyses, including TATA box, initiator, CCAAT-box and GC-box; BRE, MTE and DPE; and BANP.

44The complexity analysis of STARR-seq libraries generated from the TF motifs, genomic DNA and random enhancer libraries are described in Supplementary Methods; read counts are given in Supplementary Table 7.

45First, 150-base-long STARR-seq RNA and input DNA paired-end reads were combined using the FLASH program (version 1.2.11; options —min-overlap = 130 —max-overlap = 134 −x 0.25 −z −t 4), and only combined sequences of length 170 were chosen. To avoid including PCR duplicates of the same sequence with few mismatches due to sequencing errors, the sequences were sorted four times based on 40-bp-long nonoverlapping subsequences from base 6 to 165, and only one sequence per identical subsequence at each sorting step was taken. This ensured that from sequences that had Hamming distance less than 4, only one was taken. Only one representative sequence from the similar sequences was used for downstream analysis, so each sequence is either present or absent in the sample. The sequences are sampled from a huge input DNA library, which prohibits precise determination of initial input frequencies of individual sequences. Thus, our analysis relies on finding common features of different selected sequences instead of their counts. The numbers of preprocessed sequences used in downstream analysis are shown in Supplementary Table 7.

46The active enhancers were identified by calling the peaks from the STARR-seq–enriched RNA fragments against the plasmid input sample using MACS2 (ref. ). The preprocessing of genomic STARR-seq data, peak calling, overlap with genomic features, de novo motif mining and conservation of genomic STARR-seq elements are described in detail in Supplementary Methods.

47The STARR-seq enhancer sequences derived from RNA were mapped to corresponding promoter–enhancer pairs in the input DNA by exact matches of the first 20 bases of the 150-base-long enhancer sequences. Duplicate sequences were removed as described for random enhancers, except that three 40-bp-long subsequences from 16 to 135 were used, thus ensuring that only one of the sequences with Hamming distance less than 3 was chosen (Preprocessing of random enhancer library sequences). Then, promoter and enhancer sequences were filtered separately by removing (1) all adapter sequences that included some (partial) adapter sequence according to cutadapt version 1.9.1 (ref. ), (2) sequences that mapped to plasmid backbone sequence using bowtie2 version 2.2.4 (ref. ) and (3) outlier sequences in terms of nucleotide composition (count of any nucleotide more than three median absolute deviations higher than the median count) that removes, for example, those high-G-content sequences that are an Illumina artifact. After preprocessing, the correlation between observed dinucleotide frequencies and those expected based on mononucleotide frequencies was over 0.99 (GP5d random replicate 1 enhancers). Input DNA sequences were processed the same way. For promoter–enhancer pair analyses, the remaining promoter–enhancer pairs were collected, and pairs containing highly similar sequence as a promoter and an enhancer were removed. The numbers of sequences used in downstream analysis are shown in Supplementary Table 7.

48First, the synthetic random sequences derived from spliced transcripts were identified using the constant sequence spanning the splice site after intron removal (cutadapt program); other sequences were not processed further. Next, the UMI sequence was removed from the 5ʹ end of each sequence, and the last 20 bp of its random part was used to recognize the corresponding promoter from the input DNA. To accurately recognize the first base of the transcript and thus the position of the TSS, it was assumed that the template-switch process had added at least three and at most four guanines to the 5ʹ end of the transcript. On this basis, only the RNA sequences starting with at least three Gs were used in the analysis. Each such sequence was aligned to the corresponding input DNA promoter sequence using an exact 20-bp match starting from the sixth base to allow for the extra Gs. Finally, the Gs added by the template switch were trimmed and discordant sequences removed according to the alignment. The frequency of the four Gs instead of three was estimated from the sequences that do not have G at the fourth position in the alignment to the input. For those that did have a G also in the input sequence, removing three or four Gs was decided randomly but so that the frequency of the fourth G matched the estimate. The two GP5d template-switch libraries were processed separately and then merged so that only one transcript was kept for each input DNA promoter sequence to prevent duplicate promoter sequences. The exact positions of the TSSs at the promoters were recorded, and the flanking sequences were used for further analysis of the positioning of different sequence features relative to TSSs. The numbers of sequences obtained are listed in Supplementary Table 7, as the number of flanking sequences fitting to the random region depends on the flank sizes. The comparison to human endogenous promoters was done using TSS positions from EPD (hsEPDnew 006).

49The motif matching was done using MOODS version 1.9.3 (ref. ), and fold changes were estimated using the function PsiLFC in R package lfc. The details of motif matching and motif spacing analysis in STARR-seq random enhancers are described in Supplementary Methods.

50For RNA and randomly sampled input DNA promoter–enhancer pairs, the number of such pairs that one motif occurs in the promoter and a second one in the enhancer was counted for each motif pair (excluding heterodimers) and motif-match strand combination (++, +−, −+ and −−). The counts over the strand combinations were summed to get the total number of pairs, and the fold change of the number of pairs between input DNA and RNA was estimated using the function PsiLFC in R package lfc. If the promoter and enhancer occurrences are independent of each other, then the expected frequency of the pair of sites is the product of the individual frequencies. The expected log2 fold change assuming independent actions of the promoter and enhancer motif was thus calculated as the sum of their individual log2 fold changes. The same analysis was done using the reversed, but not complemented, control motifs.

51To estimate the significance of the number of observed motif pairs, our null hypothesis was that the probability of observing a motif-match pair in a promoter–enhancer sequence pair was the same as estimated from the individual motif-match frequencies. The two-sided binomial tests done for 528,529 motif pairs resulted in a significant P value (multiple hypothesis-corrected P value < 0.05, Holm’s method) for 253 pairs.

52For analyzing motif positioning within active promoters derived from synthetic random sequences, motifs were matched to sequences flanking the TSSs identified from the template-switch data, and for each motif, only the highest-affinity match per sequence was considered. The number of matches for each motif was then counted separately at each position and strand. To get positional activity scores for position-specific regression analysis, motif matching was done for TSS flanking sequences from position −100 to +20 in relation to TSS and for a control set generated by sampling for each TSS sequence a subsequence of same length from the same position from an input DNA promoter (background probability of a match 5 × 10−4). The log2 fold changes of the motif match counts between TSS flanking set and control set (estimated with the lfc package) were then used as a positional activity score for each position and strand.

53To study p53 motif-match positioning relative to the STARR-seq vector, the motif was matched (background probability of a match 10−5) to highly selected sequences chosen by taking only sequences observed at least twice in both GP5d enhancer replicate experiments. A histogram of match start positions was generated by counting only the highest-affinity match in each sequence. A smoothed density estimate was generated using a Gaussian kernel (R ggplot geom_density with adjust=0.5).

54The mutual information analysis was done as previously described for the aligned STARR-seq reads (60 + 60 bp surrounding the TSS from the template-switch data); details are described in Supplementary Methods.

55The datasets used in each machine learning analysis and their division into training, test and validation sets are detailed in Supplementary Table 7. To enable sequences from genomic measurements (genomic STARR-seq and ATAC-seq) to be scored on the CNNs that were trained on the random enhancer STARR-seq data and vice versa, the length of the sequences fed to these models was standardized to 170 bp. Thus, additional preprocessing specific to machine learning analyses was done for the genomic STARR-seq and ATAC-seq data as detailed in the Supplementary Methods.

56The details of machine learning analyses using logistic regression and CNNs, discussion about the optimal classification of random enhancer STARR-seq data, details of differential expression prediction, interpretation of CNN classifiers and validation of the predicted promoter mutation effects with external data are described in Supplementary Methods. Briefly, the logistic regression classifiers were implemented using the LogisticRegression function from scikit-learn (version 0.21.3) library; the CNN models were built on Keras (https://keras.io/; version 2.2.4) using TensorFlow 1.14.0 backend; DeepLIFT version 0.6.12.0 (ref. ) was used to visualize the TERT promoter mutations and study the sequence features learned by the random enhancer STARR-seq CNN model along with TF-MoDISco version 0.5.14.1 (ref. ).

57The binary STARR-seq design allows looking for relatively short-range interactions between promoters and enhancers, and the details of machine learning analysis used for testing such interactions are detailed in the Supplementary Methods.

58All the promoter models were trained on data where the TSS is 100 bp from the start of the training sequence. Thus, scoring any 120 bp sequence with these models gives a probability that the position 100 in these sequences is a TSS of a functional promoter sequence. Each possible TSS position within ±500 bp from the TSSs of the active GP5d promoters was analyzed by taking 100 bp upstream and 20 bp downstream from the candidate TSS and scoring these sequences with the promoter models. Active GP5d promoters were defined as those EPD promoters that overlapped with a GP5d CAGE peak. For each test set, active GP5d TSS and promoter model, the position obtaining the highest promoter probability from the corresponding model was taken as the predicted TSS position.

59Human promoter coordinates were obtained from the EPD version 006, hg19 (ref. ), and their preprocessing, together with analysis of the CAGE data, is described in Supplementary Methods.

60Further information on research design is available in the Nature Research Reporting Summary linked to this article.

61Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41588-021-01009-4.