Free-energy simulations reveal molecular mechanism for functional switch of a DNA helicase
Abstract
Helicases play key roles in genome maintenance, yet it remains elusive how these enzymes change conformations and how transitions between different conformational states regulate nucleic acid reshaping. Here, we developed a computational technique combining structural bioinformatics approaches and atomic-level free-energy simulations to characterize how the Escherichia coli DNA repair enzyme UvrD changes its conformation at the fork junction to switch its function from unwinding to rezipping DNA. The lowest free-energy path shows that UvrD opens the interface between two domains, allowing the bound ssDNA to escape. The simulation results predict a key metastable 'tilted' state during ssDNA strand switching. By simulating FRET distributions with fluorophores attached to UvrD, we show that the new state is supported quantitatively by single-molecule measurements. The present study deciphers key elements for the 'hyper-helicase' behavior of a mutant and provides an effective framework to characterize directly structure-function relationships in molecular machines.
https://doi.org/10.7554/eLife.34186.001Introduction
Helicases are ubiquitous motor proteins that move along nucleic acids and separate duplex DNA or RNA into its component strands. This role is critical for various aspects of DNA and RNA metabolism; defects in helicase function in humans can lead to genomic instability and a predisposition to cancer (van Brabant et al., 2000; Brosh, 2013). Characterizing the atomistic mechanism for helicase function, although challenging, is crucial to link protein structure with their function and help engineering helicases with novel activities (Arslan et al., 2015).
DNA helicases can unwind double-stranded DNA (dsDNA) into single-stranded DNA (ssDNA), which are later copied during DNA replication or modified in DNA repair processes (Wu and Spies, 2013; Lohman et al., 2008). They are classified into six superfamilies (SF), among which SF1 and SF2 helicases are the largest superfamilies and share many similar conserved motifs. The minimal functional units for SF1 and SF2 helicases are monomers that contain two RecA-like motor domains for ATP hydrolysis (Singleton et al., 2007). SF1 helicases can unwind dsDNA by translocating on a ssDNA strand as shown in Figure 1a. Such translocation happens in a stepwise manner, during which the chemical energy from ATP hydrolysis is used to break the bonds in dsDNA via conformational changes of the motor domains (Yang, 2010; Patel and Donmez, 2006). An exemplary Escherichia coli helicase, UvrD, belonging to SF1, has many cellular roles such as methyl-directed mismatch repair (Iyer et al., 2006; Spies and Fishel, 2015) and nucleotide excision repair (Sancar, 1996) by unwinding duplex DNA. UvrD can also prevent deleterious recombination by removing RecA filaments from ssDNA (Lestini and Michel, 2007). Along with its homologous proteins PcrA and Rep, UvrD has been demonstrated in experiments to translocate on ssDNA progressively 3’ to 5’ (Matson, 1986; Mechanic et al., 1999; Dillingham et al., 2000; Myong et al., 2005; Fischer et al., 2004). Structures of UvrD-like SF1 helicase solved so far share a four-subdomain tertiary arrangement (1A/2A/1B/2B) (Singleton et al., 2007), including two RecA-like domains (1A/2A) which contain the ATP binding site and are proposed to function as the translocase (Dillingham et al., 2001; Lee and Yang, 2006), and a flexible domain (2B) which is believed to play a regulatory role in helicase activity (Lohman et al., 2008; Dillingham, 2011). In particular, the 2B domain is known to adopt different conformations (Velankar et al., 1999; Brendza et al., 2005; Jia et al., 2011; Nguyen et al., 2017) and has been proposed to act as a ‘molecular switch’ controlling UvrD unwinding (Comstock et al., 2015).
Combining optical tweezers and single-molecule FRET, Comstock et al. (2015) demonstrated that UvrD can switch its activity from DNA unwinding to rezipping (measured by optical tweezers) by dramatically changing its conformation between two states (detected by FRET). The transition from unwinding to rezipping activities was proposed to occur through switching ssDNA strands, accompanied by rotation of the 2B domain (see Figure 1b). In this model, the GIG motif on 2B serves as an anchor point on dsDNA above the fork junction, such that rotation of 2B can position the 1A/2A translocase domains on either ssDNA strand, leading to 3’ to 5’ UvrD translocation either toward (unwinding) or away from (rezipping) the DNA fork. Two crystal structures seem to support this strand-switching model (see Figure 1c): one structure of UvrD (pdb code: 2IS2) (Lee and Yang, 2006) bound to a dsDNA-ssDNA junction is expected to be the ‘unwinding’ state (defined here as the ‘closed’ state) because its 1A/2A domains would translocate UvrD into the DNA fork; the other structure (pdb code: 3LFU) (Jia et al., 2011) solved without DNA is expected to represent the ‘rezipping’ state (defined here as the apo state) because the 1A/2A domains presumably would be bound to the opposing strand, translocating UvrD away from the DNA fork. The structural differences between closed and apo states mainly involve a simple rotation of the 2B domain (Figure 1c).
However, in order for the ssDNA strand-switching to happen, the rezipping state must adopt a conformation with a gap between the 1B and 2B domain that is large enough for the bound ssDNA to escape, whereas in both the closed and apo structures the four domains 1B-1A-2A-2B form a closed ring topologically. As we show here, contrary to the common assumption that the apo structure is a functional state of UvrD, the FRET signal simulated using real fluorophores attached to the apo-state structure does not match the experimentally observed signal of the rezipping state, nor the unwinding state. Furthermore, it has been reported that cross-linking the 2B and 1B domains of the SF1 helicase Rep can change it into a superhelicase (Arslan et al., 2015), capable of unwinding thousands of base pairs processively. What are the key regulatory factors for the functional switch and is it possible to design mutants with different activities?
To characterize the conformational states of UvrD at the fork junction and the transitions between those states, we use MD simulations, which are well-suited to study atomic-level mechanisms in conjunction with crystallography, single-molecule and biochemical techniques (Russel et al., 2009; Zhao et al., 2010; Arkhipov et al., 2013; Cheng et al., 2017; Latorraca et al., 2017). However, due to the very long time-scale of conformational changes, brute-force simulations are challenging in the case of large molecular motors such as UvrD. Here, we employed a novel computational approach which integrates advanced sampling simulations with bioinformatics tools that survey structural information from homologs. We were able to identify modes of motions for function switching from principal component analysis of a ‘trajectory’ derived from the alignment of various surveyed crystal structures. Using the first two principal components as reaction coordinates, the subsequent all-atom Hamiltonian replica exchange simulations (totaling 12s) predict a metastable ‘tilted’ conformation, which has significantly lower free energy than the apo state. The lowest free-energy path is determined to describe the transition between the closed state to the ‘tilted’ state. After the closed-to-tilted transition takes place, 2B and 1B domains are separated with enough distance from each other to enable strand-switching to happen. We demonstrate that ssDNA can be disengaged from the ssDNA binding domains of UvrD in the tilted state. Furthermore, the tilted UvrD structure is shown to be able to form stable interactions with the opposing strand after ssDNA strand switching has occurred. We also highlight the role of the GIG motif in assisting 2B domain diffusion along dsDNA during strand-switching. These findings suggest principles underlying mechanisms of related molecular machines beyond what we have known from existing structures.
The properties obtained from the transition pathway are consistent with the single-molecule data (Comstock et al., 2015) as well as mutagenesis studies (Meiners et al., 2014). Firstly, we carried out equilibrium simulations of UvrD site-specifically labeled with FRET dye pair AlexaF555/AlexaF647 for both the closed state and the tilted state. The calculated average FRET efficiencies for the two states are in good agreement with those for the unwinding and rezipping states measured in single-molecule experiments, respectively. These simulations also allow us to obtain key fluorophore conformations in the tilted state to explain the shape of the experimental FRET distribution. Secondly, we illustrate the molecular basis for hyper-helicase activity of a UvrD double mutant (D403A/D404A) for the first time. Finally, a physical model integrating the simulation results and the measured equilibrium constant from optical tweezers experiments is provided to explain the helicase function-switching mechanism.
Results
Structural bioinformatics analysis of conformational ensembles of UvrD-like proteins
Our goal is to characterize UvrD conformational changes that switch its function. Recently, free-energy simulation methods have been successfully applied to study transitions between two functional conformational states of complex molecular machines (Moradi and Tajkhorshid, 2013; Ma and Schulten, 2015; Czub et al., 2017). However, for UvrD all the known structures bound to the DNA fork junction belong to the closed (unwinding) state (Figure 1c). It is unclear whether the apo state of UvrD could bind to the dsDNA-ssDNA junction. By aligning the apo state to the closed state, we found geometrical clashes between the fork junction and the apo state (Figure 1—figure supplement 1b). We thus forced UvrD at the fork junction to rotate from the closed state to the apo state using targeted molecular dynamics (Schlitter et al., 1994) (see Figure 1—figure supplement 1c for details). However, such an operation experienced large resistance (DNA was free to move), and the protein returned back to the vicinity of the closed state after the external force was released. We thus need to find new conformations that can represent the rezipping state.
In order to reach the rezipping state while bound to the fork junction, UvrD must reach some hidden metastable states, which can be far away from the 2B-domain-rotation pathway around the dsDNA axis. To identify such states, we developed an approach based on surveying the pdb database (details in Materials and methods). We used protein-protein BLAST (basic local alignment search tool) to search the swissport database with the UvrD sequence as the query sequence. Then, we downloaded the pdb files of these homologs with 40% or more sequence identity. A subsequent principal component analysis (PCA) was carried out to find out the most significant degrees of structural variations among UvrD and its homologs. The coordinates of the homolog structures were then projected onto the first two principal components (PC1 and PC2) (Figure 2). Three distinguishable clusters are shown in Figure 2: one represents the canonical closed conformation, one represents the canonical apo state, and another one represents an interesting conformation (from the replication initiator protein) in which the 2B domain is tilted from the dsDNA axis. All the structures belonging to the apo state are without nucleic acids bound. The structure in the ‘tilted’ cluster only has ssDNA bound, and thus very likely it is not a functional state of UvrD because of the absence of dsDNA interactions. To characterize the functional state for rezipping, we need to carry out all-atom free-energy simulations (the next subsection).
We next calculated the so-called involvement coefficients (Lei et al., 2009) (ICs), which are often used to show the contribution of individual modes to the overall structural displacement. For the displacement between the closed structure and the tilted structure, the ICs of the first two PCs are very high (see Figure 2—figure supplement 1a), indicating that the first two PCs are sufficient to describe the protein conformational changes based on the available UvrD homolog structures. Directions of motions along the first two PCs are shown on the closed structure (Figure 2—figure supplement 1b). We noted that PC1 is in a similar direction as the rotational movement between the closed and apo states. PC2 represents a tilting motion orthogonal to the rotation. Since the closed-to-apo rotation of the 2B domain cannot bring UvrD to the rezipping state due to steric clashes, we suspect that PC2 might make a very important contribution to UvrD conformational switching when bound to the junction.
Free-energy landscape of UvrD conformational ensembles when bound to the fork junction
Based on the information revealed by the PCA analysis, we would like to find the UvrD conformation responsible for the rezipping state when bound to the dsDNA-ssDNA junction. For this purpose, extensive enhanced sampling simulations (12 s in total) were carried out to characterize the free-energy landscape of UvrD conformations and detect any interesting metastable states in it. See Materials and methods for the setup and simulation details.
We first characterized the 2D potential of mean force (PMF) using the first two PCs as coordinates (middle panel of Figure 3). We identified two conformations located in the two local minima of the 2D PMF map, respectively (right and left panels of Figure 3). These two conformations are defined as ‘closed’ and ‘tilted’ states. The tilted state has features that have not been found in any of the existing crystal structures, as we show in the following sections. The PDB file for the newly found tilted state is provided as Supplementary file 1.
The closed and tilted conformations served as the initial and final states for a transition path finding protocol, which was employed to find the lowest free-energy path between them (see Materials and methods). The most probable transition happens in two phases, during which the 2B domain undergoes coupled rotational and tilting motions. In the first phase (closedIM), 2B carries out a large-scale tilting motion along PC2, overcoming a 4.4 kcal/mol barrier before reaching an intermediate state IM. In the second phase (IMTStilted), 2B performs mostly a rotational motion along PC1, overcoming a 1 kcal/mol barrier () at the global transition state (TS) before reaching the tilted state. Thus, the rate-limiting step is the first phase, which involves mostly a tilting motion. Figure 3—figure supplement 1 provides the PMF values and intermediate conformations along the lowest free-energy path. Video 1 shows the conformational changes of UvrD during the transition.
One can notice that the region the apo structure represents has a high energy value, which is more than 8 kcal/mol higher than the initial state. This demonstrates that the apo state, which is connected to the closed state by 2B domain rotation, is very unfavorable at the dsDNA-ssDNA junction.
We took the representative protein structure in the final and initial states and measured the gap size, which is defined by the closest C atom distance between the 2B and 1B domain. The extended ssDNA has a diameter around 10 Å (Landy et al., 2013). The initial closed state has a very small gap size of 6 Å, through which the ssDNA cannot pass. The final tilted state has a gap size of 14 Å, which is open enough for ssDNA to pass through.
The overall free-energy landscape projected along a progress variable is plotted in Figure 4. is proportional to the projection on PC1 and is scaled from 0 to 1.0 between the closed state and the tilted state. The free energy for the metastable tilted state is about 2.5 kcal/mol higher than that of the closed state. The system has to overcome a 4.2 kcal/mol energy barrier at the transition state (TS) to reach the tilted state.
Validation of the predicted tilted state
We first tested if the ssDNA can escape from the tilted structure. To accelerate the process, we used targeted molecular dynamics by adding a harmonic potential to the coordination number between UvrD and ssDNA. The targeted coordination number was forced to change from an initial value of 18 to 0 in 30 ns. As shown in Video 2, the ssDNA is seen disengaged from the ssDNA binding domains of UvrD. The final interaction energy between ssDNA and the 1A/2A/1B domains of UvrD gradually drops to zeros (see Figure 5—figure supplement 4). Further below, we also show that this tilted structure can bind stably to the opposing strand to complete the strand-switching process (Figure 8d).
To quantitatively validate our simulation results against experimental data, we compared the FRET efficiency distributions predicted for the closed and tilted states computationally to those of the functional states measured experimentally. We first obtained the 1D FRET efficiency distributions for the unwinding and rezipping state based on the raw single-molecule data (see Materials and methods for details). The distributions, shown in Figure 5a, have peak positions at 0.66 and 0.29 for unwinding and rezipping, respectively. By explicitly simulating UvrD in the two states with fluorophore labels (AlexaFluor555/AlexaFluor647) as in the single-molecule experiments, we also determined FRET efficiencies for the closed and tilted states (Figure 5b). The simulations accumulated 500 ns for each state, and we considered the orientation factor of the fluorophores in determining the FRET efficiency (Materials and methods). The predicted FRET efficiency peak for the closed state is at 0.72, whereas the peak for the tilted state is around 0.31. The close agreement between experimental and simulated FRET distributions reaffirms that the tilted state should be the protein conformation responsible for rezipping. As a control, we simulated the apo-state structure with the fluorophore labels for 500 ns as well. The apo-state FRET distribution, which peaks at 0.16, is quite different from the rezipping-state distribution (Figure 5—figure supplement 1), suggesting that the apo structure is not the conformation for UvrD rezipping at the junction.
We further examined the representative fluorophore pair conformations at the local maxima (FRET Efficiency = 0.3 and 0.6) of the tilted state FRET distribution (green curve in Figure 5b). It appears that the fluorophores have different conformations at the two different FRET values (Figure 5—figure supplement 3), due to the conformational dynamics of the dyes with the long linkers. The ‘shoulder’ of the tilted-state FRET distribution curve at 0.6 efficiency is caused by a metastable conformation of AlexaFluor555 with different pair-distance and orientation comparing to the conformation at 0.3 efficiency.
UvrD diffusion along dsDNA
In the UvrD functional switching model, the 2B domain of UvrD has to maintain contact with dsDNA; otherwise the protein might disassociate from the fork junction during the ssDNA strand exchanging. It is known that the GIG motif (motif IVc) of UvrD plays a key role in interacting with dsDNA (Myong et al., 2005), and T422 (a representative residue of GIG) is important for UvrD activity (Lee and Yang, 2006). We thus monitored the changes in the interaction between GIG and dsDNA. Figure 6 shows a free-energy landscape plotted against the DNA base ID in contact with GIG and the distance between them. For each simulation frame, we calculated the distances between every DNA residue’s O2P atom and the OG1 atom of T422. Then, the minimal distance and the corresponding DNA base ID were used as the two coordinates. Note the two strands of dsDNA share the same base ID here: for residue x in strand A (indexing according to pdb), the complementary residue in strand B has the same ID x. In the present case, frames with base ID 18 only involve strand A - T422 interaction; whereas frames with base ID 14 only involve strand B - T422 interaction.
In the closed state, residue 18 of strand A contacts the GIG motif, whereas in the tilted state, residue 14 of strand B contacts the GIG motif. Thus, there is a diffusional motion along the dsDNA during the conformational change (see Figure 6). In such a way, UvrD is able to switch the binding dsDNA strand and finds an energetically favorable configuration for the ssDNA strand-switching that will happen in the next step. The diffusion happens in a way that the DNA and T422 are disengaged first, and T422 then re-engages with another DNA residue along the double strand. The base ID in contact with T422 during the transition path from the closed state to the tilted state is shown in Figure 6—figure supplement 1. One can see that the 2B diffusion happens late during the transition. Although UvrD diffuses along dsDNA during the transition, there is no base pair unwound during the closed-to-tilted transition.
Molecular mechanism for the UvrD303 mutant
Our simulations provide a molecular explanation for the hyper-activity reported for a mutant (UvrD303) that involves two important aspartic acid residues at the 2B-1B interface. Previous experimental work (Meiners et al., 2014) discovered that UvrD303 with substitution of two residues, 403 and 404 (both from Asp to Ala), in the 2B domain exhibits a ‘hyper-helicase’ unwinding activity in vitro. The authors suggested that such mutations will reduce the 1B-2B domain interactive contacts and thus yield an intermediate conformation instead of a closed conformation. Such an intermediate state they argued would result in the hyperactivity. However, this explanation is not consistent with the single-molecule measurements (Comstock et al., 2015) showing that the closed conformation is responsible for unwinding activity.
To reconcile the conflict, we estimated for the binding free-energy between the 1B and 2B domains upon mutating D403 and D404 into alanine, based on our enhanced sampling trajectory. Here , where is the binding free-energy for the mutant and is that for the wild type. calculated for the closed state is around kcal/mol, showing a stabilization effect of the double alanine mutant. On the other hand, calculated for the tilted state is around 0. This indicates that UvrD303 actually favors the closed conformation and thus will lead to better unwinding activity. The so-called MM/PBSA method (molecular mechanics Poisson-Bolzmann surface area) (Kollman et al., 2000; Homeyer and Gohlke, 2012) was used for calculating .
Figure 7a shows the configuration of D403/D404 and key residues on 1B that contribute most significantly to the binding energy change upon the mutation in the closed state. The first five residues on 1B with the largest contribution to are listed in Figure 7b (for the tilted state, all the individual residue contributions to become zero). We noted that there are not many positively charged residues on 1B that are very close to D403/D404. The maximum number of hydrogen bonds formed between D403/D404 and the 1B domain is around two pairs during the simulations. Considering that there are also negatively charged residues of 1B (E118/E117) near D403/D404, mutating the two aspartic acid residues into alanine will not decrease but rather increase the interaction strength between 1B and 2B. We also found that there are significant numbers of nonpolar residues located around residues 403 and 404 (L186, A184, L114, I113, L122). Thus, mutating the two charged residues into hydrophobic residues instead increases the interaction strength between the nonpolar groups and the two alanine residues. Overall, the stabilization of the closed state of UvrD303 leads to consistent unwinding of UvrD helicase, reconciling the biochemical measurement (Meiners et al., 2014) with the single-molecule experiment (Comstock et al., 2015).
Discussion
We have characterized the conformational dynamics and a key metastable state of UvrD at a fork junction with a hybrid computational approach. The transition pathway as well as the free-energy landscape for UvrD functional switching at the fork junction was obtained, and we found that the opening of the 2B domain involves a major tilting motion followed by a major rotational motion. Diffusion of 2B along the dsDNA happens in the late stage of the transition, during which the GIG motif switches its contact from one strand of dsDNA to the other strand. The transition leads to a gap opening between 2B and 1B, which enables the ssDNA to escape presumably allowing the motor domains to strand-switch.
A physical model for UvrD functional switching
A schematic model can be established based on the simulation results (Figure 8a,b). The corresponding molecular models are shown in Figure 8d. The UvrD functional switching happens in a two-step manner. A first step is the opening of the 2B domain, followed by a second step of the switching of the bound ssDNA strand, in which the original ssDNA disengages from the 1A/2A/1B domain binding site and the other strand fills in.
To obtain the free-energy difference between the unwinding and rezipping states, we performed a dwell time analysis based on past single-molecule measurements (Comstock et al., 2015). The dwell times of the unwinding and rezipping states of UvrD monomers are plotted in a histogram and the calculated averaged rates for both transitions are almost equal (= 6.6 s-1 and = 7.0 s-1) at 10 M ATP (see Figure 8c). Thus, the equilibrium constant is around one and the unwinding and the re-zipping conformations should have similar free energy. This is consistent with the picture that the tilted state is a little bit less favorable than the initial state but as soon as the ssDNA releases and the other ssDNA strand binds to the UvrD, the system returns to a lower free energy (the rezipping state) (Figure 8a). For the mutant UvrD303, the free energy for the closed state drops around 3 kcal/mol, whereas G and G remain the same. The relative stabilization of the closed state leads to more persistent unwinding.
The strand switching is mostly driven by Brownian motion and does not require energy from ATP hydrolysis. To address the possible effect of ATP on strand switching, we (1) analyzed additional data from the optical tweezers experiments and compared the switching rate at two different ATP concentrations and (2) also analyzed the X-ray structures of the closed state with and without ATP. The dwell time distributions at 10 M ATP and 2.5 M ATP concentration are plotted in Figure 8c. The equilibrium constants of switching measured for the two concentrations are very similar (both around 1), which suggests that strand switching is likely an ATP-independent process. Furthermore, although our simulated system is based on an ATP-free UvrD crystal structure (2IS2), our computational approach covered the structural information from ATP (or its analogs) bound structures (Figure 8—figure supplement 1). One can see that the 2B motion between the ATP-substrate bound and empty UvrD in the closed state is small relative to the large closed-tilted conformational change. Therefore, it is not very likely that ATP binding has a noticeable impact on the closed-to-tilted transition.
DNA-UvrD conformation at the rezipping state
To explore the structure of the rezipping state further, we built a rezipping structure starting from the tilted conformation after ssDNA strand switching has occurred (see Figure 8d and Materials and methods). After a 100 ns equilibration simulation, the modeled system was stable and had a rmsd around 3 Å from the tilted state (Figure 8—figure supplement 2). The newly obtained rezipping structure satisfies the following considerations: (1) The protein conformation is very similar to the tilted conformation. (2) The interaction configuration between 2B and dsDNA remains the same between the tilted state and the rezipping state. Note that during our simulations of the closed to tilted transition, the 2B domain changed its contact from one strand of dsDNA to the other (4 bp shift, from the red strand A to the blue strand B in Figure 6). (3) The ATPase domains 1A-2A are in the correct orientation along the ssDNA (3’ to 5’), pointing away from the junction. Such a conformation enables UvrD to translocate along ssDNA, allowing the duplex to rezip behind it. In examining the rezipping structure, we found a small loop forming between the dsDNA junction and the ssDNA-1A binding site. A similar feature was proposed by a translocation model of PcrA helicase in Park et al. (2010), which suggested that PcrA can extrude a ssDNA loop while it attaches to dsDNA and translocates the 5’ ssDNA tail in an open conformation.
The apo state seen in the crystal structures without DNA bound is likely not a functional state of UvrD at the fork junction. First, the simulated apo-structure FRET distribution is quite different from the rezipping-state FRET distribution. Upon completion of the ssDNA strand-switch, we expect the conformation of UvrD to stay close to the tilted structure. The FRET signal from the single-molecule experiment shows a clear two-state distribution, and our FRET distribution for the tilted state from simulations agrees very well with the experiment. Second, we aligned the apo state structure to the tilted state in Figure 8—figure supplement 3, and there are serious clashes between the apo structure and the dsDNA. Thirdly, the apo state is highly unfavorable at the fork junction according to our simulations.
Functional insights for UvrD and its homologs
Our simulations, backed by the single-molecule measurements, provide functional insights for UvrD in several important biological processes. For example, frequent strand switching of UvrD due to 2B conformational transition results in unwinding over short distances (Comstock et al., 2015), which is consistent with the small number of basepairs unwound during nucleotide excision repair (Kisker et al., 2013). On the other hand, UvrD303 is associated with a recombination-deficient phenotype (Centore et al., 2009), possibly due to lacking such a structural transition as the closed state is over-stabilized. It has been reported that UvrD can dismantle RecA filaments from the ssDNA at a stalled replication fork (Veaute et al., 2005; Lestini and Michel, 2007). As RecA has a central role in homologous recombination (Cox, 2007), a population shift toward the closed state could enhance UvrD’s ability to disrupting RecA-ssDNA filaments and impair recombinational repair.
The tilted state and related motions found here can possibly help connect structural information with function for other SF1 helicases. A highly homologous helicase, PcrA, is known to efficiently strip RecA filaments off ssDNA in an ‘open’conformation (Park et al., 2010). The low-FRET ‘open’ conformation of PcrA could be similar to the tilted conformation revealed in this study. In this case, PcrA is anchored to the dsDNA and translocates the 5’ ssDNA strand in the direction toward the junction. A different mode of PcrA is binding to the 3’ ssDNA and the dsDNA while unwinding the duplex in the closed form (Velankar et al., 1999; Niedziela-Majka et al., 2007). Another UvrD homolog RecB, by mostly tilting its 2B domain from the putative closed state, forms interactions with other subunits in the RecBCD complex (Singleton et al., 2004; Wilkinson et al., 2016), which has a key role in initiating recombinational repair (Spies and Kowalczykowski, 2005).
It may be possible to engineer UvrD-like helicases with tunable unwinding activities. Experiments have shown that cross-linking Rep and PcrA in the closed form resulted in superhelicase activity (Arslan et al., 2015). We demonstrated that mutating the two aspartic acid residues into alanine on 2B domain stabilizes the UvrD closed conformation. The analysis of the binding free-energy change upon the mutation (Figure 7b) provides potential target residues to guide future experimental designs. For example, mutating some negatively charged residues on 1B might also result in hyper-helicase behavior. Our findings for the conformational dynamics of UvrD and the related computational strategy establish a foundation for future studies to reveal principles employed by other related helicase systems.
Materials and methods
Structural bioinformatics analysis of UvrD homologs
Request a detailed protocolOur computational study starts from analyzing the structural ensemble of UvrD homologs. There are two representative conformations for UvrD: one being the so-called ‘closed’ state (e.g. 2IS2); the other one being the apo state (e.g. 3LFU). As stated in Results, the apo state is likely not a functional structure of UvrD at the DNA fork junction. To explore the conformational space of UvrD as much as possible, we performed a structural survey for possible UvrD homologue conformations using bioinformatics sequence and structure alignment tools (Altschul et al., 1997; Cock et al., 2009; Bakan and Bahar, 2009; Bakan et al., 2014). The initial sequence alignments were obtained using NCBI blastp search (Altschul et al., 1997) of Protein Data Bank database sequences, with UvrD as the query sequence. Twenty-six structures were selected from the surveyed structures with sequence identity better than 40 and query sequence coverage larger than 60. The structure alignment was generated by ProDy (Bakan and Bahar, 2009; Bakan et al., 2014) from the pairwise sequence alignments by Biopython (Cock et al., 2009). The resulting 26 structures can be interpreted in terms of a ‘trajectory’ with the coordinates , of which each frame () contains coordinates (from C atoms) of the homologous structures that were mapped onto the original UvrD chain.
We then performed principal component analysis (PCA) (García, 1992; Bakan and Bahar, 2009; Raveh et al., 2016) with ProDy to determine a number of modes for reducing the phase space of UvrD motion. Only the C coordinates of the 2B domain were used for the PCA calculations, after aligning the 1A/2A/1B domains of all the 26 structures to those of the closed state structure (2IS2). The covariance matrix for PCA is determined via where the angular brackets denote the average over (all the frames). The eigenvectors (principal components or PCs) of the matrix are determined by . These PCs, which are ranked by their corresponding eigenvalues, represent different directions of conformational motion away from the original closed state.
The homologous structures were then projected onto the first two PCs with the largest eigenvalues. As stated in Results, a 'tilted' structure based on pdb 1UAA was found as an outstanding cluster among the homologous structures. To see the contributions of different PCs to the displacement between the closed structure and tilted structure, we further calculated the involvement coefficiency (Ma and Karplus, 1997; Lei et al., 2009) of the th PC. is defined as , where is the unit vector describing the displacement from the closed structure to the tilted structure. Only the first two PCs contribute significantly to the overall motion (Figure 2—figure supplement 1A). PC1 and PC2 are used later as coordinates to compute the free-energy landscape.
MD simulation setup
Request a detailed protocolOur simulations were initiated from the closed state (pdb 2IS2) of UvrD (see Figure 1c). The protein-DNA system was solvated in a 100 Å 100 Å 130 Å water box with 55 mM NaCl (the system had 140K atoms in total). A 210-step energy minimization was carried out and the system was then heated to 310 K in 30 ps, employing harmonic constraints with 1 kcal/(mol Å) spring constant to the C atoms. Keeping the spring constant, a 1 ns equilibration in the NPT ensemble (1 atm at 310 K) was performed with a Langevin thermostat for temperature coupling. This was followed by a 1 ns NVT-ensemble simulation, during which the spring constant was gradually decreased to zero. The system was then equilibrated for 60 ns, and the resulting configuration is referred to as the closed state. All MD simulations in our study were performed using NAMD 2.10 (Phillips et al., 2005) with the CHARMM36 force field (Best et al., 2012; Hart et al., 2012).
Free-energy simulation protocol
Request a detailed protocolTo determine the free-energy profile along a reaction coordinate, we employed the Hamiltonian replica-exchange (HREX) method (Park et al., 2012; Jiang et al., 2012; Jiang et al., 2014). HREX uses a series of replicas ( 2,. .., M) of the system, which are simulated concurrently with slightly different Hamiltonians and are exchanged frequently among themselves based on the Metropolis exchange criterion (Sugita et al., 2000). HREX can be very powerful in reconstructing rugged free-energy landscapes by exchanging external biasing potentials, which, with different biasing parameters, are added to the replicas to enhance the sampling throughout the reaction coordinate (RC). The biasing potential (or the window potential) for each replica usually assumes the form of , where is the current value of the reaction coordinate for replica , ( 2,. .., M) is the index for the biasing potentials (windows), is the preassigned parameter for the center of the harmonic potential, and is the spring constant. The centers of the biasing potentials () are selected as an ordered list of values () all over the RC to cover the reaction of interest fully. Exchanges between two neighboring replicas (replicas with neighboring values) are attempted periodically during the simulations. Without the replica-exchange strategy, this protocol reduces to the conventional umbrella sampling, which often suffers from the inefficient sampling of degrees of freedom orthogonal to the reaction coordinate (Jiang et al., 2012).
The present study chooses the projection on the first PC () as the reaction coordinate and includes M = 120 biasing windows between the closed state and the tilted state. The initial configurations for the M windows were generated through a 5 ns targeted MD simulation (Schlitter et al., 1994), by driving UvrD from the closed state to the tilted state. An exchange between two neighboring replicas was attempted every 10 ps and the spring constant of the harmonic potential was set to 100 kcal/(mol Å). The production run of each replica lasted 100 ns, and the total simulation time added up to 12 s (100 ns 120). Eventually the weighted histogram analysis method (WHAM) (Kumar et al., 1992) was applied to obtain the unbiased 1D and 2D free-energy landscapes in Figures 4 and 3. We performed the Monte Carlo bootstrap error analysis (Stine, 1989; Hub et al., 2010) to estimate the uncertainty along the reaction coordinate. The basic idea of bootstrapping is to obtain several estimates (we obtained 10 trials) for the free energy based on randomly generated subpopulations from the histogram in each window. Our simulations with HREX benefitted from a scalable multiple copy algorithm (Jiang et al., 2014) which enables simulating hundreds of replicas simultaneously on a petascale supercomputer.
As stated in Results, the tilted state structure was identified as one of the most important metastable states. Based on the free-energy landscape using the projections on the first two PCs, the lowest free-energy path describing the most probable reaction mechanism was localized between closed state and the tilted state using the optimization algorithm in (Ensing et al., 2005). The path was then smoothed and 120 images were chosen uniformly along the 2D pathway applying the curve-fitting protocol in (Ma and Schulten, 2015).
FRET efficiency calculation based on simulations with dye molecules
Request a detailed protocolTo check if the simulated closed and tilted states generate the FRET signals of the respective unwinding and rezipping states measured by the single-molecule experiments, we carried out equilibrium simulations with the actual dye molecules for both states. AlexaFluor555 and AlexaFluor647 maleimides (Molecular Probes, Eugene, OR) were modeled according to (Vrljic et al., 2010) and (Gust et al., 2014) (see Figure 5—figure supplement 2). Then the two dyes were, respectively, attached to UvrD residues 473 and 100, which were mutated to cysteine from alanine. Force field parameters for the dyes linked to a cysteine residue were obtained from the CHARMM General Force Field (CGenFF) (Vanommeslaeghe et al., 2010) using the ParamChem server. The total charges were set to 0 and −3 for the two dyes respectively (Gust et al., 2014). Partial charges on the atoms were further refined by the Force Field Toolkit (ffTK) (Mayne et al., 2013) in VMD (Humphrey et al., 1996). Parameters for bonds, angles and dihedrals from CGenFF with high penalty scores were validated or refined by ffTK.
To sample dye dynamics efficiently, we launched 50 independent standard MD simulations with random initial velocity seeds for the closed, tilted and apo states. Every single simulation lasted 10 ns and a total 500 ns simulation time was accumulated for each state.
The FRET efficiency was determined by , where is the distance between the donor and acceptor, and is the Föster radius (or the 50% energy transfer distance). is given by the relationship (Wu and Brand, 1994) , where is the index of refraction, is the donor quantum yield, is the spectral overlap integral, and is the orientation factor. is determined to be 51 Å when equals 2/3, assuming that the dyes randomize their orientations by rapid diffusion prior to energy transfer. Such an assumption can be problematic, and in the present study the orientation factor is calculated using , where is the angle between the donor and acceptor transition dipole moments and and are the angles between these two dipoles and the vector connecting the donor and acceptor (Corry and Jayatilaka, 2008). The transition dipole moments for AlexaFluor555/647 or very similar dyes have been determined in (Corry and Jayatilaka, 2008) and (Graen, 2009). The simulated FRET data were integrated to four ns per point to obtain its probability distribution using the density kernel estimation method (Parzen, 1962; Comstock et al., 2015).
Analysis of single-molecule data
Request a detailed protocolTo validate our simulation results, analysis based on the raw data from single-molecule experiments was carried out. Comstock et al. (2015) combined optical tweezers (to detect UvrD unwinding activity) and single-molecule FRET (to detect UvrD conformation) measuring both simultaneously. Example raw time traces of UvrD activity and conformation are shown in Fig. 3 and Fig. S5 in Comstock et al. (2015) (at 10 M ATP concentration). Time traces from the optical tweezers were sampled at 267 Hz. Time traces for donor and acceptor intensities were integrated to 30–60 ms per data point. The time-dependent FRET efficiency was calculated by (Ha et al., 1999; Choi et al., 2010), where and are the measured donor and acceptor intensities, and is a correction factor accounting for the different detection efficiencies for the two dyes, and can be measured from photobleaching events. is determined to be 0.78 from 20 acceptor photobleaching events, where and are the acceptor and donorintensity changes upon acceptor photobleaching, respectively.
To measure the FRET efficiency distribution for the unwinding and rezipping states individually, we needed to assign each raw data point to the two states separately. Since the helicase velocity and FRET efficiency were measured concurrently, we used helicase velocity to define whether each data point in the traces belonged to unwinding or rezipping states (see Fig. S5 in [Comstock et al., 2015]). Time intervals were determined during which the helicase was either unwinding, rezipping, or paused (positive velocity indicates UvrD is in the unwinding state; negative velocity indicates rezipping; absolute unwinding velocity smaller than 20 bp/s indicates a pause). Paused states were not considered in the analysis. FRET efficiencies over each time interval were collected for the unwinding and rezipping states from 141 time intervals (13 molecules in total). We then used the density kernel estimation method to obtain the experimental FRET distribution (Parzen, 1962; Comstock et al., 2015). A density kernel plot is a summation of small Gaussians centered at each FRET data point. We used a standard deviation of 0.06 for the Gaussians.
We also analyzed the dwell times for both the unwinding state (high FRET) and the rezipping state (low FRET) of UvrD monomers. For this purpose, the duration of each time interval defined above was measured using the traces from optical tweezers measurements. We chose to select intervals and calculate the dwell time using the tweezers signal because it has a higher time resolution than the FRET signal (about one order of magnitude higher). The dwell time distribution was obtained by histogramming the collected duration values for the unwinding and rezipping state separately. In order to assess the effect of ATP concentration on UvrD functional switching, we analyzed optical tweezers data of UvrD activity at two different ATP concentrations. Figure 8c plots the distributions of dwell times at 10 M and 2.5 M ATP concentration. The rates of the transitions were estimated by the inverse of the averaged dwell times.
Modeling of the rezipping-state structure
Request a detailed protocolTo construct a structure of the rezipping state (after ssDNA strand switching) starting from the tilted conformation (Figure 8d), we consider the following constraints: (1) the 2B domain maintains contact with the dsDNA while the ssDNA binding domains (1A and 2A) disassociate from one ssDNA strand and bind to the other ssDNA strand. Otherwise, the entire protein would dissociate from DNA. (2) The interaction configuration between ssDNA and the motor domains (1A-2A) must remain the same after strand switching. The motor domains move from 3’ to 5’ on the ssDNA in both the unwinding and rezipping modes. With these considerations, we created a structural model in which we repositioned the UvrD-bound ssDNA segment from the 3’ ssDNA tail (strand A) to the 5’ ssDNA tail (strand B) of the junction in the tilted state. This was achieved by attaching the 5’ terminus of strand B to the 3’ terminus of strand A, and by cutting the ssDNA (strand A) at the junction position. We then equilibrated the modeled system in a water box with 55 mM NaCl for 100 ns.
Data availability
The PDB file of our predicted structure (tilted state) has been uploaded as Supplementary File 1. CHARMM36 force field parameters for the fluorophore molecules used in the simulations have been uploaded as Supplementary file 2.
References
-
Gapped BLAST and PSI-BLAST: a new generation of protein database search programsNucleic Acids Research 25:3389–3402.https://doi.org/10.1093/nar/25.17.3389
-
Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral anglesJournal of Chemical Theory and Computation 8:3257–3273.https://doi.org/10.1021/ct300400x
-
DNA helicases involved in DNA repair and their roles in cancerNature Reviews Cancer 13:542–558.https://doi.org/10.1038/nrc3560
-
Single-molecule FRET-derived model of the synaptotagmin 1-SNARE fusion complexNature Structural & Molecular Biology 17:318–324.https://doi.org/10.1038/nsmb.1763
-
Motoring along with the bacterial RecA proteinNature Reviews Molecular Cell Biology 8:127–138.https://doi.org/10.1038/nrm2099
-
Mechanochemical energy transduction during the main rotary step in the synthesis cycle of F1-ATPaseJournal of the American Chemical Society 139:4025–4034.https://doi.org/10.1021/jacs.6b11708
-
Superfamily I helicases as modular components of DNA-processing machinesBiochemical Society Transactions 39:413–423.https://doi.org/10.1042/BST0390413
-
A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactionsThe Journal of Physical Chemistry B 109:6676–6687.https://doi.org/10.1021/jp045571i
-
Mechanism of ATP-dependent translocation of E.coli UvrD monomers along single-stranded DNAJournal of Molecular Biology 344:1287–1309.https://doi.org/10.1016/j.jmb.2004.10.005
-
Large-amplitude nonlinear motions in proteinsPhysical Review Letters 68:2696–2699.https://doi.org/10.1103/PhysRevLett.68.2696
-
ThesisDye’ing for FRET: The HTH motif in smFRET-MD simulationGeorg-August-Universität, Göttingen.
-
Optimization of the CHARMM additive force field for DNA: Improved treatment of the BI/BII conformational equilibriumJournal of Chemical Theory and Computation 8:348–362.https://doi.org/10.1021/ct200723y
-
Free energy calculations by the molecular mechanics poisson-boltzmann surface area methodMolecular Informatics 31:114–122.https://doi.org/10.1002/minf.201100135
-
g_wham—A free weighted histogram analysis implementation including robust error and autocorrelation estimatesJournal of Chemical Theory and Computation 6:3713–3720.https://doi.org/10.1021/ct100494z
-
VMD: visual molecular dynamicsJournal of Molecular Graphics 14:33–38.https://doi.org/10.1016/0263-7855(96)00018-5
-
DNA mismatch repair: functions and mechanismsChemical Reviews 106:302–323.https://doi.org/10.1021/cr0404794
-
Rotations of the 2B sub-domain of E. coli UvrD helicase/translocase coupled to nucleotide and DNA bindingJournal of Molecular Biology 411:633–648.https://doi.org/10.1016/j.jmb.2011.06.019
-
Calculation of Free energy landscape in multi-dimensions with hamiltonian-exchange umbrella sampling on petascale supercomputerJournal of Chemical Theory and Computation 8:4672–4680.https://doi.org/10.1021/ct300468g
-
Generalized scalable multiple copy algorithms for molecular dynamics simulations in NAMDComputer Physics Communications 185:908–916.https://doi.org/10.1016/j.cpc.2013.12.014
-
Prokaryotic nucleotide excision repairCold Spring Harbor Perspectives in Biology 5:a012591.https://doi.org/10.1101/cshperspect.a012591
-
Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum modelsAccounts of Chemical Research 33:889–897.https://doi.org/10.1021/ar000033j
-
THE weighted histogram analysis method for free-energy calculations on biomolecules. I. The methodJournal of Computational Chemistry 13:1011–1021.https://doi.org/10.1002/jcc.540130812
-
Limiting law excess sum rule for polyelectrolytesPhysical Review E 88:052315.https://doi.org/10.1103/PhysRevE.88.052315
-
Segmented transition pathway of the signaling protein nitrogen regulatory protein CJournal of Molecular Biology 392:823–836.https://doi.org/10.1016/j.jmb.2009.06.065
-
Non-hexameric DNA helicases and translocases: mechanisms and regulationNature Reviews Molecular Cell Biology 9:391–401.https://doi.org/10.1038/nrm2394
-
Ligand-induced conformational changes in ras p21: a normal mode and energy minimization analysisJournal of Molecular Biology 274:114–131.https://doi.org/10.1006/jmbi.1997.1313
-
Mechanism of substrate translocation by a ring-shaped ATPase motor at millisecond resolutionJournal of the American Chemical Society 137:3031–3040.https://doi.org/10.1021/ja512605w
-
Escherichia coli helicase II (urvD gene product) translocates unidirectionally in a 3' to 5' directionThe Journal of Biological Chemistry 261:10169–10175.
-
Rapid parameterization of small molecules using the Force Field ToolkitJournal of Computational Chemistry 34:2757–2770.https://doi.org/10.1002/jcc.23422
-
Escherichia coli DNA helicase II is active as a monomerJournal of Biological Chemistry 274:12488–12498.https://doi.org/10.1074/jbc.274.18.12488
-
The UvrD303 hyper-helicase exhibits increased processivityJournal of Biological Chemistry 289:17100–17110.https://doi.org/10.1074/jbc.M114.565309
-
Repetitive shuttling of a motor protein on DNANature 437:1321–1325.https://doi.org/10.1038/nature04049
-
Bacillus stearothermophilus PcrA monomer is a single-stranded DNA translocase but not a processive helicase in vitroJournal of Biological Chemistry 282:27076–27085.https://doi.org/10.1074/jbc.M704399200
-
Transmembrane helix assembly by window exchange umbrella samplingPhysical Review Letters 108:108102.https://doi.org/10.1103/PhysRevLett.108.108102
-
On estimation of a probability density function and modeThe Annals of Mathematical Statistics 33:1065–1076.https://doi.org/10.1214/aoms/1177704472
-
Mechanisms of helicasesJournal of Biological Chemistry 281:18265–18268.https://doi.org/10.1074/jbc.R600008200
-
Scalable molecular dynamics with NAMDJournal of Computational Chemistry 26:1781–1802.https://doi.org/10.1002/jcc.20289
-
The structural dynamics of macromolecular processesCurrent Opinion in Cell Biology 21:97–108.https://doi.org/10.1016/j.ceb.2009.01.022
-
DNA excision repairAnnual Review of Biochemistry 65:43–81.https://doi.org/10.1146/annurev.bi.65.070196.000355
-
Targeted molecular dynamics: a new approach for searching pathways of conformational transitionsJournal of Molecular Graphics 12:84–89.https://doi.org/10.1016/0263-7855(94)80072-3
-
Structure and mechanism of helicases and nucleic acid translocasesAnnual Review of Biochemistry 76:23–50.https://doi.org/10.1146/annurev.biochem.76.052305.115300
-
Mismatch repair during homologous and homeologous recombinationCold Spring Harbor Perspectives in Biology 7:a022657.https://doi.org/10.1101/cshperspect.a022657
-
BookHomologous Recombination by the RecBCD and RecF PathwaysIn: Higgins N. P, editors. The Bacterial Chromosome. Washington, DC, United States: American Society of Microbiology Press. pp. 389–403.https://doi.org/10.1128/9781555817640.ch21
-
An introduction to bootstrap methodsSociological Methods & Research 18:243–291.https://doi.org/10.1177/0049124189018002003
-
Multidimensional replica-exchange method for free-energy calculationsThe Journal of Chemical Physics 113:6042–6051.https://doi.org/10.1063/1.1308516
-
DNA helicases, genomic instability, and human genetic diseaseAnnual Review of Genomics and Human Genetics 1:409–459.https://doi.org/10.1146/annurev.genom.1.1.409
-
CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fieldsJournal of Computational Chemistry 31:671–690.https://doi.org/10.1002/jcc.21367
-
Molecular mechanism of the synaptotagmin-SNARE interaction in Ca2+-triggered vesicle fusionNature Structural & Molecular Biology 17:325–331.https://doi.org/10.1038/nsmb.1764
-
Overview: what are helicases?Advances in Experimental Medicine and Biology 767:1-16.https://doi.org/10.1007/978-1-4614-5037-5_1
-
Resonance energy transfer: methods and applicationsAnalytical Biochemistry 218:1–13.https://doi.org/10.1006/abio.1994.1134
-
Lessons learned from UvrD helicase: mechanism for directional movementAnnual Review of Biophysics 39:367–385.https://doi.org/10.1146/annurev.biophys.093008.131415
Article and author information
Author details
Funding
National Institute of General Medical Sciences (9P41GM104601)
- Zaida Luthey-Schulten
- Klaus Schulten
National Science Foundation (PHY-1430124)
- Yann R Chemla
- Zaida Luthey-Schulten
- Klaus Schulten
National Institute of General Medical Sciences (R01 GM120353)
- Yann R Chemla
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The present contribution is dedicated to Klaus Schulten (1947–2016), whose visionary developments in high performance simulation tools for molecular dynamics permit characterizing biomolecular motor action at biologically relevant timescales. This work has been supported by grants from the NIH (9P41GM104601 to ZL-S and KS; R01 GM120353 to YRC) and from the NSF PHY-1430124 (Center for the Physics of Living Cells to YRC, ZL-S and KS). The authors gladly acknowledge supercomputer time provided by the Texas Advanced Computing Center via Extreme Science and Engineering Discovery Environment grant NSF-MCA93S028 and the Blue Waters sustained-petascale computing project, which is supported by NSF (OCI-0725070 and ACI-1238993) and the State of Illinois. We also thank former Chemla lab member Dr. Matthew J Comstock for scientific discussions.
Copyright
© 2018, Ma et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 4,006
- views
-
- 407
- downloads
-
- 15
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Computational and Systems Biology
- Microbiology and Infectious Disease
Variant calling is fundamental in bacterial genomics, underpinning the identification of disease transmission clusters, the construction of phylogenetic trees, and antimicrobial resistance detection. This study presents a comprehensive benchmarking of variant calling accuracy in bacterial genomes using Oxford Nanopore Technologies (ONT) sequencing data. We evaluated three ONT basecalling models and both simplex (single-strand) and duplex (dual-strand) read types across 14 diverse bacterial species. Our findings reveal that deep learning-based variant callers, particularly Clair3 and DeepVariant, significantly outperform traditional methods and even exceed the accuracy of Illumina sequencing, especially when applied to ONT’s super-high accuracy model. ONT’s superior performance is attributed to its ability to overcome Illumina’s errors, which often arise from difficulties in aligning reads in repetitive and variant-dense genomic regions. Moreover, the use of high-performing variant callers with ONT’s super-high accuracy data mitigates ONT’s traditional errors in homopolymers. We also investigated the impact of read depth on variant calling, demonstrating that 10× depth of ONT super-accuracy data can achieve precision and recall comparable to, or better than, full-depth Illumina sequencing. These results underscore the potential of ONT sequencing, combined with advanced variant calling algorithms, to replace traditional short-read sequencing methods in bacterial genomics, particularly in resource-limited settings.
-
- Cancer Biology
- Computational and Systems Biology
Assay for Transposase-Accessible Chromatin sequencing (ATAC-Seq) is a widely used technique to explore gene regulatory mechanisms. For most ATAC-Seq data from healthy and diseased tissues such as tumors, chromatin accessibility measurement represents a mixed signal from multiple cell types. In this work, we derive reliable chromatin accessibility marker peaks and reference profiles for most non-malignant cell types frequently observed in the microenvironment of human tumors. We then integrate these data into the EPIC deconvolution framework (Racle et al., 2017) to quantify cell-type heterogeneity in bulk ATAC-Seq data. Our EPIC-ATAC tool accurately predicts non-malignant and malignant cell fractions in tumor samples. When applied to a human breast cancer cohort, EPIC-ATAC accurately infers the immune contexture of the main breast cancer subtypes.