Assessment of cross-assay consistency in enhancer identification

We analyzed six distinct STARR-seq and MPRA datasets produced by laboratories within the ENCODE Consortium’s Functional Characterization Center, comprising three TilingMPRA datasets, a LentiMPRA dataset, an ATAC-STARR-seq, and a WHG-STARR-seq dataset [26,27,28]. Although all assays were performed in the human K562 cell line, they differed in experimental objectives, design strategies, and data processing methods. An overview of these experimental designs is illustrated in Fig. 1a, with detailed dataset descriptions provided in the Additional file 1.

Fig. 1

Overview of experimental designs and assay consistencies across MPRAs and STARR-seq assays. a Schematic representation of the experimental workflows for four types of MPRAs and STARR-seq assays analyzed in this study. b Heatmap displaying the number of overlapping enhancer regions between assays and their percentage relative to the total number of enhancer regions identified in each assay, based on the ≥ 1-bp overlap criterion. c Heatmap presenting the Jaccard Index for pairwise comparisons between assays using ≥ 1-bp overlap criterion, quantifying overall similarity in enhancer identification

To evaluate the consistency of enhancer identification, we compared enhancer calls reported in each dataset. Data were retrieved from either the ENCODE portal [26,27,28] or original publications and processed according to each laboratory’s guidelines. The original number of lab-reported enhancer regions is summarized in Additional file 3: Table S1. To standardize comparisons, overlapping enhancer calls within each dataset were merged into unique regions, resulting in 12,919 enhancer regions across three TilingMPRA datasets, 56,840 regions from LentiMPRA, and 46,906 and 38,671 regions from ATAC-STARR-seq and WHG-STARR-seq, respectively.

We compared enhancer calls across assays by measuring the number of overlapping enhancer regions in each pairwise comparison, applying a minimal overlap threshold of 1 base pair (bp) to ensure inclusion of partially overlapping regions (Additional file 2:Fig. S1b). The highest overlap was observed between LentiMPRA and ATAC-STARR-seq, where approximately 40% (22,780 out of 56,840) of LentiMPRA regions overlapped with 44% (20,692 out of 46,906) of ATAC-STARR-seq regions. ATAC-STARR-seq and WHG-STARR-seq showed the second-highest overlap, with around 11% (5359 out of 46,906) of ATAC-STARR-seq regions overlapping with 16% (6255 out of 38,671) of WHG-STARR-seq regions. Comparisons involving LentiMPRA and WHG-STARR-seq, as well as TilingMPRA with other assays, exhibited lower overlap, reflecting differences in enhancer calls across these datasets (Fig. 1b).

To further quantify similarity across assays, we calculated the Jaccard Index (JI) for each pairwise comparison. Overall, enhancer identification exhibited low consistency, with most JI values approaching zero (Fig. 1c). The highest JI was observed between LentiMPRA and ATAC-STARR-seq (0.28), followed by ATAC-STARR-seq and WHG-STARR-seq (0.08). Applying stricter overlap criteria further reduced similarities (Additional file 2: Fig. S2a,b), highlighting the substantial variability in enhancer identification across different assays.

Unified processing of reporter assay datasets: initial data quality check

To better understand the limited consistency in enhancer identification observed across MPRA and STARR-seq datasets, we examined potential technical and biological contributors to assay variability, as detailed in the Additional file 1. While both biological and technical factors may contribute, a previous study [29] and our comparative analyses of assay design and data processing suggest that technical variation is the predominant driver of cross-assay discrepancies. These considerations underscore the importance of a standardized analytical framework that accounts for key technical differences to enable more accurate and equitable comparisons across assays. To this end, we implemented a unified processing strategy, beginning with a comprehensive quality assessment of all datasets. In addition to the factors discussed in the Additional file 1, the baseline quality of these datasets directly impacts assay consistency and reliability.

For TilingMPRA, LentiMPRA, count data were readily available through the ENCODE portal [26,27,28], or were re-processed with guidance from the original authors [30]. For ATAC-STARR-seq and WHG-STARR-seq datasets, we applied a unified genomic binning approach, creating 100-bp genomic bins with a 10-bp step size in both forward and reverse orientations (Additional file 2: Fig. S3a). Only fragments that fully covered each genomic bin were counted, allowing for orientation-independent enhancer identification and a more accurate assessment of genome-wide coverage of tested regions.

We assessed genome-wide coverage for the WHG-STARR-seq dataset and observed substantial library complexity, with over 96% of the human genome assayed after processing using the genomic binning approach. However, a more detailed evaluation revealed that a notable subset of genomic bins exhibited low read depths in the DNA library (< 10; Methods) (Fig. 2a), raising concerns about limitations in sequencing depth and transfection efficiency. Because low-read-depth regions are typically excluded from downstream analyses, the reported genome-wide coverage likely overstates the extent of the genome that was effectively tested.

Fig. 2

Evaluation of data quality across MPRAs and STARR-seq assays. a Funnel plot showing the genome-wide coverage distribution of WHG-STARR-seq at varying read depths thresholds in DNA libraries. b Funnel plot illustrating the coverage distribution of accessible regions at different read depths thresholds in DNA libraries for ATAC-STARR-seq and LentiMPRA. Accessible regions are defined by ATAC-seq peaks from ATAC-STARR-seq DNA libraries and DNase-seq narrow peaks for LentiMPRA. c Bar plot presenting average Pearson correlation coefficients for log2-transformed DNA CPM and RNA CPM, and log2(RNA/DNA) ratios across assays. d Bar plot depicting average library recovery rates in DNA and RNA libraries across assays.

We also evaluated the coverage of accessible chromatin regions in the ATAC-STARR-seq and LentiMPRA datasets, where assayed fragments were either enriched in or selected from these regions. Both datasets demonstrated the ability to capture a substantial proportion of accessible regions with high read depths (Fig. 2b). Specifically, ATAC-STARR-seq achieved almost 100% of coverage of accessible regions characterized by ATAC-seq peaks, while LentiMPRA successfully covered 44% of DNase hypersensitive sites (DHSs) at higher read-depth threshold (≥ 10) in DNA libraries.

To assess reproducibility between replicates, we calculated Pearson correlations (ρ) for log-transformed counts per million (logCPM) of DNA and RNA counts, as well as log2(RNA/DNA) ratios, as these ratios represent the primary measurement of enhancer activity in downstream analyses. Overall, TilingMPRA and LentiMPRA demonstrated strong replicate correlations, indicating high reproducibility across libraries (Fig. 2c). Specifically, LentiMPRA showed robust correlations for both logCPM of DNA and RNA counts (0.97 < ρ < 0.99) and log2(RNA/DNA) ratios (0.72 < ρ < 0.80). Among the TilingMPRA datasets, ENCSR917SFD and ENCSR363XER displayed consistently high correlations (0.96 < ρ < 0.99 for logCPM, 0.87 < ρ < 0.90 for log2(RNA/DNA)), while ENCSR394HXI had moderately lower values (0.62 < ρ < 0.89 for logCPM, 0.47 < ρ < 0.58 for log2(RNA/DNA)), suggesting some variability within this dataset.

In contrast, ATAC-STARR-seq and WHG-STARR-seq demonstrated considerably lower fragment-level reproducibility (Fig. 2c). ATAC-STARR-seq showed weak agreement between replicates (0.001 < ρ < 0.26 for logCPM, 0.12 < ρ < 0.22 for log2(RNA/DNA)), while WHG-STARR-seq exhibited even greater variability, including negative RNA correlations (Fig. 2c). Aggregating fragments into genomic bins markedly improved replicate reproducibility for DNA and RNA counts and log2(RNA/DNA) ratios in the ATAC-STARR-seq dataset and RNA counts in WHG-STARR-seq dataset (Fig. 2c). Despite these improvements, the correlations for log2(RNA/DNA) ratios remained low in both datasets (0.18 < ρ < 0.37 for ATAC-STARR-seq, 0.42 < ρ < 0.47 for WHG-STARR-seq). Further restricting analysis to accessible genomic bins in ATAC-STARR-seq provided marginal improvements but did not reach the high reproducibility observed in MPRA datasets, highlighting persistent variability in genome-wide STARR-seq measurements.

We also evaluated library recovery rates by calculating the proportion of fragments or genomic bins with at least one read in each library. TilingMPRA and LentiMPRA had high recovery rates (89–100%), whereas ATAC-STARR-seq exhibited an average library recovery rate below 40% in DNA libraries and even lower in RNA libraries (Fig. 2d). These findings suggest that many fragments were not consistently detected in ATAC-STARR-seq, possibly due to low sequencing depth or low transfection efficiency. Further analysis of fragments overlapping ATAC-seq peaks showed similar discrepancies in recovery rates between DNA and RNA libraries, pointing to limitations in data quality (Fig. 2d). WHG-STARR-seq also had low recovery rates at the fragment level (18–38%), but most genomic bins were represented in both DNA and RNA libraries (93–98%) (Fig. 2d), indicating that issues with sequencing depth and transfection efficiency were not as severe.

These results revealed substantial variability in data quality across different datasets. While MPRA assays exhibited consistently high data quality and reproducibility, genome-wide STARR-seq datasets were more susceptible to limitations such as insufficient sequencing depth and potential low transfection efficiency. These factors likely contributed to higher variability and reduced reliability in enhancer identification, and this issue can remain significant even when genomic binning is applied. Our findings highlight the necessity of applying stringent filtering criteria to exclude low-read-depth regions in the downstream analysis while also ensuring that the final reported tested regions accurately represent sequences with sufficient read depth, rather than using all assayed regions as a proxy for measuring tested region coverage.

Uniform processing of reporter assay datasets: enhancer call pipeline

While future studies should further address experimental challenges, to address the role of data processing in contributing to the observed inconsistencies, we implemented a unified enhancer call pipeline and applied it consistently across all datasets. The workflow is illustrated in Fig. 3a, with detailed methodology provided in the Methods section.

Fig. 3

Enhancer identification using a unified pipeline. a Schematic of the uniform enhancer call pipeline. The workflow begins with a raw count matrix as input, applies dataset-specific filters to exclude low-depth regions, and normalizes library size using TMM normalization. Regulatory activity is calculated as log2(RNA/DNA) values and Z-score analysis is performed to identify regions with significantly higher regulatory activity than negative control regions as enhancer regions in an orientation-independent manner. b Bar plot showing the assayed coverage, tested coverage in either orientation and tested coverage in both orientation for open chromatin regions characterized by ATAC-seq peaks derived from DNA libraries in ATAC-STARR-seq. c Bar plot summarizing genome-wide assayed coverage, tested coverage in either orientation and tested coverage in both orientations. d Meta-plots comparing the average DNase-seq and ATAC-seq signal profiles (± 1 kb from the center) for 2000 enhancer regions randomly sampled from those identified in both orientations versus 2000 regions tested in both orientations but active in only one orientation. e Meta-plots comparing the average of H3K4me3 and H3K27ac histone modification profiles (± 1 kb from the center) for 2000 enhancer regions randomly sampled from those identified in both orientations versus 2000 regions tested in both orientations but active in only one orientation. f Meta-plots comparing the average DNase-seq and ATAC-seq signal profile (± 1 kb from the center) for 2000 randomly sampled enhancer regions from laboratory-reported enhancer calls versus those identified using the uniform enhancer call pipeline. g Meta-plots comparing the average of H3K4me3 and H3K27ac histone profiles (± 1 kb from the center) for 2000 randomly sampled enhancer regions from laboratory-reported enhancer calls versus those identified using the uniform enhancer call pipeline

The pipeline begins with a raw count matrix as input and applies dataset-specific filters to remove fragments or genomic bins with low read depth. We then adapted the Trimmed Mean of M-values (TMM) normalization [31] and linear model approach from the Limma-Voom pipeline [32] to calculate log2(RNA/DNA) as a measure of regulatory activity for each fragment or genomic bin in each orientation. For targeted assays that included negative control sequences, we modified the original TMM normalization method to rely solely on negative controls for adjusting library size and composition bias. This approach provides greater accuracy in normalization, particularly for targeted assays where the assumption that most fragments lack regulatory effects may not hold.

After computing the log2(RNA/DNA) values, we assessed the regulatory activity of each fragment or genomic bin in both orientations by comparing it to the activity levels of negative controls through a Z-score analysis rather than relying on an arbitrary log2(RNA/DNA) cutoff. This comparison allowed for the identification of regions with significantly elevated activity relative to the basal transcription level defined by the negative controls in each orientation. To mitigate orientation bias, we incorporated regulatory activity in both orientations as a criterion for determining whether a fragment or genomic bin qualifies as a potential enhancer.

For genome-wide STARR-seq datasets that lacked negative controls in the original assays, we used genomic bins within exonic regions as surrogate negative controls, as enhancers are predominantly located in non-coding regions [26, 33, 34]. To ensure a clear distinction between potential enhancer regions and those likely to exhibit basal transcription, we excluded genomic bins overlapping the 300-bp flanking regions on either side of exons. This approach minimizes the risk of using genomic bins that may have counted fragments overlapping with enhancers in intronic regions, increasing the reliability of these surrogate negative controls.

Finally, our pipeline recorded both active and inactive regions identified in an orientation-independent manner, ensuring an accurate assessment of genome-wide coverage of tested regions. This comprehensive reporting approach also enables robust cross-assay comparisons. Detailed numbers of fragments or genomic bins tested in one or both orientations, the numbers of negative controls, and the numbers of enhancer regions identified are provided in the Additional file 3: Table S2.

Improved enhancer identification through unified enhancer call pipeline

We applied the uniform enhancer call pipeline to all datasets to standardize the identification of enhancer regions. In the ATAC-STARR-seq dataset, while all accessible regions characterized by ATAC-seq peaks were initially included in the assay, 91.20% were statistically tested for regulatory activity in at least one orientation (Fig. 3b). Furthermore, the effective coverage of regions tested in both orientations within accessible chromatin was reduced to 64.72% (Fig. 3b). Similarly, for the WHG-STARR-seq dataset, 96.61% of the entire human genome was included in the assay; however, only 56.15% of regions were statistically assessed in at least one orientation, with just 44.59% tested in both orientations (Fig. 3c). These findings reveal that the effective coverage of genome-wide STARR-seq datasets is significantly lower than expected, underscoring the importance of comprehensive reporting of tested regions to accurately evaluate assay performance and coverage.

Using our unified pipeline, we identified 57 enhancer regions in TilingMPRA (ENCSR394HXI), 26,874 in LentiMPRA, 11,507 in ATAC-STARR-seq, and 25,274 in WHG-STARR-seq. Notably, these enhancer regions exhibited significant regulatory activity in both orientations. For the two TilingMPRA datasets (ENCSR817SFD and ENCSR363XER), which tested elements exclusively in one orientation, we adapted our pipeline to perform orientation-dependent analysis, identifying 2117 enhancer regions in ENCSR817SFD and 3752 in ENCSR363XER.

To evaluate the significance of making orientation-independent enhancer calls, we investigated their epigenomic features by analyzing 2000-bp windows centered on these regions. Specifically, we compared the epigenomic features of orientation-independent enhancer regions to those of regions that were tested in both orientations but exhibited significant activity in only one, leveraging ENCODE datasets for DNase-seq, ATAC-seq, and ChIP-seq (H3K4me3 and H3K27ac) in the K562 call line. Orientation-independent enhancers displayed higher chromatin accessibility, as indicated by stronger DNase-seq and ATAC-seq signal intensities compared to enhancers active in only one orientation across all datasets (Fig. 3d). Additionally, they exhibited greater enrichment of both promoter- and enhancer-associated histone modifications, with a more pronounced bimodal patter around their centers (Fig. 3e). These findings suggest that orientation-independent enhancers are more robustly marked by epigenomic features characteristic of active regulatory elements and highlight the importance of making orientation-independent enhancer calls.

We also compared the enhancer regions identified through our unified processing pipeline with the original enhancer calls reported by each laboratory. Across all datasets, uniformly processed enhancer regions exhibited higher chromatin accessibility, as evidenced by stronger DNase-seq and ATAC-seq signals (Fig. 3f). Notably, while some enhancer calls from the unified pipeline were in inaccessible regions, they were still more enriched in accessible regions compared to original lab-reported peaks in the WHG-STARR-seq dataset (Fig. 3f). Additionally, histone modification profiles confirmed that orientation-independent enhancer regions identified by the unified pipeline were more strongly marked by H3K4me3 and H3K27ac compared to lab-reported enhancer regions (Fig. 3g). These results highlight the advantages of our unified pipeline in enhancing the confidence of enhancer identification and providing a more reliable foundation for comparative and functional studies.

Enhanced consistency across assay using uniform processed enhancer calls

With both active and inactive regions recorded through our uniform enhancer call pipeline, we reassessed assay consistency by evaluating how many enhancers identified in one assay were also identified as enhancers in others. To achieve this, we conducted pairwise comparisons by assessing the overlap between enhancer regions from one assay and all tested regions in another. Because our enhancer regions were defined in an orientation-independent manner, inactive regions were also generated by merging elements or genomic bins tested in both orientations that lacked significant enhancer activity.

For each pairwise comparison, enhancer regions from assay A were evaluated against all tested regions in assay B, and vice versa, as overlaps were not necessarily symmetric. In cases where an enhancer region overlapped multiple tested regions in another assay, or multiple enhancer regions overlapped a single tested region, we assigned the best overlap based on the highest number of overlapping base pairs to minimize redundancy. We then calculated the JI and recorded both the number of enhancer regions that were also classified as enhancers in the other assay and the total number of enhancer regions tested. By restricting comparisons to commonly tested regions, this approach provided a more accurate and comprehensive assessment of cross-assay consistency.

Using the minimal overlap threshold (≥ 1 bp), we observed a statistically significant improvement in assay consistency, reflected by higher JI values (Fig. 4a; one-sided Wilcoxon paired test, p = 0.02). Applying a stricter threshold of ≥ 50% reciprocal overlap likewise yielded a significant increase in cross-assay consistency, with JI values substantially exceeding those based on lab-reported enhancer regions (one-sided Wilcoxon paired test, p = 0.001). Together, these results demonstrate that adopting a uniform enhancer calling pipeline and employing refined comparison strategies enhances cross-assay consistency, underscoring the importance of standardized processing in functional characterization studies.

Fig. 4

Enhanced consistency in cross-assay comparisons using uniformly processed enhancer calls. a Box plot showing the Jaccard Index for pairwise comparisons between assays, calculated using the minimal overlap criterion of 1-bp and the stricter criterion of ≥ 50% reciprocal overlap. Results are shown for both laboratory-reported and uniformly processed enhancer calls, illustrating the improved consistency achieved through uniform processing. b,c Heatmaps displaying the number of overlapping enhancer regions between assays under the ≥ 1-bp overlap criterion (b) and under the ≥ 50% reciprocal overlap criterion (c). Each cell shows the ratio of the number of enhancer regions in the row dataset that overlap with enhancer regions in the column dataset to the number of enhancer regions in the row dataset overlapping with tested regions in the column dataset. Diagonal cells display the total number of enhancer regions identified in each dataset. d,e Heatmaps displaying the number of overlapping enhancer regions between assays under the ≥ 50% reciprocal overlap criterion in proximal regions (d) and distal regions (e)

Assay-specific factors influence cross-assay consistency

While previous comparisons using lab-reported enhancer regions showed lower agreement across assays when a stricter overlap criterion (≥ 50% reciprocal overlap) was applied, comparisons using uniformly processed data demonstrated the opposite trend: most pairwise comparisons exhibited increased JI values under the stricter criterion compared to the ≥ 1-bp threshold (Fig. 4a). For instance, when comparing LentiMPRA enhancers to tested regions in ATAC-STARR-seq and WHG-STARR-seq, the proportion of consistently active regions rose from 17 and 19% (using a ≥ 1-bp threshold) to 74 and 78% (using ≥ 50% reciprocal overlap), respectively. A similar pattern was observed in pairwise comparisons between ATAC-STARR-seq and WHG-STARR-seq (Fig. 4b,c), indicating that enhancer identification is more consistent when comparing sequences with greater overlap.

Despite the overall increase in consistency under more stringent overlap criteria, cross-assay agreement remained largely unchanged when comparing enhancer regions identified by ATAC-STARR-seq and WHG-STARR-seq to those tested in LentiMPRA, regardless of the overlap threshold (Fig. 4b,c). This suggests that assay-specific factors, rather than sequence overlap alone, dominate cross-assay agreement with LentiMPRA. Because LentiMPRA positions candidate sequences immediately upstream of a reporter gene, we suspected that its propensity to capture promoter activity, rather than enhancer activity, may explain this difference.

To test this, we stratified comparisons by TSS proximity, excluding TilingMPRA due to limited sample size. Tested regions were defined as proximal if ≥ 90% of their sequence overlapped within 500 bp of a protein-coding TSS (GENCODE [35] annotation v45) and distal otherwise. Stratification revealed that ATAC-STARR-seq and WHG-STARR-seq showed significantly higher consistency with LentiMPRA in proximal regions than in distal regions. Specifically, using a ≥ 50% reciprocal overlap threshold, ~ 76–82% of proximal enhancer regions identified by STARR-seq assays were also active in LentiMPRA, compared to only ~ 52–64% of distal regions (Fig. 4d,e). Similar trends were observed under the ≥ 1 bp overlap threshold (Additional file 2: Fig. S4a,b). These findings indicate that LentiMPRA aligns more closely with STARR-seq assays in detecting regulatory sequences in proximal regions. However, because MPRA-based assays may preferentially capture promoter rather than enhancer activity, further investigation is required to determine whether sequences active in both LentiMPRA and STARR-seq assays reflect enhancer-driven or promoter-driven regulation.

To test this, we assessed assay consistency separately in proximal and distal regions. TilingMPRA is excluded from this analysis due to limited sample size. Tested regions were classified as proximal if ≥ 90% of their sequence overlapped within 500 bp of a protein-coding TSS (based on GENCODE [35] annotation v45) and distal otherwise. Stratifying comparisons by TSS proximity revealed that ATAC-STARR-seq and WHG-STARR-seq exhibited significantly higher consistency with LentiMPRA in proximal regions than in distal regions. Specifically, ~ 62–73% of proximal enhancer regions identified by STARR-seq assays were also active in LentiMPRA, whereas only ~ 33–47% of distal enhancer regions showed consistent activity (Fig. 4d,e and Additional file 2: Fig. S4b,c). Notably, these proportions differed only when comparing distal versus proximal regions but remained largely unchanged across different overlap thresholds (Fig. 4d,e). These findings suggest that LentiMPRA is more likely capturing promoter activity rather than enhancer activity as measured in genome-wide STARR-seq assays, emphasizing that assay-specific factors play a dominant role in determining cross-assay consistency when comparing to LentiMPRA.

Taken together, these results demonstrate that implementing a uniform enhancer calling pipeline effectively reduces technical variation, ensuring that cross-assay comparisons are not confounded by differences in data processing. The observed patterns of agreement or disagreement instead primarily reflect assay-specific biological factors intrinsic to each experimental design. This underscores the importance of standardized data processing for fair cross-assay evaluation, while also indicating that the remaining inconsistencies arise from biological properties of the assays rather than technical variation.

Evaluating functional support for enhancer-like and promoter like sequences in cCREs

Epigenomic features such as DNA accessibility and histone modifications have long been recognized as key indicators of active enhancers [5, 8, 34]. Leveraging these features, the ENCODE Consortium established a registry of cCREs [8]. To assess how well these elements are functionally validated by massively parallel reporter assays, we examined their coverage and activity in LentiMPRA, ATAC-STARR-seq, and WHG-STARR-seq datasets.

Since cCREs were not specifically designed as targeted sequences in these assays, we assessed their coverage by identifying overlaps between cCRE elements and tested regions. A cCRE was considered covered if it had at least a 1-bp overlap with a tested region. To further characterize their representation across assays, we categorized covered cCREs into three mutually exclusive groups based on their overlap extent: high (≥ 80% reciprocal overlap), moderate (50–80% reciprocal overlap), and low (all other overlap). Detailed coverage statistics are provided in Additional file 2: Fig. S5a and Additional file 3: Table S3.

To evaluate the functional relevance of cCREs, we analyzed their active rates across LentiMPRA, ATAC-STARR-seq, and WHG-STARR-seq (Additional file 2: Fig. S5b). In both genome-wide STARR-seq datasets, cCREs associated with enhancer-like and promoter-like signatures—dELS, pELS, and PLS—demonstrated the highest active rates among all cCRE subtypes, whereas other cCRE categories exhibited lower active rates. Specifically, high-overlap dELS, pELS, and PLS each showed active rates ranging from 46 to 89% in ATAC-STARR-seq and WHG-STARR-seq (Additional file 2: Fig. S5b), highlighting their strong functional relevance in both genome-wide STARR-seq datasets. In contrast, the active rates of other cCRE subtypes declined sharply, with CA-H3K4me3 and CA-TF elements exhibiting moderate active rates (24–47%), followed by CA-CTCF and CA-only elements, which showed more limited active rates (5–9%). As expected, low-DNase elements, which are generally classified as inactive cCREs, displayed the lowest active rates (2–4%), only slightly higher than regions without any cCRE overlap (0.4–0.5%).

While the overall active-rate patterns were consistent across cCRE subtypes in genome-wide STARR-seq datasets, LentiMPRA displayed a distinct trend. PLS elements showed the highest active rate (57%), whereas dELS (24%) and pELS (22%) exhibited activity levels comparable to CA-H3K4me3 (30%), CA-TF (23%), and CA-only (31%). Moreover, low-DNase elements showed a 15% active rate in LentiMPRA, substantially higher than the 2–4% observed in ATAC-STARR-seq and WHG-STARR-seq, and regions without overlap with any cCREs also exhibited elevated activity (4%) compared to the minimal levels detected in STARR-seq assays (0.4–0.5%) (Additional file 2: Fig. S5b). Together, these observations highlight assay-specific biological factors shaping LentiMPRA activity profiles and suggest multiple, non-exclusive explanations: LentiMPRA may be more sensitive to promoter-associated activity, and its in-genome, integration-based readout likely introduces chromatin-context effects that differ from plasmid-based STARR-seq assays, potentially positioning certain sequences (e.g., low-DNase elements) into accessible chromatin environments and thereby inflating their apparent activity.

Collectively, these findings highlight the predictive power of cCREs in identifying active enhancers in reporter assays, particularly for dELS, pELS, and PLS, which exhibited significantly higher activity than other cCRE categories. The near absence of enhancer activity in regions lacking biochemical features underscores the essential role of chromatin accessibility and histone modifications in defining functional enhancers. At the same time, the distinct activity patterns observed in LentiMPRA, likely reflecting its preference for promoter-associated sequences and other assay-specific influences, emphasize the need to carefully consider assay-specific factors when interpreting results and integrating data from different massively parallel reporter assays. Further in-depth investigation will be required to dissect the biological factors underlying the inconsistencies observed between LentiMPRA and STARR-seq assays.

Transcription as a critical mark of active enhancers

In addition to epigenomic features, enhancers are distinguished by their ability to generate eRNAs through divergent transcription [36, 37]. Tippens et al. demonstrated that divergent transcription serves as a more precise marker of active enhancers than histone modifications and identified a fundamental enhancer unit based on divergent transcription start sites (TSSs) [18]. Expanding on this, Yao et al. showed that GRO/PRO-cap is the most effective experimental approach to identify eRNAs and their divergent TSSs, and further compiled an enhancer compendium with a unified definition of enhancers based on divergent transcription [38].

Leveraging uniformly processed enhancer calls from large-scale reporter assays, we next examined these transcriptional characteristics of enhancers. Using the same analytical framework applied to cCREs, we assessed the coverage of GRO-cap enhancers [38] (divergent elements identified by PINTS from GRO-cap data) across the three assays. Detailed statistics are provided in Additional file 3: Table S4 and Additional file 2: Fig. S5d.

High-overlap GRO-cap enhancers exhibited strong enhancer activity, with 87 and 78% being active in ATAC-STARR-seq and WHG-STARR-seq, respectively (Additional file 2: Fig. S5e). Furthermore, GRO-cap enhancers consistently displayed significantly higher active rates compared to regions that neither overlapped with any GRO-cap elements nor exhibited GRO-cap signals (Additional file 2: Fig. S5e,f). Notably, while regions devoid of both transcriptional signals and overlap with GRO-cap elements exhibited the lowest active rates across all three assays (0.7–10%), regions that did not overlap with any annotated GRO-cap elements but still contained detectable GRO-cap signals showed slightly higher, albeit low, levels of activity (2–21%) (Additional file 2: Fig. S5f). These findings reinforce the strong functional relevance of GRO-cap enhancers in reporter assays, demonstrating that divergent transcription is a defining characteristic of active enhancers and supporting the enhancer architecture defined by previous studies [18, 38].

To further explore the functional relevance of transcriptional level, we categorized tested regions in LentiMPRA, ATAC-STARR-seq, and WHG-STARR-seq into four transcription-level classes (high, medium, low, and none) based on GRO-cap signals [39] (see Methods) and calculated the active rates within each category. Our analysis revealed a clear positive relationship: regions with higher transcription levels were significantly more likely to be classified as active across all three assays (Fig. 5a). Regions with no or low GRO-cap signals exhibited minimal active rates, particularly in ATAC-STARR-seq and WHG-STARR-seq, where values remained below 1%. Regions with medium transcription levels showed moderate active rates (6–24%), whereas highly transcribed regions reached the highest rates, with 31–50% of tested regions classified as active (Fig. 5a). These findings reinforce transcription level as a key predictor of enhancer activity across reporter assays.

Fig. 5

Impact of transcription levels on active rates and assay consistencies. a Bar plot illustrating the active rates of all tested regions in LentiMPRA, ATAC-STARR-seq, and WHG-STARR-seq, categorized by transcription levels (none, low, medium, and high) determined by GRO-cap signals. b Line plot depicting Jaccard Index values for pairwise comparisons between LentiMPRA, ATAC-STARR-seq, and WHG-STARR-seq across all tested regions with varying transcription levels, calculated using the ≥ 50% reciprocal overlap criterion. c Bar plot illustrating the active rate of transcribed and untranscribed regions with high-overlap with any types of cCREs or with high-overlap with dELS, pELS, and PLS or without any overlap with cCRE and PINTS elements. d Bar plot showing the active rate of high-overlap dELS, pELS, and PLS regions with different transcription levels (low, medium, high) determined by GRO-cap signals

Despite the low active rates observed in regions with little or no transcription, thousands of such regions were still identified as active enhancers across all three assays (Fig. 5a). This raised concerns that a subset of these enhancer calls might represent false-positive hits. To explore this possibility, we examined assay consistency across transcription classes, hypothesizing that regions with lower transcription levels would exhibit reduced cross-assay agreement, suggesting a higher prevalence of false positives. Indeed, using ≥ 50% reciprocal overlap as the comparison criterion, we observed a positive relationship between transcription levels and assay consistency (Fig. 5b, Additional file 2: Fig. S6). Regions lacking detectable transcription signals exhibited the lowest Jaccard Index values across all pairwise comparisons (Fig. 5b), indicating poor reproducibility across assays. Conversely, high-transcription regions exhibited the highest assay consistencies (Fig. 5b, Additional file 2: Fig. S6d). These results support the hypothesis that these reporter assays may yield a greater proportion of false positives in regions with lower transcription.

Transcription enhances the predictive power of biochemical features for enhancer activity

Next, we assessed whether transcription improves the ability of biochemical features to predict active enhancers. We analyzed tested cCREs with high overlap (≥ 80% reciprocal overlap) with reporter assay regions and classified them as either transcribed or untranscribed based on detectable GRO-cap signals. We then compared their active rates across assays.

Untranscribed cCREs exhibited low active rates in all three assays (~ 0.8–10%), with active rates only slightly higher than untranscribed regions that lacked cCRE or PINTS annotations (~ 0.3–3%) (Fig. 5c). Untranscribed dELS, pELS, and PLS showed slightly elevated active rates (~ 0–16%), though their sample sizes were limited.

In contrast, transcribed cCREs displayed significantly higher active rates across all assays (~ 24–33%) (Fig. 5c). This trend was particularly pronounced for transcribed dELS, pELS, and PLS, which exhibited much higher active rates (~ 28–75%) than their untranscribed counterparts (Fig. 5c). These results indicate that dELS, pELS, and PLS contain a higher proportion of functional enhancers than other cCRE categories and suggest that transcription serves as an additional predictive layer beyond traditional biochemical features such as chromatin accessibility and histone modifications (H3K4me3 and H3K27ac).

Further stratification of tested dELS, pELS, and PLS by transcription levels reinforced the strong relationship between transcription and enhancer activity across all assay types (Fig. 5d). Highly transcribed dELS, pELS, and PLS exhibited particularly high active rates, reaching 83% in ATAC-STARR-seq and 73% in WHG-STARR-seq (Fig.