方法文章

小鼠胚胎干细胞中基因间/基因内增强子RNA定量的计算分析流程

DOI：

10.3791/69400

⸱

2025年10月28日

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

* These authors contributed equally

本文内容

摘要

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

本方案提供了一条简化的计算流程，用于定量新生增强子转录本。通过整合染色质可及性、染色质特征和转录数据，该方法能够在复杂的基因内区域实现对增强子活性的准确检测和链特异性分析，同时适用于不具备深入生物信息学训练背景的研究人员。

摘要

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

核心顺式调控元件（cis-regulatory elements）被称为增强子，在实现对靶基因的精确转录调控中发挥中心作用，从而控制多种细胞功能和发育过程。这些增强子通常在两个方向上都被转录，产生被称为增强子RNA（enhancer RNAs, eRNAs）的长链非编码转录本。eRNAs的表达与活跃的染色质特征（如H3K27ac修饰和共激活因子的招募）密切相关，并在功能上促进靶基因的转录激活。然而，eRNAs的检测与定量仍然具有挑战性，尤其是在其转录区域与宿主基因重叠的情况下。为解决这一问题，我们提出了一种标准化且用户友好的计算流程，用于从新生RNA测序数据中分析增强子的转录活性。该方案指导用户完成数据预处理、读段比对和质量控制，随后进行增强子相关转录的链特异性定量，并针对信号归属复杂的基因内增强子（intragenic enhancers）提供了专门的处理步骤。可视化模块可清晰展示增强子在不同基因组背景下的活性，内置选项支持对基因间和基因内增强子的分析。本流程专为生物信息学经验有限的研究人员设计，为开展一致、可重复且可扩展的增强子转录研究提供了实用框架，有助于在多种生物系统中更广泛地应用增强子生物学研究。

引言

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

增强子是一类顺式调控DNA元件，通过介导染色质成环并招募转录机器来调控靶基因的转录1,2,3。其组织特异性活性使得在发育过程和细胞谱系决定中实现精确调控成为可能4,5,6,7,8。活跃的增强子具有特征性的染色质修饰，如H3K4me1（组蛋白H3第4位赖氨酸单甲基化）和H3K27ac（组蛋白H3第27位赖氨酸乙酰化），并且通常位于DNase I超敏感区域，这些区域标志着开放的染色质结构9,10,11,12。这些特征使转录因子和RNA聚合酶II能够结合DNA，从而在增强子位点启动新生转录13,14,15,16。

这一连续的生物学过程会产生来源于增强子的转录本，称为eRNA。eRNA具有双向性、非编码性，通常不带有poly(A)尾13,14,15,16。eRNA不仅可作为增强子活性的标志物，其本身也具有效应分子的功能16,17,18,19,20,21,22,23,24。它们可通过释放停滞的RNA聚合酶II上的负延伸因子（NELF），促进有效的转录延伸16,19，并有助于稳定增强子-启动子之间的环状结构17,18,20。此外，eRNA还可能通过m6A（N6-甲基腺苷）修饰参与转录凝聚体的形成21,22,23。

然而，由基因体内调控元件启动的基因内增强子转录功能仍存在争议。一些研究表明，来自基因内增强子的增强子RNA（eRNA）可增强宿主基因的表达25，其机制可能涉及促进NELF释放以及刺激依赖性的有效延伸过程26,27。相反，其他研究则认为，此类转录可能通过RNA聚合酶II的碰撞或转录干扰阻碍宿主基因，导致转录衰减或提前终止28,29。这些相互矛盾的观察结果，加上eRNA兼具标记物与调控因子的双重作用，凸显了对其进行精确量化和功能解析的必要性。然而，检测基因内eRNA十分困难，因为它们通常与正义链的宿主转录本存在重叠13,25,26,30。当增强子位于具有嵌套基因或双链均存在重叠转录的区域时，这一挑战进一步加剧，导致增强子特异性信号难以分辨。

为了克服这些挑战，我们开发了一种生物信息学分析流程，用于检测、定量和可视化与增强子相关的转录本，特别关注基因内区域。该流程整合了转座酶可及染色质测序（ATAC-seq）、染色质免疫沉淀测序（ChIP-seq）、全局延伸测序（GRO-seq）以及基因组注释信息，即使在复杂的基因组背景下，也能实现增强子水平的分辨率。

该流程包含四个主要步骤：（i）预处理、比对、峰识别和信号生成31；（ii）利用染色质特征识别增强子；（iii）链方向分配，特别是在基因体内部；（iv）新生增强子转录本的定量与可视化。该框架特别适用于具有高分辨率测序数据的系统，例如本研究中分析的小鼠胚胎干细胞，当具备合适的数据集时，也可扩展至其他生物体。通过在现有流程难以实现的增强子特异性定量方面提供支持，该工作流程为在不同基因组背景下对基因内eRNA转录进行基准评估和深入研究提供了一种实用工具。

访问受限。请登录或开始试用以查看此内容。

方案

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

注意：工作流程中使用的所有原始数据集均列于表1中。生物信息学工具的详细信息见材料表。本流程中使用的线程数可通过修改每个脚本顶部定义的THREADS变量进行调整。用户可根据自身的CPU资源增加线程数以加快分析速度。
每一步完成后都会生成一个日志文件。为快速检查是否失败，可使用如下命令：cat StepXX_log.txt; grep -qF "ERROR" StepXX_log.txt && echo "ERROR found. Fix before next step." || echo "OK: no ERROR markers"。若出现任何ERROR提示，应视为该步骤失败，需先解决问题后再继续下一步。

1. 从 GitHub 仓库下载完整分析流程

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

  1. 启动适用于所用操作系统的命令行界面（CLI）。
    1. Windows：使用适用于 Linux 的 Windows 子系统（WSL）设置 Linux 环境。在继续操作前，请按照官方说明安装并配置 WSL32。
    2. macOS：无需额外设置，因为 macOS 基于 Unix 系统。请参考官方指南打开终端33。
    3. Linux 用户，特别是使用 Ubuntu 的用户：按照参考说明中描述的方法打开终端34。
  2. 在终端中运行 wget https://github.com/myunggeunO/Enhancer-transcript-identification-from-read-to-visualization/archive/refs/heads/main.zip -O ~/pipeline.zip  ，以下载用于增强子识别和增强子 RNA 定量的分析流程。
  3. 在终端中输入 unzip ~/pipeline.zip -d ~/。这将把所有必要文件解压到主目录中。
  4. 运行 rm ~/pipeline.zip，然后输入 mv ~/Enhancer-transcript-identification-from-read-to-visualization-main ~/Enhancer-transcript-identification-from-read-to-visualization，以删除压缩包并重命名解压后的文件夹。
  5. 输入 cd ~/Enhancer-transcript-identification-from-read-to-visualization/，然后运行 chmod +x scripts/*，使“scripts/”目录中的所有脚本可执行。

2. 为分析流程设置 mamba/conda 环境

  1. 输入 bash scripts/Step1_conda_environment_formation.sh 以创建并运行 mamba 虚拟环境。执行过程中如出现提示，输入 Y 并按 Enter 键确认安装软件包。对于 macOS 系统，请执行步骤 2.1.1；对于已安装 Mamba 或 Conda 的系统，请执行步骤 2.1.2。
    1. macOS：打开脚本，将 miniconda 下载链接替换为 macOS 版本：
      https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh
      然后，继续执行步骤 2.1。
    2. 当提示符中出现 (enhancer-env) 时，输入 bash scripts/Step2_package_installation.sh 以安装下游分析所需软件包。安装过程中如出现提示，输入 Y 并按 Enter 键确认。
    3. 运行步骤 2.2 后，检查终端输出中是否存在错误信息。如有问题，请先解决后再重新运行步骤 2.2。
    4. （可选）运行 mamba list 以确认 材料表 中所有由 mamba 管理的软件包均已安装。HOMER 为手动安装，不会出现在 mamba 列表中。请通过检查是否存在“~/homer/”目录来验证 HOMER 是否已正确安装。

3. 从 SRA（序列读取存档）下载公开可用的 ChIP-seq、ATAC-seq 和 GRO-seq 数据集

  1. 运行 cp scripts/Step{3..12}_*.sh ./ 以复制用于原始测序读段处理的必要 shell 脚本。
  2. 输入 bash Step3_download_file_list.sh > Step3_log.txt 2>&1 并按下 Enter 键，以从 SRA 下载并处理原始测序数据。
    注意：该脚本可自动下载并准备公共测序数据用于分析。它将在“MATERIAL/”目录下创建一个标准化的文件夹结构，按实验类型和重复样本（生物学重复：rep1/rep2；技术重复：trep1/trep2）进行组织。脚本内置的 SRA 登录号列表将驱动 prefetch (v3.2.0) 进行数据获取，使用 fasterq-dump (v3.2.0) 将数据转换为 FASTQ 格式（双端测序： --split-files），并使用 pigz (v2.8) 进行压缩以减少存储空间。处理流程依次为 prefetch、fasterq-dump 和 pigz，输出文件将写入相应的“00.Rawdata/”目录中。

4. 对原始测序读段进行质量控制与修剪

  1. 输入命令 bash Step4_read_trimming_and_QC.sh > Step4_log.txt 2>&1，对原始 FASTQ 文件进行读段修剪和质量控制。
    注意：该脚本用于处理来自 GRO-seq、ATAC-seq 和 ChIP-seq（H3K27ac、H3K4me1）及其相应输入对照的原始 FASTQ 文件。首先对原始读段运行 FastQC（v0.12.1）35，然后使用特定实验参数，通过 Trim Galore（v0.6.10）36 去除接头序列。对于 GRO-seq，首先去除 NextSeq G 尾巴和极短读段（--nextseq 20, --length 20），然后使用 Cutadapt（v5.1）37 去除长的 poly-A 连续序列，同时保留长度超过 20 个核苷酸的读段（ -a A{15}, -m 20）。对于 ATAC-seq，处理双端文库，并针对 Tn5/Nextera 接头（--paired, --nextera）。对于 H3K27ac 及其输入样本，处理双端 ChIP-seq 文库（--paired）。对于 H3K4me1 及其输入样本，则执行标准的单端修剪（默认设置）。修剪完成后，再次对修剪后的读段运行 FastQC。所有输出结果写入各实验对应的“01.Clean/”目录中，中间生成的 GRO-seq 接头修剪文件将被删除。

5. 准备 Bowtie2 参考索引

  1. 选择以下两种选项之一，为 mm10（Mus musculus）参考基因组准备 Bowtie2 (v2.5.4)38 基因组索引。
    1. 运行 bash Step5_1_download_reference_index.sh > Step5_1_log.txt 2>&1，使用开发者提供的预构建 Bowtie2 索引。
    2. 运行 bash Step5_2_download_reference_make_index_with_
      bowtie2.sh > Step5_2_log.txt 2>&1      ，下载原始 mm10 基因组序列并手动构建索引。
      ​注意：两种方法均会在“reference_index/”目录中生成索引文件，且在标准比对中功能等效。

6. 将修剪后的序列比对至 mm10 参考基因组

  1. 输入命令 bash Step6_alignment_to_make_bam.sh > Step6_log.txt 2>&1，以将每个实验的测序读段比对至参考基因组，并生成 BAM 文件。
    ​注意：此步骤将每个实验的测序读段比对至先前已建立索引的 mm10 参考基因组。H3K27ac ChIP-seq 及其对应的 Input 样本采用双端比对模式，以实现准确性和灵敏度的平衡（-1, -2），并使用默认灵敏度参数。H3K4me1 ChIP-seq 及其对应的 Input 样本在默认设置下采用单端比对（-U）。ATAC-seq 采用高灵敏度比对模式，以适应长度可变且较长的 Tn5 酶切片段（--very-sensitive, -X 2000），并使用双端输入数据（-1, -2）。GRO-seq 采用高灵敏度比对，以更精确地定位经过截短处理的短读段（--very-sensitive），并使用单端输入数据（-U）。随后，使用 samtools（v1.22.1）39 的 view 命令将 SAM 文件转换为 BAM 文件，并进行过滤：ChIP-seq 和 Input 样本采用中等严格程度的比对质量值（MAPQ）过滤标准（-b, -q 10），而 ATAC-seq 和 GRO-seq 则采用更严格的过滤阈值（-b, -q 30）；最终生成的 BAM 文件将保存至各数据集的“02.Align/”目录中。

7. 合并 H3K27ac ChIP-seq 数据的技术重复

  1. 运行 bash Step7_merge_trep.sh > Step7_log.txt 2>&1 以合并 H3K27ac ChIP-seq 及相应 input 的技术重复 BAM 文件。
    ​注意：该脚本使用 sambamba（v1.0.1）40 的 sort 命令对技术重复的 BAM 文件进行排序，然后使用 sambamba merge 将 H3K27ac ChIP-seq 及相应 input 重复样本合并为整合的 BAM 文件。如果用户的数据集不包含任何技术重复，则跳过此步骤，直接使用单个 BAM 文件进行后续分析。

8. 去除重复序列和非必需染色体

  1. 对 ChIP-seq 和 GRO-seq 数据去除重复序列并排序比对后的读段。
    1. 输入命令 bash Step8_1_duplicate_removal_sorting-ChIP_GRO.sh > Step8_1_log.txt 2>&1，以去除组蛋白 ChIP-seq 数据集中的 PCR 重复序列，并对 ChIP-seq 和 GRO-seq 输出的 BAM 文件进行基于坐标排序。
      ​注意：该脚本使用 sambamba markdup（-r）去除 PCR 重复序列，并使用 sambamba sort 进行坐标排序。对于 H3K27ac 数据，处理对象为技术重复样本合并后的 BAM 文件，ChIP 与 input 样本分别处理；对于 H3K4me1 数据，每个重复样本及其配对的 input 样本均单独处理。对于 GRO-seq 数据，则跳过去重步骤，仅进行坐标排序。输出文件保存在各数据集的“02.Align/”目录中。
  2. 对 ATAC-seq 数据去除重复序列并过滤线粒体染色体比对上的读段。
    1. 输入命令 bash Step8_2_duplicate_chrM_removal_sorting_ATAC.sh > Step8_2_log.txt 2>&1，以去除 PCR 重复序列、过滤线粒体读段（chrM）并对 ATAC-seq BAM 文件进行排序。
      ​注意：该脚本参考了 ENCODE 联盟推荐的 ATAC-seq 数据处理流程。首先使用 sambamba sort（-n）进行名称排序，并使用 samtools fixmate（-m）修复双端读段的配对信息，以确保在标记重复序列前正确分配配对信息。随后使用 sambamba markdup（-r）去除 PCR 重复序列。通过 samtools idxstats 生成保留列表（排除 chrM 和 *），再使用 samtools view（-b）仅保留列表中的参考序列，从而去除线粒体读段。最后使用 sambamba sort 进行坐标排序。处理后的干净 BAM 文件写入各重复样本的“02.Align/”目录中。

9. 对每个数据集进行峰检测

  1. 输入命令 bash Step9_peak_calling.sh > Step9_log.txt 2>&1 运行第9步脚本，以对 ChIP-seq 和 ATAC-seq 数据进行峰识别。
    注意：该脚本使用 MACS3 (v3.0.3)41 对 ATAC-seq 和 ChIP-seq（H3K27ac、H3K4me1）数据执行峰识别。ATAC-seq 在无模型模式下将所有重复样本的 BAM 文件作为信号输入，并设置偏移/延伸参数（-f BAMPE, --nomodel, --shift -100, --extsize 200, -q 0.01）来生成峰。H3K27ac 使用合并的技术重复样本及其匹配的输入对照，在宽峰模式下识别宽峰（-f BAMPE, --broad）。H3K4me1 则在单端宽峰模式下（-f BAM, --broad）对每个重复样本及其匹配的输入对照单独处理，并使用 bedtools (v2.31.1)42 的 intersect 工具获取重叠峰。输出结果保存至各数据集的“peak_calling/”目录中，高置信度的 H3K4me1 重叠峰位于“peak_calling/overlapped_peak/”子目录下。

10. 合并 ChIP-seq 和 ATAC-seq BAM 文件的生物学重复样本以进行下游信号分析

  1. 运行 bash Step10_merge_rep_forMakingSignal.sh > Step10_log.txt 2>&1 ，以合并 ATAC-seq 和 H3K4me1 ChIP-seq 的生物学重复样本的 BAM 文件。
    ​注意：该脚本使用 sambamba merge 工具合并 ATAC-seq、H3K4me1 ChIP-seq 及相应的 H3K4me1 input 的重复样本 BAM 文件。合并后的 BAM 文件可用于后续分析（例如 bigWig 信号生成和标准化）。输出的 BAM 文件将保存在各样本路径下的“merge/02.Align/”目录中。

11. 从比对后的序列读段生成标签目录和信号 bigWig 文件

  1. 运行 bash Step11_make_tag_to_signal.sh > Step11_log.txt 2>&1 以创建标签目录，并为每个数据集的比对读段生成 bigWig 信号文件。
    ​注意：该脚本使用 HOMER（v5.1）43 和 ucsc-bedgraphtobigwig（v482）44 工具包中的 makeTagDirectory 命令构建 HOMER 标签目录，然后通过 makeUCSCfile 命令生成 bigWig 信号轨道。所有信号轨道均基于 UCSC 基因组浏览器提供的染色体大小文件生成（https://hgdownload.soe.ucsc.edu/goldenPath/mm10/）。 GRO-seq 生成链特异性信号轨道 （-style rnaseq, -strand + / -, -bigWig）。ATAC-seq 基于合并后的 BAM 文件生成未归一化的轨道（-bigWig）。ChIP-seq（H3K27ac, H3K4me1）生成以输入样本为参照、伪计数为 1 的归一化轨道（-bigWig, -i, -pseudo 1）。输出结果分别存放在“03.TagDir/”和“04.bigwig/”目录下。

12. 准备用于增强子鉴定的文件

  1. 在终端中输入 bash Step12_E_identification_material.sh > Step12_log.txt 2>&1，以准备增强子鉴定所需的文件和参考数据。
    注意：此步骤将增强子鉴定所需的所有必要文件收集到“01.E_identification/material/”目录下，并分别组织到“ATAC/”、“Histone/”和“Annotation/”文件夹中。该步骤会将峰值文件（ATAC-seq、H3K27ac、H3K4me1）复制到相应目录。GENCODE M23 注释文件（mm10）将自动下载，同时从预定义路径复制参考文件，包括 ENCODE 黑名单45（mm10-blacklist.v2.bed）和染色体大小文件（mm10.chrom.sizes）。

13. 从注释中鉴定启动子候选区和基因体

  1. 输入 cd 01.E_identification/ 进入工作目录，然后运行 cp ../scripts/Step{13..20}_*.sh ./ 以复制增强子鉴定脚本。
  2. 运行 bash Step13_promoter_candidates_genebody_identification.sh > Step13_log.txt 2>&1，使用 GENCODE 注释生成启动子区域、基因体和蛋白质编码基因（PCG）的 BED 文件。
    注意：该脚本处理下载的 GENCODE M23 基因转移格式（GTF）文件，以创建启动子候选区域、所有基因体以及蛋白质编码基因体的 BED 文件，并将输出结果保存至“material/Annotation/”目录。启动子定义为每个转录本转录起始位点（TSS）上下游各 2 kb 的区域，通过 bedtools slop 工具（-b 2000, -g mm10.chrom.sizes）生成。基因体和蛋白质编码基因体来源于 GTF 中标注为“gene”的条目，其中蛋白质编码基因进一步通过“gene_type = protein_coding”进行筛选。

14. 处理峰以进行增强子鉴定

  1. 运行 bash Step14_ATAC_ChIP-seq_processing.sh > Step14_log.txt 2>&1 以预处理 ATAC-seq 和组蛋白 ChIP-seq 峰文件，用于增强子鉴定。
    注意：该脚本用于增强子识别前的峰集预处理。对于 ATAC-seq，使用 bedtools subtract 命令（-A）去除与黑名单区域重叠的区域，然后再次使用 bedtools subtract 命令（-A）排除与启动子候选区重叠的峰；对于组蛋白修饰标记（H3K27ac、H3K4me1），使用 bedtools slop 命令（-b 1000 -g mm10.chrom.sizes）将每个峰在两侧对称扩展 1 kb，然后使用 bedtools subtract 命令去除与启动子区域的重叠部分。

15. 鉴定并分类增强子

  1. 运行 bash Step15_inter_intragenic_E_sets_identification.sh > Step15_log.txt 2>&1，利用染色质峰数据定义并分类增强子。
    ​注意：此步骤使用 bedtools 定义和注释增强子。通过 intersect 工具（-wa, -u）获取 ATAC-seq 峰与扩展后的 H3K4me1 峰之间的重叠区域，并进一步使用 intersect（-wa, -u）筛选出同时与侧翼 H3K27ac 峰重叠的区域，将其归类为活性增强子；非活性增强子则通过 subtract 工具从完整增强子集合中去除活性区域获得。使用 intersect（-u）收集与各增强子集合重叠的 ATAC-seq 峰顶点，然后根据基因体使用 intersect（-v 或 -u）将峰顶点划分为基因间区和基因内区。最终，基于峰关联的峰顶点，使用 intersect（-u）将增强子区间分配至基因间区或基因内区类别。所有结果保存在“01.E_identification/”目录下的“01.allE/”、“02.interE/”和“03.intraE/”文件夹中。

16. 为基因内增强子分配临时链信息

  1. 输入命令 bash Step16_assign_temp_strand_from_gene_overlap.sh > Step16_log.txt 2>&1 ，为基因内增强子 BED 文件分配临时链信息。
    注意：此步骤使用 bedtools intersect（-wa, -wb）将增强子区间与基因体重叠，从而为基因内增强子分配基因链方向标签。通过 awk 保留第 1、2、3、4、5 和 16 列，随后对记录进行排序并去除重复。输出结果保存为 "03.intraE/strand_designation/01.overlapped_gene_strand/ES_E_intragenic_strand_with_dup.bed"。

17. 对重叠双向基因的增强子优先进行链分配

  1. 输入命令 bash Step17_initial_strand_assignment_for_both_strand_enhancers_PCG_based.sh > Step17_log.txt 2>&1，以确定跨越两条链基因的基因内增强子的链方向。
    注意：该脚本通过优先考虑与蛋白编码基因（PCG）重叠的情况，解决跨越两条链基因的基因内增强子的链方向不确定性问题。首先分离出存在于两条链上的增强子（使用 awk 按染色体/起始位点/终止位点/ID/链进行分组），然后利用 bedtools intersect（-s, -wa, -u）筛选出与PCG同链重叠的情况，再用 bedtools intersect（-v）保留非PCG重叠的情况。最后将筛选出的集合与保留的集合合并并排序。输出文件："03.intraE/strand_designation/02.enhancer_with_PCG_priority/ES_E_intragenic_PCG_priority.bed"。

18. 计算与PCG优先的基因内增强子重叠基因的链特异性每百万映射读段每千碱基读段数（RPKM）值

  1. 输入命令 bash Step18_RPKM_cal_from_partially_strand_assigned_enhancers.sh > Step18_log.txt 2>&1，以计算与增强子同链重叠基因的链特异性 RPKM。
    ​注意：此步骤使用第 17 步生成的文件“ES_E_intragenic_PCG_priority.bed”，该文件包含以下三类增强子：(i) 在一条链上与蛋白编码基因（PCG）重叠并已分配链方向；(ii) 在两条链上均与 PCG 重叠因而链方向仍不明确；(iii) 未与任何 PCG 重叠而保留双链信息。通过 bedtools intersect 工具（-s, -wa, -u 参数）筛选出同链重叠的基因，使用 awk 转换为 GTF 格式，并利用 GRO-seq 数据通过 featureCounts（v2.1.1）46 的链特异性计数模式（-s 1, -t gene, -g gene_id, -O, --fraction）进行定量。总比对读段数由 sambamba flagstat 获得，RPKM 值则根据基因长度、计数结果和总读段数计算得出。输出结果保存至“03.intraE/strand_designation/03.RPKM_calculation_of_overlapped_gene/”目录中。

19. 基于重叠基因表达（RPKM）的最终链分配

  1. 输入命令 bash Step19_second_strand_assignment_by_RPKM.sh > Step19_log.txt 2>&1，以完成对基因内增强子的链方向指派。
    ​注意：此步骤利用第18步获得的基因表达支持数据，为基因内增强子指派链方向。使用 bedtools intersect（-s, -wa, -wb）找出增强子（ES_E_intragenic_PCG_priority.bed）与基因之间的同链重叠区域。通过 awk/sort/join 将基因的 RPKM 值关联至对应的基因区间，生成包含基因 RPKM 信息的 BED 文件。对于每个增强子，选择与其重叠且 RPKM 值最高的基因，并使该增强子继承该基因的链方向。结果保存至“03.intraE/strand_designation/04.enhancer_strand_designation_by_RPKM_of_gene/ES_E_intragenic_PCG_priority_strand_by_gene_RPKM.bed”。

20. 为基因内增强子及增强子峰分配链信息

  1. 输入命令 bash Step20_strand_assignment_for_intragenicE_and_summits.sh > Step20_log.txt 2>&1，为所有基因内增强子（intragenic enhancers）及各个峰顶（summit）分配确定的链信息。
    注意：此步骤使用 bedtools 最终确定基因内增强子及其对应峰顶的链方向。通过 subtract 命令（-S）去除反义链区间，使用 intersect 命令（-wa, -u）将指定链的增强子与活性和非活性增强子集合进行交集分析，并通过 intersect 命令（-wa, -wb）将峰顶文件与指定链的增强子重叠，再利用 awk 提取链信息以重新注释峰顶文件。结果将保存在“03.intraE/final_strand_IntragenicE/”目录及其“summit/”子目录中。

21. 准备用于增强子验证、eRNA定量和可视化的输入文件

  1. 输入 cd ../ 或 cd ~/Enhancer-transcript-identification-from-read-to-visualization 进入流程根目录，然后输入 cp scripts/Step21_preparing_quantification_and_visualization.sh ./，以复制下游分析准备脚本。
  2. 运行 bash Step21_preparing_quantification_and_visualization.sh > Step21_log.txt 2>&1，以准备增强子聚合、GRO-seq 信号处理和 eRNA 定量所需的所有文件。
    注意：此步骤将在“02.E_visualization_quantification/”目录下创建用于增强子验证、eRNA 定量和信号可视化的输入文件及目录结构。系统将创建用于存放 bigWig 文件、增强子 BED 文件、BAM 文件、信号矩阵、计数结果和图表的子文件夹。关键输入文件（如 summit BED 文件、bigWig 文件、增强子列表和 GRO-seq BAM 文件）将被复制到相应位置。

22. 生成增强子验证的聚合图

  1. 输入 cd 02.E_visualization_quantification/ 进入工作目录，然后输入cp ../scripts/Step{22..24}_*.* ./ 复制下游分析所需的脚本。
  2. 运行bash Step22_generation_of_aggregation_plot.sh > Step22_log.txt 2>&1 生成每类增强子峰周围染色质信号的聚合图。
    注意：此步骤使用 deepTools (v3.5.6)47 中的 computeMatrix 和 plotProfile 可视化以增强子峰为中心的平均染色质信号富集情况。对于每个定义的增强子集合，computeMatrix reference-point（--referencePoint center, -a 5000, -b 5000, --missingDataAsZero）利用 ATAC-seq、H3K27ac 和 H3K4me1 的 bigWig 文件，计算增强子峰周围 10 kb 窗口内的信号密度。输出的矩阵被传递给 plotProfile，以生成不同增强子组之间的信号聚合曲线。聚合图保存在“01.Profiling/04_1.aggregation/”目录中。

23. 定量与可视化增强子RNA表达

  1. 运行 bash Step23_quantifing_eRNA_RPKM.sh > Step23_log.txt 2>&1，使用 featureCounts 从 GRO-seq 数据中量化 eRNA 表达水平。
    注意：此步骤使用 GRO-seq 以链特异性方式量化来自已定义增强子区域的 eRNA 转录。基因间增强子采用无链特异性模式（-s 0, -t enhancer, -g gene_id, -O, --fraction）通过 featureCounts 进行计数，而基因内增强子则采用反义链模式（-s 2）进行量化，以排除与基因转录重叠的信号。在计数前，BED 区域通过 awk 转换为 GTF 格式。总比对读段数来自 sambamba flagstat，计数结果通过增强子长度、读段计数和总比对读段数归一化为 RPKM。输出结果按“inter/”和“intra/”分类，存储于“02.eRNA_quantification/03.count_normalized_with_RPKM/”目录下。
  2. 运行 Rscript Step24_visualization_of_enhancer_transcript.R > Step24_log.txt 2>&1，使用 R 对不同增强子组间的 eRNA 表达水平进行可视化和比较。
    注意：R 脚本使用 ggplot2 (v3.5.2)48 和 cowplot (v1.2.0)49 软件包，基于 RPKM 值生成小提琴图和箱线图，用于比较活跃与非活跃增强子的表达水平。为便于可视化和统计检验，RPKM 值被转换为 log2 (RPKM + 1)。统计显著性通过 Wilcoxon 秩和检验进行评估。汇总图和 p 值表格均保存至“03.eRNA_visualization/”目录，供后续分析使用。
    注意： 若本流程中任何步骤失败且重新运行后问题仍存在，请在 https://github.com/myunggeunO/Enhancer-transcript-identification-from-read-to visualization/issues 提交问题报告。明确说明失败步骤并附上日志文件，有助于获得准确的故障排查支持。

访问受限。请登录或开始试用以查看此内容。

结果

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

增强子转录本定量分析流程的示意图
公开可用的 ChIP-seq（H3K27ac、H3K4me1）、ATAC-seq 和 GRO-seq 数据集（表 1）通过一个标准化流程进行处理，主要用于验证。使用 Trim Galore 和 Cutadapt 进行接头序列修剪和质量过滤，随后利用 Bowtie2 将测序数据比对至 mm10 参考基因组（详见实验方案步骤 6）。对于 ChIP-seq 和 ATAC-seq 数据，采用 MACS3 鉴定峰值区域，并使用 HOMER 和 ucsc-bedgraphtobigwig 生成信号强度轨迹文件（bigWig 文件），用于下游可视化和定量分析（图 1A）。该一致性的预处理流程可确保数据兼容性，最大限度减少不同实验方法带来的偏差，并为下游增强子定量提供可靠的分析框架。

增强子的鉴定基于...

访问受限。请登录或开始试用以查看此内容。

讨论

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

在发现增强子来源的转录本后13,14,15，准确量化增强子RNA（eRNA）仍是一个重大挑战，尤其是在基因内环境中，eRNA常与宿主基因的转录本发生重叠。这种重叠使得链的分配和信号归属变得复杂，难以区分真正的增强子转录与背景基因表达13,25,26,30。尽管eRNA increasingly被认可为增强子活性的功能性指标，但该领域一直缺乏标准化且易于获取的增强子转录本量化框架。为满足这一需求，我们开发了一种模块化且具备链特异性的分析流程，该流程整合染色质可及性与组蛋白修饰数据以识别增强子，并利用新生RNA测序技术高精度地量化增强子来源的转录。该流程在设计时注重易用性，使不具备深厚计算背景的研究人员...

访问受限。请登录或开始试用以查看此内容。

披露

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

作者无任何利益冲突需披露。

致谢

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

本研究由韩国忠南国立大学研究基金资助［2022-0582-01（S.-K.K.）和2023-0545-01（S.-K.K.）］。图1使用BioRender（https://biorender.com/）制作。

访问受限。请登录或开始试用以查看此内容。

材料

本文使用的材料清单
姓名公司目录编号评论
bedtoolsQuinlan 实验室，犹他大学 v2.31.1用于编辑 BED 文件的工具集
bowtie2Langmead 实验室，约翰斯·霍普金斯大学v2.5.4用于将测序读段比对至参考基因组的多线程比对工具
cowplotWilke 实验室，德克萨斯大学v1.2.0用于组合和对齐基于 ggplot2 的图形的工具
cutadapt斯德哥尔摩大学，生命科学中心v5.1接头序列及 poly-A/G 尾剪切工具
deeptools马克斯·普朗克研究所，生物信息学平台v3.5.6用于量化特定基因组区域中测序读段数量的工具
fastqcBabraham 生物信息学团队，Babraham 研究所v0.12.1测序读段的质量控制分析工具
featureCounts (subread)Shi 实验室，蒙纳士大学v2.1.1针对指定基因组区域的原始读段计数工具
homerBenner 实验室，加州大学圣地亚哥分校（UCSD）v5.1用于 ChIP-seq、ATAC-seq 和新生 RNA 分析的工具包；包含标签目录构建和信号谱分析功能
macs3Chan Zuckerberg 倡议组织v3.0.3用于 ChIP-seq 和 ATAC-seq 数据集的峰识别分析
pigz.v2.8用于生成 gzip 压缩文件的多线程压缩工具
sambamba圣彼得堡州立大学v1.0.1多线程 SAM/BAM 文件处理工具包
samtools惠康信托桑格研究所v1.22.1用于处理和操作 SAM/BAM 文件的工具集
sra-tools美国国家生物技术信息中心（NCBI）v3.2.0用于从 NCBI SRA 数据库下载 SRR 文件
tidyversePosit PBCv2.0.0用于数据处理与可视化的 R 软件包集合
trim-galoreAltos Labs，剑桥科学研究所v0.6.10基于多线程的接头序列及低质量碱基剪切工具
Ubuntu 20.04用于流程的开发与测试
ucsc-bedgraphtobigwig Kent 实验室，加州大学圣克鲁兹分校v482用于生成 bigWig 信号轨道的工具

参考文献

$$\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).

访问受限。请登录或开始试用以查看此内容。

重印与许可

申请许可以重复使用本 JoVE 文章的文本或图表

申请许可

标签

GRO seq ATAC seq H3K27

相关文章