Method Article

A Computational Pipeline for Intergenic/Intragenic Enhancer RNA Quantification in Mouse Embryonic Stem Cells

DOI:

10.3791/69400

⸱

October 28th, 2025

 ,  , 
1Department of Convergent Bioscience and Informatics, and Graduate School of Biological Sciences, Chungnam National University (CNU)

* These authors contributed equally

In This Article

Summary

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

This protocol provides a streamlined computational pipeline for quantifying nascent enhancer transcripts. By integrating chromatin accessibility, chromatin feature, and transcriptional data, it enables accurate detection and strand-specific analysis of enhancer activity in complex intragenic regions, while remaining accessible to researchers without extensive bioinformatics training.

Abstract

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

Core cis-regulatory elements known as enhancers play a central role in enabling precise transcriptional regulation of target genes that control diverse cellular functions and developmental processes. These enhancers are often transcribed in both directions, producing long non-coding transcripts referred to as enhancer RNAs (eRNAs). The expression of eRNAs is closely linked to active chromatin features, such as H3K27ac and co-activator recruitment, and functionally contributes to the transcriptional activation of target genes. Nevertheless, the detection and quantification of eRNAs remain challenging, especially when they overlap with host-gene transcription. To address this, we present a standardized, user-friendly computational workflow for analyzing enhancer transcription from nascent RNA sequencing data. The protocol guides users through data preprocessing, read mapping, and quality control, followed by strand-specific quantification of enhancer-associated transcription, with dedicated procedures for intragenic enhancers where signal assignment is complex. Visualization modules enable clear inspection of enhancer activity across genomic contexts, and built-in options support analyses of both intergenic and intragenic enhancers. Designed for researchers with limited bioinformatics expertise, this workflow provides a practical framework for consistent, reproducible, and scalable studies of enhancer transcription, facilitating broader application of enhancer biology across diverse systems.

Introduction

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

Enhancers are cis-regulatory DNA elements that control target gene transcription by organizing chromatin looping and recruiting the transcriptional machinery1,2,3. Their tissue-specific activity enables precise regulation during development and lineage commitment4,5,6,7,8. Active enhancers show characteristic chromatin features such as H3K4me1 (histone H3 Lysine 4 mono-methylation) and H3K27ac (histone H3 Lysine 27 acetylation) and are usually found in DNase I hypersensitive regions that mark open chromatin9,10,11,12. These features enable transcription factors and RNA polymerase II to access the DNA, initiating nascent transcription at enhancer loci13,14,15,16.

This sequential biological process produces enhancer-derived transcripts, called eRNAs, which are bidirectional, noncoding, and typically non-polyadenylated RNAs13,14,15,16. eRNAs serve as markers of enhancer activity and function as effectors in their own right16,17,18,19,20,21,22,23,24. They promote productive elongation by releasing Negative Elongation Factor (NELF) from paused RNA polymerase II16,19, and help stabilize enhancer-promoter loops17,18,20. They also support the formation of transcriptional condensates, potentially via m6A (N6-Methyladenosine) modification21,22,23.

Still, the function of intragenic enhancer transcription, initiated from regulatory elements within gene bodies, remains controversial. Some studies report that eRNAs from intragenic enhancers augment host gene expression25, potentially by promoting NELF release and stimulus-dependent productive elongation26,27. In contrast, other work suggests that this transcription can impede host genes via RNA polymerase II collisions or transcriptional interference, leading to attenuation or premature termination28,29. These conflicting observations, together with the dual role of eRNAs as markers and regulators, highlight the need for careful quantification and functional dissection. Yet, measuring intragenic eRNAs is difficult because they often overlap sense-strand host transcripts13,25,26,30. The challenge is amplified when enhancers reside in regions with nested genes or overlapping transcription on both strands, which obscures enhancer-specific signals.

To overcome these challenges, we developed a bioinformatics pipeline to detect, quantify, and visualize enhancer-associated transcripts, with a particular focus on intragenic regions. The pipeline integrates the Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq), Chromatin Immunoprecipitation Sequencing (ChIP-seq), Global Run-on Sequencing (GRO-seq), and genomic annotations to achieve enhancer-level resolution even in complex genomic contexts.

The pipeline comprises four major steps: (i) preprocessing, alignment, peak calling, and signal generation31; (ii) enhancer identification using chromatin features; (iii) strand orientation assignment, particularly within gene bodies; and (iv) quantification and visualization of nascent enhancer transcripts. This framework is particularly useful for systems with high-resolution sequencing data, such as the mouse embryonic stem cells analyzed in this study, and it can be extended to other organisms when suitable datasets are available. By enabling enhancer-specific quantification where existing pipelines fall short, this workflow offers a practical tool for benchmarking and studying intragenic eRNA transcription across diverse genomic contexts.

Access restricted. Please log in or start a trial to view this content.

Protocol

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

NOTE: All raw datasets used in workflow are listed in Table 1. Details of bioinformatics tools are provided in the Table of Materials. The number of threads used in this pipeline can be adjusted by modifying the THREADS variable defined at the top of each script. Users may increase the number to accelerate analysis depending on the user's CPU resources.
After each step, a log file is generated. For quick failure checks, use commands like cat StepXX_log.txt; grep -qF "ERROR" StepXX_log.txt && echo "ERROR found. Fix before next step." || echo "OK: no ERROR markers". If any ERROR appears, treat the step as failed and resolve it first.

1. Downloading of full analysis pipeline from GitHub repository

(https://github.com/myunggeunO/Enhancer-transcript-identification-from-read-to-visualization/)

  1. Launch the Command-Line Interface (CLI) appropriate for the operating system used.
    1. Windows: Set up a Linux environment using Windows Subsystem for Linux (WSL). Follow the official instructions to install and configure WSL before proceeding32.
    2. macOS: Continue without additional setup, as macOS is Unix-based. Refer to the official guide to open the terminal33.
    3. Linux users, particularly those using Ubuntu: Open a terminal as described in the referenced instructions34.
  2. Run wget https://github.com/myunggeunO/Enhancer-transcript-identification-from-read-to-visualization/archive/refs/heads/main.zip -O ~/pipeline.zip   in terminal to download pipeline for enhancer identification and enhancer RNA quantification.
  3. Type unzip ~/pipeline.zip -d ~/ in terminal. This will extract all necessary files into the home directory.
  4. Run rm ~/pipeline.zip and type mv ~/Enhancer-transcript-identification-from-read-to-visualization-main ~/Enhancer-transcript-identification-from-read-to-visualization to remove archive and rename extracted folder.
  5. Type cd ~/Enhancer-transcript-identification-from-read-to-visualization/, and run chmod +x scripts/* to make all scripts in "scripts/" directory executable.

2. Setting up mamba/conda environment for analysis pipeline

  1. Type bash scripts/Step1_conda_environment_formation.sh to create and run a mamba virtual environment. If prompted during execution, type Y and press Enter to confirm package installations. For macOS, follow Step 2.1.1; for systems with Mamba or Conda already installed, follow Step 2.1.2.
    1. MacOS: Open the script and replace miniconda download link with macOS version:
      https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh
      Then, follow step 2.1.
    2. Once (enhancer-env) appears in the prompt, type bash scripts/Step2_package_installation.sh to install required packages for downstream analyses. Type Y and press Enter if prompted during installation.
    3. After running Step 2.2, check the terminal output for any error messages. Resolve any issues; then re-run Step 2.2.
    4. (OPTIONAL) Run mamba list to confirm all mamba-managed packages in the Table of Materials are installed. HOMER is installed manually and will not appear in the mamba list. Verify HOMER by checking that "~/homer/" directory exists.

3. Download publicly available ChIP-seq, ATAC-seq, and GRO-seq datasets from SRA (Sequence Read Archive)

  1. Run cp scripts/Step{3..12}_*.sh ./ to copy necessary shell scripts for raw read processing.
  2. Type bash Step3_download_file_list.sh > Step3_log.txt 2>&1 and press the Enter key to download and process raw sequencing data from SRA.
    NOTE: This script automates download and preparation of public sequencing data for analysis. It creates a standardized folder structure under "MATERIAL/", organized by assay type and replicate (biological: rep1/rep2; technical: trep1/trep2). A built-in list of SRA accession numbers drives retrieval with prefetch (v3.2.0), conversion to FASTQ with fasterq-dump (v3.2.0) (paired-end: --split-files), and compression with pigz (v2.8) to reduce storage. Processing follows prefetch, fasterq-dump, pigz, with outputs written to corresponding "00.Rawdata/" directories.

4. Perform quality control and trimming of raw reads

  1. Type bash Step4_read_trimming_and_QC.sh > Step4_log.txt 2>&1 to perform read trimming and quality control of raw FASTQ files.
    NOTE: This script processes raw FASTQ files from GRO-seq, ATAC-seq, and ChIP-seq (H3K27ac, H3K4me1) with corresponding input controls. It runs FastQC (v0.12.1)35 on raw reads, then trims adapters with Trim Galore (v0.6.10)36 using assay-specific parameters. For GRO-seq, it first removes NextSeq G tails and very short reads (--nextseq 20, --length 20), then uses Cutadapt (v5.1)37 to strip long poly-A tracts while retaining reads longer than 20 nucleotides ( -a A{15}, -m 20). For ATAC-seq, it processes paired-end libraries and targets Tn5/Nextera adapters (--paired, --nextera). For H3K27ac and corresponding input, it handles paired-end ChIP-seq libraries (--paired). For H3K4me1 and its input, it performs standard single-end trimming (default). FastQC is run again on trimmed reads. Outputs are written to each assay's "01.Clean/" directory, and intermediate GRO-seq adapter-trimmed files are removed.

5. Prepare Bowtie2 reference index

  1. Choose one of the two options below to prepare Bowtie2 (v2.5.4)38 genome index for mm10 (Mus musculus) reference genome.
    1. Run bash Step5_1_download_reference_index.sh > Step5_1_log.txt 2>&1 to use pre-built Bowtie2 index provided by developers.
    2. Run bash Step5_2_download_reference_make_index_with_
      bowtie2.sh > Step5_2_log.txt 2>&1       to download raw mm10 genome sequence and build index manually.
      ​NOTE: Both approaches generate an index file to "reference_index/" directory and are functionally equivalent for standard alignment.

6. Align trimmed reads to mm10 reference genome

  1. Type bash Step6_alignment_to_make_bam.sh > Step6_log.txt 2>&1 to align each assay's reads to reference and generate BAM files.
    ​NOTE: This step maps each assay's reads to previously indexed mm10 reference. H3K27ac ChIP-seq and corresponding input use paired-end mapping for balanced accuracy (-1, -2) with default sensitivity. H3K4me1 ChIP-seq and corresponding input use single-end mapping under default settings (-U). ATAC-seq uses high-sensitivity mapping to accommodate variable, long Tn5-derived fragments (--very-sensitive, -X 2000) with paired-end input (-1, -2). GRO-seq uses high-sensitivity mapping to better place short, trimmed reads (--very-sensitive) with single-end input (-U). SAM files are converted to BAM and filtered with samtools (v1.22.1)39 view, using moderate MAPQ for ChIP/input (-b, -q 10) and stricter thresholds for ATAC-seq and GRO-seq (-b, -q 30); final BAM files are written to each dataset's "02.Align/" directory.

7. Merge technical replicates of H3K27ac ChIP-seq data

  1. Run bash Step7_merge_trep.sh > Step7_log.txt 2>&1 to merge technical replicates of H3K27ac ChIP-seq and corresponding input BAM files.
    ​NOTE: This script sorts BAM files for technical replicates with sambamba (v1.0.1)40 sort and then merges H3K27ac ChIP-seq and corresponding input replicates into consolidated BAMs with sambamba merge. If the user's dataset does not contain any technical replicates, skip this step and proceed with individual BAM files.

8. Remove duplicates and non-essential chromosomes

  1. Remove duplicates and sort mapped reads for ChIP-seq and GRO-seq.
    1. Type bash Step8_1_duplicate_removal_sorting-ChIP_GRO.sh > Step8_1_log.txt 2>&1 to remove duplicates from histone ChIP-seq datasets and sort both ChIP-seq and GRO-seq output BAMs.
      ​NOTE: This script removes PCR duplicates with sambamba markdup (-r) and coordinate-sorts with sambamba sort. For H3K27ac, processing targets BAMs from merged technical replicates, handling ChIP and input separately. For H3K4me1, each replicate and matched input are processed individually. For GRO-seq, duplicate removal is skipped and only coordinate sorting is applied. Outputs are saved in each dataset's "02.Align/" directory.
  2. Remove duplicates and filter mitochondrial chromosome-mapped reads from ATAC-seq.
    1. Type bash Step8_2_duplicate_chrM_removal_sorting_ATAC.sh > Step8_2_log.txt 2>&1 to remove PCR duplicates, filter mitochondrial reads (chrM), and sort ATAC-seq BAMs.
      ​NOTE: This script references ENCODE pipeline recommendations for ATAC-seq data processing. It begins by performing name sorting with sambamba sort (-n), and fixing mate-pair information using samtools fixmate (-m). This step ensures that mate information is properly assigned prior to marking duplicates. Afterwards, PCR duplicates are removed using sambamba markdup (-r). Mitochondrial reads are removed by generating a keep-list from samtools idxstats (exclude chrM and *) and retaining only listed references with samtools view (-b). Final coordinate sorting is performed with sambamba sort. Cleaned BAMs are written to each replicate's "02.Align/" directory.

9. Perform peak calling for each dataset

  1. Run step 9 script by typing bash Step9_peak_calling.sh > Step9_log.txt 2>&1 to perform peak calling for ChIP-seq and ATAC-seq data.
    NOTE: This script performs peak calling with MACS3 (v3.0.3)41 for ATAC-seq and ChIP-seq (H3K27ac, H3K4me1). ATAC-seq generates peaks by supplying all replicate BAMs as signal in no-model mode with shift/extension (-f BAMPE, --nomodel, --shift -100, --extsize 200, -q 0.01). H3K27ac calls broad peaks from merged technical replicates with matched input (-f BAMPE, --broad). H3K4me1 processes each replicate individually with its matched input in single-end broad mode (-f BAM, --broad), and overlapping peaks are obtained with bedtools (v2.31.1)42 intersect. Outputs are written to each dataset's "peak_calling/" directory, with high-confidence H3K4me1 peaks under "peak_calling/overlapped_peak/".

10. Merge biological replicates of ChIP-seq and ATAC-seq BAM files for downstream signal analysis

  1. Run bash Step10_merge_rep_forMakingSignal.sh > Step10_log.txt 2>&1 to merge BAM files from biological replicates of ATAC-seq and H3K4me1 ChIP-seq.
    ​NOTE: This script combines replicate BAMs with sambamba merge for ATAC-seq, H3K4me1 ChIP-seq, and corresponding H3K4me1 input. Merged BAMs support downstream analyses (e.g., bigWig signal generation, normalization). Output BAMs are saved in "merge/02.Align/" directories under each sample path.

11. Generate tag directories and signal bigWig files from mapped reads

  1. Run bash Step11_make_tag_to_signal.sh > Step11_log.txt 2>&1 to create tag directories and generate bigWig signal files for mapped reads of each dataset.
    ​NOTE: This script builds HOMER tag directories with makeTagDirectory and then generates bigWig signal tracks with makeUCSCfile command using HOMER (v5.1)43 and ucsc-bedgraphtobigwig (v482)44 packages. All signal tracks are created using chromosome size files from UCSC Genome Browser (https://hgdownload.soe.ucsc.edu/goldenPath/mm10/). GRO-seq produces strand-specific signal tracks (-style rnaseq, -strand + / -, -bigWig). ATAC-seq produces unnormalized tracks from merged BAMs (-bigWig). ChIP-seq (H3K27ac, H3K4me1) produces input-normalized tracks with a pseudo-count of 1 (-bigWig, -i, -pseudo 1). Outputs are organized under "03.TagDir/" and "04.bigwig/".

12. Prepare files for enhancer identification

  1. Type bash Step12_E_identification_material.sh > Step12_log.txt 2>&1 in the terminal to prepare necessary files and reference for enhancer identification.
    NOTE: This step collects all necessary files for enhancer identification into "01.E_identification/material/", organized into "ATAC/", "Histone/", and "Annotation/" folders. It copies peak files (ATAC-seq, H3K27ac, H3K4me1) into appropriate directories. GENCODE M23 annotation file (mm10) is downloaded automatically, and reference files, including ENCODE blacklist45 (mm10-blacklist.v2.bed) and chromosome size file (mm10.chrom.sizes) are copied from a predefined path.

13. Identify promoter candidates and gene bodies from annotation

  1. Type cd 01.E_identification/ to enter working directory, then run cp ../scripts/Step{13..20}_*.sh ./ to copy enhancer-identification scripts.
  2. Run bash Step13_promoter_candidates_genebody_identification.sh > Step13_log.txt 2>&1 to generate BED files for promoter regions, gene bodies, and Protein-Coding Genes (PCGs) using GENCODE annotation.
    NOTE: This script processes downloaded GENCODE M23 Gene Transfer Format (GTF) to create BED files for promoter candidates, all gene bodies, and PCG bodies, saving outputs to "material/Annotation/". Promoters are defined as a 2 kb window around each transcript's TSS (Transcription Start Site) with bedtools slop (-b 2000, -g mm10.chrom.sizes). Gene and protein-coding bodies are derived from GTF entries annotated as gene, with protein-coding entries further filtered by "gene_type = protein_coding".

14. Process peaks for enhancer identification

  1. Run bash Step14_ATAC_ChIP-seq_processing.sh > Step14_log.txt 2>&1 to preprocess ATAC-seq and histone ChIP-seq peak files for enhancer identification.
    NOTE: This script preprocesses peak sets for enhancer calling. For ATAC-seq, it removes regions overlapping blacklist regions with bedtools subtract (-A) and then excludes peaks overlapping promoter candidates with bedtools subtract (-A); for histone marks (H3K27ac, H3K4me1), it symmetrically expands each peak by 1 kb on each side with bedtools slop (-b 1000 -g mm10.chrom.sizes) and then removes promoter overlaps with bedtools subtract.

15. Identify and classify enhancers

  1. Run bash Step15_inter_intragenic_E_sets_identification.sh > Step15_log.txt 2>&1 to define and classify enhancers using chromatin peak data.
    ​NOTE: This step defines and annotates enhancers using bedtools. Overlaps between ATAC-seq peaks and slop-expanded H3K4me1 peaks are obtained with intersect (-wa, -u), and regions that also overlap flanked H3K27ac peaks are classified as active enhancers with intersect (-wa, -u); non-active enhancers are derived by removing active regions from the full enhancer set with subtract. ATAC-seq summits overlapping each enhancer set are collected with intersect (-u), then summits are split into intergenic and intragenic with intersect (-v or -u) against the gene body. Enhancer intervals are finally assigned to intergenic/intragenic classes based on the peak's associated summit with intersect (-u). All results are saved in "01.E_identification/" under "01.allE/", "02.interE/", and "03.intraE/".

16. Assign temporary strand information to intragenic enhancers

  1. Type bash Step16_assign_temp_strand_from_gene_overlap.sh > Step16_log.txt 2>&1 to assign temporary strand information to intragenic-enhancer BED.
    NOTE: This step assigns gene-strand labels to intragenic enhancers by overlapping enhancer intervals with gene bodies using bedtools intersect (-wa, -wb). Columns 1, 2, 3, 4, 5, 16 are retained with awk, then the records are sorted and deduplicated. Output is saved as "03.intraE/strand_designation/01.overlapped_gene_strand/ES_E_intragenic_strand_with_dup.bed".

17. Prioritize strand assignment for enhancers overlapping both-strand genes

  1. Type bash Step17_initial_strand_assignment_for_both_strand_enhancers_PCG_based.sh > Step17_log.txt 2>&1 to resolve strand direction for intragenic enhancers overlapping genes on both strands.
    NOTE: This script resolves strand ambiguity for intragenic enhancers overlapping genes on both strands by prioritizing PCG overlaps. Enhancers present on both strands are first isolated (awk grouping by chrom/start/end/id/strand), PCG-overlapping cases on the same strand are selected with bedtools intersect (-s, -wa, -u), and non-PCG cases are retained with bedtools intersect (-v). Selected and retained sets are concatenated and ordered. Output: "03.intraE/strand_designation/02.enhancer_with_PCG_priority/ES_E_intragenic_PCG_priority.bed".

18. Calculate strand-specific Reads Per Kilobase per Million mapped reads (RPKM) values for genes overlapping PCG-prioritized intragenic enhancers

  1. Type bash Step18_RPKM_cal_from_partially_strand_assigned_enhancers.sh > Step18_log.txt 2>&1 to calculate strand-specific RPKM for genes overlapping enhancers on same strand.
    ​NOTE: This step uses "ES_E_intragenic_PCG_priority.bed" from Step 17, which contains enhancers that (i) overlapped a PCG on one strand and received a strand, (ii) overlapped PCGs on both strands and remained ambiguous, or (iii) had no PCG overlap and retained both strands. Same-strand overlapping genes are selected with bedtools intersect (-s, -wa, -u), converted to GTF via awk, and quantified from GRO-seq with featureCounts (v2.1.1)46 strand-specific count mode (-s 1, -t gene, -g gene_id, -O, --fraction). Total mapped reads are obtained with sambamba flagstat, and RPKM is computed from gene length, counts, and totals. Outputs are written to "03.intraE/strand_designation/03.RPKM_calculation_of_overlapped_gene/".

19. Final strand assignment based on gene expression (RPKM) of overlapping genes

  1. Type bash Step19_second_strand_assignment_by_RPKM.sh > Step19_log.txt 2>&1 to finalize strand assignment for intragenic enhancers.
    ​NOTE: This step assigns strands to intragenic enhancers using gene-expression support from Step 18. Same-strand overlaps between enhancers (ES_E_intragenic_PCG_priority.bed) and genes are found with bedtools intersect (-s, -wa, -wb). Gene RPKM values are joined to gene intervals via awk/sort/join, producing gene-RPKM BED. For each enhancer, the overlapping gene with the highest RPKM is selected, and the enhancer inherits that gene's strand. Results are saved to "03.intraE/strand_designation/04.enhancer_strand_designation_by_RPKM_of_gene/ES_E_intragenic_PCG_priority_strand_by_gene_RPKM.bed."

20. Assign strand information to intragenic enhancers and enhancer summits

  1. Type bash Step20_strand_assignment_for_intragenicE_and_summits.sh > Step20_log.txt 2>&1 to assign decided strand information to all intragenic enhancers and each summit.
    NOTE: This step finalizes strand assignment for intragenic enhancers and matching summits using bedtools. Opposite-strand intervals are removed with subtract (-S), strand-designated enhancers are intersected with active and non-active sets using intersect (-wa, -u), and summit files are re-annotated by overlapping summits with strand-designated enhancers via intersect (-wa, -wb) and get strand with awk. Results are saved in "03.intraE/final_strand_IntragenicE/" and its "summit/" subfolder".

21. Prepare input files for enhancer validation, eRNA quantification, and visualization

  1. Type cd ../ or cd ~/Enhancer-transcript-identification-from-read-to-visualization to move to pipeline root, then type cp scripts/Step21_preparing_quantification_and_visualization.sh ./ to copy the script for preparation of downstream analysis.
  2. Run bash Step21_preparing_quantification_and_visualization.sh > Step21_log.txt 2>&1 to prepare all necessary files for enhancer aggregation, GRO-seq signal processing, and eRNA quantification.
    NOTE: This step prepares input files and directory structure for enhancer validation, eRNA quantification, and signal visualization under "02.E_visualization_quantification/". Subfolders are created for bigWig files, enhancer BEDs, BAMs, signal matrices, counts, and plots. Key inputs such as summit BEDs, bigWigs, enhancer lists, and GRO-seq BAMs are copied to appropriate locations.

22. Generate aggregation plots for enhancer validation

  1. Type cd 02.E_visualization_quantification/ to enter working directory, then type cp ../scripts/Step{22..24}_*.* ./ to copy necessary scripts for downstream analysis.
  2. Run bash Step22_generation_of_aggregation_plot.sh > Step22_log.txt 2>&1 to generate aggregation plots of chromatin signals around each type of enhancer summit.
    NOTE: This step visualizes average chromatin signal enrichment centered on enhancer summits using computeMatrix and plotProfile from deepTools (v3.5.6)47. For each defined enhancer set, computeMatrix reference-point (--referencePoint center, -a 5000, -b 5000, --missingDataAsZero) calculates signal density within a 10 kb window around enhancer summits using bigWig files for ATAC-seq, H3K27ac, and H3K4me1. Output matrix is passed to plotProfile, which generates signal aggregation curves for comparison between enhancer groups. Aggregation plots are saved in "01.Profiling/04_1.aggregation/" directory.

23. Quantify and visualize enhancer RNA expression

  1. Run bash Step23_quantifing_eRNA_RPKM.sh > Step23_log.txt 2>&1 to quantify eRNA expression levels from GRO-seq using featureCounts.
    NOTE: This step quantifies eRNA transcription from defined enhancer regions in a strand-specific manner using GRO-seq. Intergenic enhancers are counted with featureCounts in unstranded mode (-s 0, -t enhancer, -g gene_id, -O, --fraction), and intragenic enhancers are quantified in antisense mode (-s 2) to exclude signal from overlapping gene transcription. BED regions are converted to GTF with awk before counting. Total mapped reads come from sambamba flagstat, and counts are normalized to RPKM using enhancer length, read counts, and total mapped reads. Outputs are organized under "02.eRNA_quantification/03.count_normalized_with_RPKM/" by "inter/" and "intra/".
  2. Run Rscript Step24_visualization_of_enhancer_transcript.R > Step24_log.txt 2>&1 to visualize and compare eRNA expression levels across enhancer groups using R.
    NOTE: R script uses ggplot2 (v3.5.2)48 and cowplot (v1.2.0)49 packages to generate violin and box plots comparing active and non-active enhancer expression based on RPKM values. For visualization and statistical testing, RPKM values are transformed to log2 (RPKM + 1). Statistical significance is assessed using the Wilcoxon rank-sum test. Both summary plots and p-value table are saved to "03.eRNA_visualization/" for downstream interpretation.
    NOTE: If any step in this pipeline fails and persists even after re-run, report the issue at https://github.com/myunggeunO/Enhancer-transcript-identification-from-read-to visualization/issues. Clearly stating the failed step and attaching the log file ensures accurate troubleshooting support.

Access restricted. Please log in or start a trial to view this content.

Results

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

Schematic workflow for enhancer transcript quantification pipeline
Publicly available ChIP-seq (H3K27ac, H3K4me1), ATAC-seq, and GRO-seq datasets (Table 1) were processed with a standardized pipeline designed primarily for validation. Adapter trimming and quality filtering were performed with Trim Galore and Cutadapt, followed by alignment to the mm10 reference genome using Bowtie2 (detailed in protocol Step 6). For ChIP-seq and ATAC-seq, peaks were i...

Access restricted. Please log in or start a trial to view this content.

Discussion

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

Following the discovery of enhancer-derived transcripts13,14,15, accurately quantifying eRNAs has remained a major challenge, particularly in intragenic contexts where eRNAs often overlap with host gene transcripts. This overlap complicates strand assignment and signal attribution, making it difficult to distinguish genuine enhancer transcription from background gene expression13,

Access restricted. Please log in or start a trial to view this content.

Disclosures

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

The authors have no conflicts of interest to disclose.

Acknowledgements

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

This study was supported by the research fund of Chungnam National University [2022-0582-01 (S.-K.K.) and 2023-0545-01 (S.-K.K.)], South Korea. Figure 1 was created using BioRender (https://biorender.com/).

Access restricted. Please log in or start a trial to view this content.

Materials

List of materials used in this article
NameCompanyCatalog NumberComments
bedtoolsQuinlan Lab, University of Utah v2.31.1Utilities for editing BED files
bowtie2Langmead Lab, Johns Hopkins Universityv2.5.4Multi-threaded aligner for mapping reads to a reference genome
cowplotWilke Lab, University of Texasv1.2.0Tools for combining and aligning ggplot2-based figures
cutadaptScience For Life Laboratory, Stockholm Universityv5.1Adaptor and poly-A/G tail trimmer
deeptoolsBioinformatics Facility, Max Planck Institutev3.5.6Read counting tool for quantifying reads in defined genomic regions
fastqcBabraham Bioinformatics, Babraham Institutev0.12.1Quality control for sequencing reads
featureCounts (subread)Shi Lab, Monash Universityv2.1.1Raw read counting tools for specified genomic regions
homerBenner Lab, University of California San Diego (UCSD)v5.1Toolkit for ChIP-seq, ATAC-seq, and nascent RNA analysis; includes tag directory creation and signal profiling
macs3Chan Zuckerberg Initiativev3.0.3Peak calling for ChIP-seq and ATAC-seq datasets
pigz.v2.8Multi-threaded compression tool for generating gzip-compressed files
sambambaPetersburg State Universityv1.0.1Multi-threaded SAM/BAM file processing toolkit
samtoolsWellcome Trust Sanger Institutev1.22.1Tools for processing and manipulating SAM/BAM files
sra-toolsNational Center for Biotechnology Information (NCBI)v3.2.0For downloading SRR files from the NCBI SRA database
tidyversePosit PBCv2.0.0Collection of R packages for data manipulation and visualization
trim-galoreAltos Labs, Cambridge Institute of Sciencev0.6.10Adapter and low-quality base trimming using multi-threading
Ubuntu 20.04Developing and testing the pipeline
ucsc-bedgraphtobigwig Kent Lab, University of California Santa Cruzev482Tools for generating bigWig signal tracks

References

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,
  1. Bulger, M., Groudine, M. Functional and mechanistic diversity of distal transcription enhancers. Cell. 144 (3), 327-339 (2011).
  2. Smith, E., Shilatifard, A. Enhancer biology and enhanceropathies. Nat Struct Mol Biol. 21 (3), 210-219 (2014).
  3. Li, W., Notani, D., Rosenfeld, M. G. Enhancers as non-coding RNA transcription units: recent insights and future perspectives. Nat Rev Genet. 17 (4), 207-223 (2016).
  4. Whyte, W. A., et al. transcription factors and mediator establish super-enhancers at key cell identity genes. Cell. 153 (2), 307-319 (2013).
  5. Plank, J. L., Dean, A. Enhancer function: mechanistic and genome-wide insights come together. Mol Cell. 55 (1), 5-14 (2014).
  6. Alexander, J. M., et al. Brg1 modulates enhancer activation in mesoderm lineage commitment. Development. 142 (8), 1418-1430 (2015).
  7. Huang, J., et al. Dynamic control of enhancer repertoires drives lineage and stage-specific transcription during hematopoiesis. Dev Cell. 36 (1), 9-23 (2016).
  8. Xiong, L., et al. Genome-wide identification and characterization of enhancers across 10 human tissues. Int J Biol Sci. 14 (10), 1321-1332 (2018).
  9. Heintzman, N. D., et al. Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature. 459 (7243), 108-112 (2009).
  10. Creyghton, M. P., et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci. 107 (50), 21931-21936 (2010).
  11. Calo, E., Wysocka, J. Modification of enhancer chromatin: what, how, and why. Mol Cell. 49 (5), 825-837 (2013).
  12. Barakat, T. S., et al. Functional dissection of the enhancer repertoire in human embryonic stem cells. Cell Stem Cell. 23 (2), 276-288 (2018).
  13. Kim, T. -K., et al. Widespread transcription at neuronal activity-regulated enhancers. Nature. 465 (7295), 182-187 (2010).
  14. Melgar, M. F., Collins, F. S., Sethupathy, P. Discovery of active enhancers through bidirectional expression of short transcripts. Genome Biol. 12 (11), R113(2011).
  15. Djebali, S., et al. Landscape of transcription in human cells. Nature. 489 (7414), 101-108 (2012).
  16. Gorbovytska, V., et al. Enhancer RNAs stimulate Pol II pause release by harnessing multivalent interactions to NELF. Nat Commun. 13 (1), 2429(2022).
  17. Mousavi, K., et al. eRNAs promote transcription by establishing chromatin accessibility at defined genomic loci. Mol Cell. 51 (5), 606-617 (2013).
  18. Hsieh, C. -L., et al. Enhancer RNAs participate in androgen receptor-driven looping that selectively enhances gene activation. Proc Natl Acad Sci. 111 (20), 7319-7324 (2014).
  19. Schaukowitch, K., et al. Enhancer RNA facilitates NELF release from immediate early genes. Mol Cell. 56 (1), 29-42 (2014).
  20. Pnueli, L., Rudnizky, S., Yosefzon, Y., Melamed, P. RNA transcribed from a distal enhancer is required for activating the chromatin at the promoter of the gonadotropin α-subunit gene. Proc Natl Acad Sci. 112 (14), 4369-4374 (2015).
  21. Sabari, B. R., et al. Coactivator condensation at super-enhancers links phase separation and gene control. Science. 361 (6400), eaar3958(2018).
  22. Nair, S. J., et al. Phase separation of ligand-activated enhancers licenses cooperative chromosomal enhancer assembly. Nat Struct Mol Biol. 26 (3), 193-203 (2019).
  23. Lee, J. -H., et al. Enhancer RNA m6A methylation facilitates transcriptional condensate formation and gene activation. Mol Cell. 81 (16), 3368-3385 (2021).
  24. Chen, Q., et al. Enhancer RNAs in transcriptional regulation: recent insights. Front Cell Dev Biol. 11, 1205540(2023).
  25. Moon, J., et al. Embryonic stem cell-specific intragenic enhancer RNA essential for NSUN2-mediated stem cell fate regulation. Int J Biol Macromol. 245, 470(2025).
  26. Tuvikene, J., et al. Intronic enhancer region governs transcript-specific Bdnf expression in rodent neurons. Elife. 10, e65161(2021).
  27. Cheng, F., et al. Intronic enhancers of the human SNCA gene predominantly regulate its expression in brain in vivo. Sci Adv. 8 (47), eabq6324(2022).
  28. Hobson, D. J., Wei, W., Steinmetz, L. M., Svejstrup, J. Q. RNA polymerase II collision interrupts convergent transcription. Mol Cell. 48 (3), 365-374 (2012).
  29. Cinghu, S., et al. Intragenic enhancers attenuate host gene expression. Mol Cell. 68 (1), 104-117 (2017).
  30. Bressin, A., et al. High-sensitive nascent transcript sequencing reveals BRD4-specific control of widespread enhancer and target gene transcription. Nat Commun. 14 (1), 4971(2023).
  31. Lee, J., et al. Introductory analysis and validation of CUT&RUN sequencing data. J Vis Exp. (214), e67359(2024).
  32. How to install Linux on Windows with WSL. , Microsoft. https://learn.microsoft.com/en-us/windows/wsl/install (2025).
  33. Terminal user guide. , Apple. https://support.apple.com/guide/terminal/welcome/mac (2025).
  34. How to open terminal in Linux. , GeeksforGeeks. https://www.geeksforgeeks.org/linux-unix/how-to-open-terminal-in-linux/ (2025).
  35. Simon, A. FastQC: a quality control tool for high throughput sequence data. Version 0.10.1, (2010).
  36. Krueger, F. Trim Galore!: a wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files, with extra functionality for RRBS data. Babraham Inst. , (2015).
  37. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17 (1), 3(2011).
  38. Langmead, B., Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat Methods. 9 (4), 357-359 (2012).
  39. Li, H., et al. The sequence alignment/map format and SAMtools. Bioinformatics. 25 (16), 2078-2079 (2009).
  40. Tarasov, A., Vilella, A. J., Cuppen, E., Nijman, I. J., Prins, P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 31 (12), 2032-2034 (2015).
  41. Zhang, Y., et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9 (9), R137(2008).
  42. Quinlan, A. R., Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 26 (6), 841-842 (2010).
  43. Heinz, S., et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 38 (4), 576-589 (2010).
  44. Kent, W. J., Zweig, A. S., Barber, G., Hinrichs, A. S., Karolchik, D. BigWig and BigBed: enabling browsing of large distributed datasets. Bioinformatics. 26 (17), 2204-2207 (2010).
  45. Amemiya, H. M., Kundaje, A., Boyle, A. P. The ENCODE blacklist: identification of problematic regions of the genome. Sci Rep. 9 (1), 9354(2019).
  46. Liao, Y., Smyth, G. K., Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 30 (7), 923-930 (2014).
  47. Ramírez, F., et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160(2016).
  48. Wickham, H. ggplot2: elegant graphics for data analysis. , Springer. 189-201 (2016).
  49. Wilke, C. O. cowplot: streamlined plot theme and plot annotations for ggplot2. CRAN Contrib. Packages. , (2015).
  50. Andersson, R., et al. An atlas of active enhancers across human cell types and tissues. Nature. 507 (7493), 455-461 (2014).
  51. Spicuglia, S., Vanhille, L. Chromatin signatures of active enhancers. Nucleus. 3 (2), 126-131 (2012).
  52. Zentner, G. E., Tesar, P. J., Scacheri, P. C. Epigenetic signatures distinguish multiple classes of enhancers with distinct cellular functions. Genome Res. 21 (8), 1273-1283 (2011).
  53. Blinka, S., Reimer, M. H., Pulakanti, K., Rao, S. Super-enhancers at the Nanog locus differentially regulate neighboring pluripotency-associated genes. Cell Rep. 17 (1), 19-28 (2016).
  54. Zhao, S., Ye, Z., Stanton, R. Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols. RNA. 26 (8), 903-909 (2020).

Access restricted. Please log in or start a trial to view this content.

Reprints and Permissions

Request permission to reuse the text or figures of this JoVE article

Request Permission

Tags

Intergenic EnhancersIntragenic EnhancersGRO seq AnalysisATAC seq DataH3K27 AcetylationChromatin Peak DataStrand Specific QuantificationAggregation Plots

Related Articles