Artificial data generation

Artificial force spectroscopy data were generated using a custom-written python script (supporting material). The fully folded part of FD curves was modeled using an equation for extensible worm-like chain (WLC) models (Eq. 4). The partially unfolded region was modeled using a combination of WLC and freely jointed chain (FJC) models (Eqs. 5 and 6). For a more detailed description, see the supporting material.

Force-ramp data analysis

For the identification of (un)folding events, we employed a derivative-based approach, which has been previously demonstrated to allow efficient step recognition (23). There are also other algorithms available that are based on probabilistic approaches, such as FEATHER (22). However, it must be noted that these tools are mostly developed for the analysis of atomic-force-microscopy-generated data (20, 21, 22, 23, 24, 25). Here, we aimed to combine step recognition with downstream data fitting and determination of work, based solely on recorded force and distance values. Furthermore, we aimed to keep the pipeline intuitive and adjustable to user requirements. Although this tool was initially developed for the analysis of Lumicks FD data in H5 format, in principle, POTATO can be employed to analyze any data set format independent of the type of OT instrument.

Statistical analysis

In force-ramp trajectories, an unfolding event is characterized by a simultaneous drop in force and a quick increase in distance as the secondary structure of the polymer undergoes a sudden transition from the folded to the unfolded state (Fig. 1 C). Refolding events have opposite characteristics, in which the distance decreases and the force increases upon refolding. When flipped, the refolding data cannot be distinguished from the unfolding data and the processing, therefore step identification can be performed in an identical manner. Ultimately, these (un)folding events can be identified as a local maximum in the derivative of the distance and a local minimum in the derivative of the force (Eq. 2). When plotted, the numerical derivative data of both distance and force show two populations of values. The first is a normal-like distribution representing the measurement noise, while outliers from the normal distribution represent the second population—the actual (un)folding events. To distinguish real (un)folding events from background noise, we calculate the moving median and the standard deviation (SD). These are then used to separate the normally distributed data from the extreme values outside a given Z score (i.e., number of SDs = 3 by default) (Fig. 1 D). This should include 99.73% of the normally distributed data points. As the initially calculated SD is affected by the outliers, a second SD is calculated from the data points inside the threshold, and the data are sorted again. The cycle is repeated until the difference between initial and secondary SD is <x (with x default = 5%). After the force and distance derivatives are sorted, our algorithm finds the local extrema of the derivatives, representing the saddle points of the (un)folding events in the FD curve. Then, it finds the adjacent crossing points of the derivative with the moving median, representing the start or end of the corresponding unfolding events.

Numerical approximation of the derivatives:

F is force, D is distance, t is time, x is position, and step d is a change in position.

Data fitting

Once the respective (un)folding steps are identified, this information is employed for data fitting. Data fitting is performed on the untrimmed data to model the trajectories more precisely. For the characterization of the mechanical properties of the (bio)polymer under tension, the extensible WLC model is commonly used, relating the applied force and molecular extension (Eq. 3) (32). For that, the FD curve is split into multiple parts. The fully folded part (until the first detectable unfolding step) is fitted with a WLC (32) to calculate the persistence length (dsL_P) of the tethered molecule, while the contour length (dsL_C) is fixed. In addition, baseline and offsets in both force and distance are included in the model to compensate for the experimental variability in the FD curves.

The partially and fully unfolded parts of the FD curves are subsequently fitted using a combined model comprising WLC (describing the folded double-stranded handles) and FJC (Eqs. 4 and 5) or another WLC model (representing the unfolded single-stranded parts) (Eq. 6) (Fig. 1 E) (32,33). To mathematically fit the models, we applied model polymer stretching functions from the free python package pylake (Lumicks).

Extensible WLC model:

X is an extension, L_C is contour length, F is force, L_P is persistence length, k_B is Boltzmann constant, T is thermodynamic temperature, K₀ is stretch modulus, F_offset is force offset, and d_offset is distance offset.

FJC:

$WLC+FJC:$

$WLC+WLC:$

Work calculations

Unfolding and refolding FD trajectories also yield crucial information on the thermodynamic properties of the molecule under study. Accordingly, the work applied by the OT instrument onto the system can be calculated from the area under the FD curve (AUC), here using composite Simpson’s rule (Eq. 7). First, we determine the work applied to the whole construct, including the handles (Fig. 2 A). The total work on the construct is the sum of the AUC of the folded model until the starting point of the step (W_ds) and work performed during the step transition (W_step), represented by the rectangular area of the step length times force average ((F_start + F_end)/2) (Fig. 2 A). In order to extract the amount of work applied only to the structure of interest (W_structure; Fig. 2 C), the work applied to the handles, represented by the AUC of the combined model (W_ss), is subtracted from the sum of the work on the whole construct (Eq. 8; Fig. 2 B and C). It shall be noted that the work derived from these calculations equals the Gibbs free energy of the studied structure provided the system is in thermodynamic equilibrium. However, if the (un)folding trajectories do not coincide, it indicates that the molecule is out of equilibrium. In non-equilibrium scenario, Gibbs free energy can be extracted from the work values (5,18,19,29,34, 35, 36) (Fig. S3). It should be noted that while POTATO performs work calculations, the estimations of free-energy values have to be derived by the user separately.

Open in a separate window

Figure 2

Work determination of a simple hairpin. (A–C) FD curve obtained during force-ramp experiment of a short stem loop of 30 nucleotides. Inlets: the optical tweezers construct stretched between the beads with gray regions indicating to what parts of the construct the calculated work relates. (A) Marked region (gray) corresponding to the work necessary for stretching of the whole construct including the structure of interest. (B) Marked region (gray) corresponding to the work necessary for stretching of the handles and the unfolded single-stranded RNA. (C) Marked region (gray) corresponding to the work necessary for stretching of the RNA structure of interest. See the subsequent analysis in Fig. S3. To see this figure in color, go online.

Numerical integration using composite Simpson’s rule:

where x_j = a + jh for j = 0, 1, …, n-1 with h=(b-a)/n; x₀ = a and x_n = b.

Non-equilibrium work calculation:

W_structure is work needed to unfold the structure of interest. W_ds is numerical integration of the fully folded model, W_ss is numerical integration of the unfolded model, and W_step is numerical integration of the step region between the two models.

Constant-force data analysis

In addition to force-ramp experiments, the algorithm we provide can also analyze constant-force data (Fig. S1). In this way, the dynamics of the structure at a given force can be investigated. This way, the equilibrium force at which the chance of the structure to be folded or unfolded are equal can be derived.

The constant-force analysis accepts the same input formats as the force-ramp batch analysis, and data preprocessing is performed similarly by downsampling and filtering of the data without trimming. First, it is necessary to display the constant-force data in order to optimize the preprocessing parameters and the plot’s axis (Fig. S1 B). At this step, two plots are generated for visualization. In the first plot, distance is plotted against time. Here, the difference in distance corresponds to the change in the contour length of the tethered molecule. The second plot is a histogram of the distance distribution (Fig. S1 C). From this histogram, the number of different folding states can be deduced. Afterward, the histogram is fitted with multiple Gaussian functions. According to the position distribution histograms, the user can interactively provide initial estimates for various parameters including the number, localization, width (SD, Z score), and amplitude of the fits. After the optimization, the model parameters are exported together with the percentage of each folding state as a table in csv format (comma separated values).

Artificial data sets to test the limits of detection

To test the limits of (un)folding events detectable by the POTATO pipeline, an artificial data set was generated (supporting material). In this data set, some curves can show a negative step length that would not be observed in real unfolding events. We considered these steps as non-identifiable and used them as negative controls. The phenomenon of negative steps can mainly be observed for small contour-length changes (ΔL_C) between the models, combined with high force drop (ΔF) values. To test the performance of the algorithm, we defined identifiable steps as events with a drop in force and a simultaneous increase in distance (supporting material). To evaluate if a specific parameter combination results in an identifiable curve, Eq. 9 with x = 0 was solved for all sets of parameters. Each time two parameters were fixed, and the third parameter was optimized.

Minimal step calculation:

where WLC corresponds to expression from Eq. 3, ss refers to the model corresponding to single-strand values, and ds describes the double-stranded region.

A hyperplane showing the interface of theoretically identifiable and non-identifiable steps was generated from these optimized values (Fig. 3 A). This allowed us to classify the generated data set based on a combination of parameters: one with curves where POTATO is expected to find an unfolding step (x > 0) and the other one where POTATO should not identify the steps (x ≤ 0). After analyzing the artificial data set (comprising 2520 curves) with different Z scores, the expected results, based on the input parameters when the data were generated, were compared with the steps identified by POTATO. For the default Z score of 3, the expected parameters were then plotted into the three-dimensional plot and colored based on the identification by POTATO (Fig. 3 A). For an unfolding force of 25 pN, the ΔF and ΔL_C values are shown in a two-dimensional plot, making it easier to identify and compare single unfolding events analyzed with different Z scores. It can be seen that all identified steps at this specific unfolding force are above the theoretical threshold and that more unfolding events are identified at Z score 2.5 than at 3 (Fig. 3 B). Accordingly, the effect of the Z score on the derivative of force (Fig. 3 C) and distance (Fig. 3 D) can be investigated for an individual FD trajectory. In the representative trajectory, the local maximum in the derivatives of distance is above the Z score threshold for both cases. In the derivative of force, the local minimum at the same position is only detected for the lower Z score (Fig. 3 C and D).

Open in a separate window

Figure 3

Testing the limits of POTATO. For each combination of the parameters unfolding force (F_U), force drop (ΔF), and contour-length change (L_C), two parameters were fixed, and the third one was optimized so that the Eq. 9 (supporting material) evaluates to zero. (A) A hyperplane was generated from the optimized values that separate the resolvable space above the hyperplane (parameter combinations that result in identifiable steps) from the unresolvable space below the hyperplane (parameter combinations that result in unidentifiable steps). Each analyzed curve is plotted in blue if its step was identified by POTATO or in gray if it was not recognized. (B) Slices of the three-dimensional plot at F_U = 25 pN were analyzed with different Z scores. The black line corresponds to the theoretical limit of resolvable/unresolvable parameter combinations. The black dots represent curves with identified steps, whereas the gray dots represent curves where POTATO could not identify the step. (C and D) The derivatives of force (C) and distance (D) of the curve that is marked with a red arrow in (B) are displayed at different Z scores.

Next, we calculated performance measures such as accuracy, precision, sensitivity, specificity, and F1 score to validate the performance of POTATO. For a Z score of 3.2, a precision score of 0.974 indicates that most of the positive classified steps were actual steps, and even for a Z score of 2.5, the precision was still above 0.944 (Table S2). As expected, higher precision comes with the trade-off to miss certain positive events (recall 0.870–0.939), and the optimal Z score has to be chosen depending on the application. For smaller unfolding events that are difficult to detect, lower Z score should be employed, as for distinct unfolding events, the Z score can be set to higher values. This way, the number of false-positive events detected can be minimized. Since the present data set was generated using artificial parameter combinations, those might not be found in actual OT measurements. Therefore, it is important to keep in mind that we were exploring the limits of the tool by using these strict parameter constraints. Performance measures would also vary depending on where a specific data set is located in the parameter space and which Z scores were employed.

Furthermore, we investigated how accurately POTATO estimates step parameters (F_U, ΔL_C, ΔF). For that, we compared the expected and measured values of these parameters for all curves analyzed (Fig. 4). We then calculated the linear regression of the true positive values to estimate possible biases of POTATO-estimated F_U and ΔL_C values. Our analysis shows that in the case of F_U (Fig. 4 A), the values determined by POTATO are in perfect agreement with the expected values (slope of the linear regression = 0.9912). For ΔL_C (Fig. 4 B), the comparison shows a broader distribution of the measured values, with an overall trend suggesting a minor overestimation (slope of the linear regression = 1.0282) of around 3%. Lastly, in the case of ΔF (Fig. 4 C), the trend shows a slight underestimation of the measured values (slope of the linear regression = 0.8517), resulting in a bias of 12%–15%. Taken together, our performance-measures analysis suggests that the presented tool successfully identifies most (un)folding events correctly with only few false classifications (false positives/false negatives). Accordingly, in most of the cases, performance measures were above 0.9 (Table S2). Moreover, we show that POTATO can precisely estimate the parameter values describing the (un)folding events (F_U, ΔL_C, ΔF; Fig. 4). Overall, the performance measures and the accuracy of the estimates show that POTATO represents a reliable tool for optical tweezer data analysis.

Open in a separate window

Figure 4

Evaluation of the performance of POTATO. The parameters used for the generation of the data set compared with the parameters identified by POTATO are plotted against each other. All three parameters used for the data generation are evaluated with a Z score of 3. (A–C) The values of the true positive steps (black) and the values of the false-positive steps (gray) are visualized for (A) the unfolding force (F_U), (B) the contour length change (ΔL_C), and (C) the force drop (ΔF). A dashed line represents the theoretical perfect correlation between measured and expected value.

Applicability of POTATO on real experimental data

Next, we employed POTATO to test its performance on real experimental data generated from FD measurements of the programmed ribosomal frameshifting element of the encephalomyocarditis virus and severe acute respiratory syndrome coronavirus 2 (27,28). We compared the POTATO results with manually annotated steps of a subset of our data set. The results obtained with manual step identification and data fitting were in good agreement with the automated analysis using the pipeline (Fig. S2 A). Harnessing POTATO in the data processing allowed us to speed up the analysis significantly compared with previous manual analysis. Furthermore, we saw that POTATO is not only suitable for curves with a single (un)folding event like in the artificial data set, but we successfully fit FD curves with as many as five unfolding steps, and we were able to identify even short-lived intermediate states of the unfolding process (Fig. S2 B and C). In addition to the contour-length change obtained by curve fitting, the Gibbs free energy is also an important variable to conclude on the nature of the (un)folded structure as it is dependent on the base pairing of the RNA. We were able to use the work calculated by the POTATO to estimate the Gibbs free energy of the structures and thereby distinguish between different secondary structures (27). Here, to demonstrate the energy calculation, we used a stem-loop mRNA of 30 nucleotides in length (Fig. S3) (28). First, we used mfold (37) to predict the secondary structure and its Gibbs free energy (Fig. S3 A). Then, we plotted the unfolding as well as refolding work distributions calculated by POTATO (Fig. S3 B). We then employed the results of POTATO analysis to estimate the Gibbs free energies by applying 1) Crooks fluctuation theorem and 2) Jarzynski equality with bias correction (Fig. S3 C) as described in (18,34, 35, 36).

To evaluate the performance of POTATO on other published data sets generated using a self-built OT instrument, we analyzed the severe acute respiratory syndrome coronavirus 2 pseudoknot RNA FD data by Neupane et al. (29). Since the data set provided had a lower data frequency, resulting in less than 250 data points per FD curve, we first had to artificially augment the datapoints (see supporting material). Despite that, we could still successfully assign the steps and reproduce the unfolding force distribution (Fig S2) as well as the contour-length estimate (Table S3). We were also able to detect the refolding steps’ force distribution and detected steps as low as 6 pN (Fig. S2). In conclusion, regardless of the system used, we demonstrate that the pipeline output matched well with manual data analysis on real-experiment data sets and that POTATO performed analysis of FD trajectories with multiple steps or even short-live intermediates in a reliable way. Therefore, POTATO represents a versatile tool for high-throughput OT data analysis for many upcoming studies.

Limitations of the study

Processing automation comes with trade-offs (38,39). First, the statistical analysis applied in the pipeline might be prone to false-positive event discoveries due to external causes, such as vibration that might induce step-like events in the FD profile of gathered data. We split the FD data and analyze the derivatives of force and distance separately to minimize this effect. Only the events found by both approaches are considered real (un)folding events. Therefore, the robustness of the analysis is increased.

Second, the pipeline output strongly depends on parameters and threshold values that are applied throughout the analysis. The default values were set empirically to suit our needs. Therefore, it might require optimization to fit specific needs and reach an analysis output consistent with the manual data analysis. User input is required despite the user-friendly GUI environment, and an understanding of the analysis workflow is necessary to adjust the parameters rationally.

The current algorithm does not annotate the repeated folding and unfolding of a structure during force-ramp measurements and identifies this oscillation as independent steps. Nevertheless, this mainly occurs at slow loading rates and does not affect the contour-length estimates. To overcome any unexpected issues with the automated analysis, POTATO also includes a tab that allows full manual analysis of the force-ram data files. This should help to eliminate bias caused by omission of certain files from the analysis during the automated analysis.

We thank Vojtech Vrba for helpful python discussions. We thank Dr. Anke Sparmann for critically reviewing the manuscript. The work in our laboratory is supported by the Helmholtz Association and European Research Council (ERC) grant no. 948636.