Quantifying how post-transcriptional noise and gene copy number variation bias transcriptional parameter inference from mRNA distributions

We have four independent datasets from which we determined mRNA count and nascent RNA distributions. Figure 3a shows an example cell with mature single RNAs in the cytoplasm, and a bright nuclear spot representing the site of nascent transcription. Spots and cell outlines were identified using automated pipelines. Importantly, to obtain an accurate estimation of transcriptional parameters, the experimental input distributions of mRNA count and nascent RNAs require high accuracy. We therefore first determined how technical artifacts in the analysis affects the inference estimates.

Figure 3

Download asset Open asset

First, if the number of mRNA transcripts per cell is high, accurate determination of the number of transcripts may be challenging, as transcripts may overlap. To determine if this occurred in our datasets, we analyzed the distributions of intensities of the cytoplasmic spots, which revealed unimodal distributions where ∼90% of the detected spots fell in the range 0.5× median – 1.5× median (Figure 4a). We therefore concluded that overlapping spots are not a large confounder in our data. In fact, in our experiments, the number of detected mature mRNA transcripts per cell was lower than expected, based on the number of nascent transcripts (compare Figure 3 with Figure 4). This discrepancy between nascent and mature transcripts likely arises because the addition of the PP7 loops to the GAL10 RNA destabilizes the RNA, resulting in faster mRNA turnover compared to most endogenous RNAs (Miller et al., 2011; Wang et al., 2002; Holstege et al., 1998; Geisberg et al., 2014). Previously, both shorter and longer mRNA half-lives from the addition of stem loops have been observed, which may be caused because changes in the 5’ UTR length or sequence affect its recognition by the mRNA degradation machinery (Heinrich et al., 2017; Tutucci et al., 2018; Garcia and Parker, 2015). In our case, we note that such high turnover should aid transcriptional parameter estimates, as it closely reflects transcriptional activity.

Figure 4

Download asset Open asset

A second possible source of error is cell segmentation. To test how cell segmentation errors contribute to the mature mRNA distribution and the transcriptional bursting estimates, we compared two independent segmentation tools, where segmentation 1 often resulted in missed spots (Figure 3b), resulting in an underestimation of the mean mRNA count and of the variance (compare Figure 3b and c). We inferred the transcriptional parameters using the algorithm described in Methods Section Steps of the algorithm to estimate parameters from mature mRNA data. In the absence of an experimental measurement of the degradation rate, we could only estimate the three transcriptional parameters normalised by $d$. The best fits of dataset 1 are shown in (Figure 3b and c) and the transcriptional parameters (for all four datasets) are summarized in (Figure 3e). Note that the estimated parameters for all four datasets, using both segmentations, are shown in Appendix 3—table 1 and the associated best fit distributions in Appendix 3—figure 1a. Notably the segmentation algorithms led to similar estimates for the burst frequency but considerably different estimates for the rest of the parameters. In particular segmentation 1 suggested that burst expression is infrequent (≈20% of the time) whereas segmentation 2 was consistent with burst expression occurring half of the time. Given that accurate cell segmentation remains challenging, this analysis illustrates that parameter estimation from mature mRNA counts may be affected by technical errors. For the remainder of the mature mRNA analysis, we have used only segmentation 2 data.

Lastly, it may be challenging to distinguish the nascent transcription site from a mature RNA, especially if few nascent RNAs are being produced. Either one can decide to include all cellular spots in the total mRNA count, including the transcription site, with the result that the number of mature transcripts is overestimated with one RNA for cells which show a transcription site. Or conversely, one can decide to exclude the transcription site by subtracting one spot from each cell, with the result that the number of mature mRNAs may be underestimated by one RNA for cells that are transcriptionally silent. To understand how this choice affects the accuracy of parameter inference, we compared both options in (Figure 3c, d and e), where seg2 included all spots, and seg2-TS excluded transcription sites (by subtracting 1 from each cell). The estimated parameters for all four datasets are shown in Appendix 3—table 1 and the associated best fit distributions in Appendix 3—figure 1a. Although the mean was lower when transcription sites were excluded, all the parameters except the burst frequency ${\sigma }_{\mathrm{on}}$ were within the error, indicating that the choice of whether or not to include the transcription site in the mature mRNA count had a small influence on parameter estimation. For the remainder of the analysis, we included all spots, and counted the transcription site as one RNA.

The above analysis was performed using the merged data from all cells, irrespective of their position in the cell cycle. The inferred parameters of all four datasets are shown in Figure 3g (grey). To understand the effect of the cell cycle on these parameter estimates, we compared this inference with cell-cycle-specific data. We used the integrated nuclear DAPI intensity as a measure for DNA content to classify cells into G1 or G2 cells (Figure 3f (left)) to obtain separate mature mRNA distributions for G1 and G2 cells.

To infer the transcriptional parameters from mature mRNA data of cells in G1, the inference protocol remained the same. However for cells in the G2 stage, this protocol needed to be altered since G2 cells have two gene copies, whereas the solution of the telegraph model assumes one gene copy. Assuming the transcriptional activities of the two gene copies are independent, the distribution of the total molecule number is the convolution of the molecule number (obtained from the telegraph model) with itself for mature mRNA data. This convolved distribution was used in steps (ii) and (iii) of the inference algorithm in Methods Section Steps of the algorithm to estimate parameters from mature mRNA data. A difference between our method of estimating parameters in G2 from that in the literature (Skinner et al., 2016) is that we do not assume that the burst frequency is the only parameter that changes upon replication, and we estimated all transcription parameters simultaneously.

Note that the independence of gene copy transcription has been verified for genes in some eukaryotic cells (Skinner et al., 2016) where the two copies can be easily resolved. For yeast data, as we are analyzing in this paper, it is generally not possible to resolve the two copies of the allele in G2 because they are within the diffraction limit. However, in the absence of experimental evidence, the independence assumption is the simplest reasonable assumption that we could make (see later for a relaxation of this assumption).

For both G1 and G2 cells, we performed inference for cell-cycle specific mature mRNA data, the results of which are shown in Figure 3f (centre and right) and Figure 3g – see Appendix 3—table 2 for the confidence intervals of the estimates calculated using profile likelihood. As expected, the mean number of mRNAs in G2 cells was larger than that in G1 cells. For both merged and cell-cycle specific data, the parameters ordered by increasing variability of the estimates from independent samples (the standard deviation divided by the mean) were: $\rho$, $f_{\mathrm{ON}}$, ${\sigma }_{\mathrm{ON}}$, burst size and ${\sigma }_{\mathrm{OFF}}$, and the same order was predicted by the relative error (from ground truth values) from our synthetic experiments (compare with $f_{\mathrm{ON}}=0.50$ and $f_{\mathrm{ON}}=0.80$ in the middle and right panels of Figure 2d) and by sample variability (Appendix 1). In Appendix 3 and Appendix 3—table 3 we show that the relaxation of the assumption of independence between the allele copies in G2 (by instead assuming perfect state correlation of the two alleles) had practically no influence on the inference of the two best estimated parameters ($\rho$, $f_{\mathrm{ON}}$).

A comparison of the two types of data predicted different behaviour (Figure 3g bottom): merged data indicated behaviour consistent with the gene being ON half of the time and small burst sizes, while cell-cycle-specific data implied the gene is ON ≈80% of the time with large burst sizes. We note that the burst sizes have considerable sample variability, exemplifying burst size estimates of transcriptional parameters from mature mRNA distributions have to be treated with caution. Nevertheless, in line with this high fraction ON and large burst size, which start to approach constitutive expression, the variation introduced by the transcription kinetics is relatively modest with Fano factors not far from one: $2.43\pm 0.21$ for merged data and $1.75\pm 0.45$ for cell-cycle data (the slightly higher value for merged data likely was due to heterogeneity stemming from varying gene copy numbers per cell).

Comparing the mean rates between the G1 and G2 phases, we found that ${\sigma }_{\mathrm{off}}$, ${\sigma }_{\mathrm{on}}$, $\rho$ decreased while $f_{\mathrm{ON}}$ and the burst size increased upon replication. However, taking into account the variability in estimates across the four datasets, the only two parameters which were well-separated between the two phases were ${\sigma }_{\mathrm{on}}$ and $\rho$. These two parameters decreased by 65% and 21%, respectively, which suggests that upon replication, there are mechanisms at play which reduce the expression of each copy to partially compensate for the doubling of the gene copy number (gene dosage compensation) (Skinner et al., 2016).

In conclusion, what is particularly surprising in our analysis is the differences in the inference results using merged and cell-cycle specific data: the former suggests the gene spends only half of its time in the ON state while the latter implies the gene is mostly in its ON state.

To determine the number of nascent transcripts at the transcription site, we selected the brightest spot from each nucleus and normalized its intensity to the median intensity of the cytoplasmic spots. As the distribution of intensities of the cytoplasmic mRNAs followed a narrow unimodal distribution, its median likely represents the intensity of a single RNA (orange distribution in the central panel of Figure 4a). The inference of transcriptional parameters using the merged data was done using the algorithm described in Methods Section Steps of the algorithm to estimate parameters from nascent mRNA data.

Similar to above, to account for two gene copies in G2 cells, we assumed that the transcriptional activities of the two gene copies are independent. The distribution of the total fluorescent signal from both gene copies was the convolution of the signal distribution (obtained from the extended delay telegraph model, i.e. Equation (1)) with itself. This convolved distribution was then used in steps (ii) and (iii) of the inference algorithm.

The inference of transcriptional parameters from nascent RNA data was done using a fixed elongation time, which was measured previously at a related galactose-responsive gene (GAL3) at $65\mathrm{bp}/\mathrm{s}$ (Donovan et al., 2019). Since the total transcript length is $3062\mathrm{bp}$ (see Figure 1c), the elongation time ($\tau$ in our model) is $\approx 47.1s\approx 0.785\text{min}$. The fixed elongation rate enabled us to infer the absolute values of the three transcriptional parameters ${\sigma }_{\text{off}},{\sigma }_{\text{on}}$ and $\rho$.

Best fits of the extended delay telegraph model to the distribution of signal intensity of nascent mRNAs at the transcription site are shown in Figure 4a and b for dataset 1; for the other datasets see Appendix 4—figure 1. The corresponding estimates of the transcriptional parameters are shown in Appendix 4—table 1 and also illustrated by bar charts in Figure 4c. The confidence intervals of the transcriptional parameters (computed using the profile likelihood method) are shown in Appendix 4—table 2.

Comparing this estimation with that from mature mRNA, we observed that in both cases $f_{\mathrm{ON}}\approx 0.5$ for merged data and in the range $0.7-0.8$ for cell-cycle-specific data. Also in both cases, the Fano factors of merged data were larger than those of cell-cycle-specific data. Hence, we are confident that not accounting for the cell cycle phase leads to an over-estimation of the time spent in the OFF state and of the Fano factor. In addition, comparing the burst sizes in Figure 3g and Appendix 4—table 1, we found that not taking into account post-transcriptional noise (by using mature mRNA data) led to an lower estimation of the burst size (2.6-fold, 2.6-fold, and 1.1-fold lower for inference from merged, G1 and G2 data, respectively). We note that it would be useful to directly compare the absolute estimates of the other transcriptional parameters from mature and nascent mRNA data. However, this was not possible because the telegraph model only estimates the switching rates and the initiation rate scaled by the degradation rate, and the latter is unknown. On the other hand, the estimates from nascent data were rates multiplied by the average elongation time, which is known and hence the absolute rates can be estimated from nascent mRNA data only. The only quantities that could be directly compared were the burst size and the fraction of ON time, since these are both non-dimensional.

Comparing the variability of the parameter estimates, we found that $\rho$ and $f_{\mathrm{ON}}$ were the parameters with the smallest variability across samples for the nascent data, as for inference from mature data. However, the inferred parameter variability across samples was on average about 2.5-fold lower for nascent data compared to mature mRNA data (this was obtained by computing the standard deviation divided by the mean for each parameter and then averaging over all parameters and over merged, G1 and G2 data). Likely this is because nascent data does not suffer from post-transcriptional noise. Indeed, synthetic experiments suggested that the errors in parameter inference using nascent data are often less than those in mature data when $f_{\mathrm{ON}}\approx 0.80$ (Figure 2d). In summary, we have more confidence in the parameter estimates from nascent data, in particular those from cell-cycle separated data.

To further investigate the hypothesis that estimates from cell-cycle-specific data are more accurate than merged data, we compared the estimates from merged and cell-cycle-specific data to previous live-cell transcription measurements of the same gene (Donovan et al., 2019). Because live-cell traces and simulated traces with the estimated transcriptional parameters are difficult to compare directly, we instead compared their normalized autocorrelation functions (ACFs). Specifically, the live-cell traces displayed cell-to-cell variation in overall fluorescent intensities arising from differences in the PP7 coat protein expression level, precluding a direct comparison of the live-cell intensities with the smFISH distributions. The normalized ACFs are normalized per trace and thus can be used to directly compare the kinetics. For this, we fed the parameter estimates to the SSA to generate synthetic live-cell data and then calculated the corresponding ACF (Appendix 5). We found that the estimates from cell-cycle-specific data produced ACFs that match the live-cell data closer than that from the merged data (Figure 4d). This was also clear from the sum of squared residuals which for each dataset was smaller for the ACF computed using the cell-cycle-specific estimates rather than those from merged data (Figure 4e).

Using nascent data, we also reinvestigated the hypothesis that the gene exhibits dosage compensation. Comparing the mean rates between the G1 and G2 phases, we found that ${\sigma }_{\mathrm{off}}$, ${\sigma }_{\mathrm{on}}$, $\rho$, $f_{\mathrm{ON}}$ decreased while the burst size increased upon replication. However, taking into account the variability in estimates across the four datasets, the only two parameters which were cleanly separated between the two phases were ${\sigma }_{\mathrm{on}}$ and $f_{\mathrm{ON}}$. These two decreased by 41% and 5%, respectively. These results had some similarity to those deduced from cell-cycle separated mature mRNA data (the decrease of ${\sigma }_{\mathrm{on}}$) but they also displayed differences. Namely, from mature mRNA data it was predicted that $\rho$ decreased upon replication while from nascent data we predicted that $\rho$ did not change and it was rather $f_{\mathrm{ON}}$ that decreased by a small degree. The decrease of the burst frequency ${\sigma }_{\text{on}}$ after replication has also been reported for some genes in mammalian cells (Skinner et al., 2016; Padovan-Merhar et al., 2015), indicating that this could be a general mechanism for gene dosage compensation. Our results are consistent with a population-based ChIP-seq study (Voichek et al., 2016) that showed DNA dosage compensation after replication in budding yeast. We note that our single-cell analysis only revealed partial dosage compensation, where the mean signal intensity of nascent mRNAs in G2 is not the same as in G1, but 1.7-fold higher in G2 than in G1 (Figure 4c).

Although inference on cell cycle separated data outperformed inference on merged data, we noticed that the corresponding best fit distributions did not match well to the experimental signal distributions in the lower bins (Figure 4b and Appendix 4—figure 1). In all cases, the experimental distributions showed high intensities in bins 1, 2, and 3, which was likely an artifact of the experimental data acquisition system. Since we defined the transcription site as the brightest spot, this implies that in the absence of a transcription site, a mature transcript can be misclassified as a nascent transcript. We therefore investigated two methods to correct for this, the ‘rejection’ method and the ‘fusion’ method.

The rejection method removed all data associated with the first $k$ bins of the experimentally obtained histogram of fluorescent intensities (Figure 5a shows the fits for dataset 1; for the other datasets see Appendix 4—figure 2). We found that the parameter estimates varied strongly when the number of bins from which data was rejected ($k$) was changed (Figure 5b; see also Appendix 4—table 3). Although the distributions fit well to the experimental histograms (Appendix 4—figure 1), comparison with the live-cell normalized ACF indicated that the estimates actually became worse than non-curated estimates, with a higher sum of squared residuals (Figure 5c and d). The rejection method therefore does not produce reliable estimates.

Figure 5

Download asset Open asset

Next, we considered another data curation method which we call the fusion method. This works by setting to zero all fluorescent intensities in a cell population which were below a certain threshold. In other words, we fused or combined the first $k$ bins of the experimentally obtained histogram of fluorescent intensities, thereby taking into account that the true intensity of bin 0 was artificially distributed over some of the first bins.

Figure 6 and Appendix 4—table 4 show that the fusion method led to estimates that varied little with $k$ which enhanced our degree of confidence in them (note that $k=1$ is the same as the uncurated data). The peak at the zero bin for both G1 and G2 was better captured using the fusion method than using non-curated data (compare Figure 4b and Appendix 4—figure 1, with Figure 6b). Comparison to the autocorrelation function of the live-cell data shows that correction with the fusion method also led to improved transcriptional estimates, as indicated by a reduction in the sum of the squared residuals for all four data sets (Figure 6c).

Figure 6

Download asset Open asset

Overall, we conclude that for inferring parameters from the smFISH data, the optimal method is to use nascent cell-cycle-specific data, corrected by the fusion method. The optimally inferred parameters for the four data sets in our study are those given in Appendix 4—figure 2d. The profile likelihood estimates of the 95% confidence intervals of these parameters are shown in Appendix 4—table 5. Note that in line with our synthetic data study in Figure 2, the parameters suffering from the least sample variability were $f_{\mathrm{ON}}$ and $\rho$. The rest of the parameters (${\sigma }_{\mathrm{off}},{\sigma }_{\mathrm{on}}$ and burst size) suffered more sample variability because the fraction of ON time was high; however since their standard deviation divided by the mean (computed over the four datasets) was not high (in the range of 10-20%), they still can be regarded as useful estimates. Note also that the previous prediction that gene dosage compensation involves regulation of the burst frequency did not change upon correction of the nascent data using the fusion method. All these results were deduced assuming that the two copies in G2 are independent from each other. Inferring rates under the opposite assumption of perfectly synchronized copies (Appendix 4—table 6) gave very similar estimates for $\rho$ and $f_{\mathrm{ON}}$ (to be expected since according to the synthetic data study, these two are the most robustly estimated parameters for genes spending most of their time in the active state) but different estimates for the rest of the parameters. While such perfect synchronization of alleles is unlikely, some degree of synchronization is plausible and further improvement of the transcriptional parameters in the G2 phase will require its precise experimental quantification.

The steady-state solution of the telegraph model of gene expression (Peccoud and Ycart, 1995) gives mature mRNA distributions. The reaction steps in this model are illustrated in Figure 1a. Next we describe the generation of synthetic mature mRNA data and the algorithm used to infer parameters from this data.

We generate parameter sets on an equidistant mesh grid laid over the space:

where the units are inverse minute. Furthermore we apply a constraint on the effective transcription rate

In each of the three dimensions of the parameter space, we take 10 points that are equidistant, leading to a total of 1000 parameter sets which reduce to 789 after the effective transcription rate constraint is enforced.

We additionally fix the degradation rate $d=1$ min^-1. Note that we choose not to vary the degradation rate (as we did for the other three parameters) since it is not possible to infer all four rates simultaneously – this is because the steady-state solution of the telegraph model is a function of the non-dimensional parameter ratios $\rho /d,{\sigma }_{\text{off}}/d$ and ${\sigma }_{\text{on}}/d$ (Raj et al., 2006).

Once a set of parameters is chosen, we use the stochastic simulation algorithm (SSA Gillespie, 2007) to simulate the telegraph model reactions in Figure 1a and generate 10⁴ samples of synthetic data. Note that each sample mimicks a single cell measurement of mature mRNA.

The inference procedure consists of the following steps: (i) select a set of random transcriptional parameters; (ii) use the solution of the telegraph model to calculate the probability of observing the number of mature mRNA measured for each cell; (iii) evaluate the likelihood function for the observed data; (iv) iterate the procedure until the negative log-likelihood is minimized; (v) the set of parameters that accomplishes the latter provides the best point-estimate of the parameters of the telegraph model that describes the measured mature mRNA fluctuations.

For step (i), we restrict the search for optimal parameters in the following region of parameter space

The degradation rate is fixed to $d=1$ min^-1.

Step (ii) can be obtained either by computing the distribution from the analytical solution (Peccoud and Ycart, 1995 or by using the finite state projection (FSP) method Munsky and Khammash, 2006). Here, for the sake of computational efficiency, we use the FSP method to compute the probability distribution of mature mRNA numbers.

For step (iii) we calculate the likelihood of observing the data given a chosen parameter set $\theta$

where $P(n_j;\theta )$ is the probability distribution of mature mRNA numbers obtained from step (ii) given a parameter set $\theta$, n_j is the total number of mature mRNA from cell $j$ and $N_{\text{cell}}$ is the total number of cells.

Steps (i) and (iv) involve an optimization problem. Specifically we use a gradient-free optimization algorithm, namely adaptive differential evolution optimizer (ADE optimizer) using BlackBoxOptim.jl (https://github.com/robertfeldt/BlackBoxOptim.jl; Feldt and Stukalov, 2022) within the Julia programming language to find the optimal parameters

The minimization of the negative log-likelihood is equivalent to maximizing the likelihood. Note the optimization algorithm is terminated when the number of iterations is larger than 10⁴; this number is chosen because we have found that invariably after this number of iterations, the likelihood has converged to some maximal value. Note that the inference algorithm is particularly low cost computationally, with the optimal parameter values estimated in at most a few minutes.

Once the best parameter set ${\theta }^{*}$ is found, we calculate the mean relative error (MRE) which is defined as

where ${\theta }_i^{*}$ and ${\theta }_{\text{true},i}$ represent the $i$-th estimated and true parameters respectively, and $M$ denotes the number of the estimated parameters. Thus, the mean relative error reflects the deviation of the estimated parameters from the true parameters.

The steady-state solution of the delay telegraph model (Xu et al., 2016) gives the distribution of the number of bound Pol II. In Appendix 6, we present an alternative approach to derive the steady-state solution. The reaction steps are illustrated in Figure 1a.

The position of a Pol II molecule on the gene determines the fluorescence intensity of the mRNA attached to it. In particular for fluorescence data acquired from smFISH PP7-GAL10, the fluorescence intensity of a single mRNA on the DNA locus looks like a trapezoidal pulse (see Figure 1b for an illustration). This presents a problem because although we can predict the distribution of the number of bound Pol II using the delay telegraph model, we do not have any specific information on their spatial distribution along the gene. However, since the delay telegraph model implicitly assumes that a Pol II molecule has fixed velocity and that Pol II molecules do not interact with each other (via volume exclusion), it is reasonable to assume that in steady-state, the bound Pol II molecules are uniformly distributed along the gene. This hypothesis is confirmed by stochastic simulations of the delay telegraph model where the position of a Pol II molecule is calculated as the product of the constant Pol II velocity and the time since its production.

By the uniform distribution assumption and the measured trapezoidal fluorescence intensity profile, it follows that the signal intensity of each bound Pol II has the density function $g$ defined by

where $L_1=862\text{bp}(\text{base pairs}),L_2=2200\text{bp},L=L_1+L_2$ as defined in Figure 1b. The indicator function ${\chi }_{[0,1]}(s)=1$ if and only if $s\in [0,1]$ and ${\delta }_1(s)$ is the Dirac function at 1. The probability of the signal $s$ being between 0 and 1 is due to the first part of the trapezoid function and hence is multiplied by $L_1/L$ which is the probability of being in this region if Pol II is uniformly distributed. Similarly, the probability of $s$ being 1 is due to the L₂ part of the trapezoid and hence the probability is $L_2/L$ by the uniform distribution assumption. Note that the signal $s$ from each Pol II is at most 1 because in practice, the signal intensity from the transcription site is normalized by the median intensity of single cytoplasmic mRNAs (Zenklusen et al., 2008).

The total signal is the sum of the signals from each bound Pol II. Hence, the density function of the sum is given by the convolution of the signal densities from each bound Pol II. Defining $p(s|k)$ as the density function of the signal given there are $k$ bound Pol II molecules, we have that $p(s|k)$ is the $k$–th convolution power of $g$, that is

where ${\displaystyle {\delta }_0(s)}$ is the Dirac function at. Finally we can write the total fluorescent signal density function as

where $P(k;\theta )$ is the steady-state solution of the delay telegraph model giving the probability of observing $k$ bound Pol II molecules for the parameter set $\theta$. Hence Equation (8) represents the extension of the delay telegraph model to predict the smFISH fluorescent signal of the transcription site.

Comparison to the algorithm in Xu et al., 2016. Both algorithms take into account the fact that the signal intensity depends on the position of Pol II on the gene, albeit this is done in different ways. In Xu et al., 2016 a master equation is written for the joint distribution of gene state and the number of nascent mRNA. In this case the number of nascent RNAs can have non-integer values since it represents the experimentally measured signal from the (incomplete) nascent RNA. Solution of this master equation proceeds by (a) a discretization of the continuous nascent mRNA signal into bins which are much smaller than one; (b) solution using finite state projection (FSP). This approach can lead to a large state space which incurs a large computational cost. In contrast, in our method, we use FSP to solve for the delay telegraph model, i.e. the distribution of the discrete number of bound Pol II from which we construct (using convolution) the approximate distribution of the continuous nascent mRNA signal by assuming the Pol II is uniformly distributed on the gene. Since the state space of bound Pol II is typically not large, our method will typically be more computationally efficient than the one described in Xu et al., 2016.

We generated synthetic smFISH signal data by using the SSA, modified to include delay to simulate the delay telegraph model (Fu et al., 2022). Specifically, we use Algorithm 2 described in Barrio et al., 2006. One run of the algorithm simulates the fluctuating number of bound Pol II molecules in a single cell.

The total fluorescence intensity (mimicking smFISH) is obtained as follows. When a particular bound Pol II is produced by a firing of the transcription reaction $G\to G+N$, we record this production time; since the elongation rate is assumed to be constant, given the production time we can calculate the position of the Pol II molecule on the gene at any later time and hence using Figure 1b we can deduce the fluorescent signal due to this Pol II molecule.

Specifically we normalize each transcribing Pol II’s position to $[0,1]$ and map the position to its normalized signal by

where $x$ is the normalized position on the gene. Thus at a given time, the total fluorescent signal from the $n$-th cell (the $n$-th realization of the SSA) equals

where $J_n$ is the number of bound Pol II molecules in the $n$-th cell, and $\{x_j\}$ with $j=1,\mathrm{\dots },J_n$ is the vector of all Pol II positions on the gene. The total signal from each cell is a real number but it is discretized into an integer.

The kinetic parameters are chosen from the same region of parameter space as in (2), on the same equidistant mesh grid and with the same constraint on the effective transcription rate. Unlike the mature mRNA case, here there is no degradation rate; instead we have the elongation time, which we fix to $\tau =0.5\text{(min)}$. Note that fixing this time is necessary since it is not possible to infer the three transcriptional parameters rates and the elongation time simultaneously because the steady-state solution of the delay telegraph model is a function of the non-dimensional parameter ratios $\rho \tau ,{\sigma }_{\text{off}}\tau$ and ${\sigma }_{\text{on}}\tau$.

Once a set of parameters is chosen, we use the modified SSA (as described above) to simulate the signal intensity in each of 10⁴ cells.

The inference procedure is essentially the same as steps (i)-(v) described in mature mRNA inference except for the following points.

In step (ii), the probability of observing a total signal of intensity $i$ from a single cell is obtained by integrating $p(s;\theta )$ in Equation (8) on an interval $[i-1,i]$ for $i\in \mathbb{N}$ which, in our numerical scheme, means

Note that the integration over the interval of length 1 is to match the discretization of the synthetic data and $\theta \in \mathrm{\Theta }$. Intuitively, one can always choose a positive integer $K$ such that $P(k)=0$ for any $k\ge K$. The computation of the solution of the delay telegraph model $P(k)$ can be done either using the analytical solution (evaluated using high precision) or using the finite state projection algorithm (FSP) Munsky and Khammash, 2006. In Appendix 6—figure 1 and Appendix 6—table 1, we show that the two methods yield comparable accuracy and CPU time.

For step (iii) we calculate the likelihood of observing the data given a chosen parameter set $\theta$

where q_j is the discretized total signal intensity from cell $j$ and $N_{\text{cell}}$ is the total number of cells. In the optimization, we aim to find

The whole procedure (for both mature and nascent mRNA inference) is summarized by a flow-chart in Figure 1c.

In the main text, Figure 2, we showed how the addition of 5% external noise to synthetic mature mRNA data degrades the inference accuracy. In Appendix 1—figure 1 we show how the addition of a larger amount of external noise (10%) causes an even larger loss of accuracy. In particular for 91% of the parameters, the inference accuracy is higher when using nascent mRNA data (Appendix 1—figure 1a) and the median relative errors become very high for most parameters, especially for $\rho$ and ${\sigma }_{\mathrm{off}}$ in the limit of large $f_{\text{ON}}$ (Appendix 1—figure 1b).

Appendix 1—figure 1

Download asset Open asset

Appendix 1—figure 2

Download asset Open asset

In the main text, Figure 2, we showed how the accuracy of parameter estimation is not uniform across parameter space. Here we investigate if there is any relationship between this accuracy and how well is a distribution of mature mRNA numbers / signal intensity fit by the inference algorithm. For 12 parameter sets (6 for the telegraph model – Appendix 1—figure 2a) and (6 for the delay telegraph model – Appendix 1—figure 2b), we evaluate the fits to the distribution of synthetic data in Appendix 1—figure 2c–d. The results show that independent of the accuracy of the parameters estimated by the inference algorithm, the fits of the delay telegraph and telegraph model distributions to the distributions generated from synthetic data are generally excellent.

To obtain a better understanding of the variability of the inference procedure, for each of the six parameter sets in Appendix 1—figure 2b and i.e. for the inference using the delay telegraph model, we generated 10 independent sets of synthetic data and used maximum likelihood to infer the parameters for each dataset. The mean and standard deviation of the parameters (computed over all 10 datasets) are shown in Appendix 1—table 1. The means are close to the true parameter values (in Appendix 1—figure 2b); this shows that the inference procedure is working correctly and the deviations from ground truth values are mostly due to noise in the synthetic datasets. We can define the sample variability of a parameter as the standard deviation divided by the mean. Ordering parameters by this quantity, we find that for all six parameter sets the error is largest for ${\sigma }_{\mathrm{off}}$. For small $f_{\mathrm{ON}}$, the parameters ordered by increasing error are: ${\sigma }_{\mathrm{on}}$, burst size, $\rho$, $f_{\mathrm{ON}}$ and ${\sigma }_{\mathrm{off}}$. While for large $f_{\mathrm{ON}}$, the order is: $\rho$, $f_{\mathrm{ON}}$, ${\sigma }_{\mathrm{on}}$, burst size and ${\sigma }_{\mathrm{off}}$. Note that this order is the same as we determined using the relative errors (from the ground truth) which is shown in Figure 2 of the main text.

We perform a profile likelihood study (Kreutz et al., 2013) on the 12 parameter sets of synthetic mature/nascent mRNA data described in Appendix 1—figure 2a-b. We obtain the 95% confidence interval for each parameter. The results are shown in Appendix 1—table 2. We also compare the relative errors for each parameter (computed as (estimated value - ground truth)/ground truth using the data in Appendix 1—figure 2a and b) with the profile likelihood error (computed as (upper bound - lower bound)/optimal estimate using the bounds in Appendix 1—table 2 and the optimal values in Appendix 1—figure 2a and b). The results are shown in Appendix 1—table 3. Note that in most cases, the parameters ordered by relative error are in agreement with the parameters ordered by profile likelihood error.

To assess the reliability of the inference results due to errors in spot counting, we redid the inference with synthetic mRNA data (and using the telegraph model) perturbed randomly by minus 1/plus 1/unchanged with probability $1/3$. The results are shown in Appendix 1—table 4.

We found that when the fraction of ON time was very small (Set 1), there is a considerable effect of the perturbations on the values of the inferred parameters – this is because in this case, the mean number of mRNA is very small and hence a perturbation of one molecule is very significant. However as expected, the inference results are quite robust when the fraction of ON time is not too small (Sets 2–6).

Given the synthetic nascent mRNA signal data ${\{S_i\}}_{i=1}^N$ (where $N$ represents the number of samples) generated by delay telegraph model, for each $S_i$, we perform a random perturbation $\mathrm{\Phi }$ under the following conditions

where $\mathrm{\Phi }$ is a stochastic perturbation satisfying the following constraints

${\displaystyle \mathrm{\Phi }(S_i)}$ is a random variable sampled from the distribution Log-normal ${\displaystyle (\alpha ,\beta )}$ whose mean equals ${\displaystyle S_i}$, and the standard deviation equals ${\displaystyle 0.1\ast S_i}$.

This means the random perturbation keeps the mean value of the signal $S_i$ unchanged but adds noise with a coefficient of variation equal to 0.1.

In Appendix 1—table 5 we compare the results of inference using synthetic nascent mRNA data with the aforementioned stochastic perturbation and without. As for mature mRNA data, we find that the perturbation has only a significant impact when the fraction of ON time is very small.

In main text Section Inference from nascent mRNA data: cell cycle effects, experimental artifacts and comparison with mature mRNA inference, we used independent segmentation tools to study segmentation artifacts. The inference results under segmentation 1, segmentation 2 and segmentation 2 without counting the transcriptional site are summarised in Appendix 3—table 1.

Appendix 3—figure 1

Download asset Open asset

We obtain 95% confidence intervals for the parameters estimated from mature mRNA data under segmentation 2 (shown in Figure 3f of the main text). The confidence intervals are shown for ${\sigma }_{\text{off}},{\sigma }_{\text{on}},\rho$ in Appendix 3—table 2.

In the main text, we assumed that the transcription of two allele copies in the G2 phase are independent. Here we consider the opposite scenario where the two allele copies in G2 phase are perfectly synchronized with each other, i.e. when one copy switches from on to off, the other copy also switches from on to off. Note that the time at which mRNA is transcribed from each copy is however not the same. Using this modified simulation algorithm, we re-perform the inference using the experimental data under segmentation 2. The results are shown in Appendix 3—table 3. Comparing these values with those estimated for phase G2 under the assumption of allele independence (main text Figure 3f right panel), we find that $\rho$ and $f_{\mathrm{ON}}$ are very similar but the other parameters vary considerably.

We show the inference results of the experimental non-curated nascent mRNA data (with and without taking into consideration the cell-cycle) for the four datasets (Appendix 4—table 1) and their 95% confidence interval using the profile likelihood estimate (Appendix 4—table 2). The best fit distributions for merged and cell-cycle-specific data are shown in Appendix 4—figure 1.

Appendix 4—figure 1

Download asset Open asset

We show the inference results using nascent mRNA data curated with the rejection method in Appendix 4—table 3 and Appendix 4—figure 2, and with the fusion method in Appendix 4—table 4. We also used the profile likelihood method to obtain the 95% confidence intervals for the fusion method with $k=4$; the results are shown in Appendix 4—table 5. We also reinferred the parameters in the G2 phase for fusion corrected data by assuming that the allele states are perfectly synchronised; these are shown in Appendix 4—table 6.

Appendix 4—figure 2

Download asset Open asset

The delay telegraph model includes four reactions

where $G$ and $G^{\star }$ stand for the active (ON) and inactive (OFF) gene state, respectively, and ${\sigma }_{\text{off}}$ and ${\sigma }_{\text{on}}$ are the activation and inactivation rates, respectively. Once a nascent mRNA is produced, it will be removed from the system after a deterministic time $\tau$. We aim to find the quantitative description of the kinetics of nascent mRNA ($N$). We proceed by deriving the delay chemical master equation (CME) of Equation (15). The same delay models have also been studied by other authors .

Let $P(0,n,t)$ and $P(1,n,t)$ be the probability of observing $n$ nascent mRNAs while the gene state is inactive and active at time $t$, respectively. Consequently, we can state

Note that by instant reactions, we mean all reactions except the delayed reaction that removes nascent mRNA. Since Part A includes all Markovian reactions, it follows from standard arguments Gillespie, 2007 that

The contribution of Part B requires a careful consideration of the history of the process: (i) $(1,n^{\prime },t-\tau )$ leads to $(1,n^{\prime }+1,t-\tau +\mathrm{\Delta }t)$ when the gene was ON; (ii) $(1,n^{\prime }+1,t-\tau +\mathrm{\Delta }t)$ leads to $(i,n+1,t)$ where the gene state $i$ can be either OFF or ON; (iii) $(i,n+1,t)$ leads to $(i,n,t)$. Indeed, since the removal of nascent mRNA is a delayed reaction, (iii) is due to a production reaction in (i) a time interval $\tau$ earlier. The probability of (i) occurring is $\rho \mathrm{\Delta }tP(1,n^{\prime },t-\tau )$. The probability of (iii) occurring is 1 since every production event is followed by a removal event a time $\tau$ later. The probability of (ii) occurring is:

As for Equation (18), the two ‘+1’s means that the new nascent mRNA was produced during the time interval $(t-\tau ,t-\tau +\mathrm{\Delta }t)$, and it did not participate in any other reactions, thereby not influencing the probability at time $t+\mathrm{\Delta }t$. In short,

Furthermore, all $n^{\prime }$ nascent mRNAs must leave the system prior to time $t$ as they were all born before time $t-\tau$. This implies that

Hence the probability that event B occurs is given by the product of the probability of events (i), (ii) and (iii) and summing over all values of $n^{\prime }$

where we used ${\sum }_{n^{\prime }}P(1,n^{\prime },t-\tau )=P(1,t-\tau )$.

Finally, the contribution of Part C is obtained from simple probability arguments

where we used the same argument as in Equation (21). Using Equations 16; 17; 21; 22 and taking the limit of small $\mathrm{\Delta }t$, we finally obtain the set of delay CMEs

where $E^{k}P(i,n,t)=P(i,n+k,t)$ is the step operator. Summing over all possible $n$ for the two equations in Equation (23), we obtain

initiated with the inactive state, in which $P(0,t)$ and $P(1,t)$ are the probabilities of finding a cell at the inactivated and activated state at time $t$, and their solutions are