pulseTD: RNA life cycle dynamics analysis based on pulse model of 4sU-seq time course sequencing data

Description of the GEO dataset

The 4sU-seq experimental datasets were obtained from GEO database (https://www.ncbi.nlm.nih.gov/geo/). The first dataset was an RNA-seq dataset from mouse DC (GSE56977). RNA was sampled from mouse DC every 15 min for the first 3 h (total 13 samples) of their response to LPS, followed by a short (10 min) metabolic marker pulse(4sU) before the sampling time point (Rabani et al., 2014). The samples of GSE59717 dataset were the infection of primary human foreskin fibroblasts with a wild-type simplex virus strain 17 with a multiplicity of infection of 10 (Rutkowski et al., 2015).

The algorithm flow of pulseTD model

The rates of RLCD within 24 h of the pulse state were evaluated, which represented a complete RNA life cycle, and the mode of oscillation was discharged (La Manno et al., 2018). We assumed that the rates of transcription, processing and degradation conforms to the functional form of the pulse model in the pulse state. Expression value data of pre-mRNA and mRNA at any time point could be obtained by 4su-seq and RNA-seq. The parameters were estimated by using an optimization algorithm. Next, we introduced the principles of the model and the use of the software. After that we created simulation data and evaluate performance.

First, we used the R package (GenomicAlignments, GenomicFeatures, etc.) to analyze RNA-seq and 4sU-seq bam files, in order to quantify intronic (reads alignment to the intron region) and exonic (reads alignment to the exon region) in Reads Per Kilobase of exon model per Million mapped reads (RPKM), Transcripts Per Kilobase of exon model per Million mapped reads (TPM) or Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) of each gene. Some methods (min–max normalization, log normalization) had been added to standardize expression profile data before evaluating RLCD. In the cell life cycle, transcription continuously generates new pre-mRNA expressed by $P\left(t\right)$. Processing converts newly generated pre-RNA into mature mRNA expressed by $M\left(t\right)$. Finally, mRNA is targeted for degradation (Moore & Proudfoot, 2009). The equation according with the process of RLCD as follows: (1.1) $$\{\begin{array}{c}\stackrel{\cdot }{P}=\mathrm{\alpha }-\mathrm{\gamma }\\ \stackrel{\cdot }{M}=\mathrm{\gamma }-\mathrm{\beta }\end{array}$$

Among them, $\mathrm{\alpha },\mathrm{\beta },\mathrm{\gamma }$ represent the expression level of synthesized, degraded and processed RNA per unit time, respectively. Since the labeling time ($t_L$) is very short during the 4sU experimental, we assumed that RNA was not degraded during experimental labeling time. The total labeled RNA expression level $T_L\left(t\right)$ can be expressed as: (1.2) $$\{\begin{array}{c}{\stackrel{\cdot }{T}}_L=\mathrm{\alpha }\\ \mathrm{\beta }=0\end{array}$$

The transcription level of intracellular RNA was in a steady state without being stimulated by outside world. When some factors (transient pulse stimulation) disrupted the equilibrium, the rates of genes at different stages were changed to cushion the stimulus. After a period of time, most cells returned to a steady state due to some factors such as cell morphology and stress. As mentioned above, the change in gene expression values was mainly the result of the combination of transcription, processing and degradation. Therefore, external stimulus conditions, pulse stimuli directly affected RLCD. To this end, we hypothesized that the three processes of gene expression were transcription $\mathrm{\alpha }$, processing $\mathrm{\gamma }$, degradation $\mathrm{\beta }$ and $\mathrm{\alpha }=\mathrm{\beta }$, $\mathrm{\gamma }=\mathrm{\alpha }$ in steady state, which were broken when the external stimulus pulse was stimulated. After a period of time, it reached steady state again, and ${\mathrm{\alpha }}^{{}^{\prime }}={\mathrm{\beta }}^{{}^{\prime }}$, ${\mathrm{\gamma }}^{{}^{\prime }}={\mathrm{\alpha }}^{{}^{\prime }}$. The functional form of $\mathrm{\alpha },\mathrm{\beta },\mathrm{\gamma }$ needed to be determined.

In this case, the rates of transcription, processing and degradation changed from $\mathrm{\gamma }~\mathrm{\alpha }~\mathrm{\beta }$ to ${\mathrm{\gamma }}^{\mathrm{\prime }}~{\mathrm{\alpha }}^{\mathrm{\prime }}~{\mathrm{\beta }}^{\mathrm{\prime }}$, which could be considered as a process of RNA from a steady state to a new steady state under pulsed stimulation. So, we used the pulse model $f\left(x\right)$ (Chechik & Koller, 2009) to represent functional form of $\mathrm{\alpha },\mathrm{\beta },\mathrm{\gamma }$: (1.3) $$f_{\mathrm{\theta }}\left(x\right)={\displaystyle \frac{1}{h_1}\left(h_0+\left(h_1-h_0\right){\displaystyle \frac{1}{1+e^{-\mathrm{\beta }\left(x-t_1\right)}}}\right)\left(h_2+\left(h_1-h_2\right){\displaystyle \frac{1}{1+e^{\mathrm{\beta }\left(x-t_2\right)}}}\right)}$$

The $\mathrm{\theta }$ is the parameter vector represented by $\left(h_0,h_1,h_2,t_1,t_2,\mathrm{\beta }\right)$. $h_0,h_1,h_2$ represent the initial state rates value, the peak time rates value and the new steady state rates value. $t_1,t_2$ are the maximum time for the first and second rise or fall changes, and $\mathrm{\beta }$ is the slope of the both transitions.

Taking into account the effects of existing RNA, we first estimated the global normalization factors in the model. We assumed that real mRNA $R$ levels were composed of the labeled mRNA $S$ and pre-existing mRNA $N$: (1.4) $$R=S+N$$

Make the total mRNA observations proportional to $S$ and $N$, and the scale factors were both$w$. At the same time, we assumed that the labeled mRNA observations included two parts, labeled $S$ and unlabeled $N$, with scale factors of $w_1$ and $w_2$, respectively. Then: (1.5) $$\{\begin{array}{c}w\left(S+N\right)=O_T=T(t)\hfill \\ w_{1}S+w_{2}N=O_L=T_L(t)\hfill \end{array}$$

Three scale factors are constant in one sample, and $O_T,O_L$ represent observations for total RNA and labeled RNA, respectively. According to formula (1.5): (1.6) $$\left(w_1-w_2\right)S+{\displaystyle \frac{w_2}{w}{O}_T=O_L}$$

Here $P(t)$, $T(t)$ and $T_L(t)$ are represented as a linear combination of pulse function integrals. According to formula (1.1~1.6): (1.7) $$\{\begin{array}{l}P\left(t\right)=P\left(0\right)+{\displaystyle \underset{0}{\overset{t}{\int }}{f}_{{\mathrm{\theta }}_{\mathrm{\alpha }}}}\left(t\right)dt-{\displaystyle \underset{0}{\overset{t}{\int }}{f}_{{\mathrm{\theta }}_{\mathrm{\gamma }}}}\left(t\right)dt\\ T\left(t\right)=T\left(0\right)+{\displaystyle \underset{0}{\overset{t}{\int }}{f}_{{\mathrm{\theta }}_{\mathrm{\alpha }}}}\left(t\right)dt-{\displaystyle \underset{0}{\overset{t}{\int }}{f}_{{\mathrm{\theta }}_{\mathrm{\beta }}}}\left(t\right)dt\\ T_L\left(t\right)=\frac{w_2}{w}T\left(t\right)+\left(w_1-w_2\right){\displaystyle \underset{t-t_L}{\overset{t}{\int }}{f}_{{\mathrm{\theta }}_{\mathrm{\alpha }}}}\left(t\right)dt\end{array}$$

Where $P\left(t\right),T\left(t\right),T_L\left(t\right)$ are the observation data of the measurement time node, $t_L$ is the labeled time, ${\mathrm{\theta }}_x=\left(h_0^x,h_1^x,h_2^x,t_1^x,t_2^x,{\mathrm{\beta }}^x\right)$, so the parameter vector $X=\left({\mathrm{\theta }}_{\mathrm{\alpha }},{\mathrm{\theta }}_{\mathrm{\gamma }},{\mathrm{\theta }}_{\mathrm{\beta }}\right)$. The total objective function is (1.8) $$J\left(X\right)={\displaystyle \frac{1}{2}\sum _{j=0}^{t_{end}}{\left({\hat{P}}_j-N{P}_j\right)}^2+{\left({\hat{T}}_j-N{T}_j\right)}^2+{\left({{\hat{T}}_L}_j-N{T_L}_j\right)}^2}$$

Among them $\hat{P},\hat{T},{\hat{T}}_L$ are the model prediction value and $NP,NT,N{T}_L$ are the standardized observation data. Therefore, $\mathrm{\alpha },\mathrm{\beta },\mathrm{\gamma }$ can be solved by the following constrained optimization problem: $$X=\arg minJ\left(X\right)$$ $$\mathrm{s}.\mathrm{t}.X>0,$$ $${\mathrm{h}}_0^{\mathrm{\alpha }}~{\mathrm{h}}_0^{\mathrm{\beta }}~{\mathrm{h}}_0^{\mathrm{\gamma }},$$ $${\mathrm{h}}_2^{\mathrm{\alpha }}~{\mathrm{h}}_2^{\mathrm{\beta }}~{\mathrm{h}}_2^{\mathrm{\gamma }},$$

First, we evaluated the global normalization factor using a gradient descent method with random initialization parameters. We minimized the objective function by using the nlminb method, which was available in the stats R packages. The best fit result was chosen from the results of 100 (default) random initial values. At the same time, the chi-square test of goodness of fit was used to verify the statistical validity of pulseTD.

The interpretability of the model output

Here, the rates of transcription, processing and degradation were defined as the expression level of RNA transcribed, processed or degraded per unit time (unit: RNA/min or RNA/hour). In order to increase the interpretability of pulseTD and make it easier to compare with other tools, a conversion method was applied to pulseTD output. Let $\mathrm{\gamma }={\mathrm{\gamma }}^{k}P$ and $\mathrm{\beta }={\mathrm{\beta }}^k\left(T-P\right)$ at any time t. We assumed that the processing rates were directly proportional to the pre-mRNA expression and the ratio was ${\mathrm{\gamma }}^k$. Similarly, degradation rates were directly proportional to mature mRNA expression value, the ratio was ${\mathrm{\beta }}^k$. According to formula (1.1): (1.9) $$\{\begin{array}{c}\stackrel{\cdot }{P}=\mathrm{\alpha }-{\mathrm{\gamma }}^{k}P\\ \stackrel{\cdot }{T}=\mathrm{\alpha }-{\mathrm{\beta }}^k\left(T-P\right)\end{array}$$

This was similar to previous researches (De Pretis et al., 2015; Rabani et al., 2011, 2014; Sun et al., 2012). $\mathrm{\alpha }$ represents the transcription rates in units of mRNA/min. Where $\stackrel{\cdot }{P}$ and $\stackrel{\cdot }{T}$ represent the derivatives of the functions P(t) and T(t) at time t, respectively. And they are functions of time t. Here we used $\mathrm{\alpha },{\mathrm{\gamma }}^k,{\mathrm{\beta }}^k$ to compare with other methods, and it was easy to explain.

Generation of simulation data

To evaluate the effectiveness of the model, we generated simulation data for 1,000 genes by randomly drawing from a specific distribution. Transcription, processing, and degradation rates were first generated, as well as randomly generated scale factors. The simulation expression value was then created based on the rates using Runge–Kutta method.

First of all, based on previous researches, we knew that the half-life of RNA was considered to follow a normal distribution (Friedel et al., 2009). Next, we determined the distribution of the transcription rates based on the mean ($\mathrm{\mu }$) and variance ($\mathrm{\sigma }$) of the observations. And the transcription rates were randomly extracted t from the normal distribution $N\left(\mathrm{\mu },{\mathrm{\sigma }}^2\right)$. We had expected to determine the degradation rates in the same way. However, the degradation rates and the transcription rate were dependent. To simulate this dependency, we evaluated the correlation between transcription and degradation rates at any time based on the pearson correlation coefficient, which was $k_t$. Therefore, the degradation rates approximately obeyed the $N\left(\mathrm{\mu }/{\mathrm{k}}_{\mathrm{t}},{\left(\mathrm{\sigma }/{\mathrm{k}}_{\mathrm{t}}\right)}^2\right)$ distribution. Similarly, processing rates were randomly drawn in the same way. Finally, we also randomly generated global scale factors, which were used to simulate existing and newly generated RNA.

After determining all rates, we used the fourth-order Runge–Kutta method of the R package (deSolve) to evaluate the expression levels of pre-mRNA, mature mRNA and labeled RNA as a simulation data set. Among them, the initial value of the expression value was randomly selected from the distribution of the observed data.

In general, we determined the rate distribution based on the experimental data and generated simulated data parameters based on the mean and variance of the rate. Here, the time range of the experimental analysis was 0–180 min, sampling was performed every 15 min, and the 4sU marking was performed 10 min before sampling. Finally, we got a simulation data set of 1,000 genes.

Software framework and description

The pulseTD combines the pulse model to predict the steady state of the RLCD. We defined the rates of transcription processing and degradation as a pulse function which had 6 parameters, a total of 18 parameters. To standardize RNA expression levels, additional global scale factors needed to be evaluated. For the evaluation of each gene, 100 random initializations were required. We used a multi-threaded approach to reduce runtime. The software supports different ways to evaluate expression levels such as counts, RPKM, TPM, FPKM. The workflow of pulseTD software is as follows (Fig. 1): (i) The RNA-seq and 4sU-seq data are aligned to the reference genome, and the results are used as the input files of the software. Expression values of pre-RNA, mRNA and labeled-RNA can be calculated in the form of RPKM, corresponding to the R function named extimateExpression. (ii) Subsequently, the expression profile is used to optimize the parameters of the pulse model. This is an optimization problem with six parameters, which are determined by MSE (minimizing the mean square error). The extimateParams function can estimate the pulse parameters of the transcription, processing and degradation rates. (iii) The parameters can be re-estimated using the correctionParams function because of the influence of random initial values. (iv) Next, transcriptional dynamic rates are solved flexibly, including the rates of transcription, processing and degradation, or the steady state rates are predicted. Complete guide reference documentation.

Compared with other methods

Several published studies had revealed RNA RLCD from different perspectives. We compared the following common software tools: INSPEcT (De Pretis et al., 2015), DRiLL (Rabani et al., 2014), DTA (Schwalb et al., 2012) and pulseR (Uvarovskii & Dieterich, 2017). pulseTD, INSPEcT and DRiLL provide evaluation functions for gene expression and processing rates, except for DTA and pulseR. The encoding of DRiLL is MATLAB, which requires a linux system and has a large operational limit. In terms of model, pulseTD directly applies the pulse model to different stages of RLCD, making the calculation of RLCD continuous. Although the pulse model was also used in INSPEcT and DRiLL, the purpose was to fit the expression pattern. Most importantly, pulseTD has the ability to predict RLCD (Table 1). Overall, pulseTD improves accuracy between biological repetitions and enhances the performance of evaluating low expression data. At the same time, it adds the ability to predict steady state. The assumption that the transcription rate is constant within the 4sU labeling time is abandoned in the solution process, and dynamic rates and gene expression values can be predicted at some future time points.

Performance analysis

The efficacy of pulseTD was evaluated on the simulation dataset. We calculated the Pearson correlation coefficient (PCC) for pulseTD between real and simulated transcriptional dynamic rates, and found that the PCC values for transcription, processing and degradation rates were 0.95, 0.95 and 0.77, respectively (Figs. 2A–2C). The PCCs for the transcription, processing, and degradation rates of INSPEcT were 0.94, 0.87 and 0.42, respectively. However, the rates of RCLD evaluated by INSPEcT were less than the true rates (Figs. 2D–2F). The MSE between real and simulated RCLD rates were calculated using pulseTD, INSPEcT and curve-fitting methods. The results showed that the MSE value of pulseTD was <0.1 (Fig. 3). These results suggested that pulseTD was accurate in assessing the RCLD rates.

To eliminate the bias of simulation data, simulation data produced by INSPEcT was used to evaluate the performance of pulseTD. We used the INSPEcT tool to generate two biological replications dataset, which contained 10 time nodes. The correlation of total RNA expression dataset was 0.98. Then, we used pulseTD to evaluate the RCLD of the biological replications, where the correlations of transcription, processing and degradation rates were 0.97, 0.90 and 0.92 (Figs. 4A–4C). As a comparison, we used INSPEcT to evaluate the RCLD of the simulation data. The correlations of transcription, processing, and degradation rates were 0.98, 0.86 and 0.84 (Figs. 4D–4F). In general, pulseTD had shown good performance in different simulation data.

The robustness of pulseTD was tested using two biologically repeated gene expression datasets from GEO (GSE59717). The PCC values of transcription, processing, and degradation rates between replications were 0.95, and the RCLD rates distribution also showed high consistency (Figs. 5A–5C). This reflected high robustness of pulseTD. For comparison, the PCC values of INSPEcT were calculated using the same data, they were 0.92, 0.93 and 0.84 (Figs. 5D–5F). The PCC for its degradation rates was low, which might be due to more outliers, during the evaluation process. And its density distribution produced a large deviation. pulseTD showed higher stability in biological replicate data. At the same time, in order to judge the impact of expression levels of gene on software performance, we calculated the mean expression levels of the total exons, total introns and labeled exons of each gene during the detection time and ranked the mean. The first 1/3 genes were considered to be low expression, and the last 1/3 genes were considered to be high expression. Both of them were used to evaluate RCLD. The PCC value of degradation rates and processing rates obtained by INSPECT were 0.58 and 0.68, but the PCC value of pulseTD remained stable. This showed that pulseTD had a higher tolerance for poor quality data and the evaluation results were more stable. We also explored the effects of pre-RNA, total RNA and labeled RNA expression levels on model optimization. The correlation of the transcription, processing and degradation rates of each gene was calculated between the two biological replicate data sets. At the same time, we sorted them according to the mean expression levels of the total exons, total introns and labeled exons, respectively. We found that the number of genes with an average correlation greater than 0.5 was 96.21% and overall showed a high correlation (Figs. 6A–6C). This indicated the expression levels of total RNA, pre-RNA and labeled RNA had little effect on the optimization of the model.

Predict RLCD at unknown time nodes

The pulseTD can predict dynamic transcription rates and gene expression because of its pulsed model characteristics. To evaluate prediction effectiveness, we divided the GEO (GSE56977) dataset into two parts. We selected the first 5, 7, 9 and 11 time points from 13 measurement time points as training samples (Assuming they are experimental measurement data) to estimate the rates of RCLD and used the remaining time points (test samples) for prediction. We estimated and predicted the rates of RCLD based on the training samples by using pulseTD and some other methods. Then, the MSE values between the predicted results and the test samples were calculated to evaluate the prediction performance. The results showed that pulseTD had the lower average MES values of RLCD rates than other methods (Figs. 7A–7C). These results showed that pulseTD had good prediction capabilities. We also found that as the training time points increase, the MSE values gradually decreased. When estimating RLCD, we recommend increasing the number of measurement time points in the experiment in order to have more accurate predictions.