
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae246
bbae246
Problem Solving Protocol
AcademicSubjects/SCI01060
Predicting RNA polymerase II transcriptional elongation pausing and associated histone code
https://orcid.org/0000-0001-9054-0809
Ren Lixin School of Mathematics and Physics, University of Science and Technology Beijing, 30 Xueyuan Road, Haidian District, Beijing 100083, China

Ma Wanbiao School of Mathematics and Physics, University of Science and Technology Beijing, 30 Xueyuan Road, Haidian District, Beijing 100083, China

Wang Yong CEMS, NCMIS, MDIS, Academy of Mathematics and Systems Science, National Center for Mathematics and Interdisciplinary Sciences, Chinese Academy of Sciences, 55 Zhongguancun East Road, Haidian District, Beijing 100190, China
School of Mathematical Sciences, University of Chinese Academy of Sciences, 19A Yuquan Road, Shijingshan District, Beijing 100049, China
Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, 32 Jiaochang Donglu, Wuhua District, Kunming 650223, China
Key Laboratory of Systems Biology, Hangzhou Institute for Advanced Study, University of Chinese Academy of Sciences, Chinese Academy of Sciences, 1 Xiangshan Zhi Nong, West Lake District, Hangzhou 330106, China

Corresponding author. Wanbiao Ma, School of Mathematics and Physics, University of Science and Technology Beijing, 30 Xueyuan Road, Haidian District, Beijing 100083, China. E-mail: wanbiao_ma@ustb.edu.cn; Yong Wang, CEMS, NCMIS, MDIS, Academy of Mathematics and Systems Science, National Center for Mathematics and Interdisciplinary Sciences, Chinese Academy of Sciences, 55 Zhongguancun East Road, Haidian District, Beijing 100190, China. E-mail: ywang@amss.ac.cn
7 2024
23 5 2024
23 5 2024
25 4 bbae24622 1 2024
12 4 2024
07 5 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

RNA Polymerase II (Pol II) transcriptional elongation pausing is an integral part of the dynamic regulation of gene transcription in the genome of metazoans. It plays a pivotal role in many vital biological processes and disease progression. However, experimentally measuring genome-wide Pol II pausing is technically challenging and the precise governing mechanism underlying this process is not fully understood. Here, we develop RP3 (RNA Polymerase II Pausing Prediction), a network regularized logistic regression machine learning method, to predict Pol II pausing events by integrating genome sequence, histone modification, gene expression, chromatin accessibility, and protein–protein interaction data. RP3 can accurately predict Pol II pausing in diverse cellular contexts and unveil the transcription factors that are associated with the Pol II pausing machinery. Furthermore, we utilize a forward feature selection framework to systematically identify the combination of histone modification signals associated with Pol II pausing. RP3 is freely available at https://github.com/AMSSwanglab/RP3.

Pol II pausing
machine learning
histone modification
multi-omics data integration
National Key Research and Development Program of China 10.13039/501100012166 2022YFA1004800 2020YFA0712402 CAS Project for Young Scientists in Basic Research YSBR-077 National Natural Science Foundation of China 10.13039/501100001809 12025107 12371481
==== Body
pmcIntroduction

RNA polymerase II (Pol II) serves as the principal enzyme responsible for the transcription of all nuclear protein-coding genes and a substantial proportion of non-coding genes in eukaryotic cells [1–3]. Transcription can be categorized into three stages: initiation, elongation, and termination [4]. During elongation, Pol II translocates along the gene body and incorporates nucleotides into the nascent RNA [5]. However, this process is not continuous and uniform [6], as Pol II can pause at specific sites and regulate the rate and timing of transcription [7, 8]. Pol II pausing is widespread in eukaryotic genomes [3] that up to 40% of protein-coding genes undergo this regulatory mechanism during transcription [9, 10].

Pol II pausing has been implicated in various biological processes and pathologies, such as, differentiation and development [11, 12], inflammation [13], cardiac hypertrophy [14] and cancers [15]. Recently, the native elongating transcript sequencing (NET-seq) approach based on high-throughput sequencing of 3′-ends of the nascent RNA providing a global and strand-specific quantitative measure of Pol II abundance in cells [7, 16]. The powerful NET-seq technique allows for the identification of the full spectrum of Pol II pausing events at single-nucleotide resolution [7]. However experimental profiling of Pol II pausing remains technically challenging and costly, resulting in only a limited fraction of cell types have been analyzed [7, 17]. Furthermore, despite the advancement of high-throughput DNA sequencing techniques enabling the characterization of Pol II genomic location distribution and abundance, the underlying molecular mechanisms and regulatory factors that govern Pol II pausing are not fully understood. Therefore, developing computational prediction methods is a desirable strategy for interrogating the information of Pol II pausing and exploring regulatory mechanisms for cell types of interest.

Recently, PEPMAN was developed as a deep learning framework for modeling Pol II pausing by utilizing cellular context non-specific raw DNA sequences as input features [17]. However, since Pol II pausing events can be cell type-specific, reliance on context non-specific DNA sequences means that this method may not accurately capture or reflect the cell-specific epigenetic factors and nuances. Gajos et al. developed a model of Pol II pausing based on random forests by incorporating purely sequence-dependent features (such as Z-DNA etc.) and DNA methylation levels [3]. Although the model includes DNA methylation, it overlooks other important epigenetic factors like histone modifications and chromatin accessibility. Recent studies have shown that histone modifications are involved in the regulation of Pol II pausing. For example, H3K4me2/3 can modulate the stability of Pol II pausing [18], H4K20me3 enforces Pol II promoter-proximal pausing by antagonizing hMOF-mediated H4K16Ac [19], and H3K9ac promotes the release of Pol II pausing by directly recruiting the super elongation complex (SEC) to chromatin [20]. Moreover, trans-acting transcription factors (TFs) also contribute to Pol II pausing. For example, CTCF, YY1, and Reb1 induce Pol II pausing by creating obstacles on the DNA template [21–23], TF-binding sites for ESR1 and PAX5 are significantly enriched in gene body pauses [3]. Collectively, increasing evidence suggests pivotal roles for histone modifications and TFs in Pol II pausing regulation. However, a systematic way to leverage the widely available transcriptomic and epigenomic data to identify key histone modifications and TFs orchestrating Pol II pausing is still lacking.

Here, our aim is to mine widely available multi-omics data to interrogate the information of transcript elongation Pol II pausing, and to identify the potential regulators involved in this process. Specifically, we develop RP3, a machine-learning framework based on network regularized logistic regression to predict transcriptional elongation Pol II pausing by integrating widely available multilayer omics data. This model is distinct in its ability to integrate diverse omics data layers, including genome sequence, histone modifications, gene expression, chromatin accessibility, and protein–protein interaction (PPI). This comprehensive approach surpasses existing models that typically rely on a narrower range of data types, thereby providing a more holistic perspective on Pol II pausing. RP3 exhibits exceptional performance in accurately predicting Pol II pausing events not only within individual cell lines but also across various cellular contexts. Our unique implementation of the network regularization strategy incorporates topological information from TF PPI networks and co-expression networks across cell lines for feature selection, effectively identifying biologically significant TFs associated with Pol II pausing [24]. Furthermore, we employed a forward feature selection framework to systematically investigate the combination of histone modification signals, i.e. histone code, that are associated with Pol II pausing. Overall, our work provides a robust way to predict Pol II pausing events on human genome by using widely accessible omics data in different cell lines and offers biological insights into the regulatory roles of histone codes and TFs in governing the mechanism of Pol II pausing.

Materials and methods

Data collection

We collected NET-seq and RNA-seq data for HeLa-S3 and HEK293T cell lines from [22]. The corresponding DNase-seq was obtained from the ENCODE project. ChIP-seq data for 10 histone modifications (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me2, H3K4me3, H3K79me2, H3K9ac, H3K9me3, H4K20me1) in HeLa-S3 were downloaded from ENCODE. For HEK293T, ChIP-seq data for five histone modifications (H3K27ac, H3K36me3, H3K4me1, H3K4me2, H3K4me3) were sourced from GEO (accession numbers GSE51633 and GSE99407). RNA-seq data for 832 human samples were collected from ENCODE and ROADMAP. PPI data were gathered from the BIOGRID database.

Construction of gold-standard data

To extract the gold-standard positive (GSP) samples, i.e. Pol II pausing sites, from the NET-seq data, we followed the same guidelines as described in [16]. Specifically, we used the following criteria: (i) the read density was at least three standard deviations above the mean over a surrounding 201 bp window, (ii) in cases where the genomic distance between two positive samples was < 201 bp, we kept only the sample with the higher read density, (iii) the read count was larger than four, regardless of sequencing coverage. Furthermore, to focus exclusively on Pol II pausing events in expressed genes, we excluded those genes without any RNA-seq reads mapped to their gene bodies. In total, we obtained 79 572 and 58 116 positive samples for HeLa-S3 and HEK293T, respectively. Then, we generated gold-standard negative (GSN) samples by randomly selecting regions from the whole genome that did not contain any Pol II pausing events. The sample size of both GSP and GSN was kept the same.

Feature extraction and combination

Genomic features including histone modification signal strength (H), chromatin openness (O), TF expression (TF), TF binding strength (B), and TF expression specificity (TFS) were extracted from multi-omics data include ChIP-seq, DNase-seq, RNA-seq, and sequence data (see details in Supplementary Materials). These features are not independent, as TF binding strength, TF expression, and TF expression specificity each provide unique information about the activity of TFs. To better capture the overall activity of TFs, we combine these features and define a genomic feature called TF activity (TFA) as the product of TF binding strength, TF expression, and TF expression specificity:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ TFA= TF\cdot B\cdot TFS $$\end{document}

TF activity indicates that TF has a significant motif match on genomic regions, high expression, and specific expression across cell types, which facilitates Pol II pausing.

After the extraction and combination of the functional genomic features, a set of genomic features was obtained as model input, including histone modification signal strength (H), chromatin openness (O), and TF activity (TFA).

Logistic regression algorithm

We utilized logistic regression (LR), a robust discriminative method commonly utilized for binary classification problems, to predict the occurrence of Pol II pausing. Let us consider gold-standard dataset that comprises both GSP and GSN, and has a size of n, with the predictor dimensionality being p. The response variable value for the ith sample, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_i$\end{document}, can take two values [0,1], where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_i=1$\end{document} indicates the genomic regions of the ith sample having pol II pausing and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_i=0$\end{document} indicates the absence of Pol II pausing. The model is fitted using genomic features extracted from the input data, including histone modification signal strength (H), chromatin openness (O), and TF activity (TFA). The mathematical expression of the model can be formulated as follows:

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{align*} \mathit{\log}\frac{P\left({Y}_i=1\left| TF,{O}_i\right.,H\right)}{1-P\left({Y}_i=1\left| TF,{O}_i,H\ \right.\right)}&={\alpha}_0+{\alpha}_1{O}_i+{\sum}_{k\in S}{\beta}_k{B}_{i,k}\cdot T{F}_k\cdot TF{S}_k\\\nonumber&+{\sum}_{j\in P}{\gamma}_j{H}_{i,j} \end{align*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $S= PPI\left(\mathrm{Pol}\ \mathrm{II}\right)\cap MB,$\end{document}  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P$\end{document} is the set of histone modifications.

(2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{align*} & {u}_i=P\left({Y}_i=1\left| TF,{O}_i\right.,H\right)\nonumber\\&\qquad =\frac{\exp \left({\alpha}_0+{\alpha}_1{O}_i+{\sum}_{k\in S}{\beta}_k{B}_{i,k}\cdot T{F}_k\cdot TF{S}_k+\sum_{j\in P}{\gamma}_j{H}_{i,j}\right)}{1+\exp \left({\alpha}_0+{\alpha}_1{O}_i+{\sum}_{k\in S}{\beta}_k{B}_{i,k}\cdot T{F}_k\cdot TF{S}_k+\sum_{j\in P}{\gamma}_j{H}_{i,j}\right)} \end{align*}\end{document}

The log-likelihood function of Equation (1) is formally defined as:

(3) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} \ell \left(\alpha, \beta, \gamma \right)=\frac{1}{n}{\sum}_{i=1}^n{y}_i\mathit{\log}{\mu}_i+\left(1-{y}_i\right)\mathit{\log}\left(1-{\mu}_i\right) \end{equation*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\alpha =\left({a}_0,{\alpha}_1\right)$\end{document} a vector of coefficients that are independent of histone modifications and TFs, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta =\left({\beta}_1,{\beta}_2,{\beta}_3,\dots \dots \right)$\end{document} represents the coefficients that are dependent on TFs, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma =\left({\gamma}_1,{\gamma}_2,{\gamma}_3,\dots \dots \right)$\end{document} represents the coefficients that are dependent on histone modifications, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_i$\end{document} represents the training label of sample i.

The parameters are estimated by minimizing the negative log-likelihood function:

(4) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} \underset{\alpha, \beta, \gamma }{\min }-\left[\frac{1}{n}{\sum}_{i=1}^n{y}_i\mathit{\log}{u}_i+\left(1-{y}_i\right)\log \left(1-{u}_i\right)\right] \end{equation*}\end{document}

Network regularization strategy for logistic regression

The network regularization strategy employed in our study was taken from our previous work [24], which focused on reconstructing the network topology among predictors for regularized logistic regression. Here, we applied this strategy with the assumption that TFs involved in assisting Pol II to pause act in a cooperative rather than independent manner. The objective of our network regularization approach is to unveil the biological relationships among these TFs by exploiting the network structural information within the TFs. Specifically, it mainly involves three steps:

(i) Constructing TF networks

TF set related with Pol II pausing was constructed as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $S= PPI\left(\mathrm{Pol}\ \mathrm{II}\right)\cap MB$\end{document} by considering PPI with Pol II and motif occurrence. Let S contain q TFs. To explore the network structure relationship among TFs, we initially analyzed the physical PPIs among TFs and extracted the TF PPI network from available PPI data. The TF PPI network was represented as an unweighted and undirected graph, denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_1\left(S,{E}_1\right)$\end{document}, described by the binary adjacency matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A={\left({a}_{ij}\right)}_{q\times q}$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_{ij}=1$\end{document} if there exists the PPI between the ith TF and the jth TF; otherwise, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_{ij}=0$\end{document}. Secondly, we investigated the co-expression relationships among TFs by analyzing a large amount of gene expression data in diverse contexts, and generated TF co-expression network. The TF co-expression network was represented as a complete graph, denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_2\left(S,{E}_2\right)$\end{document}, described by the adjacency matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $B={\left({b}_{ij}\right)}_{q\times q}$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${b}_{ij}$\end{document} is the weight of the edge between the ith and jth TF, and determined based on their correlation in expression levels across diverse contexts.

(ii) Integrating TF networks

Next, we combined the TF PPI and co-expression networks by computing the Hadamard product of their respective adjacency matrices:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ C=A\circ B={\left({C}_{ij}\right)}_{q\times q} $$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${c}_{ij}={a}_{ij}\cdot{b}_{ij}$\end{document}.

The resulting matrix C represents the integrated TF network, where each element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${c}_{ij}$\end{document} in C is calculated by taking the Hadamard product of the corresponding elements in the PPI and co-expression adjacency matrices, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${c}_{ij}={a}_{ij}\cdot{b}_{ij}$\end{document}. The resulting integrated network represents the combined information of both PPI and co-expression relationships among the TFs.

(iii) Constructing the candidate TF network (CTN)

Define the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_Q(x)$\end{document} as follow:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ {f}_Q(x)=\left\{\begin{array}{c} \mid x \mid \mid x\mid \ge Q\\{}0 \mid x\mid <Q\end{array}\right. $$\end{document}

Then we update the integrated TF network represented by matrix C by defining a new matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $R=$\end{document}  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\left({r}_{ij}\right)}_{q\times q}$\end{document}, where each element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${r}_{ij}={f}_Q\left({c}_{ij}\right)$\end{document} is obtained by applying the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_Q(x)$\end{document} to the corresponding element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${c}_{ij}$\end{document} in C. This function emphasizes the strong co-expression relationships between TFs and preserves the PPI relationships. Thus, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $R$\end{document} represents a network that captures both the PPI and strong co-expression relationships among TFs.

Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${s}_i={\sum}_{j\ne i}{r}_{ij}$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal{I}=\left\{i,{s}_i\ne 0\ \right]$\end{document}, then we extracted a sub-matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${R}^{\prime }=R\left[\mathcal{I},\mathcal{I}\right]$\end{document} of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $R$\end{document} to characterize the CTN and constructed feature weights \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${w}_i=1/{s}_i$\end{document}. Using the CTN and derived weights, we regularized the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta$\end{document} in the logistic regression model (3) as follow:

(5) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} \underset{\alpha, \beta, \gamma }{\min }-\left[\frac{1}{n}{\sum}_{i=1}^n{y}_i\mathit{\log}{u}_i+\left(1-{y}_i\right)\log \left(1-{u}_i\right)\right]+\lambda{\sum}_{k\in S}{w}_k\left|{\beta}_k\right| \end{equation*}\end{document}

In practice, we utilize the glmnet R package to fit the network-regularized model by configuring the vector of weights as ‘penalty factors’ (see details in Supplementary Materials). The selection of the value Q in the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_Q(x)$\end{document} is determined by striking a balance between optimizing prediction accuracy and enhancing feature selection regularization. We have set a default Q value of 0.75, while also allowing users the flexibility to choose their own preferred Q value by providing it as a parameter in the software.

Histone code identification

Given a set of M histone modifications \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\left({h}_1,{h}_2,\dots \dots, {h}_M\right)$\end{document}, we construct a feature library, denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\theta (h)=({h}_1,{h}_2,\dots \dots, {h}_1{h}_2,\ \dots\\ \dots, {h}_1{h}_2{h}_3,\dots \dots, {h}_m^M)$\end{document}, comprising nonlinear combinations of each histone modification. This library includes all possible combinations of individual modifications and their higher-order interactions (pairs, triples, and up to combinations raised to the Mth power). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\theta (h)$\end{document} enables modeling complex interactions between histone modifications for determining the histone code of Pol II pausing. We utilized a forward feature search strategy based on ordinary logistic regression (OLR) to find the optimal subset of features. This strategy starts with a single feature and sequentially adds features, evaluating the model’s performance at each step. During each iteration, the feature that results in the highest increase in the performance metric is selected for inclusion in the model. The process continues until no further improvement is observed or a predefined stopping criterion is met. The resulting optimal feature subset contains nonlinear combinations of histone modifications considered as Pol II pausing-related histone codes. In our study, we collected a total of five histone modification ChIP-seq data that were common to both HeLa-S3 and HEK293T cell lines. Therefore, we set M = 5.

Statistical analysis

Comparison of histone modification signal intensity and chromatin openness between Pol II pausing and non-pausing sites were evaluated using the Student’s t-test. RP3 was compared with other machine learning methods also using the Student’s t-test. The statistical significance of the overlap between motif enriched TF set at Pol II pausing sites and Pol II PPI TF set were determined using hypergeometric test. P-value <.05 was considered statistically significant.

Results

Overview of RP3 - a novel machine learning method to predict Pol II pausing

To overcome the technical limitations of genome-wide Pol II pausing measurements and gain insights into its mechanisms, we introduced RP3, a novel machine learning model for predicting Pol II pausing events using multi-omics data. We formalized the modeling of Pol II pausing prediction as a classification task for the RP3 algorithm. Figure 1A illustrates the RP3 framework, where gold-standard positive samples were derived from NET-seq data, and negative samples were randomly selected non-pausing regions (details in Materials and Methods). We adopted stringent criteria in constructing a gold standard set to minimize potential biases, such as those caused by GC content (See online supplementary material for a colour version of Supplementary Materials and Supplementary Figs 14-16).

Figure 1 Overview of RP3. (A) Schematic of the RP3 model. RP3 consists of three main components: input, model and output. Both context-specific and non-context-specific data for model input, encompassing paired chromatin accessibility and gene expression, histone modifications, across-cell-type gene expression, sequence data, protein–protein interactions, and TF binding motif data, which were employed to extract genomic features. TF expression, TF binding strength, and TF expression specificity were subsequently integrated into a single feature referred to as TF activity. Genomic features include openness, TF activity, and histone mark signal strength was utilized to generate the training dataset. Gold-standard data were derived from NET-seq. The reconstruction of the TF PPI and co-expression network were used to generate the CTN for the purpose of regularization. The TFs with weights of the CTN was employed to train a regularized LR model. Refer to Table 1 for model components and mathematical notations and ‘Materials and Methods’ section for details of RP3. Model output include the pausing probability of Pol II in a genomic region and TFs potentially involved in regulating Pol II pausing. (B) Illustration of the Network regularization strategy.

RP3 integrates both context-specific (histone modification, chromatin accessibility, and gene expression for a specific cell line) and non-context-specific data (gene expression from various cell types, sequence, PPI, and TF binding motifs) as model inputs. Genomic features include histone modification signal strength (H), chromatin openness (O), TF expression (TF), TF binding strength (B), and TF expression specificity (TFS) were extracted from input dataset (Supplementary Materials). Then we combined these TF-related features and defined a genomic signature called TF activity (TFA), which was calculated as the product of B, TF, and TFS, to capture overall information of TF and reduce the feature dimension. Recognizing the cooperative nature of metazoan TFs in achieving specific functions [25], we employed our previously developed network regularization (NR) approach [24] to gain a deeper understanding of the mechanisms underlying Pol II pausing and explore biological relationships among TFs. NR facilitates feature selection by utilizing network topology information among predictors, enhancing model interpretability. Specifically, we checked two types of network relationships among TFs: PPI and co-expression networks. NR comprises the following three parts: (i) Constructing TF networks; (ii) Integrating TF networks; (iii) Constructing the candidate TF network (CTN) (Fig. 1B) (see Materials and Methods for details).

Finally, the genomic features (H, O, TFA) and TF-related feature weights derived from CTN were utilized to train a regularized logistic regression model. The model output includes the probability of the Pol II pausing event occurring in a given genomic region and identifies biologically significant TFs contributing to Pol II pausing. See Table 1 for the model components and corresponding notations.

Table 1 Model component and notations of RP3

Description of data and variables	Notation	Example	
Context-specific data			
Expression of TF	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${TF}_k:$\end{document} expression of TF k	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${TF}_{BRAC1}=20.75$\end{document} in HeLa-S3	
Accessibility of a genomic region	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${O}_i$\end{document} : degree of openness of genomic region i	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${O}_{chr6:36561836-36562037}=17.06$\end{document} in HeLa-S3	
Histone modification signal strength of a genomic region	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${H}_{i,j}:$\end{document} signal intensity of histone modifications \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} in the genomic region i	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${H}_{chr6:36561836-36562037,H3K27 ac}=25.78$\end{document} in HeLa-S3	
Pol II pausing status in a genomic region	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Y}_i:$\end{document} indicator for whether Pol II pausing in genomic region i	Pol II pausing in the genomic region chr6: 36561836–36,562,037 in HeLa-S3	
Non-context-specific data			
Interacting TFs for Pol II	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $PPI\left( Pol\ II\right)$\end{document} : set of TFs known to interact with Pol II	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $PPI\left( Pol\ II\right)$\end{document} contains BRCA1	
TFs with motif match in a genomic region	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $MB:$\end{document} set of TFs with significant motif match in the genomic region i	BRAC1 has motif match in genomic regions chr6: 36561836–36,562,037	
Motif matching strength of TF in a genomic region	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${B}_{i,k}:$\end{document} binding strength of TF k in the genomic region i	\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${B}_{chr6:36561836-36562037, BRCA1}=17.04$\end{document}	

Histone modifications, TF activity, and chromatin openness have the ability to predict Pol II pausing

Our comparative analyses between Pol II pausing and non-pausing sites revealed that histone modifications, TFs, and chromatin accessibility may potentially play a regulatory role in Pol II pausing (See online supplementary material for a colour version of Supplementary materials and Supplementary Figs 1-4). Using OLR for each feature, we assessed their predictive ability and found that histone modifications, TF activity, and chromatin openness were effective predictors of Pol II pausing in both HeLa-S3 and HEK293T cell lines. Each genomic feature demonstrated strong predictive ability, with AUROC and AUPR values exceeding 0.73 (Fig. 2). Interestingly, the predictive capacity of these features appeared consistent across cell lines, with minimal differences in AUROC values for chromatin accessibility and TF activity (< 0.01) between HeLa-S3 and HEK293T (Fig. 2A and C). The ranking of features’ ability to predict Pol II pausing remained consistent, with histone modifications showing the highest predictive power, followed by TF activity and chromatin accessibility. This consistent trend highlights the potential robustness and significance of these regulatory factors in modulating Pol II pausing.

Figure 2 Performance evaluation of RP3 within individual cell lines. (A-B) ROC and PR curves of single genomic features and RP3 in predicting Pol II pausing for HeLa-S3, respectively. (C-D) ROC and PR curves of single genomic features and RP3 in predicting Pol II pausing for HEK293T respectively.

In addition, in univariate LR, AUROC values for individual histone signals ranged from 0.5 to 0.77 in HeLa-S3 and 0.5 to 0.71 in HEK293T (See online supplementary material for a colour version of Supplementary Figs 5-6). However, the predictive accuracy of all histone signals was significantly improved when combined, with an AUROC of 0.84 in HeLa-S3 and 0.82 in HEK293T. This finding underscores the potential of combining multiple histone signals to improve the accuracy of predicting Pol II pausing, as we will discuss in detail later.

Overall, our findings provide compelling evidence for the utility of these genomic features in predicting Pol II pausing and suggest that a combination of histone modifications, TF activity, and chromatin accessibility could provide a more comprehensive understanding of Pol II pausing regulation.

Performance assessment of RP3

The RP3 model, designed for Pol II pausing prediction by integrating multiple genomic features, exhibited superior performance in both HeLa-S3 and HEK293T cell lines, as shown in Fig. 2. The ROC and PR curves for RP3 predictions demonstrated an AUROC of 0.8827 and AUPR of 0.8902 in HeLa-S3, and an AUROC of 0.8365 and AUPR of 0.8620 in HEK293T. Notably, RP3 outperformed single-type feature predictions, underscoring the significance of integrating multiple layers of omics data for enhanced accuracy.

To further validate the effectiveness of our method, we first compared RP3, which utilizes network regularization logistic regression, with ordinary logistic regression (OLR). RP3 demonstrated comparable accuracy to OLR while utilizing fewer features for prediction (See online supplementary material for a colour version of Supplementary Figs 7-8). We proceeded to compare RP3 with other three regularization methods: LASSO [26], elastic-net [27], and the adaptive LASSO. RP3 consistently outperformed LASSO, and elastic-net in selecting TF-related features. Specifically, RP3, LASSO, and elastic-net identified 13.41%, 71.86%, and 76.26% of all TFs in HeLa-S3, and 12.41%, 51.86%, and 76.16% of all TFs in HEK293T, respectively (See online supplementary material for a colour version of Supplementary Table 1). Compared to RP3, adaptive LASSO selects a much smaller number of TFs, but the prediction accuracy significantly decreases. These results suggest that RP3 enhances interpretability by reducing predictor dimensions while maintaining high predictive accuracy.

Moreover, we compared RP3 with four other machine learning classification methods, namely, Naive Bayes, LinearSVM, Decision Tree and Random Forest. Consistently, RP3 outperformed Naive Bayes, LinearSVM, and Decision Tree in predictive performance in both cell lines (See online supplementary material for a colour version of Supplementary Figs 7-8). Moreover, RP3 achieved prediction accuracy similar to Random Forest while utilizing significantly fewer TFs. We also included statistical validation through Student’s t-test, further confirming these results (See online supplementary material for a colour version of Supplementary Figs 17-18). These findings collectively underscore RP3’s superior performance over other machine learning classification methods.

In addition, independent validations with POLR2A ChIP-seq and GRO-seq data in HeLa-S3 exhibited notably higher signals at RP3-predicted pausing sites (See online supplementary material for a colour version of Supplementary Fig. 19). Extensive cross-validation across both cell lines consistently affirmed RP3’s accuracy (See online supplementary material for a colour version of Supplementary Fig. 20). These validations, combining POLR2A ChIP-seq, GRO-seq, and extensive cross-validation, substantially strengthen our computational approach’s credibility.

Cross cell line prediction

Following the observation of RP3’s excellent performance within individual cell lines, we further tested the model’s predictive ability in a different cell line by utilizing the model trained on one cell line to make predictions and assessments on another cell line. This cross-cell line prediction is a more realistic and challenging scenario, as it tests the generalizability of the model. In terms of the genomic features of histone modifications, we utilized five histone marks (H3K27ac, H3K36me3, H3K4me1, H3K4me2, and H3K4me3) that are shared by both HeLa-S3 and HEK293T cell lines for the purpose of cross-cell line prediction. Our results showed that RP3 trained with data from one cell type can accurately predict Pol II pausing specific to another cell type (Fig. 3). Notably, the performance of RP3 in cross-cell line prediction was only slightly lower than within individual cell types, with differences in AUROC and AUPR values remaining <0.02 for both HeLa-S3 and HEK293T (Fig. 3).

Figure 3 Performance evaluation of RP3 across cell lines. (A-B) ROC and PR curves for across cell line prediction of RP3 between HeLa-S3 and HEK293T, respectively. HeLa-S3 and HEK293T represent the performance of RP3 within individual cell lines. HeLa-S3_HEK293T and HEK293T_HeLa-S3 depict the across cell line performance of models trained on data from the former and evaluated on data from the latter. (C-D) Comparison of POLR2A ChIP-seq signal intensity at regions predicted as Pol II pausing sites by the model trained on HeLa-S3, with regions predicted as non-pausing sites in GM12878 and MCF7 cell lines, respectively. P-values were determined using Student’s t-tests. *** indicted P-value <.001.

To further validate the generalizability of RP3, we applied the model trained on HeLa-S3 and HEK293T to predict open regions in two other cell lines, specifically GM12878 and MCF7. As GM12878 and MCF7 lacked corresponding NET-seq data, we employed ChIP-seq data for POLR2A, the largest subunit of Pol II, as an alternative strategy to validate Pol II pausing predictions. As illustrated in Fig. 3C and D, the regions predicted as Pol II pausing sites by the model trained on HeLa-S3 showed significantly higher POLR2A ChIP-seq signal intensity compared to the predicted non-pausing sites. Similar results were observed with the model trained on HEK293T (See online supplementary material for a colour version of Supplementary Fig. 21). These findings demonstrating a good level of generalizability of our model.

RP3 uncovers biologically meaningful TFs associated with Pol II pausing

We hypothesize that TFs participate in the regulation of pausing through physical PPI with Pol II. As depicted in Fig. 4A, the original PPI network reveals the presence of 455 TFs that interact with Pol II through both first-order and second-order PPIs (Fig. 4C). To identify biologically meaningful TFs associated with Pol II pausing, PR3 takes advantage of the network regularization strategy for TF-related genomic feature selection. The results demonstrate that RP3 effectively outputs 13.4% and 14.7% of the TFs from the original PPI network in HeLa-S3 and HEK293T cell lines, respectively, leading to improved biological interpretability (Fig. 4B and C). Moreover, the comparative analysis with OLR and other regularization methods confirms RP3’s outstanding proficiency in feature selection, providing compelling evidence for its ability to identify key TFs associated with the regulatory mechanism of Pol II pausing. Notably, the majority of the identified TFs by RP3 in both HeLa-S3 and HEK293T cell lines are shared, indicating these TFs may possess conserved regulatory roles or functional relevance in the context of Pol II pausing (Fig. 4B).

Figure 4 Identification of TFs potentially implicated in the regulation of Pol II pausing. (A) The second-order TF protein–protein interaction network of Pol II. (B) The TF network identified by the RP3 for Pol II pausing in HeLa-S3 and HEK293T cell lines represents a subset of the Pol II TF protein–protein interaction network. TAL1 is specific to the HeLa-S3 cell line, while PAX7, MYOD1, GATA1, LMO2, ESR1, GRHL2, and EGR1 are specific to the HEK293T cell line. The remaining TFs reported by RP3 are shared between both HeLa-S3 and HEK293T. The two TFs connected by a solid line edge demonstrate both protein–protein interactions and co-expression across various cell types, with a Pearson correlation coefficient of gene expression >0.75. Two TFs connected by a dotted edge represent only exhibit protein–protein interactions. (C) The bar plot displays the number of TFs present in Pol II TF PPI network and RP3 output in HeLa-S3 and HEK293T. (D) The pie chart represents the proportion of TFs exported by RP3 that are motif enriched in Pol II pause sites. (E) The number of Pol II pausing sites (Positive) and non-pausing sites (Negative) that overlap with the transcription factor’s ChIP-seq peak. Ratio represents the fold change between Positive and Negative.

To validate the effectiveness of TFs identified by RP3, we conducted motif enrichment analysis using Homer software. The analysis revealed that 65.0% and 61.2% of TFs identified by RP3 exhibited significant enrichment on Pol II pausing sites in HeLa-S3 and HEK293T cell lines, respectively (Fig. 4D). Subsequently, we collected independent TF ChIP-seq data from the data-rich HeLa-S3 cell line and conducted a comparative analysis to assess the overlap between TF binding and both Pol II pausing and non-pausing positions. For the 13 TFs identified by RP3 with available ChIP-seq data, there was a significantly higher overlap between Pol II pausing positions and TF ChIP-seq peaks compared to non-pausing positions. The fold change ranged from 4.76 to 163.45, indicating a strong association between TF binding and Pol II pausing (Fig. 4E). Recent studies suggest that promoter-proximal paused Pol II plays a crucial role in enhancing enhancer-promoter 3D interactions [28], with TFs acting as key regulators in 3D genome organization [29–31]. We thus checked the Hi-C data of HeLa-S3 and observed that multiple TFs identified by RP3 were concurrently bound to both the left and right anchors of certain chromatin loops associated with the Pol II pausing (See online supplementary material for a colour version of Supplementary Figs 9-11). For example, in the loop from chr11:9775000–9800000 to chr11:10325000–10350000, the binding peaks of 11 out of the 13 TFs from the collected ChIP-seq datasets were simultaneously detected on both the left and right anchors (See online supplementary material for a colour version of Supplementary Fig. 9). These findings suggest that the TFs identified by RP3 may exert a critical influence over 3D genome architecture, particularly in the context of Pol II pausing. Further research may be necessary to fully understand the functional implications of these TFs in shaping the 3D genome and their role in gene regulation.

Additionally, we conducted an extensive literature search and found multiple TFs outputted by RP3 that have been previously reported to be associated with the regulation of Pol II pausing. For instance, TRIM28 is implicated in regulating Pol II pausing and transcriptional elongation in numerous mammalian genes, and CTCF promotes Pol II pausing, linking DNA methylation to splicing processes. A total of 27 TFs identified by RP3, including BRCA1, MYC, EGR1, etc., have literature evidence supporting their association with the regulation of Pol II pausing (See online supplementary material for a colour version of Supplementary Table 2). In summary, RP3 unveils biologically significant TFs associated with Pol II pausing, offering valuable insights into the underlying mechanism of Pol II pausing. Pol II pausing at distinct positions may be governed by diverse TFs and combinations of TFs, highlighting the intricate nature of its regulatory mechanism (See online supplementary material for a colour version of Supplementary Fig. 11).

Identification of histone code for pol II pausing

Multivalent histone modifications provide a powerful mechanism for regulating gene expression and Pol II pausing. Combinations of multiple modifications on the same histone tail create a complex histone code, exerting precise regulatory functions [32]. To decode the histone code associated with Pol II pausing, we constructed a feature library encompassing all nonlinear combinations of five histone marks (H3K27ac, H3K36me3, H3K4me1, H3K4me2, and H3K4me3) shared between HeLa-S3 and HEK293T cell lines (Fig. 5A). Using a forward feature search strategy based on OLR (Materials and Methods), we identified the optimal subset of features that exhibit the strongest association with the histone code linked to Pol II pausing (Fig. 5A). As illustrated in Fig. 5B and C, the top three features demonstrated high predictive power, close to the maximum possible, in both HeLa-S3 and HEK293T cell lines (See online supplementary material for a colour version of Supplementary Figs 6 and 12). Two of the top three histone modification combinations, (H3K27ac, H3K4me2) and the quadratic of H3K36me3, were shared between the cell lines, indicating their critical and conserved role in regulating Pol II pausing across cellular contexts.

Figure 5 Deciphering the histone code associated with Pol II pausing. (A) Schematic illustration of the identification of histone code for Pol II pausing. (B-C) Plot representing histone combination features versus AUROC values derived from forward feature selection algorithm in HeLa-S3 and HEK293T, respectively. (D) Genome browser snapshot illustrates examples of various histone modification combinations at Pol II pausing sites. Different-colored bars represent distinct combinations of histone modifications at various Pol II pausing sites. This snapshot was generated using the WashU Epigenome Browser (http://epigenomegateway.wustl.edu/browser/).

H3K27ac and H3K4me2 can co-occur in the same region of chromatin, resulting in a ‘bivalent’ modification state. This bivalent state is thought to play a role in gene regulation during development and differentiation [33]. Studies have suggested that H3K27ac modification can facilitate Pol II release from the promoter and promote productive elongation of transcription. On the other hand, the presence of H3K4me2 modification has been linked to the maintenance of Pol II pausing stability, as it promotes the recruitment of negative elongation factors that pause Pol II [18]. Thus, the co-occurrence of H3K27ac and H3K4me2 modifications in the same chromatin region can create a bivalent state that regulates the balance between Pol II pausing and transcriptional elongation. The quadratic relationship observed for H3K36me3 indicates a more intricate and nonlinear association with Pol II pausing. The precise mechanism by which H3K36me3 affects the regulation of Pol II pausing remains unclear. However, H3K36me3 is a well-known epigenetic mark for transcription elongation [34]. This suggests that H3K36me3 may be involved in promoting productive elongation and reducing Pol II pausing. In HeLa-S3, H3K4me3 also exhibits a quadratic relationship with Pol II pausing, similar to H3K36me3. Recent studies have provided evidence indicating that the acute depletion of H3K4me3 leads to a widespread decrease in transcriptional output, an increase in Pol II pausing, and slower elongation. H3K4me3 participates in the regulation of Pol II pause release by recruiting the integrator complex subunit 11 (INTS11) [35]. Furthermore, the bivalent modification of H3K4me2 and H3K36me3 in HEK293T may possess analogous functions to the bivalent modification of H3K4me2 and H3K27ac regulating the delicate balance between Pol II pausing maintenance and release.

Pol II pausing at distinct genomic locations may be regulated by a variety of histone combinations (Fig. 5D). The intricate interplay of histone modifications, including bivalent states and quadratic relationships, elucidates the complexity of the regulatory mechanisms governing Pol II pausing. These findings provide valuable insights into the regulatory role of multivalent histone modifications in Pol II pausing, advancing our understanding of the complicated interplay between histone modifications and the regulation of gene transcription.

RP3 predicts Pol II promoter-proximal pausing more accurately

Pol II promoter-proximal pausing is a prevalent mechanism that governs the dynamics of gene expression, and current research in metazoans primarily focused on this phenomenon [3, 11]. Hence, we proceeded to assess the efficacy of RP3 in predicting Pol II promoter-proximal pausing. Specifically, when a Pol II pausing site was positioned within 300 nucleotides downstream of the transcription start site (TSS), it was categorized as ‘promoter-proximal’ [3]. We detected 11% of all identified pause sites as ‘promoter-proximal’ in both HeLa-S3 and HEK293T cell lines (See online supplementary material for a colour version of Supplementary Fig. 13). As depicted in Fig. 6 and Supplementary Fig. 25, RP3 exhibits superior accuracy in predicting promoter-proximal pausing compared with non-promoter-proximal pausing. In both HeLa-S3 and HEK293T, the AUROC and AUPR values for RP3’s predictions of promoter-proximal pausing were consistently close to 1, showing a substantial contrast with the predictions of non-promoter-proximal pausing, which differed by > 0.1. It is noteworthy that histone modifications exhibit the highest predictive ability for both promoter-proximal and non-proximal pausing. In contrast, the predictive abilities of chromatin openness and TF activity were notably reduced for non-proximal promoter pausing when compared to proximal-promoter pausing. This suggests that the regulatory mechanisms for pausing events in different genomic regions may vary, and histone modifications are particularly relevant and robust across both types of pausing events.

Figure 6 RP3 predicts promoter-proximal and non-proximal Pol II pausing. (A-B) ROC and PR curves of single genomic features and RP3 in predicting promoter-proximal Pol II pausing for HeLa-S3, respectively. (C-D) ROC and PR curves of single genomic features and RP3 in predicting non-promoter-proximal Pol II pausing for HeLa-S3, respectively.

Discussions and conclusions

Po II pausing is key regulatory step for implementing precise gene expression programs and holds critical significance in a multitude of essential biological processes [11, 36]. However, the genome-wide profiling of Pol II pausing is currently accessible only in a limited number of cell lines, as experimental measurement techniques remain laborious and costly. Therefore, we developed RP3, leverages a network-regularized logistic regression approach to predict Pol II transcriptional elongation pausing. Our model advances the field by integrating multi-omics data, including genome sequence, histone modifications, gene expression, chromatin accessibility, and protein–protein interactions, offering a comprehensive perspective on Pol II pausing. Feature ablation analysis indicates the effectiveness of each input feature (See online supplementary material for a colour version of Supplementary Materials and Supplementary Figs 22-24). This comprehensive integration represents a significant advancement over existing models, which often rely on fewer data types. RP3 demonstrates outstanding predictive performance for Pol II pausing in diverse cellular contexts. Uniquely, by utilizing the network regularization strategy to incorporate topological information from TF PPI and across cell lines co-expression networks for feature selection, RP3 can effectively identify biologically meaningful TFs associated with Pol II pausing. RP3 outperforms regularization methods such as LASSO, elastic-net, and adaptive LASSO in feature selection, maintaining high predictive accuracy while reducing predictor dimensions for improved interpretability. It also exhibits superior performance compared to other classification methods, including Naive Bayes, LinearSVM, Decision Trees, and Random Forests. Additionally, we systematically investigated how interplay between different histone modification signals in orchestrating the regulation of Pol II pausing. Our extensive examination of the ENCODE and Roadmap databases reveals 37 distinct cell lines with datasets that meet the input requirements for RP3, indicating its broad applicability in various scenarios (See online supplementary material for a colour version of Supplementary Table 3).

Our work not only offers valuable insights into the regulatory mechanisms of transcriptional pausing but also provides a robust predictive model in computational biology. However, several limitations should be acknowledged. First, due to the limited availability of other histone ChIP-seq data, our study concentrated on a subset of histone modifications, specifically focusing on five key modifications. As more data become accessible, we hope to explore a more extensive range of histone modifications to unravel a comprehensive histone code associated with Pol II pausing. Second, although our model is effective in identifying Pol II pausing sites, it does not fully capture the dynamic and temporal aspects of gene regulation, which are crucial for a comprehensive understanding of transcriptional processes. Several aspects of our research are worthy of further exploration in future work. Firstly, our study elucidates the presence of multiple regulatory factors governing Pol II pausing, underscoring the intricate and sophisticated nature of its regulation. Further exploration of the interplay between these regulators can help us reveal more precise underlying mechanisms that regulate Pol II pausing. Secondly, recent studies have demonstrated the involvement of additional regulatory factors, such as R-loops formed by nascent RNA [37] and enhancer-promoter interaction [38, 39], in the inducing RNA polymerase pausing. Incorporating information about other regulatory factors promising for enhancing the predictive accuracy and interpretability of the model. Lastly, comprehensively elucidating the molecular mechanisms underlying Pol II pausing and its impact on important biological processes, such as differentiation [12], development [11], inflammation [13], and tumorigenesis [15], can greatly advance our understanding of gene transcription regulation and facilitate the identification of therapeutic targets for disease.

Key Points

We developed RP3, a machine learning framework for predicting RNA Polymerase II (Pol II) transcriptional elongation pausing, by integrating genomic sequence, histone modification, gene expression, chromatin accessibility, and protein–protein interaction data.

By leveraging network regularization to embed topological information from TF PPI and co-expression networks for feature selection, RP3 enables effective identify biologically relevant transcription factors involved in Pol II pausing.

We systematically identified the combination of histone modification signals associated with Pol II pausing.

Supplementary Material

Supplementary_Figures_bbae246

Supplementary_Materials_bbae246

Funding

The National Key Research and Development Program of China (2022YFA1004800 and 2020YFA0712402), CAS Project for Young Scientists in Basic Research (YSBR-077), and the National Natural Science Foundation of China (12025107 and 12371481).

Data and code availability

All the datasets used in this study are publicly available and accessible. Detailed data sources are described in the ‘Data collection’ section of ‘Materials and methods.’ Source codes are available at https://github.com/AMSSwanglab/RP3.

Author Biographies

Lixin Ren is a Ph.D candidate at University of Science and Technology Beijing. His research interests include bioinformatics, machine learning, and computational biology.

Wanbiao Ma is a professor at University of Science and Technology Beijing. His research interests are mathematical biology.

Yong Wang is a professor at Chinese Academy of Sciences. His research interests are computational systems biology.
==== Refs
References

1. Cramer  P . Organization and regulation of gene transcription. Nature  2019;573 :45–54.31462772
2. Djebali  S, Davis  CA, Merkel  A, et al.  Landscape of transcription in human cells. Nature  2012;489 :101–8.22955620
3. Gajos  M, Jasnovidova  O, van  Bömmel  A, et al.  Conserved DNA sequence features underlie pervasive RNA polymerase pausing. Nucleic Acids Res  2021;49 :4402–20.33788942
4. Shandilya  J, Roberts  SG. The transcription cycle in eukaryotes: from productive initiation to RNA polymerase II recycling. Biochim Biophys Acta  2012;1819 :391–400.22306664
5. Chen  FX, Smith  ER, Shilatifard  A. Born to run: control of transcription elongation by RNA polymerase II. Nat Rev Mol Cell Biol  2018;19 :464–78.29740129
6. Jonkers  I, Lis  JT. Getting up to speed with transcription elongation by RNA polymerase II. Nat Rev Mol Cell Biol  2015;16 :167–77.25693130
7. Mayer  A, Churchman  LS. Genome-wide profiling of RNA polymerase transcription at nucleotide resolution in human cells with native elongating transcript sequencing. Nat Protoc  2016;11 :813–33.27010758
8. Adelman  K, Lis  JT. Promoter-proximal pausing of RNA polymerase II: emerging roles in metazoans. Nat Rev Genet  2012;13 :720–31.22986266
9. Williams  LH, Fromm  G, Gokey  NG, et al.  Pausing of RNA polymerase II regulates mammalian developmental potential through control of signaling networks. Mol Cell  2015;58 :311–22.25773599
10. Chen  F, Gao  X, Shilatifard  A. Stably paused genes revealed through inhibition of transcription initiation by the TFIIH inhibitor triptolide. Genes Dev  2015;29 :39–47.25561494
11. Abuhashem  A, Garg  V, Hadjantonakis  A-K. RNA polymerase II pausing in development: orchestrating transcription. Open Biol  2022;12 :210220.34982944
12. Abuhashem  A, Chivu  AG, Zhao  Y, et al.  RNA pol II pausing facilitates phased pluripotency transitions by buffering transcription. Genes Dev  2022;36 :770–89.35981753
13. Yu  L, Zhang  B, Deochand  D, et al.  Negative elongation factor complex enables macrophage inflammatory responses by controlling anti-inflammatory gene expression. Nat Commun  2020;11 :2286.32385332
14. Anand  P, Brown  JD, Lin  CY, et al.  BET bromodomains mediate transcriptional pause release in heart failure. Cell  2013;154 :569–82.23911322
15. Rahl  PB, Lin  CY, Seila  AC, et al.  c-Myc regulates transcriptional pause release. Cell  2010;141 :432–45.20434984
16. Churchman  LS, Weissman  JS. Nascent transcript sequencing visualizes transcription at nucleotide resolution. Nature  2011;469 :368–73.21248844
17. Feng  P, Xiao  A, Fang  M, et al.  A machine learning-based framework for modeling transcription elongation. Proc Natl Acad Sci  2021;118 :118.
18. Hu  S, Song  A, Peng  L, et al.  H3K4me2/3 modulate the stability of RNA polymerase II pausing. Cell Res  2023;33 :403–6.
19. Kapoor-Vazirani  P, Kagey  JD, Vertino  PM. SUV420H2-mediated H4K20 trimethylation enforces RNA polymerase II promoter-proximal pausing by blocking hMOF-dependent H4K16 acetylation. Mol Cell Biol  2011;31 :1594–609.21321083
20. Gates  LA, Shi  J, Rohira  AD, et al.  Acetylation on histone H3 lysine 9 mediates a switch from transcription initiation to elongation. J Biol Chem  2017;292 :14456–72.28717009
21. Shukla  S, Kavak  E, Gregory  M, et al.  CTCF-promoted RNA polymerase II pausing links DNA methylation to splicing. Nature  2011;479 :74–9.21964334
22. Mayer  A, Di Iulio  J, Maleri  S, et al.  Native elongating transcript sequencing reveals human transcriptional activity at nucleotide resolution. Cell  2015;161 :541–54.25910208
23. Colin  J, Candelli  T, Porrua  O, et al.  Roadblock termination by reb1p restricts cryptic and readthrough transcription. Mol Cell  2014;56 :667–80.25479637
24. Ren  L, Gao  C, Duren  Z, et al.  GuidingNet: revealing transcriptional cofactor and predicting binding for DNA methyltransferase by network regularization. Brief Bioinform  2021;22 :bbaa245.33048117
25. Lambert  SA, Jolma  A, Campitelli  LF, et al.  The human transcription factors. Cell  2018;172 :650–65.29425488
26. Tibshirani  R . Regression shrinkage and selection via the lasso. J R Stat Soc B Methodol  1996;58 :267–88.
27. Zou  H, Hastie  T. Regularization and variable selection via the elastic net. J R Stat Soc Series B Stat Methodology  2005;67 :301–20.
28. Barshad  G, Lewis  JJ, Chivu  AG, et al.  RNA polymerase II dynamics shape enhancer–promoter interactions. Nat Genet  2023;55 :1370–80.37430091
29. Stadhouders  R, Filion  GJ, Graf  T. Transcription factors and 3D genome conformation in cell-fate decisions. Nature  2019;569 :345–54.31092938
30. Kim  S, Shendure  J. Mechanisms of interplay between transcription factors and the 3D genome. Mol Cell  2019;76 :306–19.31521504
31. Ren  L, Ma  W, Wang Y. SpecLoop predicts cell type-specific chromatin loop via transcription factor cooperation, Comput Biol Med  2024;171 :108182.
32. Arita  K, Isogai  S, Oda  T, et al.  Recognition of modification status on a histone H3 tail by linked histone reader modules of the epigenetic regulator UHRF1. Proc Natl Acad Sci  2012;109 :12950–5.22837395
33. Fox  S, Myers  JA, Davidson  C, et al.  Hyperacetylated chromatin domains mark cell type-specific genes and suggest distinct modes of enhancer function. Nat Commun  2020;11 :4544.32917861
34. Wagner  EJ, Carpenter  PB. Understanding the language of Lys36 methylation at histone H3. Nat Rev Mol Cell Biol  2012;13 :115–26.22266761
35. Wang  H, Fan  Z, Shliaha  PV, et al.  H3K4me3 regulates RNA polymerase II promoter-proximal pause-release. Nature  2023;615 :339–48.
36. Mayer  A, Landry  HM, Churchman  LS. Pause & go: from the discovery of RNA polymerase pausing to its functional implications. Curr Opin Cell Biol  2017;46 :72–80.28363125
37. Skourti-Stathaki  K, Proudfoot  NJ, Gromak  N. Human senataxin resolves RNA/DNA hybrids formed at transcriptional pause sites to promote Xrn2-dependent termination. Mol Cell  2011;42 :794–805.21700224
38. Schaaf  CA, Kwak  H, Koenig  A, et al.  Genome-wide control of RNA polymerase II activity by cohesin. PLoS Genet  2013;9 :e1003382.23555293
39. Fay  A, Misulovin  Z, Li  J, et al.  Cohesin selectively binds and regulates genes with paused RNA polymerase. Curr Biol  2011;21 :1624–34.21962715
