Inferring DNA sequences from mechanical unzipping data: the large-bandwidth case

V. Baldazzi, S. Bradde, S. Cocco, E. Marinari, and R. Monasson

Phys. Rev. E 75, 011904 – Published 8 January 2007

Abstract

The complementary strands of DNA molecules can be separated when stretched apart by a force; the unzipping signal is correlated to the base content of the sequence but is affected by thermal and instrumental noise. We consider here the ideal case where opening events are known to a very good time resolution (very large bandwidth), and study how the sequence can be reconstructed from the unzipping data. Our approach relies on the use of statistical Bayesian inference and of Viterbi decoding algorithm. Performances are studied numerically on Monte Carlo generated data, and analytically. We show how multiple unzippings of the same molecule may be exploited to improve the quality of the prediction, and calculate analytically the number of required unzippings as a function of the bandwidth, the sequence content, and the elasticity parameters of the unzipped strands.

24 More

Received 19 May 2006

DOI:https://doi.org/10.1103/PhysRevE.75.011904

Authors & Affiliations

V. Baldazzi^1,2,3, S. Bradde^2,4, S. Cocco², E. Marinari⁴, and R. Monasson³

¹Dipartimento di Fisica, Università di Roma Tor Vergata, Roma, Italy
²CNRS-Laboratoire de Physique Statistique de l’ENS, 24 rue Lhomond, 75005 Paris, France
³CNRS-Laboratoire de Physique Théorique de l’ENS, 24 rue Lhomond, 75005 Paris, France
⁴Dipartimento di Fisica and INFN, Università di Roma La Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy

Article Text (Subscription Required)

Click to Expand

References (Subscription Required)

Click to Expand

Issue

Vol. 75, Iss. 1 — January 2007

Reuse & Permissions

Access Options

Author publication services for translation and copyediting assistance advertisement

Images

Figure 1
Sketch of a fixed-force unzipping experiment: the adjacent $5^{'}$ and $3^{'}$ extremities of a DNA molecule are submitted to a constant force $f$ . The distance between the extremities, $x$ , is measured as a function of time. $x$ is proportional to the number $n$ of open base pairs (bp) up to some fluctuations due to the floppiness of the unzipped strands. The number $n$ of open bp increases or decreases by one with rates $r_{o}$ and $r_{c}$ , respectively, see dynamical model in Sec. 2A.Reuse & Permissions
Figure 2
Free energy $G$ (units of $k_{B} T$ ) to open the first $n$ base pairs, for the first 50 bases of the DNA $λ$ phage at forces 15.9 (dashed curve) and $16.4 pN$ (full curve). For $f = 15.9 pN$ the two minima at bp 1 and bp 50 are separated by a barrier of 12 $k_{B} T$ . Inset: additional barrier representing the dynamical rates (3) to go from base 10 to 9 (barrier equal to $g_{s} = 2.5 k_{B} T$ ), and from base 9 to 10 [barrier equal to $g_{0} (b_{9}, b_{10}) = 3 k_{B} T$ ], see text.Reuse & Permissions
Figure 3
Number of open base pairs as a function of the time for various forces (shown on figure). Data show one numerical unzipping (for each force) obtained from a Monte Carlo simulation of the random walk motion of the fork with rates (3).Reuse & Permissions
Figure 4
Fork position $n$ as a function of time $t = i \times Δ t$ with $i$ integer valued; the sojourn times on each base are given. We call $t_{i}$ the total time spent on base $i$ , and $u_{i}, d_{i}$ the numbers of $i \to i + 1, i \to i - 1$ transitions, respectively. Assuming the fork does not come back to $n = 1$ or 2 at later times, we have: $t_{1} ∕ Δ t = 9$ , $u_{1} = 2$ , $d_{1} = 0$ , and $t_{2} ∕ Δ t = 5$ , $u_{2} = 1$ , $d_{2} = 1$ .Reuse & Permissions
Figure 5
Number of open bases as a function of applied force, and for $5 \times 10^{6}, 10^{7}, 10^{8}$ Monte Carlo steps. Data are averaged over 100 samples. The durations of the unzippings are, respectively, of 7, 15, and $140 seconds$ . The DNA $λ$ phage includes 48 502 bp. In the inset we report the theoretical estimate of the number of open base pairs, for $10^{7}$ and $10^{8}$ MC steps, of Sec. 4C1.Reuse & Permissions
Figure 6
Fraction $ϵ$ (14) of mispredicted bases as a function of the force for the $λ$ -phage sequence. Data are averaged over 100 samples and shown with standard deviations. The dotted line $ϵ = 0.75$ shows the failure rate for a random choice of one base among the four base values.Reuse & Permissions
Figure 7
(a) Error $ϵ$ as a function of the number of unzippings for the phage. (b) Same as a but without distinguishing $A$ from $T$ and $G$ from $C$ , see text.Reuse & Permissions
Figure 8
Probability $ϵ_{i}$ (a) that bp $i$ is not correctly predicted and Shannon entropy $σ_{i}$ (b) for the first $450 bp$ of the DNA $λ$ phage. Inference is made from $R = 1$ unzipping (dashed line) and $R = 40$ unzippings (full line). The force is $f = 16.4 pN$ , and data are averaged over 1000 MC samples.Reuse & Permissions
Figure 9
(Top) free energy landscape for unzipping at force $f = 16.4 pN$ . Local minima correspond to the portion of the sequence that are best predicted. (Bottom) pairing free energy as a function of the base pair index, with and without window average (Gaussian weight over 20 base pairs).Reuse & Permissions
Figure 10
Error rate $ϵ_{i}$ (semilog scale) as a function of the number of repeated unzippings for base pairs $i = 6$ (a) and $i = 27$ (b) arbitrarily selected, for forces $f = 17.4$ and $40 pN$ . Numerical data are averaged over 25 000 to $10^{7}$ samples, see error bars.Reuse & Permissions
Figure 11
Top: error $ϵ_{i}$ for the first 50 bases of the $λ$ DNA for $R = 1, 50, 200$ unzippings. Bottom: connected correlation $χ_{j, i}$ for bases $j = 6$ and $j = 27$ (black dots) for $R = 50$ unzippings. $χ_{27, i}$ is multiplied by 10 to be more visible; data correspond to $f = 40 pN$ (large force).Reuse & Permissions
Figure 12
Probability distribution $P_{R}$ of the sojourn time $t$ spent on a weak [ $g_{0} (W) = - 1.06$ , ${⟨ t ⟩}_{W} = 0.8 μ s$ , dashed line] and strong [ $g_{0} (S) = - 3.9$ , ${⟨ t ⟩}_{S} = 13.7 μ s$ , full line] bases. Time is rescaled by $1 ∕ R$ (see horizontal axis). The number of unzippings is $R = 1$ (left), $R = 2$ (middle), and $R = 10$ (right). The probability $ϵ$ (12) that a $W$ (respectively, $S$ ) base is not correctly predicted is the area under the dashed (respectively, full) curve right (respectively, left) to the crossing point. As $R$ increases time distributions are more and more concentrated, and the error gets smaller and smaller.Reuse & Permissions
Figure 13
Errors on sequences of, respectively, strong (full line) and weak (dashed line) bases as a function of the number $R$ of unzippings in the infinite force limit and without stacking interaction. The difference of pairing free energies $Δ$ is, from bottom to top, 0.5, 1, and 2.8. We show the results of numerical simulations for $ϵ_{R}^{W}, ϵ_{R}^{S}$ with the error bars for $Δ = 0.5, 2.8$ (full dots: $S$ sequence, empty dots: $W$ sequence).Reuse & Permissions
Figure 14
Probability of misprediction for repeated $W W$ (full line) and $S S$ (dashed line) sequences as a function of the number $R$ of unzippings in the infinite force limit and in the presence of stacking interactions. Here, $g_{0} (W, W) = - 1.5, g_{0} (S, W) = g_{0} (W, S) = - 2.5, g_{0} (S, S) = - 3.5$ . The strong and weak sequences are repeated $S S$ and $W W$ sequences, respectively. Simulation results are shown with the error bars. Notice that the slope of $\ln ϵ$ is about twice the one for the nonstacking case with $Δ = 1$ (Fig. 13), see Eq. (56) and attached discussion.Reuse & Permissions
Figure 15
Probabilities $ϵ_{R}^{S W, S}$ and $ϵ_{R}^{S W, W}$ of mispredicting, respectively, a $S$ (black dots, full curve) and $W$ (empty dots, dashed curve) base in a repeated $S W$ sequence as a function of the number $R$ of unzippings in the infinite force limit. The stacking interactions are $g_{0} (W, W) = - 1.5, g_{0} (S, W) = g_{0} (W, S) = - 2.5, g_{0} (S, S) = - 3.5$ . Simulation results are shown with the error bars, while continuous curves correspond to the theoretical expression (52). As $R$ grows the prediction on a single base becomes essentially random ( $(ϵ \to \frac{1}{2})$ since $S W S W \dots$ and $W S W S \dots$ sequences cannot be distinguished from one another).Reuse & Permissions
Figure 16
Connected correlation function $χ_{1}$ at distances $d = 1$ for, respectively, repeated $S S$ (a) and $W W$ (b) sequences as a function of the number $R$ of unzippings in the infinite force limit ( $f = 40 pN$ in simulations). For comparison we show the $d = 0$ correlation function, $χ_{0} = ϵ (1 - ϵ)$ .Reuse & Permissions
Figure 17
Average number $⟨ u_{i} ⟩$ of openings of bp $i$ for the $λ$ -phage sequence during one unzipping for forces 16.4 and 17.4 pN. (a) Theoretical values in the limit of infinite time. (b) Numerical values from MC simulations with $M = 10^{7}$ steps. Note that the infinite time theoretical values coincide with the numerical values up to some base index $i_{\max}$ such that $\sum_{i < i_{\max}} ⟨ u_{i} ⟩ ⪡ M$ , e.g., $i_{\max} ≃ 200$ for $f = 16.4 pN$ and $M = 10^{8}$ steps.Reuse & Permissions
Figure 18
Probability $ϵ$ of misprediction on repeated sequences of $W$ (empty dots, dashed lines) and $S$ (black dots, full lines) bases for pairing free-energy differences $Δ = 2.8$ (a) and $Δ = 0.5$ (b) in the absence of stacking. For each case we show the error as a function of the number $R$ of unzippings for forces above the critical force by 0.5, 2, and 10 pN. The decay constants have for $Δ = 2.8$ the following values: $R_{c} (f = \infty) = 32$ ; for $f = f_{c} + 10$ pN, $R_{c} = 28.5$ ; for $f = f_{c} + 2$ pN, $R_{c} = 10.9$ ; for $f = f_{c} + 0.5$ pN, $R_{c} = 3.4$ .Reuse & Permissions
Figure 19
Phase diagram in the number of unzippings vs force plane. Efficient prediction is possible above the critical line $R_{c} (f)$ (72). Here $g_{0} (W) = - 1.06$ , $g_{0} (S) = - 1.55$ . The full line indicates the repeated $W$ sequence, the dashed line corresponds to the repeated $S$ sequence. For forces smaller than the critical value $f_{c} ≃ 9 pN$ for the W sequence, $f_{c} ≃ 12 pN$ for the S sequence (vertical lines) the molecule remains closed. At large force the number of required unzippings reaches a common value $R_{c} ≃ 30$ .Reuse & Permissions
Figure 20
Theoretical values for the number $R_{c} (f, i)$ of unzippings necessary for a good prediction of base $i$ at force $f = 17.4$ (full line) and $f ⩾ 40 pN$ (dashed line) for the first 400 bases of the $λ$ phage sequence obtained from formula (70).Reuse & Permissions
Figure 21
Decay constant $R_{c}^{+ & -} (f, i)$ of the prediction error on base $i$ for the first 450 base pairs of the $λ$ -phage DNA, at force $f = 17.4 pN$ from the two-way unzipping numerical unzipping.Reuse & Permissions
Figure 22
Probability distribution of a $j$ -base jump for $Δ t$ between $10^{- 5} s$ and $10^{- 3} s$ and for forces $f = 16.4$ and $17.4 pN$ . Notice that the probability is not necessarily a monotonously decreasing function of $∣ j ∣$ , see extra humps in the right column, due to sequence effects.Reuse & Permissions
Figure 23
Error $ϵ$ as a function of the delay $Δ t$ between measures for various ranges (shown on figure). (a) Case of one unzipping $(R = 1)$ of a $λ$ -phage DNA molecule at $f = 16.4 pN$ . (b) Case of $R = 20$ unzippings of a uniform sequence of weak bases at $f = 11.8 pN$ . Results are averaged over 50 samples in both panels.Reuse & Permissions
Figure 24
(a) Fraction of mispredicted bases $ϵ$ as a function of the number of unzippings for different temporal resolutions $Δ t$ . The value of the range is $J = 4$ . (b) Same as right but we only discriminate among strong and weak bases. Data refer to the opening of a $λ$ -phage sequence at $f = 16.4 pN$ and they are averaged over 50 samples.Reuse & Permissions
Figure 25
Distribution $A (x ∣ i)$ of the extension $x$ of the open ssDNA at fixed position of the opening fork, $i = 1$ and $i = 10$ . The r.m.s. of the distribution (at a force of $16 pN$ ) increases as $\sqrt{i}$ . The apparent value of the number of opened bases corresponding to a given $x$ , $i^{a}$ (88), is shown on the top axis.Reuse & Permissions
Figure 26
Value of the number of unzippings controlling the decay of the error in predicting a base, $R_{c} (i)$ , as a function of the base index $i$ . The sequence is made of bases $S$ with a single $W$ base at position $i$ . The dotted line shows the value of $R_{c}$ in the absence of ssDNA fluctuation, for a difference of free energy between $S$ and $W$ bases equal to $Δ = 0.5$ . (a) Case $Δ t = ⟨ t ⟩ ∕ 10$ . The decay constant $R_{c} (i)$ for the Bayesian error (92) grows as $\sqrt{i}$ (dashed line). (b) Case $Δ t \to 0$ . The full line shows $R_{c} (i)$ for the Viterbi procedure without deconvolution; for $i ⩾ 7 R_{c} (i)$ is infinite, meaning that the $W$ base is almost surely predicted to be of $S$ type. With appropriate deconvolution the dotted line value for $R_{c}$ is recovered.Reuse & Permissions
Figure 27
Probability that a base is correctly predicted, $1 - ϵ_{i}^{W}$ , as a function of its location $i$ in the case of: a repeated sequence of $S$ bases with a single $W$ base at position $i$ (black dots), an alternate sequence $S W S W \dots$ (empty dots). In both cases the rate of success decreases exponentially with the square root of $i$ . The difference of free energies between $S$ and $W$ bases is $Δ = 2.8$ .Reuse & Permissions
Figure 28
Symbolic representation of the calculation of the distribution of likelihood at site $i$ (a), and the joint distribution at sites $i$ and $i + 1$ (b). The left part of the sequence induces a left-to-right likelihood $h_{i}^{\to}$ on base $i$ , while the right part contribution is $h_{i}^{\leftarrow}$ (a) or $h_{i + 1}^{\leftarrow}$ (b).Reuse & Permissions
Figure 29
Patterns A, B, and C present in the transition trace around base pair $i$ . See text for definition.Reuse & Permissions
Figure 30
(a) Rate function $ω (\hat{u})$ governing the large deviations of the number $\hat{u}$ of openings of a base per unzipping. $ω$ vanishes when û equals its average value, $⟨ \hat{u} ⟩$ , and is strictly negative otherwise. (b) $R_{c}$ vs logarithm of the number of samples, $μ = \ln M ∕ R$ , for the ninth base of the $λ$ -phage sequence.Reuse & Permissions
Figure 31
Logarithm $μ_{1}$ of the number of samples (divided by $R$ ) to obtain a good estimate of $R_{c} (i)$ vs base pair index $i$ . $μ_{1}$ strongly depends on the base index $i$ , e.g., we need to sample over $M \sim e^{8 \times R}$ to accurately estimate $R_{c}$ for all bases. The force is $f = 17.4 pN$ .Reuse & Permissions

Physical Review E

covering statistical, nonlinear, biological, and soft matter physics