RNA structure calculation applying topographic restraint
Given that duplexes are well-conserved and they are the predominant building blocks in folded RNA structures by far36, they can be considered semi-rigid bodies within a folded RNA structure. Since they are covalently connected, these duplexes can be treated as kinematic chains. Adding kinematic constraints between rigid bodies will significantly decrease the degrees of freedom of a rigid body system57, and imposing the topographic constraints in addition to the kinematic constraints further reduces the degrees of freedom of sampling space.
A high-resolution AFM image is more than just a ‘frame’ of a molecule. The width and pitch of an A-form RNA duplex are ~25 and 30 Å, respectively, which are on a similar scale to a sharp AFM probe and sensitive to detection. Thus, given an achievable imaging resolution of 10–15 Å34 (Fig. 2a,b), major structural features such as grooves and pitches of long duplexes, along with molecular shapes and topological folds of larger structured RNAs in solution, are discernable in high-resolution AFM images58,59,60,61. Thus, a high-resolution AFM image of a molecule is a 3D frame with details about topographic information on individual molecules. The explicit expression of the physical relationship between a molecular structure and the topographic molecular surface is defined. A detailed description of the implementation is provided in Supplementary Information.
Unsupervised machine learning
Our UML approach assumes that the classical molecular dynamics simulation guided by topographic information can sample the real native conformational space of the RNA, and that the correct models can be identified based on the established hierarchical folding principle62, energetics63 and agreement with topographic restraints. Our UML algorithm is able to decipher the underlying correlation of the dataset, resulting in the recognition of generalizable models without pre-training or data labelling. Each analysed dataset (trajectory) is unique, and the machine does not have any expected pre-labelled output from a given input. Our UML algorithm consists of three main steps: (1) energy filtering; (2) PCA and clustering; and (3) cohort model selection. A detailed description of the UML is provided in Supplementary Information.
Supervised DNN
Based on the question of whether the most fundamental characteristics of models such as their energetics and known topology of a structure contained in the AFM experimental data would provide enough information to consistently determine the r.m.s.d. between the structural model and an unknown ground-truth structure, we designed a DNN64 to explore how these fundamental characteristics could be intrinsically correlated. A detailed description of the supervised DNN is provided in Supplementary Information.
Underfitting and overfitting
To avoid overfitting and to be able to keep increasing the complexity of our ANN, we added regularization penalties to the training. Within the known regularizers, we evaluated training using ridge regression (L2 regularizer) and the dropout technique. Ridge regression adds a penalty to the loss function term for all the weights squared, preventing the weights from assuming excessively high values. In the dropout technique, in each step of the training optimization, some neurons have a given (set) chance to be turned off. We also tested increasing the size and variety of the dataset by adding more data (trajectories) (Supplementary Table 3).
Optimized architecture
To train the DNN and assess its performance, we split dataset BM0 (which contains only one kind of RNA) into two parts: the training, and the initial training test set (called ‘training validation’), where we could check the regularization effect over the same trajectory and assure that the regularization was blocking the train from overfitting the trajectory over the split, thus providing similar loss on both sets. The training set had 80% of the 3.5 million trajectory models, while the training-validation set had 20%. The optimized training dataset that yielded the best performance was built using all data from the BM0, with an additional 5% of data from each of the trajectories of BM1 and BM2 (Extended Data Fig. 5).
The validation set was created by using a different RNA trajectory simulation, BM5, so that the best loss on the validation set would point to the place where the training and learning with a given RNA was still generalized to another RNA trajectory, applying early stopping on the evaluated loss considering r.m.s.d. values up to 10 Å to weight a better performance on smaller r.m.s.d. values than larger ones. Hence, the validation set was used for both tuning the hyperparameters and for selecting the best-trained model, while further tests over the benchmarks address if our model can generalize its findings and learnings to other RNAs not contained in the training data, with different RNA sizes and folds, assessing what would be the real performance of our model to other unknown RNAs and trajectories than the one used for training.
We optimized the architecture for this work by many step-by-step random searches and subsequent fine-tuning of the hyperparameters, which include the number of layers, the number of neurons per layer, weight initialization, neuron activations, regularization penalties and types, the optimizer algorithm as well as the learning rate. Additionally, more than 50 different compositions (data, kappa and noise) of the training dataset were used for training the models (Supplementary Table 3). The number of hidden layers tested (also by a random search) was between one and ten hidden layers. The number of neurons in each layer, on the other hand, was tested basically in 3 types: (1) starting with a high number of neurons in the first layer and decreasing this as the number of layers increases; (2) starting with a medium number of neurons in the first layer, and increasing the number of neurons on the next layers until reaching the middle layer, then decreasing as we continue to the last layer; and (3) through a random search, where the number of neurons per layer was picked randomly as a multiple of 8, being able to assume values from 8 to 256 neurons per layer. For architectures with five or more layers, we included batch normalization within layers.
The non-linear activations tested were relu, leaky-relu, elu and gelu for each layer separately, or a selu65 activation set for all layers. For regularization, each layer could use either the ridge regression and/or a dropout66 chance (for selu the alpha dropout65 was used instead of Dropout to keep the self-normalizing properties). Our optimized architecture has only 3 hidden layers with a decreasing number of neurons, 128 in the first layer, 64 in the second, and 16 in the third, using elu activation with a common dropout rate of around 20% as the regularizing agent. Deeper networks also had a good performance, but with the cost of many weights to train without clear improvement. The total number of trainable parameters with the current architecture is around 11k. Within initializations, we tested Glorot uniform, Lecun normal and He normal, with the latter achieving the best performance as the weight initializer and using Adam as the optimizer algorithm with a standard learning rate of 0.001, with the mini-batch size of 128 and using Huber loss.
The models were trained using NIH-HPC (Biowulf) k80/k100x nodes: K100x node: 36 × 2.3 GHz (Intel Gold 6140), hyperthreading, 25 MB secondary cache, 4 x NVIDIA V100-SXM2 GPUs (32 GB VRAM, 5120 cores, 640 Tensor cores); K80 node: 28 × 2.4 GHz (Intel E5-2680v4), hyperthreading, 35 MB secondary cache, 2 x NVIDIA K80 GPUs with 2 x GK210 GPUs each (24 GB VRAM, 4992 cores).
RNA sample preparation
RPR
The RPR was prepared as described33. In brief, the RPR was transcribed in vitro with recombinant T7 phage RNA polymerase from a double-strand DNA template that was amplified by PCR from linearized DNA plasmid, which encodes a full-length RPR from B. stearothermophilus with an upstream T7 RNA polymerase promoter. Transcribed RNA was purified by denaturing polyacrylamide gel electrophoresis containing tris-borate with EDTA (TBE) and 8 M urea. The RNA was excised and eluted from the gel in RNA elution buffer (300 mM Sodium acetate pH 5.3, 0.1 mM EDTA) for 12 h at 4 °C. The eluted RNA was filtered using a 0.2-μm Ultrafree-MC centrifugal filter device (Millipore). Purified RNA was subjected to several buffer exchanges using a Centricon unit (Millipore) with 30 kDa molecular weight cut-off membrane against refolding buffer (50 mM MES buffer pH 6.8, 100 mM KCl, 1 mM MgCl2), then concentrated to 2 μM, aliquoted, and stored at −80 °C before utilization.
For AFM experiments, the RNA sample at 2 μM concentration was annealed in the refolding buffer (50 mM MES buffer pH 6.8, 100 mM KCl, 10 mM MgCl2) at 65 °C for 2 min followed by stepwise cooling to 37 °C over 30 min, and then kept at 4 °C before AFM measurements. To dilute the RNA sample to the required concentration (20 nM) for AFM, 1:100 volume of low-salt buffer (50 mM MES buffer pH 6.8, 10 mM KCl, 1 mM MgCl2 (preequilibrated at 4 °C) was used, and the diluted sample was immediately deposited onto mica pre-treated with 1-(3-aminopropyl) silatrane (APS) for immobilization26. The functionalization of mica with APS is widely used for the nondisruptive immobilization of nucleic acids primarily through the electrostatic interactions between protonated amino groups of the APS-mica substrate and the negatively charged nucleic acid backbone.
RRE RNA
RRE sample was prepared following the same protocol described previously in detail52. The fresh sample was used for the AFM experiments with a concentration of 20 nM in 50 mM MES buffer pH 6.8, 10 mM KCl, 1 mM MgCl2. The sample was loaded on a freshly cleaved mica pre-treated with APS and incubated for 30 min before imaging.
Peptide design, synthesis of P46 and modelling
The two ARMs (H2N-RRRDRRLRQRARRRAAAA-COOH) are joined by the amino groups of a lysine main and side chains via chemical synthesis (Shengnuo). This compound is patented under US Patent Number 10,464,970.
The monomeric ARM structural model was built using Pep-Fold67. Then, two ARM structural models were linked in parallel using the bond build function of PyMol (PyMol Molecular Graphics System, version 2.0 Schrödinger). A ~3.0-μs coarse-grained molecular dynamics simulation using CafeMol31,32 was performed to obtain the distance distribution between the 2 ARMs. The molecular dynamics trajectory was generated applying constant temperature simulation of 300 K Langevin dynamics and Go model for a total of 60 × 106 steps.
Electrophoretic mobility shift assay
RRE (1 μM) was mixed with various ratios of P46 and Rev protein in a buffer containing 10 mM HEPES (pH7.5), 300 mM KCl, 1 mM MgCl2, 0.5 mM EDTA, 0.1 μg μl−1 BSA, 0.2 μl SUPERase•In RNase Inhibitor (Thermo Fisher Scientific). The total reaction volume was 10 μl. The reactions were incubated at room temperature for 30 min, then 4 μl of each reaction was loaded into a Novex 6% TBE gel (Thermo Fisher Scientific). The gel was run for 80 min at 120 V, and the image was taken using a Gel Doc EZ Imager (Bio-Rad) after staining with SYBR Gold Nucleic Acid Gel Stain (Thermo Fisher Scientific). Adenine riboswitch (RibA71) RNA (150 μM), consisting of three helices, was used for competitive non-specific binding via peptide–major groove interactions. Samples of the Rev protein and RibA71 were prepared as reported previously53,68. The uncropped gel image file is provided in Supplementary Fig. 1.
AFM experiments and image processing
Experimental AFM image acquisition
The detailed procedure for the AFM image acquisition is described elsewhere2. Full-length RPR particle images, P1, P2 and P3 (Extended Data Fig. 1), were recorded under the solution conditions described above using a Cypher VRS AFM (Asylum Research, Oxford Instrument) at 4 °C with amplitude-modulated AC mode at a scan rate of 1 Hz (commonly known as tapping mode) using FASTSCAN-D-SS probes (Bruker). For RNA immobilization, 50 mM APS stock was freshly diluted 300-fold in deionized water right before use and then used to coat a freshly cleaved muscovite mica (Grade V1) (Ted Pella) and incubated for 30 min, followed by rinsing the mica surface with deionized water and drying gently with filtered nitrogen gas.
AFM noise estimation
For quantification of the noise present in the z coordinates, we used the cropped single molecule from the full recorded AFM topography as input, and the z values were collected for defined x and y coordinates of the ‘empty’ area around the molecule (Extended Data Fig. 2a–d). The z-coordinate values of the empty horizontal and vertical spaces can be described by a normal function, where the mean value of this distribution represents the mean noise and the uncertainty as the standard deviation (sigma). The mean noise value and uncertainty were evaluated for P1, P2 and P3 before and after image processing (Extended Data Fig. 2a–d). In this analysis, we are considering all the experimental sources that result in noise randomly distributed over all recorded data as a background signal.
AFM image resolution estimation
The topography resolution assessment was performed using an ACV approach34. There are two principal steps to be performed in this method. First, using the processed image, we calculate the 2D FFT of the AFM image and a defined ring-size (in pixels) cut-off filter is applied to select a portion of the image in Fourier space (Extended Data Fig. 2e–g). Afterwards, the image is back-calculated to real space; this step is described as a low-pass filtered Fourier ring. In the second step, we calculate the ACV between the original image (R) and the resulting one from the inverse fast Fourier transform for each of the low-pass filtered rings (R′). The comparison between the original image with its resulting image from the low-pass filter is performed using the auto-correlation equation (Supplementary Information). A loop interaction was applied starting from low to high frequency, where the ACV starts from low correlation values up to values near to 1.0 where the low-pass cut-off is close to the particle dimension in real space. In Fig. 2, we demonstrate some intermediate steps of the Fourier ring filter applied to P1, P2 and P3 particles. The formula for estimation of ACV value is provided in Supplementary Information.
AFM image processing
The detailed procedure for AFM image processing is described elsewhere2. In brief, raw images were first processed using SPIP (Scanning Probe Image Processor) software (https://www.imagemet.com/products_/spip/): plane levelling to the particle-free region by applying third-order polynomial, followed by spike filtering to remove artefact streaks, and fast Fourier transform to remove high-frequency noise (Extended Data Fig. 2). The final image resolution was increased to 4,096 × 4,096 pixels by doubling the number of pixels twice. Single-particle images were cropped from the processed images and converted to pseudo-AFM with a digital resolution of 5-Å per pixel in MountainsSPIP software (https://www.imagemet.com/products_/spip/) for structure calculation.
Benchmark information and design
BM0
BM0 was designed to provide the bulk of the training data by using a representative RNA with a known structure at an acceptable resolution. For this purpose, we chose the crystal structure of the RPR catalytic domain (PDB: 3DHS)40. The residues that were missing in 3DHS due to insufficient electron density were modelled using SimRNA69 and further refined using Coot70. This model, representing the ground-truth structure for BM0, was subjected to coarse-grained dynamic fitting in CafeMol31,32 to an experimental AFM image of this RNA, and the trajectory of models was scored using ARES28. The model with the best ARES score of 9.9 from this pool of models, named k158597 (Extended Data Fig. 6a,b), was then used as the initial model for BM0. The r.m.s.d. between k158597 and the ground-truth structure is 21.4 Å. AFM images of the ground-truth model were calculated using a resolution of 5.0 Å per pixel, with 7 different simulated Gaussian noise levels added—that is, 5, 10, 15, 20, 30, 40 and 50% of the maximum z height (Fig. 2f). The dynamic fitting using k158597 as the initial model and the AFM topography of the ground-truth structure was performed for a total of 20 × 106 steps (~0.9 µs) for all noise levels (Fig. 2h).
BM1 and BM2
BM1 and BM2 were designed to tackle cases in which only the primary sequence and secondary structure information may be known, and the starting model must be predicted. For this task, we first applied the FARFAR246 – rna_denovo application, generating 10000 structure models of RPR catalytic domain, using the primary sequence, secondary structure, and atom pair distance constraints of the well-known loop interaction L15.1–L5.1 described in detail previously71 (Supplementary Table 10). For structure refinement the minimize_rna function was applied as a potential during the FARFAR2 structure prediction run, using parallel jobs on a 28-core 2.3 GHz x2695 processor.
The FARFAR2 scoring function was calculated for all the predicted models and analysed as a function of the main energy terms. The sampled refined structures show a range of r.m.s.d. with a maximum of 46 Å and a minimum of 14 Å from the crystal model (PDB: 3DHS). We selected three models from all predicted structures from FARFAR2 using the following criteria: one model with the best ARES score (S257), one being located in the region of both ARES and FARFAR2 lowest scores and an r.m.s.d. from the ground-truth structure of at least 20 Å (S1076), and one model with the lowest r.m.s.d. (S142), (Extended Data Fig. 6).
ARES selected model S257 as the best model from the FARFAR2 ensemble, a model that presents dramatically different folds from the crystal model, with an r.m.s.d. of ~30 Å (Extended Data Fig. 6d). Using an r.m.s.d. threshold of 20 Å and scoring the models using the energetic scoring function of FARFAR2 and final score of ARES, the best model was S1076 with an r.m.s.d. of 22.0 Å, where this model shows a folding similar to the crystal model (Extended Data Fig. 6e).
BM3
We applied our method to the adenosylcobalamin riboswitch (Cbl) using the crystal structure (PDB: 4GMA), which has a folding and size (210 nt) different from the RNA used in training (BM0–B2). The structure calculation was performed with a total of 20 × 106 steps (~0.9 µs) using an AFM topography generated with 5 Å per pixel (Extended Data Fig. 8a). The final trajectory consisted of ~6.6 million models, which were analysed using UML and DNN (Figs. 3, 4).
BM4 and BM5
These two benchmarks were designed to test our method using initial models determined using low-resolution experimental data. In this case, we used the topological structures of RPR (298 nt) and group II intron (387 nt), generated by RS3D48,72 using secondary structure information and SAXS data simulated from their respective crystal structures, PDB: 2A6433 and PDB: 4E8K49. The ground-truth AFM images were calculated using the crystal models for the RNAs with a resolution of 5 Å per pixel (Extended Data Fig. 8b,c). The structure determination in each case was performed using a trajectory with a total of 60 × 106 steps (~2.7 µs). The final trajectories consisted of ~13.4 million models, which were analysed using UML and DNN (Figs. 3 and 4).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
