Designing molecular RNA switches with Restricted Boltzmann machines

Domain: physics
Published: 2025-12-18
Editorial outcome: major_revision

Abstract
Riboswitches are structured allosteric RNA molecules that change conformation upon metabolite binding, triggering a regulatory response. Here we focus on the de novo design of riboswitch-like aptamers, the core part of the riboswitch undergoing structural changes. We use Restricted Boltzmann machines (RBM) to learn generative models from homologous sequence data. We first verify, on four different riboswitch families, that RBM-generated sequences correctly capture the conservation, covariation and diversity of natural aptamers. The RBM model is then used to design new SAM-I riboswitch aptamers. To experimentally validate the properties of the structural switch in designed molecules, we resort to chemical probing (SHAPE and DMS), and develop a tailored analysis pipeline adequate for high-throughput tests of diverse sequences. We probe a total of 476 RBM-designed and 201 natural sequences. Designed molecules with high RBM scores, with 20% to 40% divergence from any natural sequence, display ≈30% success rate of responding to SAM with a structural switch similar to their natural counterparts. We show how the capability of the designed molecules to switch conformation is connected to fine energetic features of their structural components.

Introduction
Riboswitches are regulatory RNA elements found mostly in bacterial and in some eukaryotic messenger RNAs. Usually located upstream of coding sequences, they modulate the expression of the downstream gene at the transcriptional or translation level in the presence of a specific metabolite1–4; some riboswitches placed within genes even regulate alternative splicing5. In order to perform their function, these RNA motifs switch between two stable conformations in response to binding of their cognate metabolite to the aptamer domain of the riboswitch (Fig. 1). This change of conformation, in turn, affects the expression platform, where the regulation signals are located. Understanding how the aptamer domain by itself is able to implement a structural switch in response to the ligand, and how this is encoded in the sequence, is an important step towards the characterization of the full riboswitch regulation.

The sequence-to-function mapping of structured RNAs is a complex problem. In the course of evolution, sequence patterns necessary for function are conserved, suggesting that large sequence datasets can shed light on this mapping. Comparative analysis of homologous RNA sequences collected in Multiple Sequence Alignments (MSA)6 have been successful to predict secondary RNA structures, tertiary structural motifs, and even the entire three dimensional architecture of complex RNA7–12. Covariation analysis has also been used to predict pseudoknots and other tertiary contacts from statistical couplings inferred from conservation and covariation across the MSA columns13,14, or by including positive and negative evolutionary information such as in the Cascade covariation Folding Algorithm (CacoFold)15. Machine learning approaches have recently shown promising results in RNA structure prediction. Among them Rosetta FARFAR216 uses Monte-Carlo-based fragment assembly methods and can be aided by geometric deep learning approaches such as ARES17 to score putative structures. DeepFoldRNA18 significantly outperformed the state-of-the-art tertiary structure prediction from sequence only. Although these approaches look promising, AlphaFold-level accuracies19 (for proteins) are not yet reached in RNA structure prediction20,21.

The mirroring problem of designing RNA sequences capable of folding in a particular target structure or of performing a desired function has also long been investigated. One successful approach is based on directed evolution (SELEX). RNA sequences are selected from an initial random library to optimize a target function, such as the switching dynamics for bistable aptamers22. Models trained on such data are capable of classifying sequences according to their functionality and of extracting key sequence-features for the desired function23–27. Classifiers have been used downstream of random mutagenesis to filter out good sequences, but this approach only works if the libraries already contain good candidates. In parallel, much effort has been devoted to the rational design of secondary structures, in particular with minimum free energy approaches28,29. However, due to algorithmic complexity9,30, those approaches often ignore pseudoknots and other tertiary contacts known to be essential for the function of some RNAs, such as riboswitches or ribozymes.

To date, building generative models effective in designing RNA sequences with tertiary structural targets remains a challenging problem. From this point of view, riboswitches, in addition to their fundamental interest in biology and relevance for the RNA world hypothesis5, offer a difficult design problem, as their sequences encode not only two conformational structures but also a metabolite-mediated switching mechanism between them. In the present work, we address this challenging issue and show how to design functional RNA switches (albeit devoid of expression platform) from natural sequence data.

One of the largest identified groups of riboswitches recognize S-adenosyl-methionine (SAM) as their effector metabolite31,32. While six different SAM binding structural motifs have been identified, this study focuses on those harboring type I SAM aptamers (SAM-I)33. Figure 1A shows the secondary structure of the aptamer domain in absence of SAM, where transcription is allowed (ON state), while panel B depicts the structure when SAM is bound and transcription continuation is prevented (OFF state). Upon SAM binding, the aptamer cooperatively folds into the closed structure characterized by the stabilization of P1, three triple base pairs and a pseudoknot (red in the figure)34. The closed state of the aptamer is stabilized by direct tertiary contacts between SAM and specific nucleotides forming the SAM binding pocket31,35.

Hereafter, we employ Restricted Boltzmann machines (RBM), a two-layer generative neural network to design new SAM aptamers (Fig. 2A). RBMs have recently been shown to provide interpretable models of proteins in various contexts36–39, with application to design26,40. By learning the sequence statistics of the SAM-I riboswitch family, the RBM models the constraints that enable aptamers to adopt the correct secondary structure, form tertiary contacts and effect a conformational switch in response to SAM presence.

The RBM model was used to design 476 sequences, which we experimentally tested with SHAPE-MaP and DMS, two chemical probing methods giving information about paired and unpaired residues in the structures. Comparison of the reactivity profiles in the presence or absence of SAM allows us to assess the effectiveness of the structural switch for each tested molecule. This high-throughput analysis is made possible by the introduction of an automated Bayesian analysis of the SHAPE and DMS reactivity profiles. Our results for RBM-generated aptamers are compared to experiments on 201 natural sequences, and 58 sequences designed by RFAM Covariance Models, another generative model capturing local conservation and secondary-structure covariation only.

Results
Our pipeline is described in Fig. 2 and includes: sequence data acquisition from Rfam41, training and sampling the RBM to design artificial SAM-I aptamers, experimental characterization of SAM-induced conformational switch in natural and designed sequences by chemical probing (SHAPE42,43 and DMS44), and statistical analysis of the measured reactivities.

Generative models of SAM riboswitch aptamers
We train an RBM (Fig. 2A, B) on a MSA of natural homologs of the aptamer domain of SAM-I riboswitches, gathered from the Rfam41 database (RfamID: RF00162). RBM are energy-based generative models, that once trained, define a score, −Eeff(v), over all possible sequences v. Sequences with high scores (equivalently, low energies) are then “good” fits to the family, according to the model. Artificial sequences of high score can be generated by sampling the resulting Boltzmann measure, PRBM(v) ∝ e−Eeff(v), see Fig. 2C and “Methods” for details.

The weights between visible units, carrying the RNA sequence, and hidden units, extract latent factors of variation in the data (Fig. 2D). After marginalization over those latent variables, effective interactions between pairs of residues can be computed36, defining epistatic scores between sites (Fig. 2E and Supplementary Eq. S12 for precise definition). Pairs of sites with large epistatic scores correspond to major secondary and tertiary contacts in folded aptamers, see heatmap in Fig. 3A. Interestingly, epistatic scores at P1 are weaker than in other helices, reflecting the flexibility of P1, which is able to open or close in concert with SAM binding (Fig. 1). The pseudoknot is also correctly identified (red in Fig. 3A), proving the capability of RBM to identify tertiary motifs. Besides structural contacts, the RBM hidden units capture extended motifs, most likely relevant for tertiary structure formation and SAM binding, see weights in Fig. 3B, C.

We then evaluate the sequences designed by the RBM by comparing their scores to the ones of natural sequences and sequences designed by Covariance Models (CM). CM capture the conservation of residues along the sequence, as well as correlations due to the complementarity of base pairs in the secondary structure45, but are unable to model tertiary motifs (such as pseudoknots). As Rfam sequence alignments46 are based on CM6, our first baseline model for RF00162 was directly downloaded from Rfam (Methods) and will be referred to as Rfam CM (rCM) in the following.

In Fig. 4A, we show a scatter plot of rCM vs RBM scores for natural, RBM- and CM-generated sequences. RBM-generated sequences have rCM scores comparable to the natural ones, indicating that RBM samples satisfy the constraints imposed by the rCM model to the same extent as natural sequences. Moreover, RBM samples have RBM scores comparable to natural sequences, while rCM samples have significantly smaller scores, suggesting that the RBM impose further constraints beyond those captured by rCM, such as tertiary contacts (e.g., pseudoknot), which could be important for the aptamer function. We also check that R-scape9 supports significant covariation across pseudoknot sites for RBM samples, contrary to rCM samples as expected (see Supplementary Section E for details). In addition, RBM recapitulates several statistical properties of natural sequences in the MSA, including conservation, covariation, distribution of lengths, and distributions of Hamming distances between sequences (see Supplementary Section B and Supplementary Figs. S2 and S3).

Next, we carry out principal component analysis (PCA) of the natural MSA. The top principal component (PC) captures a mode of variation associated to deletion of the P4 helix, as can be seen from the large number of gaps in this region (Supplementary Fig. S5). Figure 4B shows the projections of the natural sequences, annotated by their taxonomic class, onto the top two PCs. The PCs appreciably separate taxonomic clusters of natural sequences. In particular, a group of Actinomycetota, in the top left corner, have very short or no P4 helix segments. SAM aptamers can function in the absence of P447, although the affinity for SAM decays with decreasing length of P448.

RBM-generated sequences also span the PC space, covering all the taxonomic clusters (Fig. 4C and Supplementary Fig. S5). In contrast, rCM-generated sequences, shown in Fig. 4D remain confined to a central region. The capability of RBM to capture complex constraints in the sequence distribution allows them to model the full variability present in homologs.

We then select a fraction of the generated sequences for experimental validation, see “Methods” for details about the selection criteria. Their PCA projections are shown in Fig. 4E, colored by their origin (Natural, rCM, RBM), and span a wide range of the natural variability.

Reactivity profiles of natural and generated aptamers with SHAPE and DMS
We resort to high-throughput chemical probing to characterize the structure of generated aptamers and their possible changes upon SAM addition. DMS mainly focuses on single-stranded A and C nucleotides, while SHAPE is sensitive to the conformational flexibility of individual nucleotides42. Generally speaking, paired nucleotides tend to show lower reactivities than residues left single stranded. Similarly, aptamer nucleotides bound to SAM are expected to be less reactive. SHAPE and DMS probing are routinely used to monitor aptamer structure, complexion with their ligand and structural rearrangment49–58.

The general result of an experiment for an aptamer is two profiles of site-dependent reactivities, one in the absence and the other in the presence of SAM (Fig. 2F). Changes in reactivities between the two conditions are expected to be informative about sites involved in interactions with SAM and in the structural switch, see Fig. 2G. However, because of the delicate nature of reactivity measurements, it is useful to benchmark the approach with natural aptamers, before turning to the analysis of the generated aptamers.

We probe a set of 208 natural sequences with SHAPE and a subset of 152 sequences with DMS in the presence or absence of SAM. These sequences are representative of Rfam ID RF00162 (Methods) and are shown by black crosses in Fig. 4E. We first present our approach and results for SHAPE-MaP. After standard processing42, we obtain the reactivity values ri,n,c assigned to each site i, for each aptamer n, and in each condition tested c (with or without SAM). We can then compute the difference in reactivities with and without SAM, Δri,n = ri,n,SAM − ri,n,no SAM. Figure 5 shows reactivity profiles from our experiments for two selected aptamers. Panel A displays the profiles obtained for yitJ aptamer from B. subtilis, for which a ligand-bound crystal structure was reported in ref. 51 (PDB id: 4KQY). Interaction with SAM is confirmed by strong reactivity changes (Fig. 5B) due to the ligand at various key sites, such as SAM contacts, and sites involved in a base-triple (Fig. 1B). The T. tengcongensis aptamer35 (PDB id: 2GIS) shows a similar behavior (Supplementary Fig. S26). In both cases, reactivity is low along the pseudoknot in absence of SAM, consistent with previous studies59 that report this element is already stable in the apo form (requiring only Mg+ for its formation). Figure 5C, D show another aptamer (from Deltaproteobacteria), where SAM response is evidenced by reactivity drops at SAM contacts, the base-triple and also the kink-turn and the pseudoknot. Our data may thus reveal the existence of variable responses to SAM across aptamers, in terms of which sites (e.g., the pseudoknot) become more protected when SAM is present or not.

The difference in reactivities with and without SAM, Δri,n, once averaged over all probed natural sequences n, to better extract functional sites at the level of the family42, defines a site-dependent Δ-reactivity template, 〈Δri〉nat, shown in Fig. 5E. We observe reactivity decreases (also called protection) for the pseudoknot (sites 25–28, 77, 79), sites involved in base triples (24, 76, 100, 73, 74) or flanking them (75), and for some of the sites directly in contact with SAM (10, 11, 46, 103). These hallmark sites, listed in Supplementary Table S2, were previously recognized for their relevance to the structural switch by previous studies using crystal structures, chemical probing, and mutagenesis experiments35,51,52, see Fig. 1B. Supplementary Section Q summarizes the literature supporting these choices.

Results for DMS probing are compatible with the above findings. We report in Supplementary Fig. S34 the reactivity profiles ri,n,c of the same natural sequences as in Fig. 5A–D obtained with DMS. The profiles are sparser due to the generally low reactivities of sites carrying G or U nucleotides. Figure 5F shows the site-dependent differential reactivity profile, 〈Δri〉nat, averaged over all 152 probed natural sequences. Contrary to its SHAPE counterpart (Fig. 5E), this differential profile vanishes on most sites along the sequence. This is expected from the fact that sites may often be occupied by G or U nucleotides (Fig. 5G) and therefore weakly sensitive to DMS probing. As a result, DMS data are often less informative about SAM-induced changes than their SHAPE counterparts. However, we also observe that the few sites on which DMS differential reactivities are non zeros show finer spatial resolution, e.g., on site i=100, and lower sequence-to-sequence variability around the average profile (gray band around the average DMS signal), see for instance site i=28 and its neighborhood. Interestingly, this latter site, which carries mostly G’s and U’s, is sensitive to DMS probing, as it is located at the junction of a stem and a loop58.

In summary, both SHAPE and DMS average differential profiles confirm that the natural sequences probed in our experiments are mostly SAM binders and, moreover, recapitulate expected structural changes upon binding. Sequences in the seed alignment (a manually curated subset41) show the same average reactivity responses (Supplementary Fig. S16).

The reactivity profiles of two representative RBM generated sequences are reported in Fig. 6A–D. Panels A, B show an example of a RBM-generated sequence for which the differential reactivity profiles are compatible with a global structural switch, as evidenced by reactivity changes (highlighted by arrows) in most of the hallmark sites (Supplementary Table S2), including sites in direct contact with SAM, but also the pseudoknot, the kink-turn and a base-triple motif that are known to be stabilized by the presence of SAM. Figure 6C, D show another RBM-generated aptamer for which the differential reactivity is localized to fewer hallmark sites. In contrast to the previous example, sites at the kink-turn and pseudoknot do not exhibit significant reactivity changes in response to SAM. Reactivity changes in the base-triple and SAM contact sites strongly suggest a ligand-binding event, and are compatible with a global structural switch from an open to a closed conformation.

We emphasize that the variety in the patterns of response to SAM seen across generated aptamers is reminiscent of what is observed in natural ones. Manual inspection of all experimentally tested 201 natural aptamers, reveals that some molecules rearrange structurally upon binding SAM, others bind without significant conformational shift, and some showing no evidence of binding (no reactivity change). Examples are shown in Supplementary Figs. S11 and S12. Global results of this manual inspection are summarized below.

We report in Fig. 6E the average differential reactivity profile of RBM-generated sequences having high scores (>300). An excellent match with the differential reactivity profile of natural sequences is observed. In particular, protections compatible with SAM binding and the expected structural switch are found at hallmark sites. We also check that these RBM-generated sequences reproduce the reactivity response to magnesium of natural sequences (Supplementary Fig. S17). In contrast, RBM sequences with lower scores (<300) show clear discrepancies (Supplementary Fig. S18) with the average profile of natural sequences.

For the sake of comparison, we show in Fig. 6F the average differential reactivities of sequences sampled from rCM (in red). Contrary to high-score RBM-generated sequences, this group of sequences shows an appreciable lack of protection at key sites, such as 10–11 (SAM contact), 25–28 (pseudoknot), 73–76 (base triples), and 103 (SAM contact in P1). Differential reactivity profiles for DMS are shown in Supplementary Fig. S35.

In summary, RBM-generated sequences with high scores exhibit, on average, the same structural response to SAM as natural aptamers. In contrast, aptamers generated by the rCM and RBM sequences with lower scores do not reproduce the characteristic features associated with structural switch (Supplementary Fig. S18).

Statistical evaluation and properties of generated aptamers
Reactivity profiles are notoriously variable at the single-site level, with small differences between the distributions of reactivities expected for paired and unpaired sites. This variability can be ignored when looking at average effects over a large class of many molecules, e.g., natural or generated sequences, as done above. However, predictions for single sequences require the introduction of a proper statistical framework that integrates reactivities over a set of multiple hallmark sites and enhances the statistical signal.

SHAPE and DMS reactivities are intrinsically stochastic, and the distinction between closed and open bases should be understood in probabilistic terms. We show in Fig. 7A the histogram of SHAPE reactivities of sites expected to be base-paired (teal) or unpaired (gold) in presence of SAM according to the consensus secondary structure. Unpaired sites are characterized by a different distribution of reactivities with a longer tail on high values than base-paired sites; further validating the consensus secondary structure33 obtained by the covariation in the alignment and the large epistatic scores in Fig. 3 for secondary contacts. This picture also holds for DMS reactivity distributions, see histograms for base-paired and unpaired nucleotides in Fig. 7B.

A clear confirmation that structural information can be extracted at the distribution level is presented in Figs. 7C and D corresponding to, respectively, SHAPE and DMS data. The histogram of the reactivities of the sites associated with the pseudoknot (black) in the absence of SAM is compatible with the histogram of unpaired sites, consistently with the expected conformation of most aptamers in this condition (Fig. 1). In the presence of SAM, the histogram of pseudoknot reactivities shifts toward the distribution of paired sites. This is consistent with the occurrence of a conformational switch in most aptamers, leading to formation of the pseudoknot upon SAM addition. Similar observations can be made for the P1 helix (Supplementary Fig. 15).

Based on the findings above, we introduce a statistical approach to capture the information about structural changes present at the distribution-level in reactivity data. Let M be the set of hallmark sites showing significant reactivity changes in natural aptamers in response to SAM (Fig. 5E). This set includes the pseudoknot, SAM contacts, a kink turn and sites involved in base triples (see Supplementary Table S2).

We then define, for each aptamer and each condition (with or without SAM), a Protection Score S for the propensity that sites in M are paired. Formally, S is a log-likelihood ratio between these sites being all paired and all unpaired11,12 computed from the histograms of paired and unpaired sites in Figs. 7A (SHAPE) & B (DMS). The score S also accounts for sampling noise arising from limitations on the sequencing depth42, which may strongly impact some experiments, see “Methods”. We emphasize that aggregating multiple sites in the score is crucial to reduce the statistical noise intrinsic to chemical probing measurements (see Fig. 2G and Supplementary Figs. S25, S39). Furthermore, when SHAPE and DMS data are available for the same aptamer, the two protection scores can be summed up to obtain a more robust predictor, which we refer to as DMS+SHAPE below. Figure 7E reports the SHAPE protection scores without (left) and with (right) SAM for natural (top) and generated (bottom) aptamers. For aptamers switching in response to SAM, we observe that S shifts from negative values in the absence of SAM (indicating the hallmark sites are likely to be unpaired) to positive values in the presence of SAM (indicating that these sites are involved in an interaction). Hereafter, we will call
• responsive every aptamer, whose protection score S is lower than −S0 in the absence of SAM and larger than +S0 in the presence of SAM;
• non-responsive every aptamer, whose protection score S is larger than −S0 in the absence of SAM or lower than +S0 in the presence of SAM;
• inconclusive if either score (with or without SAM) is smaller than S0 in absolute value.

We adopt a 5-fold significance threshold S0 = ln(5), see “Methods”.

As shown in Fig. 7G, aptamers responsive according to SHAPE protection scores (both natural and generated) tend to have high RBM scores. In particular, 35% of RBM-designed aptamers with RBM score >300 structurally switch in response to SAM, exhibiting significant responses in the hallmark sites. These sequences differ by 10 to 30 residues from the closest natural sequences (Supplementary Fig. S4).

In the case of failing RBM-generated sequences, the structural motifs (pseudoknot, P1, etc.) remain either protected even in the absence of SAM, or reactive in the presence of SAM. We find that most of the 45 RBM non-responsive sequences fail in the second manner: they do not have the necessary contacts even in presence of SAM. Non-responsive natural sequences can fail in both ways. None of the sequences generated with rCM is functional, possibly due to the inability of rCM to model tertiary motifs6,45. Let us stress that the number of inconclusive sequences is deeply affected by the read depth of the experiment, with lower depth leading to more inconclusive sequences, see Methods section “Statistical analysis of reactivities” for a detailed analysis of this effect.

The outcomes of the manual and automated analysis based on protection scores are compared in Fig. 7H. The two analyses are in agreement for 110 out of the 142 (77.5%) aptamers where they are both conclusive. Out of the 32 disagreements, 27 (19% of conclusives for both) are responsive in the manual analysis but not in the automated analysis. The protection-score-based analysis is more stringent, requiring a global response compatible with a structural switch across most hallmark sites. The automated analysis can also detect local responses, by focusing on smaller subsets of the hallmark sites (see last two rows of Fig. 7H, and Supplementary Section N).

To provide evidence for the reproducibility of our results, we perform two replicates of the experiment, the first one on the total set of 301 natural and artificial sequences and the second one on the 201 natural sequences only, see Supplementary Section K for a detailed description. Although some aptamers in the first replicate exhibit an overall lower response to SAM (natural and artificial), the fractions of responsive sequences in each group are consistent with the results reported in Fig. 7. Moreover, 80% of identified responders in the replicates were also responsive in the first experiment, confirming the robustness of the automated analysis (Supplementary Fig. S19).

The results above, obtained from SHAPE data, are corroborated by chemical probing with DMS. Using Eq. (11), we compute protection scores combining SHAPE and DMS reactivity data for enhanced discrimination. Figure 7I compares the results from SHAPE alone and combined DMS+SHAPE. Let us focus on natural sequences first. SHAPE and DMS+SHAPE provide the same classification (responsive, non-responsive, or inconclusive) for about 86% of the aptamers. Among the remaining 14%, more than 12% are inconclusive for one of the two approaches, and SHAPE and DMS+SHAPE disagree on less than 2% of the aptamers only.

Similar patterns are observed for RBM-generated aptamers. For RBM score >300, we obtain consistent responsive rates (ratio of the numbers of responsive and conclusive sequences) of 35%, whether estimated from SHAPE or DMS+SHAPE data. Interestingly, 48% of RBM sequences that were inconclusive with SHAPE alone can be classified with DMS+SHAPE, with one quarter responding and three quarters not responding. No rCM-generated aptamer is considered as responsive by either SHAPE nor DMS+SHAPE. A complete comparison of the analysis of the SHAPE and DMS data is reported in Supplementary Fig. S37.

Inspired by previous experimental observations for other riboswitches60,61 and Sabatier’s principle for enzymes, which require intermediate substrate binding energies for proper function62, we compute the thermodynamic energies brought by P1 helix formation using the Turner energy model as implemented in the ViennaRNA package63 (Methods). Figure 8A shows that the sequences that respond to SAM through P1 helix stabilization are confined to a thermodynamic energy window ranging from −10 to 0 kcal/mol. Similarly, pseudoknot (Pk) formation in response to SAM tend to occur for aptamers having a Pk pairing energy comprised between −8 and −3 kcal/mol (Fig. 8B). As P1 and Pk consists of, respectively, 8 and 4 base pairs, the flexible energetic window spans a range of 1.25 kcal/mol per base pair in both cases, close to a weak base-pairing energy63,64. The leftmost panels in Fig. 8 show that RBM samples preferentially have pairing energies in this intermediate band for both P1 and Pk, and are thus compatible with the structural switch required for riboswitch function.

The tables in Fig. 8 give a summary of these results. Interestingly, 25 out of the 27 aptamers that stabilize P1 in response to SAM are also responsive, in the sense of Fig. 7E, F, and show broad structural responses in other Hallmark sites (Supplementary Table S2). Similarly, 29 out of 32 aptamers that stabilize Pk are also responsive. On the other hand, out of 112 identified responsive aptamers, in natural and artificial sequences, only 19 do not stabilize P1 significantly after binding SAM. These aptamers must exhibit significant compensatory stabilization of other structural Hallmark motifs from Supplementary Table S2. It is important to note that P1 can have a more flexible behavior in the full riboswitch due to competitive interaction with the expression platform, compared with the aptamer only. As shown in Supplementary Fig. S41 the P1 helix can be destabilized in the full riboswitch context, whereas other helices like P2 or P4 are not affected, see Supplementary Figs. S42 and S43. Taken together, these results are consistent with the known importance of the pseudoknot and P1 in the response of the aptamer.

Notice that, in the central panels of Fig. 8, we show only aptamers for which the the statistical analysis yields a conclusive response for P1 or Pk. Inconclusive aptamers also tend to have intermediate pairing energies for P1 and Pk, consistent with structural flexibility (e.g., breathing).

Further explorations of RNA switch diversity through design
We then perform a second batch of design and experimental validation to further assess the limits of our generative models. We probe a total of 450 generated aptamers, whose sequences are projected onto the MSA PCs in Fig. 9A.

First, we sample sequences with the RBM model exhibiting higher distances from their closest natural counterpart, focusing on RBM scores >300. In addition, as some natural sequences lack P4, we retain a subset of RBM generated sequences having severely diminished P4 lengths. These are clearly seen in Fig. 9A, clustered at the top-left corner of the plot (recall the top PC1 represents P4 deletion). We also sample more RBM sequences of high scores (>300 and >310) to obtain better statistics on the fractions of working aptamers.

Second, we consider two variations of rCM, which is over-regularized to capture distant sequences in Rfam alignments46. We rebuild a non-regularized CM trained on the same MSA, which we call Denoised CM, or dCM for short (Supplementary Fig. S8 and “Methods”). Furthermore, as CM are unable to model pseudoknots, we devise a permutation of the MSA columns that undoes the pseudoknot, see Fig. 9B. We trained a new CM variant on the permuted MSA, that we call Unknotted CM (uCM), properly taking into account covariations in the pseudoknot. We generate sequences with such model and permute back the pseudoknot columns (Methods).

Interestingly, both dCM and uCM share some of the properties of rCM noted previously. First, CM-generated sequences from all variants have predominantly low RBM scores <300, see Supplementary Fig. S9. Second, CM generated sequences exhibit restricted diversity, concentrating in a central region of the PCA plot, as in Fig. 4D. In particular, all CM are unable to generate sequences without the P4 helix. Sequences sampled from uCM have better complementarity and Turner energies favorable for base-pairing along the pseudoknot.

We then perform SHAPE-MaP experiments and analysis. Results are summarized in Fig. 9C, and show the RBM scores of the probed aptamers against the Hamming distances to the closest natural sequence.

Out of the 248 conclusive RBM sequences in the second batch, 22% switch in response to SAM (Table in Fig. 9D). The percentage of responsive among the sequences closer to the natural ones is higher and compatible to what we obtained in replicate 1 considering the error bars, see Fig. 7G.

Moreover, 25% of the RBM aptamers having P4 length ≤1, respond to SAM; an example reactivity profile is shown in Fig. 9E. We also find a few switching aptamers differing by 30–50 sites from any natural sequence. An example reactivity profile for such sequence is shown in Fig. 9F. The reactivity profile is compatible with the consensus secondary structure, with most reactivity peaks tending to occur in unpaired loops (except a portion of P3 that remains reactive), and an overall protection in response to SAM compatible with binding and stabilization of the aptamer. Notice that RBM generate diversity not only in highly variable parts of the sequence, but also in more conserved sites (Supplementary Fig. S4).

These results support the generalization ability of the RBM. In contrast, only 3 out of 20 conclusive dCM samples switch in response to SAM (15%), and only 1 out of 16 from uCM (≈6%). Thus the dCM and uCM perform better than rCM, but not as good as RBM.

Methods
MultiplesequencealignmentofSAM-Iriboswitches
The RF00162 family from the Rfam database41 groups sequence homologs of the aptamer domain of the SAM-I riboswitch. We downloaded a manually curated seed alignment from Rfam (version 14.7), containing 457 aptamer sequences supported by literature evidence. These seed sequences are aligned to a consensus secondary structure (shown in Fig. 1B) that has been informed by the holo-form of SAM-I riboswitch crystal structures35,51. After removing extended stems and variable loops, labeled as insertions in the alignment, we obtain 108 matched positions (including gaps that mark deletions) spanning four helices that interleave around a central four-way junction. We trained a covariance model (CM)45 on this seed alignment using Infernal6 with default settings. Following standard protocols46, we acquired 6161 additional sequences from Rfam, collected from genome databases and filtered for significant matches to the CM. We constructed a MSA with these sequences, that we refer to as the full MSA, to distinguish it from the seed MSA consisting only of the 457 manually curated seed sequences. The sequence conservation logo of the full MSA is shown in Fig. 1C.

Infernalpipeline
Infernal6 is a set of computational tools to facilitate modeling RNA sequence families under a profile stochastic context-free grammar formalism, also known as covariance models (CM)45. A CM is capable of modeling the conservation profile of important sites along the sequence, as well as correlations between distant sites required by the complementarity of base-pairs in a given secondary structure. Infernal is routinely used in the maintenance of alignments in the Rfam database41,46. We employed Infernal to construct the RF00162 full MSA, that we use to train the RBM. By restricting to covariations in the secondary-structure, CM can be efficiently implemented with dynamical programming algorithms45. However, these assumptions also imply that CM is unable to include additional constraints in the probabilistic sequence model, such as pseudoknots and other tertiary contacts in the 3-dimensional fold of the RNA molecule. The CM model downloaded from Rfam is significantly regularized so that it is more effective in fetching far homologs of a family in deep genome searches76. We will refer to this CM model as Rfam CM, or rCM for short.

DenoisedCM. Since rCM is strongly regularized, in this work, we also trained a CM model variant on the full MSA, with no regularization, which we call Denoised CM, or dCM for short. This model reproduces more closely some statistics of the full MSA (conservation and covariances associated with the secondary structure).

UnknottedCM. A CM model cannot model pseudoknots and other tertiary contacts. Based on our knowledge of the consensus secondary structure of the SAM-I riboswitch aptamer (Fig. 1B), we devised a third CM model able to account for sequence covariation in pseudoknot sites constructed as follows. Columns 77–80 of the MSA, corresponding to the sites on the 3′-end part of the pseudoknot, were moved and inserted after site 28, right next to the the sites at the 5′-end of the pseudoknot. In this way, the pseudoknot is “unknotted”, and is now representable in the CM model as part of a pseudo-secondary structure corresponding to the permuted MSA. Accordingly, we proceeded to train a CM model on the rearranged full MSA. We call the resulting model Unknotted CM, or uCM for short.

SamplingtheCM. To better understand the limitations of CM models and the advantages of RBM, we sampled 10,000 sequences from each of the three CM described above. For the uCM, the rearranged columns are permuted back to their original positions after sampling. We used Infernal’s cmemit program with default parameters, and without insertions. Infernal computes a score of sequences aligned to the CM, related to the likelihood of the CM to emit a given sequence (also called bit-scores). We computed this score using cmalign, with -g (global) option to avoid local approximations76.

RestrictedBoltzmannmachines
RBM77 are bipartite graphical models over N visible variables v={v1,v2,...,vN} and M hidden (or latent) variables h={h1,h2,...,hM}, see Fig. 2A. Here N=108 corresponds to the sequence length of the RF00162 alignment, and vi encodes the nucleotide present at position i of a sequence. For RNA, vi can take one of q=5 possible values, corresponding to the nucleotides A, C, G, U, and the alignment gap symbol (⊟). The hidden variables hμ are here real-valued. The two layers are connected through the interaction weights wiμ. An RBM defines a joint probability distribution over v and h through

P(v,h) = 1/Z e−E(v,h), (1)
where Z is a normalization factor, known as the partition function, and the energy E(v,h) is given by

E(v,h) = ∑i=1^N Vi(vi) + ∑μ=1^M Uμ(hμ) − ∑i=1^N ∑μ=1^M wiμ(vi)hμ (2)
The functions Vi(vi), Uμ(hμ) are potentials biasing the distributions of single units. The visible units vi can take a finite number of possible values, and therefore the quantities Vi(vi), also called ‘fields’, can be stored as a q×N matrix. Similarly, the weights wiμ(vi) can be stored as a q×N×M three-dimensional tensor. The hidden variables, on the other hand, are continuous, and we chose to parameterize their potentials with the double Rectified Linear Units (dReLU) form proposed in ref. 36,

Uμ(hμ) = γμ+ hμ + γμ− hμ + 1/2 (hμ)^2 − θμ+ hμ − θμ− hμ (3)
with real parameters γμ±, θμ±, satisfying γμ±>0. The dReLU is an attractive choice because it is expressive enough to cover several interesting settings. When γμ+ = γμ− and θμ+ = θμ−, Eq. (3) becomes a quadratic (i.e., Gaussian) potential, and is closely related to Direct-Coupling Analysis models popular in protein sequence modeling13,78–82. However, the Gaussian choice is unable to parameterize more than two-body interactions, which can be a limitation in RNA structure where some interactions are known to involve more than two sites (e.g., stacking interactions64,83), as well as functional interactions that can span complex, extended structural and sequence motifs. dReLU can also adopt a bimodal form when θμ+ > 0 > θμ−, which is helpful for clustering.

The likelihood of visible configurations under the RBM can be obtained by marginalizing over the states of the hidden units:

P(v) = 1/Z ∫ e−E(v,h) dh = 1/Z e−Eeff(v) (4)
where −Eeff(v) is the resulting RBM score that incorporates effective interactions arising from the marginalized latent variables (see Fig. 2C):

Eeff(v) = ∑i=1^N Vi(vi) − ∑μ=1^M ln ∫ e^(∑i wiμ(vi)hμ) − Uμ(hμ) dhμ (5)
Although evaluating P(v) is computationally difficult (because the partition function Z is intractable), Eq. (5) shows that the score −Eeff(v) can be computed efficiently.

The computation of epistatic scores follows36. Further details about our RBM implementation for training and sampling are given in Supplementary Section A.

Biophysicalenergycalculations
We computed biophysical pairing energy predictions for the formation of P1 and the pseudoknot of various sequences using the Turner energy model, as implemented in the ViennaRNA package63 (Methods), with the RNAeval program.
• For the P1 helix, we computed the energy difference of each sequence in the consensus secondary structure of the aptamer domain, where P1 is paired (Fig. 1B), and in a conformation where P1 is unpaired (Fig. 1A).
• To estimate the pairing energy associated to pseudoknot formation, we used RNAeval on a virtual secondary structure where only the pseudoknot sites are base-paired, and all other sites are unpaired. We then considered only interior loop contributions to the resulting folding energy.

Note that, in both cases, intrinsic limitations of the ViennaRNA algorithmic implementation imply that we cannot model the pseudoknot together with other structural elements (and other tertiary contacts).

Selectionofsequencesforfirstbatch
We probed a total of 306 sequences, breaking down as follows.

RBM sequences. We generated sequences from the RBM by Gibbs sampling. Equilibration was assessed by monitoring the average score of the sample. We found that 5000 steps were more than sufficient. We then sorted these sequences by their RBM score (−Eeff), and selected 70 sequences at random, uniformly spanning the range of scores observed in the sample. The table of sequences and their associated RBM scores is reported in the Supplementary Code listing84, see Section N.

Infernal sequences. We then sampled sequences from the rCM of the RF00162 family, downloaded from Rfam. We used the Infernal cmemit program (see “Methods”) to sample a large batch of sequences. We selected 30 sequences uniformly spanning the range of bit-scores of the samples.

Natural sequences. We selected 151 sequences members of the seed MSA and 55 sequences members of the full MSA, as described in section “Multiple sequence alignment of SAM-I riboswitches”. The selected natural sequences are diverse, spanning various taxonomic classes (see Fig. 4B). A listing of probed sequences can be found in Supplementary Data 2.

Selectionofsequencesforsecondbatch
In the second experiment, we generated a total of 450 sequences to be probed, of different origins. We considered:
• 58 CM sequences, with 29 from uCM and 29 from dCM (see section “Infernal pipeline” for definitions of these CM variants).
• 392 sequences sampled from the RBM, filtered to have RBM scores >300. In particular, 49 of them were selected because they had no P4 helix, while 100 of them were selected because they had larger Hamming distances from any natural sequences.

The full list of designed sequences is provided as part of the Supplementary Code listing, see Section N.

SelectionofsequencesforDMSprobing
We selected a subset of aptamers from batches 1 and 2 for DMS probing. From batch 1,
• 84 sequences generated by RBM;
• 16 sequences generated by rCM;
• 152 natural sequences.

From batch 2,
• 102 sequences generated by RBM;
• 10 sequences generated by uCM or dCM.

The full list of sequences probed by DMS is provided as part of the Supplementary Code listing, see Section N.

Chemicalprobingexperiments
RNA preparation. DNA oligonucleotides representing the 206 SAM-I natural sequences, and the two batches (100 and 450) of artificial sequences, preceded by the T7 promoter (5′CGGCGAATCTAATACGACTCACTATAGG3′) and followed by a tag sequence representing a 10 nucleotide barcode unique for each aptamer and a primer binding site, were purchased as an oligonucleotide pool (Twist bioscience). The Tag sequence was designed to avoid interference with the aptamer secondary structure using RNAFold63 (see ref. 85 for the tag design method). The oligo pool was PCR amplified using the T7 promoter as forward primer and five different reverse primers (5′ GGAAGGAGGCGGGCAGACG3′, 5′CGTATTACCGCGGCTGCTGG3′, 5′CGACGAGATAGGCGGACACTGG3′, 5′CGACGAGATAGGCGGACACTGG3′, 5′GAAGTCGTAACAAGGTAGCCGAT3′), provided in Supplementary Data 1. RNA was transcribed, prepared, and checked for the absence of aberrant products on a 1% agarose gel86. See Supplementary Sec. S for details.

Read depths vary with the choice of the primer, see Supplementary Figs. S27, S28, S29. As explained in Methods section “Statistical analysis of reactivities” we have verified that our statistical analysis give consistent rates of responsive aptamers, even for primers with lower coverage.

SHAPEandDMSprobing. SHAPE chemical probing was performed as described previously42. Briefly, 10 pmol of RNA were diluted in 12 μL of water and denatured for 3 min at 85∘C. Then, 6 μL of 3X pre-warmed folding buffer with or without magnesium (0.3 M HEPES pH 7.5, 0.3 M KCl, 15 mM MgCl2) were added and the solution was allowed to cool down to room temperature. Samples were then incubated at 30∘C for 5 min. SAM was added at final concentrations of 0, 0.1 or 1 mM and samples were incubated 15 min at 30∘C. 9 μL (corresponding to 5 pmoles) were aliquoted and 2 μL of 50 mM 1M7 (1-Methyl-7-nitroisatoic anhydride) or DMSO (Mock reaction) was added and allowed to react for 6 min at 30∘C. For dimethyl-sulfate (DMS) probing, 0.9 μL of 600 mM DMS stock solution (or 0.9 μL of ethanol for mock reactions) was added and allowed to react for 10 min at 30∘C. DMS probing reaction was then quenched by adding Tris pH 8.0 at 400 mM. RNA were then reverse transcribed with the Superscript III reverse transcriptase (Invitrogen®) and NGS libraries were prepared using NEBNext Ultra II DNA Library Prep Kit (New England Biolabs®). Final products were sequenced by using the Illumina technology (NextSeq 500/500 Mid 2×150 flow cell). Sequencing data were analyzed and reactivity maps were derived using ShapeMapper287. In the end, the 306 selected sequences were probed in the following conditions:
• 30∘C, without Mg2+ and without SAM.
• 30∘C, with magnesium (Mg2+).
• 30∘C, with magnesium and two concentrations (0.1 and 1 mM) of SAM.

Each probing reaction was repeated in triplicate. The two SAM concentrations were analyzed together to improve statistics, since we found no significant effect of varying the SAM concentration in the reactivity responses of the aptamers (see Supplementary Fig. S21). The reading efficiency per site (read depths reported by Shapemapper) is plotted for the tested aptamers as grouped by primers in Supplementary Figs. S27, S28 and S29.

Manual inspection of reactivity profiles. IPANEMAP88 was used to generate RNA secondary structure models for each sequences. For manual inspection, we considered the reactivity of the nucleotides known to be directly involved in SAM binding (U7, G11, A46, U69, G70, U103) and of those known to be protected from shape reactivity in the closed state, i.e., nucleotides in P1 (1–8; 101–108), in the pseudoknot (25–28; 77–80), those involved in the three base triple interactions (24, 73, 74, 76, 100). Nucleotide numbering follows the profile shown in Fig. 1C. An aptamer was considered to bind SAM if at least three of these elements are noticeably less reactive upon SAM addition, and if none of the binding determinant remain highly reactive. Note that P1 and the Pk are each considered as one element, and that some of the elements may be unreactive even in absence of SAM.

Statisticalanalysisofreactivities
Reactivity definition. SHAPE-MaP experiments result in measurements of sequencing error rates at each site of the RNA sequence, that correlate to the locations where the SHAPE probe has reacted with the RNA. For each site i=1,…,N of a sequence n, the reactivity is defined by ref. 42:

rin = (min − uin) / din (6)
where min is the mutation rate in presence of the reagent, uin is the mutation rate in its absence accounting for mutational background of the experiment, and din is the mutation rate in a denaturating condition where the RNA is expected to be unfolded, intended to cancel sequence-dependent biases. Working with rin is usually better since this form should cancel site-dependent biases in the raw SHAPE mutation rates, min. The basis of the SHAPE-MaP procedure relies on differences in the distribution of reactivities in base-paired and unpaired sites42. We have confirmed such differences are observed in our data in Fig. 7 (and also Supplementary Fig. S13).

Statistical analysis. The finite number of sequencing reads collected at a site implies a statistical error in the reactivity computed by Eq. (6). Therefore, we cannot directly access the true reactivity rin at a site, but rather an experimental measurement ~rin that fluctuates according to the number of reads taken at the site. To model this uncertainty, we make the simplifying assumption that the ideal reactivity of a site, rin, depends only on whether the site is base-paired (bp) or not (np). Under this assumption, we can write:

P(~rin|bp) = ∫ P(rin|bp) P(~rin|rin) drin (7)
where:
• P(~rin|bp) is the probability of measuring reactivity ~rin at site i of sequence n, given that the site is base-paired and conditioned on the finite number of reads taken at this position.
• P(~rin|rin) is the probability of measuring reactivity ~rin at site i of sequence n, on account of fluctuations due to a finite number of reads, conditioned on this site having a real reactivity of rin.
• P(rin|bp) is the probability distribution of reactivities of base-paired sites, at infinite read-depth, assumed to be homogeneous across sites.
• P(~rin|np) and P(rin|np) are defined in a similar manner for non-paired sites.

We approximate the distributions P(rin|bp) and P(rin|np) by kernel density estimators fit on the corresponding empirical histograms (shown in Fig. 7A for the first experiment). The kernel function used corresponds to standard normal, with a bandwidth set according to the Silverman rule89. To better estimate the histograms, we use the experimental conditions with SAM, where the secondary structure of the aptamer is expected to be more stable. We also find that these histograms can depend on the particular experiment, and therefore we fitted P(rin|bp), P(rin|np) for each replicate.

Applying Bayes theorem90 in Eq. (7), we can write:

P(~rin|bp) / P(~rin|np) = ∫ (P(rin|bp)/P(rin|np)) P(rin|~rin) drin (8)
where P(rin) is the histogram of real reactivities, regardless of whether a site is paired or not. The posterior P(rin|~rin) quantifies the uncertainty of the real reactivity rin at site i of sequence n, conditioned on our information of the measurement taken at this site. This uncertainty arises from the finite sequencing reads available, which induce an experimental error in our estimate of the quantities m, u, d appearing in Eq. (6). Since the mutation count at a site can be modeled by a Poisson distribution42, the posteriors of the mutation rates m, u, d are Gamma distributions, with a convenient choice of conjugate prior90. Then, we can produce a Monte-Carlo estimate of P(rin|~rin) by sampling the posterior Gamma distributions of m, u, d, and computing the reactivity through Eq. (6). If the sampled reactivities fall predominantly far in the tails of the histograms P(rin|bp) or P(rin|np), respectively, the reactivity measurement is discarded as an outlier. In practice, we find that 1000 samples for each site are sufficient. These samples can then be used to approximate the numerator and denominator of the right-hand side of Eq. (8). In this way, we produce estimates of the ratios P(~rin|bp)/P(~rin|np), quantifying the odds that a site is paired. Supplementary Fig. S23B shows a scatter plot of reactivities in our dataset, with the standard-error estimated by the standard SHAPE-Mapper pipeline42 (which does a first-order error propagation through the Poisson count statistics), with each point colored according to the value of the log-odds-ratio Eq. (8). Dashed lines are approximate contours separating points that are over twice more likely to be paired (blue) or unpaired (red). The fact that these contours are not straight vertical lines indicates that, using Eq. (8), we are considering both the reactivity value and its uncertainty in assessing the plausibility that a site is paired or not. A similar approach has been proposed by 11,12. See also Supplementary Section I for further discussion and tests.

Protection scores. We can exploit the likelihood ratios P(~rin|bp)/P(~rin|np) computed above to estimate the probability of the presence of a structural motif in a sequence. We define a motif of length 2L as a set of base-paired sites, M={i1,j1,...,iL,jL}. For example, the P1 helix motif corresponds to {1,108,2,107,…,8,101}. We then probabilistically assess the presence or absence of the motif M in molecule n by comparing the value of the protection score

Sn(M) = ∑i∈M ln (P(~rin|bp)/P(~rin|np)) (9)
to some thresholds ±S0, see section “Restricted Boltzmann machines”. This approach allows us to combine multiple reactivity measurements into a robust probabilistic measure, achieving more statistical power than when site reactivities are analyzed one by one.

This approach can be applied to SHAPE or DMS reactivity data. As DMS probing is efficient in detecting interactions involving nucleotides A or C predominantly, we only consider DMS reactivities estimated using only reactivities measured at sites with A or C nucleotides.

CombiningSHAPEandDMSdata. When both SHAPE and DMS data are available for the same aptamer, we can combine them to obtain better predictions about the base-pairing status of a site. Since the SHAPE and DMS reactivities are obtained in independent experiments,

P(~rin,SHAPE, ~rin,DMS|bp) = P(~rin,SHAPE|bp) × P(~rin,DMS|bp) (10)
where ~rin,SHAPE and ~rin,DMS denote SHAPE and DMS reactivity data at the same site i of aptamer n. This independence implies that the log-odds ratio of the pairing status of a site or a structural motif (as in Eq. (9)), in presence of both kinds of data, can be computed by simply adding the protection scores obtained from each kind of probing alone:

Stot(M) = SSHAPE(M) + SDMS(M) (11)
where SSHAPE is the protection score obtained from SHAPE data, and SDMS the protection score obtained from DMS data.

Error bars on the rates of responsive aptamers. Given Nconc. conclusive probed sequences, Nresp. of which are found to be globally responsive, we estimate the response rate by p=Nresp./Nconc.. The uncertainty of p is, according to the binomial law, err = sqrt(p(1−p)/Nconc.). The response rates in Figs. 7, 9 are reported as (p±err)×100%.

We have investigated the dependence of these uncertainties on the SHAPE-Mapper read depths, which varies with the primers. Supplementary Fig. S32 shows that the inconclusive rate is strongly anti-correlated with the read depth. Both rates of responsive and non-responsive sequences increase with the read depth, a consequence of the decrease of the statistical noise. To mitigate this statistical effect, throughout this work, the response rate is computed as the ratio of responsive molecules over the number of conclusive ones, compare top and bottom panels in Supplementary Fig. S32. The dispersion due to this statistical noise are accounted for by the error bars in the results shown in Fig. 7G as explained above. We have also investigated the dependence of DMS results on the read depth (Supplementary Fig. S33). As with SHAPE, the inconclusive rate increases with the read depth.

SHAPEprotectionscoresareinagreementwithconsensussecondarystructure
Sequence homologs in the RF00162 family are collected based on similarity to a group of manually curated sequences in the seed. Overall, for many of these sequences (both in the seed and in the full alignment), direct experimental evidence of their actual behavior and structure is limited, except for specific cases, such as the Thermoanaerobacter tengcongensis and the Bacillus subtilis yitJ SAM riboswitches, which have been extensively studied in the literature fueled by detailed knowledge of their published crystalized structures35,51. For many other sequences in the MSA, their actual behavior is at most hypothesized based on indirect evidence.

We have here obtained detailed SHAPE data of Bseed =151 sequences of the seed alignment. Our data show that, in average, these sequences are compatible with the consensus secondary structure of the RF00162 family, shown in Fig. 1B. Indeed, we have computed the average protection scores 〈S(i)〉seed for each site i, over the sequences in the seed alignment probed in our experiments,

〈S(i)〉seed = 1/Bseed ∑n∈seed ln (P(~rin|bp)/P(~rin|np)) (12)

Figure 10B plots 〈S(i)〉seed in the conditions with SAM and without SAM. Overall, the averaged protection scores are in detailed agreement with the consensus secondary structure of the aptamer, depicted in Fig. 10A. Helices P2, P3, P4 are seen to be base-paired in average in all conditions, with a mild overall increase in the values of S with the addition of magnesium and then SAM, indicating overall structural stabilization. The central junction loop (CL), and the loops on the second helix L2, the third helix L3, and the fourth helix L4, are consistently measured as reactive when SAM is not present, indicating that these sites are unpaired, as expected. Besides these major structural motifs, we also appreciate finer details such as the reactivity of single isolated bulge sites in positions 46 and 65 in absence of SAM.

Next, comparing the behavior across different conditions, we appreciate the effect of magnesium and SAM on the structure. We highlight (in green) sites that change significantly in response to SAM. These include sites in direct contact with SAM (as known from the crystal structure35), and other tertiary motifs known to form in response to SAM. We discuss these next.

Selectionofhallmarksites
We selected 24 hallmark sites across the aptamer sequence, for which we could rationalize observed reactivity changes in response to SAM binding, and which are consistent with expectations from previous chemical probing studies on SAM-I riboswitches and previous structural data. These sites also exhibit significant reactivity responses across natural sequences in our data, see Fig. 10. They are listed in Supplementary Table S2. In Supplementary Section Q we include further discussion and references to several previous literature reports justifying the choices of each of these sites.

Our results are robust to minor variations in the selection of Hallmark sites used to evaluate the response of aptamers to SAM. For example, although we could not find previous reports of reactivity responses in J4/1, we find in some cases that sites 98 and 99 exhibit protection upon SAM binding (see Fig. 10). We tried adding few selected sites (such as 98, 99), or excluding some, and confirmed that our main results (such as numbers of responsive sequences) remain unchanged. Additional results are reported in Supplementary Section Q.

Principalcomponentanalysis
We carried out a PCA of the natural MSA. First, we one-hot encode the natural sequences in a q×N×B binary tensor D, where B=6161 is the number of sequences in the full MSA collected above. The tensor has Din,a = 1 if sequence n of the alignment has symbol a∈{1, . . . , 5} at position i, and otherwise Din,a = 0. We then compute a covariance tensor, defined as follows

Cijab = 1/B ∑n Din,a Djn,b − (1/B ∑n Din,a) (1/B ∑n Djn,b) (13)

We flatten the tensor Cijab into a qN×qN matrix, and then perform a standard eigenvalue decomposition on it. Individual sequences are then projected along the two top components (with largest eigenvalue) of the decomposition.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors propose and test the implementation of Restricted Boltzmann Machines (RBM) as a suitable neural network architecture for functional prediction of ON-OFF behavior in SAM-I riboswitches. The trained RBM model is then used to screen for artificial allosteric SAM-I aptamers with desired properties, followed by experimental validation of aptamer functionality using 476 RBM-suggested sequences.\nIn principle, the generation of improved machine learning models to accelerate the design of aptamers for biological research and related applications such as the one presented in this work can be both relevant and timely. However, in my opinion, this group has failed to provide a comprehensive enough landscaping of their work to compare their method with comparable published work (2-6 years ago), such as the following:\nGroher, A.C., Jager, S., Schneider, C., Groher, F., Hamacher, K. and Suess, B., 2018. Tuning the performance of synthetic riboswitches using machine learning. ACS synthetic biology, 8(1), pp.34-44.\nAngenent-Mari, N.M., Garruss, A.S., Soenksen, L.R., Church, G. and Collins, J.J., 2020. A deep learning approach to programmable RNA switches. Nature communications, 11(1), p.5057.\nIwano, Natsuki, Tatsuo Adachi, Kazuteru Aoki, Yoshikazu Nakamura, and Michiaki Hamada. \"Generative aptamer discovery using RaptGen.\" Nature Computational Science 2, no. 6 (2022): 378-386.\nWhile I believe the use of Restricted Boltzmann Machines for the task of unsupervised learning of RNA switch functionality to aid construct design is somewhat novel, the authors do not provide a clear rationale of why RBMs are chosen as the primary modeling architecture, as opposed to simple standard machine learning models like random forests or even other types of neural network architectures that can learn probability distributions over its set of sequence inputs such as beta variational autoencoders, where latent features can also be investigated for interpretability. Instead, the authors cite 4 references (three of which come from the same research group and share a set of authors as this work in review ) where RBMs are used to aid in protein design as the sole justification for using RBMs in the task of interpretable unsupervised RNA switch functionality prediction and design. If this were one of the first-ever demonstrations of RNA switch functionality prediction and design using machine learning techniques, then this work would be acceptable for publication as is, but given the fact that various other demonstrations of such concept have been shown in the literature already in the last years (see 3 above references), then I believe the authors should aim to improve by providing the following:\n1) Conduct an in-silico comparative analysis of their RBM technique with a range of other unsupervised learning architectures (standard and emergent) to demonstrate that this technique is indeed advancing the field of RNA switch design using ML techniques and is not just “another” technique of comparable predictive capacity to others published in the past when used on their particular use case. Furthermore, the authors should generate a baseline for their predictive method through a comparison with the use of traditional non-machine learning energy-based calculations for secondary structure prediction using thermodynamic principles predictive capacity of using traditional energy-based calculations such as through the use of the features automatically provided by software packages like NUPACK and ViennaRNA that can calculate minimum free energy and ensemble defect of RNA sequences.\n2) The authors provide a Git with their source code and a Pluto.Jl notebook implemented in Julia. This is great, but unfortunately, this choice requires the reader to host their code locally to reproduce the results of this work. I believe this is a big barrier to reproducibility, and the authors should aim to reduce this barrier even through either of the two paths. A) Hosting a live demo of their technique for RNA riboswitch functionality prediction and sequence design as an online tool or B) Generate a version of their code that runs in a free cloud-hosted service for scientific computing such as a Google Colaboratory, which would allows users to easily run their code with all necessary imports, including the data provided from the paper. This Google Colab notebook should run training and inference of their RBM architecture along with any comparison analytics and result generation, which the audience of nature comm would greatly appreciate. I understand this second option would require the authors to make an implementation of their source code using Python instead of Julia, which is why the first option may be a preferable path for the authors to improve their work. In my opinion, the second option would be easier to maintain and would allow readers to access a working version of this code with ease and without relying on the author's servers.\nBeyond the previous comments, this work's methodology appears sound, and its main conclusions and claims seem to be justified, particularly through its experimental data with 476 RBM-designed sequences. As an experimentalist myself, I appreciate that and congratulate the authors for doing so.", "reviewer_2": "This manuscript describes the use of Restricted Boltzmann machines (RBM) for the design of RNA switches that respond to the metabolite S-adenosyl methionine (SAM). The authors train their model based on riboswitches in the SAM-I family and apply it for generative design of artificial SAM-responsive RNAs. These artificial RNAs are then tested using SHAPE-MaP probing and compared to riboswitches from nature and those generated using covariance models. While they find that ~67% of the natural riboswitches show SAM-drive switching in SHAPE-MaP, they find ~22-39% show switching depending on the set of RBM-generated designs tested. While the results are interesting and could constitute an important advance for RNA design, the manuscript should address the issues below to make it suitable for Nature Communications:\n1. At present, the paper focuses entirely on the SAM-I family of riboswitches, so the broader utility of the approach for generating new functional RNAs is unclear. Can the RBM approach also be applied to riboswitches with other ligands or other functional transcripts?\n2. The manuscript is focused on riboswitch design but the characterization of the RNAs is limited to SHAPE-MaP assays, which interrogate RNA structure but do not test the riboswitches in their normal regulatory context. Testing some of the designs to determine if they regulate gene expression or transcription in response to SAM would be a more convincing test of their design approach and confirm true riboswitch function.\n3. While the authors do make use of Hamming distances to compare their RBM designs to natural riboswitches, it is unclear from the manuscript exactly where and what sequences the RBM model is changing for the designs that it generates. Are the modifications typically occuring in stem regions and loops, or are they being made in more sensitive regions necessary for riboswitch functon(e.g. P1, Pk, triple bases)?\n4. For Fig. 7A, the P1 pairing energy can take on positive values. It is unclear from the Methods description how this energy can take on positive values. Is the anti-terminator being included in the structure (Fig. 1A)? A more detailed description of this calculation should be provided.\n5. The descriptions of the sequences actually tested is very lacking. I can only find a single supplementary spreadsheet listing 201 RNA sequences (presumably the natural sequences from Fig. 6). There does not appear to be any information on the other sequences, although an unnumbered table is mentioned on line 949. Full sequence information and status of the testing results for each sequence should be provided. The RBS score should be included with all the riboswitch sequences.\n6. In the Supplemental Material, multiple places have \"??\" in place of figure numbers. This should be corrected.\n6. In Fig. 1A-B and elsewhere, the diagram that the authors use for indicating addition and removal of SAM does not clearly show in which direction SAM is added or subtracted. This diagram should be updated.", "reviewer_3": "This manuscript uses restricted Boltzman machines (RBM) to produce synthetic sequences that resemble SAM-I riboswitch aptamer molecules. Using information obtained from alignments of naturally occurring SAM-I riboswitch aptamers, the RBM trains a probabilistic model specific for this molecule, which can then be used to sample novel sequences. Chemical probing (SHAPE-MaP) is used to characterize 476 RBM designed synthetic sequences with respect to naturally occurring SAM-I riboswitch aptamers.\nThe manuscript reports a \"30% success rate of correctly structured aptamers that undergo a structural switch in response to SAM.\" My main criticism to this manuscript is that the SHAPE-MaP experimental data cannot and does not give evidence of a \"correctly structured aptamer\" and whether the sequence undergoes or not \"a structural switch in response to SAM\".\nPerhaps the more concise positive-control experimental result in the paper is given in Figure 6E that provides the reactivity profile of a naturally occurring SAM-I aptamer, and Figure 5CD in gray that provide the change in reactivity (with/without SAM-I) in a naturally occurring aptamer. But. It would be impossible to infer the helices of the aptamer by looking at Figure 6E, and while the gray pattern in Figure 5CD looks distinct, that first impression gets questioned once one realizes that the shading correspond to a +/- half standard deviation, that is, more than 60% (assuming a normal distribution) of the actual data falls outside the shading.\nIt would be advisable to add other forms of evidence besides SHAPE when possible.\nSome specific comments and suggestions centered around the figures are,\n* Figure 1\n- Which is the actual natural sequence used to depict the SAM-I aptamer structure?\n- It would be more informative to draw the 2D structure in a way that reflects more of the 3D structure where P1/P4 coaxially stack and P2/P3 also coaxially stack and the two columns connect by the pseudoknot, and helped by the kink turn. See for instance Figure 1B in Stoddard...Batey, Structure 2010,\nhttps://www.sciencedirect.com/science/article/pii/S0969212610001851\n- At the very least, please mark in Figure 1A-1B the elements that you later in the manuscript: the 3-way contact, the nucleotides involved in binding SAM and the kink-turn in P2.\n- Figure 1C is ok, but not too informative with respect to the structure, in my opinion. I suggest replacing that with structural covariation information also obtained directly from the alignment, as Rfam does using R-scape, for instance.\n* Figure 4\nThe comparison to sequences generated with an Infernal model for the aptamer is ok, but not sure really necessary. If you keep it, some additional analysis would be recommended. The PCA analysis results in very nice figures but hard to interpret unless there is some clear association to those two first components. Maybe the Infernal sequences are longer, shorter, more conserved? Some stats about some standard variables may help explain. Also for instance, using Infernal you could align all RBM-generated sequences and also all Infernal-generated sequences amd use R-scape to see if there are any systematic difference between those alignments, and how well the two sets of synthetic sequences comparatively supported the aptamer structure.\n* Figure 5A\n- Please fix the frequency range to be [0,1]\n- In discussing Figure 5A, it is said:\n\"As expected, base-paired sites tend to have lower reactivities, and unpaired sites tend to have higher reactivities\" (page 6, 226)\nThis is an overstatement. While the two histograms have a somewhat different tail, trying to predict which residues are pair and which are unpaired from those two histograms results in a very poor prediction.\nIt is also said:\n\"This also confirms that most natural aptamers are compatible with the expected secondary structure [1]\" (page 6, 229)\nthis is also an over statement, There are several different ways to confirm the X-ray structure of the aptamer, one is the covariation observed in the Rfam alignment. However, the histogram of Figure 5A does not.\n* Figure 5CD\n- I have already discussed above the poor confidence interval used to compare natural sequences from synthetic sequences.\n- It would be informative to know as one reads the Figure which is the number of synthetic sequences averaged in those two Figures\n- The selection of RBM sequences with high scores is interesting. (A supplemental figure shows that RBM sequences with lower scores does not show significant reactivity differences with or without SAM). Would it be possible that an split of Infernal-generated sequences would produce a similar result? The figure presents results for all of tehm together.\n- A complementary analysis to asses how well a set of synthetic sequences preserve the structural properties of the aptamer could be to make an alignment of all generated sequences using the Infernal model and assessing the covariation observed in that alignment in comparison to that observed in the input alignment of natural sequences.\n* Figure 6AB\n- in these figures I am missing an assessment of false positives. A null model could be generated by repeating the same protection scores for hallmark sites after the reactivities of the naturally occuring aptamers are shuffled along the sequence.\n- Please clarify what are the different definitions: responsive vs non-responsive and switcher vs non-responsive. What is the difference between responsive and switcher? are those necessary? The reasoning is hard to follow."}, "author_response": "We would like to clarify that the focus of our work is not on the prediction of function. Such prediction tasks (e.g., does this sequence have the expected functionality?) can often be achieved with supervised machine learning methods, such as random forest models successfully used by Groher et al, or deep forward neural networks as in Angenent-Mari et al. In particular, Groher et al also exploit their models to filter the best sequences out of in-lab selection experiments (SELEX or directed evolution), starting from an initial library, and improving the target performances. These models, however, cannot be directly used to generate de novo functional sequences.\nIn contrast to these works, our RBM is a generative model, capable of sampling new, completely de novo, functional sequences. The RBM is trained on a dataset of natural sequences in an unsupervised fashion. Intuitively, the model learns to extract frequent statistical patterns (at the level of conservation or correlations) that are then combined in non-trivial ways for the generation of novel aptamers. We have then performed high-throughput generation of sequences and SHAPE-driven experimental verification of their functionality.\nIn addition, we believe it is fair to state that the RNA functionalities considered in the works cited above, such as\n(1) Optimizing the P1 stem of tetracycline-aptamers (6-10 nucleotides) starting from a random pool of sequences, to improve the dynamic range in the switching behavior by Goher et al.,\n(2) Realizing an one toehold switcher structure, with a stem, a bulge and a loop (12 nucleotides in the toehold and 18 in the stem) by Angenent-Mari et al, with a ~20-40 nt variable region designed from SELEX sequence data,\nare somewhat simpler than ours. In our manuscript, we design new RNA sequences coding for two alternative structures (including about 100 nucleotides, 4 helices with a complex pattern of secondary and tertiary contacts) controlled by the presence/absence of a metabolite.\nLast of all, in our work, the target structure/function is not unique: we do not constrain the generated sequences to correspond to one particular aptamer structure in the SAM-riboswitch family. We design new sequences that reproduce the variability of structures and functions of natural SAM riboswitches. Not all the aptamers among the natural ones which we have probed are switching (~60% responsive in 150 conclusive tested aptamers ), see also Fig. 5 in the new manuscript showing SHAPE-MaP reactivities on two natural aptamers. We thus expect that in the SAM-riboswitch family each aptamer has a slightly different structure and response to SAM, and that the designed sequences reproduce this rich variability.\nWe now give a more detailed description of the differences/similarities of our method and the works cited by Reviewer 1:\n● Groher et al. 2018.\nIn this paper a pipeline combining experiments and multiple ML methods (secondary structure predictions, graph modeling, features classification, random forest, convolutional neural networks) is introduced to select good P1-stem (with 8 to 10 nucleotides) candidates to improve the dynamic range (quantitative ratio of the on and off states) for tetracycline aptamers. These models exploit carefully manually chosen features (free energy, GC content, hydrogen bonds, …) on which a model is trained, to predict the aptamer performance and optimize a subset of the sequence (P1 stem).\nOur approach differs in several key aspects: (1) our model is trained on sequence data only in an unsupervised fashion: features are thus extracted in an automatic way and not chosen a priori; (2) it is simpler, requiring a single neural network model (RBM) for the generation of novel functional sequences: RBM is a generative model, not a supervised classifier like the neural networks used in the cited work; (3) we tackle the difficult task of de novo designing a complete riboswitch aptamer domain with >100 nts, involving several helices, a pseudoknot, and other structural elements, compared to the single stem (~10 nts) in Groher et al.\n● Angenent-Mari et al. 2020.\nThis work focuses on functional prediction with supervised neural networks, for toeholds having a one stem-loop structure, with an on-off regulatory effect on protein translation, and a total of 75 variable nucleotides in the sequence (see Fig. S1A in the paper). In lab-selection from an initial library of natural sequences is performed, and ML is used to classify sequences and predict on/off dynamical ratio, with predictions validated by experimental scans.\nAs a summary, the main differences between this work and ours are: (1) SAM riboswitches have more complex and diverse possible structures (with multiple helices and a complex network of secondary & tertiary contacts) than a single hairpin-loop structure, making the design problem harder; (2) we introduce ML methods to design novel sequences, not to classify them.\n● Iwano et al. 2022\nThis paper implements the RaptGen model that focuses on aptamer design in SELEX. These aptamers have much shorter sequences (~30 nts) than the ones considered in our work (~108nts). To answer the Reviewer’s point 1, we have considered below the similar RfamGen model (from Sumi et al. 2023), which, similarly to RaptGen, is a variational autoencoder incorporating a covariance model grammar, and compared with the RBM model on the SAM-I riboswitches studied in our paper.\nPlease notice that, in a previous work (Di Gioacchino et al. PLOS Comp. Bio. 2021), we have trained an RBM model from SELEX data on small (40 nts) aptamers selected to bind thrombin, with a setting more similar to the works above than the present manuscript.\nWe agree with Reviewer 1 that the introduction of our manuscript was not properly pointing out the state-of-the-art on the use of machine learning for the prediction of functionality in RNA aptamers and the Discussion Section was not clearly comparing our method with previous works. We have therefore expanded our literature review in the introduction (starting at line 53 of new version) and in the final discussion of our paper (line 676), with an explicit in silico comparison with other ML methods as asked in the point 1 below. Following point 1 of Reviewer 2 we have also extended the computational experiments to other riboswitch families (line 214 of the new version).\n1) Conduct an in-silico comparative analysis of their RBM technique with a range of other unsupervised learning architectures (standard and emergent) to demonstrate that this technique is indeed advancing the field of RNA switch design using ML techniques and is not just “another” technique of comparable predictive capacity to others published in the past when used on their particular use case. Furthermore, the authors should generate a baseline for their predictive method through a comparison with the use of traditional non-machine learning energy-based calculations for secondary structure prediction using thermodynamic principles predictive capacity of using traditional energy-based calculations such as through the use of the features automatically provided by software packages like NUPACK and ViennaRNA that can calculate minimum free energy and ensemble defect of RNA sequences.\nAs emphasized above, we believe our approach is new, as we are not aware of previous works using ML that (1) design RNA sequences able to switch structures with generative models learned from natural riboswitch sequences of the same family only, (2) validate the structural switch by a novel high-throughput experimental SHAPE test of the designed sequences (as kindly noticed by Reviewer 1 in his/her report).\nWe agree with Reviewer 1 that base-line models are essential for the assessment of performances. In our work, we have used Covariation Models (Infernal) as a base-line framework reproducing only the constraints coming from the secondary structure (with the corresponding 1- and 2-point statistics); we have also considered an unknotted version of CM, which allowed us to include the pseudoknot. As shown in the Table of Fig. 9, we found that only 3 out of 20 sequences designed with denoised CM and 1 out of 16 sequences generated by the unknotted CM could switch. These results confirm that generating functional riboswitches (with reasonable success probability) requires going beyond secondary-structure constraints.\nRegarding the comparisons with secondary-structure folding energies: In Fig. 8 (previously numbered as Fig. 7), we compare the RBM sores with the ViennaRNA pairing energies of the P1 helix and the pseudoknot PK for all tested sequences. Interestingly, we find that the sequences able to switch are not the most stable for P1 and Pk, but have intermediate ranges of binding energies, likely to facilitate the switching behavior thanks to a larger flexibility.\nLast of all, we compare in the new manuscript our generative model with the recent RfamGen model (from Sumi et al 2023), which like RaptGen, is a variational autoencoder (VAE) using a covariance model grammar for the encoding RNA sequences. In the new Supplementary Materials, we have added Supplementary Section L “Comparison to RfamGen” and Fig. S25 (reproduced below) to compare the generation of SAM-riboswitches by the RFamGen and RBM models. We find that RBM-generated sequences (blue) reproduce the distribution of natural sequences (gray) even according to the assigned by the VAE (panel B) of designed sequences.\nThis suggests that RfamGen does not distinguish between RBM-generated and natural sequences. On the contrary, RfamGen sequences are shifted to lower RBM scores compared to natural sequences (panel A). Taken together, these results suggest that the scoring function of the RBM model is more accurate than the one of RfamGen.\nAs a last point we would like to underline that in contrast to VAEs, RBM produces interpretable sequence motifs, informative about key nucleotides and positions for the riboswitch functionality through the position weight matrices (weights of RBM) connecting the sequence and the hidden (latent) units. See Figure 3B and C. In contrast, deep VAEs are far less interpretable.\n2) The authors provide a Git with their source code and a Pluto.jl notebook implemented in Julia. This is great, but unfortunately, this choice requires the reader to host their code locally to reproduce the results of this work. I believe this is a big barrier to reproducibility, and the authors should aim to reduce this barrier even through either of the two paths. A) Hosting a live demo of their technique for RNA riboswitch functionality prediction and sequence design as an online tool or B) Generate a version of their code that runs in a free cloud-hosted service for scientific computing such as a Google Colaboratory, which would allows users to easily run their code with all necessary imports, including the data provided from the paper. This Google Colab notebook should run training and inference of their RBM architecture along with any comparison analytics and result generation, which the audience of nature comm would greatly appreciate. I understand this second option would require the authors to make an implementation of their source code using Python instead of Julia, which is why the first option may be a preferable path for the authors to improve their work. In my opinion, the second option would be easier to maintain and would allow readers to access a working version of this code with ease and without relying on the author's servers.\nWe believe that maintaining the codebase in Julia, instead of Python, facilitates reproducibility. Julia’s package manager automatically downloads and installs all necessary dependencies, with compatible versions, and is usually much more robust than Python’s dependency resolvers.\nWe have worked on the documentation of our package to facilitate local installation for new users. We simplified the setup and wrote simple instructions in the README.md page. We have also provided Jupyter notebooks (besides the Pluto notebooks), which might be more popular for some users.\nWe will also consider hosting a server in the future. However, due to the computational effort to train a RBM (requiring a few hours of compute on a standard CPU), we cannot offer such a service on-demand at the moment.\nLastly, we created a new Repository in Python: https://github.com/cossio/SamApp2024Py, that allows to train and sample RBM on RNA sequence families. This repository is now linked from our main repository, and should help users who prefer Python over Julia. We note that both Julia and Python notebooks can be executed in Google Colab, if remote execution is preferred (this also mentioned in our repository README now).\nBeyond the previous comments, this work's methodology appears sound, and its main conclusions and claims seem to be justified, particularly through its experimental data with 476 RBM-designed sequences. As an experimentalist myself, I appreciate that and congratulate the authors for doing so.\nWe appreciate the comment, and thank Reviewer 1 for her/his efforts to assess our work."}, {"round": 2, "reviewer_comments": {"reviewer_1": "While I agree that previously reported supervised learning models cannot be directly used to generate de novo functional sequences, they can be used downstream of a random mutagenesis loop to brute force verify de novo functional sequences, which has been previously reported in at least Angenent-Mari et al. In that sense, the task of de novo functional sequence generation can be achieved in different ways, which is why I fundamentally believe those precedents should be carefully discussed and not dismissed.\nWhile the authors did not fundamentally change many things in their original manuscript to address my concerns, I think they have done a good job of justifying their methodology and novelty.\nI have no further comments.", "reviewer_2": "The authors have generally done a good job at revising the manuscript. However, the lack of experimental data verifying that the designed molecules can really operate as riboswitches remains a major weakness. While it is good that they previously showed data for 476 sequences using SHAPE-MaP, these tests were done using oligo pools so it is not a major experimental feat and these experiments only showed changes in structure rather than modulation of transcription or translation. Testing a few of these sequences for actual riboswitch function using something as simple as in vitro transcription is not that challenging and really should be included in the paper.", "reviewer_3": "This revised version of the manuscript includes some changes and clarifications, but it does not address in my opinion the main issue with this manuscript which is that not providing sufficient evidence that their designed riboswitches would actually experience a conformational change, while claiming that their SHAPE-MaP analysis shows when a riboswitch synthetic sequences \"correctly undergoes a structural change in response to SAM\" (abstract).\nThe SHAPE-MaP data is at best inconclusive. For example, consider the two natural riboswitches examples described in Figures 5C and 5E. I would have a really hard time deciding that those two sets of reactivities (or delta_reactivities) correspond to the same functional molecule. In addition, if the SHAPE-MaP signal is supposed to distinguish paired from unpaired residues, there is nothing in the two profiles for 5C and %E what would make one guess that the structural change involves the opening of the P1 stem. I would find that bothersome, towards using SHAPE-MaP as a signal describing structural conformational changes for these riboswitch sequences. Yet, the subsequent analyses proceed by comparing these signals averaged over many different molecules.\nMoreover, how are those delta_reactivates calculated? line 263 describe those just as subtractions, but By visual inspection of 5C and 5E, it looks they might have been subjected to some \"re-scaling\" so that they do not have a magnitude larger than 1.5, which would make the obvious differences observed less pronounced. For example, for site 11 (SAM), in 5C the delta_reactivity is about 1 unit, while in 5E it is about 4 units, yet, the green delta_reactivity bars for cases are about identical.\nFurthermore, based on the SHAPE-MaP data, Figure 7 introduces the concepts of \"responsive\" and \"non-responsive\" riboswitches, which are not defined with precision anywhere. And while the application of that \"responsiveness\" concept to naturally occurring riboswitches gives a sense of false negatives, I don't think that there is an estimation anywhere of false positives. The selected synthesized riboswitches are valuable candidates, but further experimentation would be required to test whether they actually structurally behave like SAM riboswitches."}, "author_response": "We thank Reviewer 1 for her/his comments. In our new manuscript, we cite and discuss previous works (including Angenent-Mari et al.) with additional explanations compared to the previous version.\nWe have now added a Python repository, see https://github.com/cossio/SamApp2024Py. The README and Jupyter notebooks have enough comments to be easily usable and show how to train and sample RBMs on a RNA sequence family. This is now mentioned in the “Accompanying software” section.\nWe remark that training and sampling a generative model such as RBM demands computational power, and unfortunately we are unable to provide a server hosting a live demo to users due to resource constraints. As an alternative, we have decided to publish accessible examples of our approach as interactive Google Colab notebooks:\n● One in Python:\nhttps://colab.research.google.com/drive/1nOfFLWCwLy7a0aZ52cFHKUfF7erAMp5f?usp=sharing\n● One in Julia (which also uses GPU):\nhttps://colab.research.google.com/drive/1lfY5t6m-j8n19EXHLnV-lRBBfJ_jLk8y?usp=sharing\nThese examples do not require any local setup (they run on the web). We believe this will make it even easier for users to try our approach. This is now mentioned in the “Accompanying software” section of the current version of the manuscript.\nReviewer #2 (Remarks to the Author):\nThe authors have generally done a good job at revising the manuscript. However, the lack of experimental data verifying that the designed molecules can really operate as riboswitches remains a major weakness. While it is good that they previously showed data for 476 sequences using SHAPE-MaP, these tests were done using oligo pools so it is not a major experimental feat and these experiments only showed changes in structure rather than modulation of transcription or translation. Testing a few of these sequences for actual riboswitch function using something as simple as in vitro transcription is not that challenging and really should be included in the paper.\nWe certainly agree with Reviewer 2 that designing functional RNA riboswitches capable of regulating transcription/translation is an important objective. Our work is an intermediary step towards this objective: we are designing “only” the aptamer part of the riboswitch, responsible for the structural change of conformation upon SAM addition. However, this aptamer alone is not functional, it must act with the expression platform to acquire its regulation capability. Designing full riboswitches would require us to change our model and incorporate the expression platform. The P1 helix of the aptamer and the expression platform is expected to strongly interact (see answer to Reviewer 3 and Figure 1 below for additional analysis on P1 flexibility when considering the expression platform), so combining a generated aptamer with an unrelated expression platform is unlikely to result in a fully functional RNA. While we plan to model full riboswitches and then test their regulatory functions in future, we believe this is out of scope for the present work, as this would entail overcoming significant new modelling challenges, such as:\n● We need to collect a large enough training set of full riboswitch sequences to train the model. We have so far managed to collect 1306 sequences from the RiboD database (that we use in Supplementary Information, section “Stability of P1 helix in the context of the full riboswitch sequence”), which is insufficient to train a model, but plan to collect more in the future.\n● Full riboswitch sequences exhibit significant length variability, with hard-to-align regions. This poses challenges for our current RBM model, which is trained on multiple sequence alignments of homologous sequences. We believe that a different architecture, such as an RBM with a convolutional layer, could be used to model unaligned regions. But this would be a significant modification of the model that we leave for future work.\nWe emphasize that designing the aptamer alone is already a non-trivial task, as we are not aware of another study proposing a generative model capable of embedding two conformations plus the switching mechanism in a single RNA sequence.\nThis point is now made very clear at different points in the manuscript, see for instance the discussion. Let us also add that we have done more experiments in the new submission, with complementary chemical probing studies, based on DMS, which provide additional support for the ability of our designed RNA aptamers to respond to SAM in a similar way to natural aptamers.\nReviewer #3 (Remarks to the Author):\nThis revised version of the manuscript includes some changes and clarifications, but it does not address in my opinion the main issue with this manuscript which is that not providing sufficient evidence that their designed riboswitches would actually experience a conformational change, while claiming that their SHAPE-MaP analysis shows when a riboswitch synthetic sequences \"correctly undergoes a structural change in response to SAM\" (abstract).\nThe SHAPE-MaP data is at best inconclusive. Reviewer 3 is concerned by the fact that SHAPE-MaP measurements are intrinsically noisy, and might not provide enough evidence to support the existence of structural changes induced by SAM. As a general answer to this concern:\n- We fully agree that SHAPE reactivities are highly variable. In fact, we report in the manuscript (now Fig. 7) the histograms of reactivities at sites which are base-paired and unpaired according to the consensus secondary structure. The two histograms largely overlap, with slight but visible discrepancies in the (left and right) tails. This observation underlines that interpreting SHAPE results must be done with care.\n- To increase the signal-to-noise ratio we first resort to a very simple approach: we average the SHAPE signal over many molecules belonging to large classes. This is what we do for the classes of natural, RBM-generated and rCM-generated aptamers. Averaging over tens to hundreds aptamers allows us to reduce the statistical uncertainty over the reactivity profile, in such a way that the error bars become smaller than the differences between the mean profiles with and without SAM, see new Fig. 6. In addition, this approach shows that these differential profiles are quite similar in natural aptamers and RBM-generated ones (when the RBM score is high enough), with an average pattern of protection along the sequences compatible with the expected structural changes.\n- Obviously, this average approach cannot be used to study a single, given aptamer. We therefore introduce a probabilistic Bayesian analysis of the reactivity profile to detect significant response to SAM. Importantly and differently from previous analysis (e.g. Ipanema), we do not assume any model of the secondary structure in this analysis, but compute, for each molecule, a protection score that estimates the odd ratio of a set of hallmark sites (known to be relevant for structural changes from previous studies) to be all paired or all unpaired. Though such an odd ratio for a single site would be too much affected by noise and meaningless, considering a set with 24 sites substantially amplifies the signal-to-noise ratio and offers a much more reliable estimator of response to SAM for a single aptamer.\n- We have also carried out another series of experiments, based on DMS, an alternative chemical probing method. These extra data allow us to further reduce the experimental noise and to confirm our findings, see new Fig. 7.\nIn conclusion, while we agree with Reviewer 3 that the analysis of SHAPE data is delicate, we disagree with her/his statement that “the SHAPE-MaP data is at best inconclusive”. Detailed answers to the questions with additional information about the treatment of noise can be found below."}, {"round": 3, "reviewer_comments": {"reviewer_3": "The new version of this manuscript includes new data in the form of DMS chemical probing, and the authors had taken steps to address the questions raised by the reviewers. The manuscript has an important shortcoming by not having any experimental verification of riboswitch function for some of the proposed designs. I agree with reviewer #2 on that. The manuscript could be presented as providing a list of strong candidates and that is acceptable and has value. However the writing of the manuscript goes forward by implying that their analyses show evidence of functionality. For instance as read in the abstract: \"Designed molecules [...] display ~30% success rate of correctly undergoing a structural switch in response to SAM\".\nThis new version of the manuscript has added DMS data to support their claims that they have experimental evidence of functionality for some of the synthetic SAM riboswitches. That is not enough, as DMS has similar issues as SHAPE.\nThe manuscript includes assertions such as \"DMS confirms SHAPE results\" which is based on the comparison of Fig 5E to 5F. I am sorry but I do not see it.\nI do not see any report of their protection scores based on DMS or a comparison to SHAPE. Overall, there is no information either added or confirmed by DMS in this manuscript, at least as it is presented.\nI would like to justify my skepticism for their data by proposing a challenge: take Fig 5/Fig 10B-10C, select any cutoff for reactivity, or delta reactivity, or protection score, average or not, and report sen and ppv in predicting the correct structure of the riboswitch.\nWhile the author do not question that data, they are eager to question the value of alignments in detecting the conserved secondary structure. For instant in\n1099 \"I. SHAPE protection scores are in agreement with consensus secondary structure\", the evidence for such claim in Fig 10B, but they do question the Rfam RF00162 alignment with the sentence,\n1110 \"For many other sequences in the MSA, their actual behavior is at most hypothesized based on indirect evidence.\"\nThe RF00162 alignment is quite strong and accurate, and comprehensive in identifying all the base pairs (the pk and A-U lone pair included) and to single out the k-turn. This manuscript could also improve by using alignment information in the analysis of their candidates instead of dismissing it and relying solely on chemical probing. I suggested some approached along those lines in my very first review, that were not followed by the revisions.\nI am intrigued by Fig 2G which has an overview cartoon showing the typical histograms of SHAPE reactivities for paired vs unpaired residues. Importantly, it also has a small insert with similar histograms for \"total protection scores\" which looks really good. Unfortunately, there is not any place, that I have found (maybe I missed them, maybe Fig S25?), where those histograms are presented in detail for real data. Section I, mentioned in the caption for Figure 2G, describes the overall method to calculate protection scores, but not histograms of protection scores for actual data. Besides protection scores are calculated for \"hallmarks sites\" of motifs M? or total? There seem to be different uses throughout the paper. Table 1 describes as hallmark independent positions. In the definition of the score (page 22), the argument is \"a motif\" M that is a whole helix.\nThe protection scores still have the issue that the selection of motifs or hallmarks is arbitrary, and could be completely deceiving. Nevertheless, if you can show at least for the SAM riboswitch that the distributions of those protection scores for real ribo sequences increase the separation between paired/unpaired residues as much as Fig 2G reports, that should be one important point to make and use. (Maybe it is related to Fig S25? what does pair/unpaired for 1 site or 8 sites means?)"}, "author_response": "There is a confusion about the meaning of functional/functionality here. As explained in extenso in our answer to the second round of reviews, our objective is to design aptamers undergoing a structural switch in response to SAM, not full riboswitches. Our experimental and statistical analysis provide solid evidence for a ~30% success rate of the designed aptamers to undergo a structural switch in response to SAM. This reported success rate is supported by a careful statistical analysis, in particular of the significance of the protection-score analysis (Results Section C and Fig. 7E), which is affected by limited sequence coverage, see Fig. 10 in the manuscript. We have also estimated error bars on the success rate, see Table in Fig. 7G.\nAs discussed in the new version of the manuscript, both methods are standardly used to detect closed/open nucleotides in the RNA structures, as extensively documented in literature. SHAPE and DMS differ as DMS mainly detects single-strand A and C nucleotides, while SHAPE is sensitive to conformational flexibilities of all nucleotides (see description in Section B of Results). The reactivity profiles are of course ‘noisy’, as all experimental data, and their quality depends in particular on the depth of sequencing coverage as mentioned above. This is why we have introduced a sophisticated and automated Bayesian analysis of the SHAPE and DMS reactivity profiles.\nThese statements are erroneous. We would be grateful to Reviewer #3 to carefully read the revised version of the manuscript (and the detailed answers to the reviews):\n● In the main Fig. 7I we compare the outcomes of our protection-score statistical method to assess structural response to SAM using SHAPE experimental probing only to the one using both SHAPE and DMS on all tested molecules in the first round of experiments (natural and designed sequences).\n● We write on lines 530-534: ”SHAPE and DMS+SHAPE provide the same classification (responsive, non-responsive, or inconclusive) for about 86% of the aptamers. Among the remaining 14%, more than 12% are inconclusive for one of the two approaches, and SHAPE and DMS+SHAPE disagree on less than 2% of the aptamers only.”\n● We write line 299: “We report in Supplementary Fig. S33 the reactivity profiles r of the i,n,c same natural sequences as in Fig. 5A-D obtained with DMS.”\n● We write line 544: “A complete comparison of the analysis of the SHAPE and DMS data is reported in Supplementary Fig. S36.”\n● We agree with Reviewer #3 that we did not plot the protection DMS scores for the molecules as shown in main Fig. 7B for SHAPE. Here is this figure, which we have now included in Supplementary Information, Fig. S37. We have also added the counterpart of Fig. S25 for DMS data, see Fig. S38. See also the caption of main Fig. 7B for explanations with DMS instead of SHAPE.\nWe have shown, over the multiple of reviews, that we are willing to do further computational studies to improve the quality of the manuscript. But we do not understand what this “challenge” means. Our statistical pipeline does not rely on de novo prediction of structural conformations. What we do is to consider the two known structural conformations of the riboswitch in the presence and in the absence of SAM, identify a set of hallmark residues that are paired in one conformation and unpaired in the other, and use the SHAPE/DMS reactivity profile of each molecule to statistically decide if the tested aptamers (natural and designed) undergo or not structural switches upon SAM addiction. Determining the structure in one of the two conditions as suggested by Reviewer #3 would merely test the performance of folding softwares. This is not what we focus on in this manuscript.\nWe do not understand how Reviewer #3 can think we are “eager” to question the quality of the RFAM alignment. We are actually convinced that the RF00162 alignment is great! Our RBM model, trained on this MSA, is able to generalize and generate novel functional aptamer sequences. This would have not been possible if the RF00162 alignment were of poor quality. The reviewer has misunderstood the sentence in line 1110: we are simply writing that most natural sequences in the alignment had not been tested experimentally prior to our work and were assumed to undergo a structural switch only because they could be aligned along the RFAM profile (this is what we call indirect evidence). In our work, we test 201 natural sequences, among which 97 undergo a structural switch according to our protection-score analysis and 48 do not (the rest being statistically inconclusive), see Table in Fig. 7C.\nWhile the sentence above is not fully explicit, we presume that Reviewer #3 is alluding to the use of R-scape here, which she/he recommended in the first round of reviews. We did follow the suggestion in our first-round reply and carried out analyses with R-scape. Appendix E shows that RBM generated sequences exhibit statistical covariation compatible with natural sequences as reported by R-scape. Again, we would appreciate that Reviewer #3 takes time to carefully read our replies, in particular to her/his requests.\nReviewer #3 missed the important point that Figure 2G is not a cartoon, but is the histogram for real data. This was explained in detail in our answer to the second round of reviews.\nThroughout the paper, protection scores are computed either for all Hallmark sites, or for subsets of sites localized on specific motifs for local analyses (e.g., P1 helix only or Pk only in Fig. 8). The choice of the Hallmark sites is based on a careful comparison of the aptamer structures with and without SAM. We remove sites that are close to the extremities of the aptamers due to insufficient coverage with SHAPE. Moreover, SAM response at these Hallmark sites is also supported by extensive literature, detailed in Supporting Information, Appendix Q.\nFig. S25 is precisely doing what Reviewer #3 is suggesting: the figure reports the protection scores computed from experimental reactivity data collected on our natural aptamers. The panels of Fig. S25 illustrate the principle that combining data from several sites enhances statistical discriminability. Fig. S25B shows 8 sites (the number of sites in the pseudoknot) as an example. For completeness, we now add a panel with 24 sites (the number of Hallmark sites), which corresponds precisely to the inset in Fig. 2G. We also added the counterpart of Fig. S25 for DMS data, see Fig. S38."}]

Decision letter