Energy-based scheme for reconstruction of piecewise constant signals observed in the movement of molecular machines

Analyzing the physical and chemical properties of single DNA-based molecular machines such as polymerases and helicases requires to track stepping motion on the length scale of base pairs. Although high-resolution instruments have been developed that are capable of reaching that limit, individual steps are oftentimes hidden by experimental noise which complicates data processing. Here we present an effective two-step algorithm which detects steps in a high-bandwidth signal by minimizing an energy-based model (energy-based step finder, EBS). First, an efficient convex denoising scheme is applied which allows compression to tuples of amplitudes and plateau lengths. Second, a combinatorial clustering algorithm formulated on a graph is used to assign steps to the tuple data while accounting for prior information. Performance of the algorithm was tested on Poissonian stepping data simulated based on published kinetics data of RNA polymerase II (pol II). Comparison to existing step-finding methods shows that EBS is superior in speed while providing competitive step-detection results, especially in challenging situations. Moreover, the capability to detect backtracked intervals in experimental data of pol II as well as to detect stepping behavior of the Phi29 DNA packaging motor is demonstrated.

Figure 1

Flowchart of the two-stage process for finding steps in noisy stepping data. The input data is first denoised via solving the convex TVDN problem. This process requires no intervention, as the regularization parameter $\lambda$ is determined automatically. This results in a lower-dimensional, compressed representation of tuples (amplitude, length). This discrete data is then handed to a graph cut algorithm which solves a combinatorial clustering problem on a graph. The graph cut allows further customization by the use of regularization parameters ${\rho }_i$ and a predefined level set.

Figure 2

Breakdown of the TVDN model dependent on $\lambda /{\lambda }_{\max }$. The number of produced steps and plateaus vs different $\lambda$. For small $\lambda /{\lambda }_{\max }$ solving the TVDN problem reproduces the input signal and the number of steps equals the number of data points. For big $\lambda /{\lambda }_{\max }$ the number of steps is significantly lower. The point ${\lambda }_h/{\lambda }_{\max }$ (magenta) before the number of steps increases suddenly marks the value of the TVDN regularization parameter that we choose in our heuristic. Plotted is a constant signal with added Gaussian white noise (red), a signal with exponentially distributed dwell times and Gaussian white noise (blue), and the same signal without white noise (cyan).

Figure 3

Illustration of the graph cut algorithm. (a) The initial graph structure we use to model the step signal. The nodes $i\in \{1...M\}$ represent the variables ${\xi }_i$. (b) In a second step, energies are mapped to capacities of edges. (c) The Boykov-Kolmogorov graph cut algorithm [36] finds the max flow and cuts the graph into two subgraphs $\mathcal{G}=\mathcal{S}\cup \mathcal{T}$, where $\mathcal{S}$ is the part connected to $s$ and $\mathcal{T}$ is the remaining part connected to $t$.

Figure 4

EBS algorithm correctly detects steps in presence of high noise. (a) Noise reduction after application of TVDN. (b) Step detection using combinatorial clustering. Shown is a magnified interval of the simulated noisy data (gray points), boxcar averaged noisy data (20 times reduced, green), simulated steps (blue), denoised signal from TVDN [magenta (a)], and detected steps after application of combinatorial clustering by graph cut [red (b)].

Figure 5

Overfitting transition depends on sampling frequency. Each step signal has 1016 Poissonian distributed steps sampled with different frequencies and covered by noise. The signals are 100 s long. The color bar shows the value of $ln(n/N)$, where $N$ is the total number of data points and $n$ the number of denoised steps.

Figure 6

Influence of SNR on over fitting transition. Each step signal has 980 Poissonian distributed steps with different noise amplitudes at 6 kHz. Every signal is 100 s long and consists of $6\times {10}^5$ data points. The scale of the color bar is chosen analog to Fig. 5.

Figure 7

Detection of steps in noisy signals corrupted by drift. (a) Shown is the noisy signal with drift (gray dots), the drift (black), simulated steps + drift (cyan), and the detected steps (red). The drift is modeled by an Ornstein-Uhlenbeck process ($1/k=9\text{min}, D=10{\mathrm{nt}}^2/\mathrm{s}$) and mimics low-frequency instrumental fluctuations around an equilibrium position. Here, over the simulated time interval of 50 s the signal has a drift of $16.6\text{nt}$ (peak to peak, indicated by green arrow dotted lines). (b) TVDN heuristic of noisy signals with simulated drift is preserved. In the intermediate regime where ${\lambda }_h$ is located, the $\lambda$-heuristic curve of the noisy signal with drift (red) is slightly above the curve of the same signal. (c) With increasing drift steps are detected less precisely. Shown is the precision (blue) and recall (red) of peak-to-peak difference of the drift [mean of 25 signals, simulated according to the slow scenario; SEM bars (not shown) are smaller than the size of the squares].

Figure 8

Performance of step-detection algorithms with respect to slow scenario (blue bars), intermediate scenario (green bars), and fast scenario (yellow bars). Panel (a) shows, in percent of the total number of simulated steps, the number of detected steps. (b) Percentage of correct steps among the simulated steps. Error bars are SEM. Panel (c) shows average step size histograms with 1 bp binning of the detected steps of the fast scenario for $t$ test, HMM, K & V, and EBS (from upper to lower histogram).

Figure 9

Dwell time distribution of detected steps in the fast scenario. Displayed are dwell time histograms (blue bars) of the $t$ test, HMM, K & V, and EBS algorithm (from top to bottom panel). For better comparison, the histogram of simulated dwell times is depicted in the lowest panel. Each dwell time histogram is fitted by a double exponential decay (red line), which yields the following rates ($k_{\mathrm{pause}}/k_{\mathrm{elong}}$ in Hz): $t$ test ($0.017/8.9$), HMM ($0.088/3.9$), K & V ($0.01/2.8$), EBS ($1.3/13$), simulation ($2.8/22$). Note that in case of the detected dwell times the first two bars are not taken into account since steps with short dwell times are likely to be skipped by step-detection algorithms.

Figure 10

The Kullback-Leibler divergence of dwell time histograms of the detected dwell time distribution with respect to the simulated one.

Figure 11

(a) 2.5-bp substeps in $\phi 29$ bacteriophage data (circles) measured in an optical tweezers experiment, EBS (blue) and $t$ test (red). (b) Paused regions in experimental Pol II transcription elongation data. Shaded regions indicate paused intervals found by the SG-filtering method with a velocity threshold of two standard deviations of the pause peak (green) and EBS (red). 1 kHz sampled transcription data (gray), SG-filtered data of polynomial order 3 and frame size: 2.5 s (cyan) and step-detection result of our method (blue). (c) Step size histogram of detected steps by EBS in the Pol II data shown in (b). (d) Magnification of a detected pause. Shown is EBS signal (blue), SG-filtered signal (black), measured data at 1 kHz (black circles), and paused regions (SG, green; EBS, red).

Figure 12

Situation in an Markov random field concerning a single variable $v_i$ and edges $A, B$ to special variables $s$ and $t$ relevant for the data term.

Figure 13

Two neighboring variables in a Markov random field and edges relevant for the pairwise term. Here the two $s$ vertices represent the same vertex in the graph and are just drawn separated to make the diagram clearer. The same is true for the $t$ vertices.

Figure 14

Mean run time performance of combinatorial clustering stage for 10 simulated noisy step signals. (a) Graph cut computation time versus number of tuples for a label grid of 800 levels. Second order polynomial fit to computation times (red curve). (b) Graph cut computation time versus label grid size for a fixed system size of 750 tuples. Linear scaling of performance with increasing number of levels in the label grid set. The error bars are SEM.

Figure 15

Graph cut MCMC comparison. Energies of optimal solutions with increasing system size for MCMC and graph cut algorithms.

Figure 16

Step and noise simulation procedure. State transition model and corresponding stepping rates determine the probability distribution from which a piecewise constant signals (first inset) is sampled. In a second step noise is simulated with the piecewise constant signals as the mean (second inset).

Figure 17

Simplified stepping model of pol II with an elongation and backtracked states.

Figure 18

Histogram of step sizes of the slow scenario (a) and the intermediate scenario (b) for $t$ test, HMM, K & V, and EBS (from upper to lower histogram).

Figure 19

Cumulative fraction of found steps for the EBS in the intermediate scenario. Plotted is the cumulative number of found steps/simulated steps of the signal plotted in Fig. 4 for different dwell times. The number of detected steps compared to simulated ones is smaller for short dwell time steps. Dwell time histograms with a binning of 87.3 ms were determined for the detected and simulated steps, respectively. The number of detected steps for each dwell time was divided by the corresponding number of simulated steps and cumulatively summed up. In total, $60.5\%$ of the number of simulated steps were found.

Figure 20

Prior terms and level grid regularize combinatorial clustering. (a) Relative frequencies of correct steps among the number of detected steps (precision (a), recall (b)) as a function of prior potential strength $1/{\rho }_S$ (${\rho }_S/{\rho }_P=0.08$ is kept constant) for simulated data using the intermediate scenario. Shown is the precision for clustering with a level grid of 1/4 bp spacing (black triangles) and with a spacing of 1 bp (gray squares). (c) Step size histograms of detected steps with a label grid of 1/4 bp without prior terms (blue), with a spacing of 1/4 bp and prior terms (green), and with a spacing of 1 bp and prior terms of the same strength (red). The computed precision (a) and recall (b) corresponding to the three histograms is encircled with the respective colors (a).

Figure 21

Detected steps and intermediate results of EBS application to the $\phi 29$ example (Results and Discussion). Shown is the noisy data [gray circles (a)], the corresponding $\lambda$-characteristic curve [green curve (b)], TVDN with respect to optimal ${\lambda }_h$ [black box (b) and red signal (a)], level grid for CC [black horizontal lines (a)], and steps detected by EBS [blue signal (a)]. Panel (c) shows the sum of detected step sizes of a set of 40 $\phi 29$ measurements.

Figure 22

Backtracked pause detection in simulated data. Shown are the noisy input signal (circles), the simulated step signal (blue), the Savitzky-Golay filtered signal (black), and the detected step signal from EBS (red). Pauses in simulated data are highlighted in green and paused regions in step detected data are indicated by the blue shaded areas.