Matan Alexander Churkin

Drory Retwitzer

Vladimir Reinharz

Yann Ponty

Jérôme Waldispühl

Danny Barash
email: dbarash@cs.bgu.ac.il

Matan Drory Retwitzer

Design of RNAs: comparing programs for inverse RNA folding

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The inverse RNA folding problem for designing sequences that fold into a given RNA secondary structure was introduced in the early 1990's in Vienna [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF]. Mathematically, much like the typical situation with inverse problems, it is not a well-posed problem by the standard definition of Hadamard, which makes it even more challenging to solve. As the well-known mathematician Andrey Tikhonov once noted, the class of ill-posed problems includes many classical mathematical problems and, most significantly, that such problems have important applications. Indeed, new emerging sub-fields that are of significant importance to a variety of functional RNAs [START_REF] Taft | Non-coding RNAs: regulators of disease[END_REF][START_REF] Hammann | Searching genomes for ribozymes and riboswitches[END_REF][START_REF] Strobel | RNA catalysis: ribozymes, ribosomes, and riboswitches[END_REF][START_REF] Breaker | Prospects for riboswitch discovery and analysis[END_REF][START_REF] Serganov | A decade of riboswitches[END_REF], such as RNA synthetic biology [START_REF] Isaacs | RNA synthetic biology[END_REF][START_REF] Chappell | A renaissance in RNA synthetic biology: new mechanisms, applications and tools for the future[END_REF] and RNA nanostructure [START_REF] Jaeger | TectoRNA: modular assembly units for the construction of RNA nano-objects[END_REF][START_REF] Guo | The emerging field of RNA nanotechnology[END_REF], are fast developing and are utilizing in their arsenal the methods for solving inverse RNA folding [START_REF] Mueller | Live attenuated influenza virus vaccines by computeraided rational design[END_REF][START_REF] Bindewald | Multi-stranded RNA secondary structure prediction and nanostructure design including pseudoknots[END_REF][START_REF] Afonin | In silico design and enzymatic synthesis of functional RNA nanoparticles[END_REF][START_REF] Fernandez-Chamorro | Designing synthetic RNAs to determine the relevance of structural motifs in picornavirus IRES elements[END_REF][START_REF] Dotu | Complete RNA inverse folding: computational design of functional hammerhead ribozymes[END_REF][START_REF] Findeiß | Design of transcription regulating riboswitches[END_REF][START_REF] Wachsmuth | Design criteria for synthetic riboswitches acting on transcription[END_REF][START_REF] Ben-Yehezkel | Rationally designed, heterologous S. cerevisiae transcripts expose novel expression determinants[END_REF].

A brute force approach that searches all the possible sequences is not a viable option because the number of sequences grows exponentially as 4n, where n is the length of the sequence, while the number of valid designs can be arbitrarily small. This upper bound can be refined by noting that paired positions have to form valid base pairs under the standard A-U, C-G, G-U base pairing scheme. This implies that 6 p/2 4 u sequences are compatible with a secondary structure having respectively u unpaired and p paired nucleotides. As a practical consequence, a typical 74 nts-long tRNA, including p = 40 paired and u = 34 unpaired ones, would require investigating ∼ 10 36 compatible sequences.

RNA inverse folding also has deep connections with theoretical evolutionary studies, where the sequence/structure relationship in RNA is a popular model for studying genotype/phenotype maps [START_REF] Schuster | From sequences to shapes and back: a case study in RNA secondary structures[END_REF][START_REF] Greenbury | Genetic correlations greatly increase mutational robustness and can both reduce and enhance evolvability[END_REF]. For instance, the identification of undesignable motifs [START_REF] Aguirre-Hernández | Computational RNA secondary structure design: empirical complexity and improved methods[END_REF] in empirical design studies immediately implies that only a negligible, exponentially decreasing on the length, proportion of secondary structures can be designed. Conversely, neutral evolution provides theoretical foundations for the practice of RNA design, and studies of the neutral network confirm a highly-variable numbers of admissible designs within structures of the same length [START_REF] Jörg | Neutral network sizes of biological RNA molecules can be computed and are not atypically small[END_REF]. It is interesting to note that the distribution of the neutral network [START_REF] Schuster | From sequences to shapes and back: a case study in RNA secondary structures[END_REF] could help us understand how to further develop efficient local search strategies to reach the target structure.

The approach to solve the inverse RNA folding problem by stochastic optimization relies on the solution of the direct problem using software available in RNA folding prediction web servers, e.g. the RNAfold server [START_REF] Hofacker | Vienna RNA secondary structure server[END_REF][START_REF] Lorenz | ViennaRNA Package 2.0[END_REF] or mfold/UNAFold [START_REF] Zuker | Mfold web server for nucleic acid folding and hybridization prediction[END_REF][START_REF] Markham | UNAFold: software for nucleic acid folding and hybridization[END_REF] as well as RNAStructure [START_REF] Mathews | RNA secondary structure analysis using RNAStructure[END_REF], by performing energy minimization with thermodynamic parameters [START_REF] Mathews | Expanded sequence dependence of thermodynamic parameters provides improved prediction of RNA secondary structure[END_REF][START_REF] Mathews | Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure[END_REF]. It should be noted that in principle, although far less popular in practice in the context of inverse RNA folding, other programs based on probabilistic models and posterior decoding that have been benchmarked in [START_REF] Puton | CompaRNA: a server for continuous benchmarking of automated methods for RNA secondary structure prediction[END_REF], e.g. PFold [START_REF] Knudsen | PFold: RNA secondary structure prediction using stochastic context free grammars[END_REF] and CentroidFold [START_REF] Sato | CentroidFold: a webserver for RNA secondary structure prediction[END_REF], can also be used to solve the direct problem. Initially, a seed sequence is chosen, after which a local search strategy is used to mutate the seed and apply repeatedly the direct problem of RNA folding prediction by energy minimization. Then, in the vicinity of the seed sequence, a designed sequence is found with desired folding Figure 1: The standard inverse RNA folding problem and the generalized inverse RNA folding problem that is shape aware and fragment-based (i.e., fragment selection enabled) are illustrated on the purine riboswitch aptamer in the middle (A). The predicted structure of an output designed sequence is shown on the right (C) for the standard inverse folding problem and on the left (B) for the generalized inverse RNA folding problem.

properties according to the objective function in the optimization problem formulation. Chronologically, algorithmic improvements that relate to this approach which was pioneered in Vienna's RNAinverse [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF] have been worked out in INFO-RNA [START_REF] Busch | INFO-RNA-a fast approach to inverse RNA folding[END_REF], RNA-SSD [START_REF] Andronescu | A new algorithm for RNA secondary structure design[END_REF] and NUPACK:Design [START_REF] Dirks | Paradigms for computational nucleic acid design[END_REF][START_REF] Zadeh | Nucleic acid sequence design via efficient ensemble defect optimization[END_REF]. Alternative methods to the adaptive random walk [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF] and the stochastic local-search [START_REF] Busch | INFO-RNA-a fast approach to inverse RNA folding[END_REF][START_REF] Andronescu | A new algorithm for RNA secondary structure design[END_REF] include genetic algorithms (belonging to the class of evolutionary algorithms) [START_REF] Dromi | Reconstruction of natural RNA sequences from RNA shape, thermodynamic stability, mutational robustness, and linguistic complexity by evolutionary computation[END_REF][START_REF] Taneda | Multi-objective genetic algorithm for pseudoknotted RNA sequence design[END_REF][START_REF] Lyngsø | Frnakenstein: multiple target inverse RNA folding[END_REF], an improved evolutionary algorithm [START_REF] Esmaili-Taheri | Evolutionary solution for the RNA design problem[END_REF][START_REF] Esmaili-Taheri | ERD: a fast and reliable tool for RNA design including constraints[END_REF], constraint programming [START_REF] Garcia-Martin | RNAiFold: A constraint programming algorithm for RNA inverse folding and molecular design[END_REF][START_REF] Garcia-Martin | RNAiFold 2.0: A web server and software to design custom and Rfam-based RNA molecules[END_REF] and ant-colony optimization [START_REF] Kleinkauf | antaRNA: Ant colony-based RNA sequence design[END_REF][START_REF] Kleinkauf | antaRNA-Multi-objective inverse folding of pseudoknot RNA using ant-colony optimization[END_REF]. On the methodological side, two separate ideas that were gradually developed and investigated were to sample the sequence space more efficiently and to include user-selected fragments in the design. The first idea started in [START_REF] Busch | INFO-RNA-a fast approach to inverse RNA folding[END_REF] and then by presenting a global sampling approach named RNA-ensign [START_REF] Levin | A global sampling approach to designing and reengineering RNA secondary structures[END_REF], followed by a weighted sampling algorithm called incaRNAtion [START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution[END_REF] that is a glocal methodology combining the global sampling approach with local search strategies. The second idea started by generalizing the inverse RNA folding problem to include RNA designed sequences that are predicted to fold into a prescribed shape [START_REF] Dromi | Reconstruction of natural RNA sequences from RNA shape, thermodynamic stability, mutational robustness, and linguistic complexity by evolutionary computation[END_REF][START_REF] Avihoo | RNAexinv: An extended inverse RNA folding from shape and physical attributes to sequences[END_REF], a utility named RNAexinv that considers the coarse-grain tree graph [START_REF] Shapiro | An algorithm for comparing multiple RNA secondary structures[END_REF] for shape representation (or potentially abstract shapes [START_REF] Giegerich | Abstract shapes of RNA[END_REF]). It culminated with a method called RNAfbinv [START_REF] Weinbrand | RNAfbinv: An interactive Java application for fragment-based design of RNA sequences[END_REF] that allows a user-selected prescribed fragment to be preserved exactly by secondary structure (in addition to its shape) whereas the rest of the structure is designed by the generalized shape-based approach. These two ideas were recently merged into an RNA design webserver called IncaRNAfbinv [START_REF] Retwitzer | IncaRNAfbinv: A webserver for the fragment-based design of RNA sequences[END_REF] that provides more flexibility in the design as compared to the aforementioned methods. An example of the benefit of IncaRNAfbinv is illustrated in Figure 1 where it is shown that the designed sequence on the left (Figure 1(B)) is a feasible purine riboswitch candidate as much as the designed sequence on the right (Figure 1(C)), both containing the essential nucleotides for either guanine or adenine binding, but the Figure 1(B) solution cannot be reached to-date by other programs besides IncaRNAfbinv because its secondary structure is different than the input structure depicted in Figure 1(A), although its tree graph shape is the same.

As was mentioned above and exemplified in the shape aware capability, different computational frameworks for the inverse RNA folding have been implemented in the various programs. In addition to the strict definition that restricts solutions to those sequences whose mean-free-energy structure is exactly the target structure, more relaxed frameworks like ensemble defect optimization in NUPACK have been introduced and developed [START_REF] Dirks | Paradigms for computational nucleic acid design[END_REF], in addition to minimizing the classical cost function given by the "structure distance" between the structure of the test sequence and the target structure [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF]. For a candidate sequence and a given target secondary structure, the ensemble defect is the average number of incorrectly paired nucleotides at equilibrium evaluated over the ensemble of unpseudoknotted secondary structures [START_REF] Zadeh | Nucleic acid sequence design via efficient ensemble defect optimization[END_REF]. It could also be possible, although not currently done by any available software, to add prescribed kinetic properties (fast folding, absence of kinetic traps) to the objectives of design, as explored in earlier studies [START_REF] Dirks | Paradigms for computational nucleic acid design[END_REF]. As will be noticed in the quantitative comparison performed in the continuation, these different computational frameworks will also make it impossible to draw a conclusion as to which program is better for use based on a benchmark. Instead, when facing an application that requires inverse RNA folding, the goal is to try and identify which is the most suitable program for each case. Some inverse RNA folding programs are geared towards more specific biological problems, compared to the ones mentioned up until now that are general in their application scope. For example, nanostructure design including pseudoknots was performed with Nanofolder for the case of multi-stranded RNA secondary structure [START_REF] Bindewald | Multi-stranded RNA secondary structure prediction and nanostructure design including pseudoknots[END_REF]. In another example, for designing the most stable and unstable mRNA sequences which code for a target protein, an algorithm was developed in [START_REF] Cohen | Natural selection and algorithmic design of mRNA[END_REF]. On a similar topic, in [START_REF] Terai | CDSfold: an algorithm for designing a protein-coding sequence with the most stable secondary structure[END_REF], an algorithm for designing a protein-coding sequence with the most stable secondary structure called CDSfold is provided. For designing multiple-target artificial miRNAs, a tabu search was used in conjunction with biochemical considerations in [START_REF] De Guire | Designing small multiple-target artificial RNAs[END_REF]. For designing RNA sequences that fold into multiple target structures, which makes it possible to efficiently design multi-stable RNA sequences, a program called RNAdesign was introduced in [START_REF] Zu Siederdissen | Computational design of RNAs with complex energy landscapes[END_REF]. Another specialized application, given here for completeness, is to allow game players to propose a set of rules for RNA design, as part of the Eterna project. Based on experimental results, Eterna players came up with a set of design rules and EternaBot was developed to design a sequence based on those rules [START_REF] Lee | RNA design rules from a massive open laboratory[END_REF]. Finally, for the problem of fixed backbone 3D design of RNA, the webserver RNA-redesign was put forth [START_REF] Yesselman | RNA-Redesign: a web-server for fixed backbone 3D design of RNA[END_REF]. To the best of our knowledge, this is the first implementation of an RNA design program that considers tertiary structure, although it is a local design. In this respect, in secondary structure, sequences that are generated by point mutations performed on an input sequence to optimize a certain objective function may also be considered a design procedure in the weaker (local) sense. Example of such programs are RNAmutants [START_REF] Waldispühl | Efficient algorithms for probing the RNA mutation landscape[END_REF] that employs an efficient sampling procedure based on the Boltzmann-weighted ensembles of mutants and RNAmute [START_REF] Churkin | An efficient method for the prediction of deleterious multiple-point mutations in the secondary structure of RNAs using suboptimal folding solutions[END_REF] that utilizes suboptimal solutions of the RNA free energy minimization prediction [START_REF] Zuker | On finding all suboptimal folding of an RNA molecule[END_REF][START_REF] Wuchty | Complete suboptimal folding of RNA and the stability of secondary structures[END_REF] for simulating only selective mutations from all possible ones. These procedures could potentially also be integrated into RNA design tools in the future, as was already done by incorporating RNAmutants ingredients into incaRNAtion [START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution[END_REF].

Inverse RNA folding programs are not only used for the design of artificial sequences to mimic natural ones. Recently, they have been used for the detection of novel naturally occurring RNAs as a pre-processing step before sequence-based searches [START_REF] Dotu | Using RNA inverse folding to identify IRES-like structural subdomains[END_REF][START_REF] Retwitzer | An efficient minimum free energy structurebased search method for riboswitch identification based on inverse RNA folding[END_REF]. In both of these separate works, they have shown to find attractive candidates for naturally occurring RNAs that are not available in RFAM [START_REF] Nawrocki | Rfam 12: Updates to the RNA families database[END_REF] and are missed by standard methods for RNA detection. This approach can well contribute to ongoing efforts aiming for de novo discovery of structured non-coding RNAs in genomic sequences [START_REF] Ruzzo | De novo discovery of structured ncRNA motifs in genomic sequences[END_REF].

Details of use

Installing the different programs or accessing the webservers is not a difficult task. They all contain a ReadMe file or a manual that is easy to follow with no prior knowledge assumed. However, not all programs provide both webserver and source code. Table 1 lists the various RNA design programs according to four categories: general purpose programs, shape aware programs, adaptive sampling programs and specialized programs. It then indicates the availability of webserver or source code for each program, including extended features and general remarks that relate to their capability, use, or strategy without providing more details on their specific methodology.

From all programs listed in Table 1, we picked five programs for further description on their details of use based on the following selection criteria: the programs have both a webserver and a source code available, and they were already used in the literature in a biological meaningful way by either a "wet-lab" experiment or the identification of a putative new non-coding RNAs. These criteria are of interest to practitioners who are considering the use of RNA design programs. The selection yielded the programs RNAinverse, RNAiFold, NUPACK and IncaRNAfbinv. In addition, the antaRNA program was added because it was published quite recently without a chance yet for practical use but it is considered promising as can also be observed by the program comparisons provided in next section and its overall strategy that allows much flexibility. Finally, although Nanofolder could not be added because of a lack of source code and our recommendation for the interested user would be to contact its authors [START_REF] Bindewald | Multi-stranded RNA secondary structure prediction and nanostructure design including pseudoknots[END_REF], it is worthwhile noting the significant practical experience that has been accumulated by Nanofolder as a specialized program for RNA nanostructure design. There are sequence design rules implemented in Nanofolder that have been formulated based on the concept of samelength sequence fragments called "critons" [START_REF] Bindewald | Multi-stranded RNA secondary structure prediction and nanostructure design including pseudoknots[END_REF], which have been successfully applied beforehand to the design of DNA nanostructures. These special rules with employed penalty scoring terms have been formulated to avoid unstable RNA designs and optimize the designed sequences. User experience with the five selected programs was performed on two example input secondary structures in dot-bracket notation (the first is a toy-problem, an artificial structure; the second is the structure of the guanine-binding riboswitch aptamer, a natural structure):

• • • RNAiFold [42, 43] • • • Experience in biology "wet-lab" RNAinverse [1] • • First program developed; Experience in biology "wet-lab" NUPACK [35, 36] • • Optional multi-stranded target struc- tures; experience in biology "wet-lab" INFO-RNA [33] • • RNA-SSD [34] • Frnakenstein [39] • • ERD [40, 41] • • • MODENA [38] • • • Shape aware IncaRNAfbinv [52] • • Fragment selection enabled; Experience in RNA detection RNAfbinv [51]
• Fragment selection enabled RNAexinv [START_REF] Avihoo | RNAexinv: An extended inverse RNA folding from shape and physical attributes to sequences[END_REF] • No user-selected fragment

Adaptive sampling

IncaRNAfbinv [START_REF] Retwitzer | IncaRNAfbinv: A webserver for the fragment-based design of RNA sequences[END_REF] • • Global-local approach; Experience in RNA detection IncaRNAtion [START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution[END_REF] • Global-local approach RNA-ensign [START_REF] Levin | A global sampling approach to designing and reengineering RNA secondary structures[END_REF] • Global approach

Specialized Nanofolder [12] • • Nanostructures; multi-stranded RNA; experience in biology "wet-lab" CDSfold [54] • • Design of protein-coding seq RNAdesign [56] • • EternaBot [57]
• Design rules set by Eterna players RNA-redesign [START_REF] Yesselman | RNA-Redesign: a web-server for fixed backbone 3D design of RNA[END_REF] • 3D: fixed backbone Table 1: A tabular overview with some basic information about the various RNA design programs.

RNAinverse in detail

RNAinverse [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF] was the first program developed for RNA design. It is available as a webserver at http://rna.tbi.univie.ac.at/cgi-bin/RNAinverse.cgi and as a standalone version in the Vienna RNA package [START_REF] Lorenz | ViennaRNA Package 2.0[END_REF]. The algorithm employs an adaptive random walk to minimize base pair distance. The distance is calculated by comparing the minimal energy folding of a mutated sequence (its predicted structure) to the target structure. To avoid folding the entire sequence, small substructures are optimized and then elongated. The algorithm also supports designing sequences which are more probable based on the partition function. Those sequences may be more stable but mostly differ from the target structure.

The server receives a secondary structure in dot bracket notation. An optional start sequence can be inserted; any lower case letter in the sequence will be conserved in the final result. The server also supports multiple energy models, folding temperature and number of sequences to generate. Once the form is submitted the result page will appear with a list of designed sequences and the calculated minimum energy for them. It may also show designed sequences that did not match the exact structure and their base pair distance away. Another set of results that will appear below that are designed sequences using the partition function if selected. The webserver also shows links to RNAfold [START_REF] Lorenz | ViennaRNA Package 2.0[END_REF] for extensive information on a specific result. The command line used to run the design in the stand-alone version is also written.

The stand-alone version of RNAinverse is part of the Vienna RNA package. The package is a C code library that includes several stand-alone programs. This means RNAinverse can be accessed both from command line and through the package API. The simple API allows the user to design and test the sequences for individual purposes directly through C. The stand-alone version is simple and examples can be generated by using the webserver as discussed above.

RNAiFold in detail

RNAiFold [START_REF] Garcia-Martin | RNAiFold: A constraint programming algorithm for RNA inverse folding and molecular design[END_REF][START_REF] Garcia-Martin | RNAiFold 2.0: A web server and software to design custom and Rfam-based RNA molecules[END_REF] is a well-established program available as a webserver at http://bioinformatics.bc.edu/clotelab/RNAiFold/ It employs constraint programming and is available in two modes. RNA-CP Design take as simplest input a secondary structure and returns up to 50 sequences or the maximum found in 2 hours, optimized for one of three different criteria: yielding the target as the minimum free energy (MFE) structure, the free energy or the ensemble defect. RNAiFold ensures the optimality of the solutions. Additional constraints can be provided such as a target amino acid sequence, and limits on the amount of GC, AU and GU base pairs as a list of admissible and forbidden base pair. The nucleotide distribution and energy model are customizable. A special mode, RNA Synthetic Design, leverages RFAM [START_REF] Nawrocki | Rfam 12: Updates to the RNA families database[END_REF] to add constraints from sequence conservation. A drop down menu allows the selection of any family and automatic extraction of the constraints. The consensus structure can be automatically selected as target. The constraints from RNA-CP Design are also available.

The output presents for each resulting sequence a number of different statistics such as its GC content, energy and entropy, and its amino acid sequence. It also provides a link to MPMATH [START_REF] Altschul | Gapped MPMATH and PSI-MPMATH: a new generation of protein database search programs[END_REF] the resulting sequence. Structure (2) was tested on the webserver. It took a short time before the query was allowed to run and returned 23 sequences after 2 hours. A version can be downloaded in source code or binaries for Linux or Mac OS X, the source requires the Vienna Package as the open source Google optimization library OR-Tools. The binary was tested on Ubuntu 12.04 while the Mac binaries still need to be updated for OS 10.11.5. It provides a simple input by command line argument or through a file in a custom format to fine-tune a few more arguments of the algorithm itself, and provides a similar output.

antaRNA in detail

antaRNA [START_REF] Kleinkauf | antaRNA: Ant colony-based RNA sequence design[END_REF][START_REF] Kleinkauf | antaRNA-Multi-objective inverse folding of pseudoknot RNA using ant-colony optimization[END_REF] is a very recent program available since 2015, as a webserver at http://rna.informatik.uni-freiburg.de/antaRNA/ It employs ant colony optimization to allow the design of structures, allowing pseudoknotted structures using hardcoded energy parameters. A sequence constraint can be provided in the IUPAC format. A visualization of the secondary structure with the sequence constraint is dynamically shown. Additionally a target GC constraint can be set. The parameters of the ant colony algorithm can be modified through advanced options.

Each output sequence, up to a 100, is shown with its MFE structure and its distance from the target GC, target structure and target sequence. A click on a sequence will show a visualization of its secondary structure with the sequence embedded. Testing for structure (2) while requesting 100 sequences took a few minutes on the webserver. All the results contained all the requested base pairs, a few additional base pairs were sometimes present. A convenient link to download all the sequences in FASTA format with or without their MFE is provided.

The python script requires the Vienna RNA Package and provides the same options as the webserver. For pseudoknots prediction a finer control is provided requiring the user to have installed one of the programs RNAShapes studio [START_REF] Janssen | The RNA shapes studio[END_REF] or HotKnots [START_REF] Ren | HotKnots: heuristic prediction of RNA secondary structures including pseudoknots[END_REF] or IPKnot [START_REF] Sato | IPKnot: fast and accurate prediction of RNA secondary structures with pseudoknots using integer programming[END_REF]. All the options from the webserver are present and must be given through command line arguments. The script removes the limit of sequences sampled.

NUPACK in detail

NUPACK [START_REF] Dirks | Paradigms for computational nucleic acid design[END_REF][START_REF] Zadeh | Nucleic acid sequence design via efficient ensemble defect optimization[END_REF] is a recent program developed in 2011. It is available as a webserver at http://www.nupack.org/design/new Its objective is to minimize the ensemble defect for a pseudoknot free structure. Its interface allows the user to specify a target secondary structure as a preference for DNA or RNA. An interesting feature of NUPACK is the ability to define unwanted motif by providing a list of forbidden sequence motifs, in IUPAC format. The webserver allows to choose between two energy models, Turner95 and Mathews99, dangles, and setting the temperature. A maximum of 10 sequences can be designed concurrently. The output presents the designed sequences. Tested on structure (2), in a few minutes 10 sequences were generated. An analysis of the sequences can be launched immediately to compute the MFE and base pair probabilities. A range of temperature can additionally be provided for this step. A link is provided to download the MFE secondary structure representation and the base pair probabilities. NUPACK is also available as command line software and was tested on a MAC OS X 10.11.5. Some options are not available on the webserver, such as setting a seed sequence or specifying the concentration of salt and magnesium. A few parameters of the algorithm can also be modified from the command line as its random seed, which is necessary to generate different sequences.

IncaRNAfbinv in detail

IncaRNAfbinv [START_REF] Retwitzer | IncaRNAfbinv: A webserver for the fragment-based design of RNA sequences[END_REF] is a very recent program for fragment-based design. It is available as a webserver at https://www.cs.bgu.ac.il/IncaRNAfbinv/ The webserver combines two base applications: incaRNAtion [START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution[END_REF] and RNAfbinv [START_REF] Weinbrand | RNAfbinv: An interactive Java application for fragment-based design of RNA sequences[END_REF]. Both applications are available as a stand-alone client. RNAfbinv uses simulated annealing with a 4-nt look ahead local search function. The function includes biologically meaningful constraints such as sequence constraints, fragment based design, and a variety of optional features. The resulting sequences fit a coarse-grained tree graph shape of the original target structure, thus allowing for flexibility. incaRNAtion augments the local search method. It uses a global sampling approach and weighted sampling techniques. The sequences generated by incaRNAtion are used as seed sequences for the local search. The incorporation of those seeds forces highly distributed results and better control of GC content.

The webserver takes as input the target secondary structure in dot bracket representation. It then converts it to coarse-grained tree graph shape for future comparison. After inserting a structure, an image will appear, as well as a list of structural motifs from which the user can select a desired motif. Additional optional parameters are: sequence constrains, target fold energy, mutational robustness, and GC content. Submission of the query leads to a web page showing design progress. Once all the results are ready, a list will appear with the designed sequences as well as their predicted structure, folding energy, mutational robustness, GC content, and distance to the original structure. The distances are calculated for both base pairs and structural motifs. Finally, a link is available for an image of the designed structure.

The stand-alone versions of both tools can be run locally, links can be found in the webserver. incaRNAtion is a simple to run Python script, thus requiring Python installed and recommends adding the MPMATH library for long sequences. It receives as input a file containing the target structure and an optional multiple sequence alignment (MSA). The user is also required to enter a value between 0 and 1, where 1 means only to regard the structure and 0 the MSA. Additional variables are available for GC content control and specific sequence constraints. A single run will generate a large amount of sequences; a minimal value for the number of outputs can also be set as input. The output is a list of sequences separated by lines.

RNAfbinv is a C application but it is also available wrapped in a Java interface. Once the Java application is running, the user must first insert the secondary structure. The user then has an option to select a specific motif to preserve as is. After inserting the structure, the user arrives to a new screen where additional control variables are available such as target folding energy and mutational robustness. The results appear in a new screen as a list including the base pair distance from the input structure.

Using the programs

The main task for using these programs is to insert an RNA secondary structure into one of them and generate sequences accepting this target structure as the MFE structure, with possible generalizations that are closely related to this framework. All programs offer the possibility of additional parameters to be chosen by the user, with default values displayed at the beginning. Some programs offer more flexibility in their constraints than others, which is indicated in Table 1 for some selected features that are shared by several programs and are general in scope. As final output, all programs offer description of parts of the analysis of designed sequences as well. RNAiFold displays the results in one webpage per solution, seemingly in the order generated and selectable through a drop down menu. It presents the MFE structure of each sequence, as an ensemble of statistics, and provides an option to MPMATH them. antaRNA directly displays each generated sequence with its MFE structure. A click on each of the solutions creates a figure of the secondary structure with the sequence. Any constraint violation is represented in red in the figure. The list of sequences can be downloaded in FASTA. NUPACK presents each sequence by increasing ensemble defect. Each sequence has a link to the analysis tool, which computes the MFE structure and base pair probabilities. In addition, the textual output provided in all these programs is substantially contributing to the analysis, as an essential step before the graphical output. In general, most of the programs are user friendly for the novice. Especially the programs that have a webserver capability (see Table 1) can also be worked out by a non-specialist user along with the corresponding manuals and instructions that are available in the websites.

Discussion

The programs listed in Table 1 have been developed in the past several years, following the first program named RNAinverse that was introduced more than twenty years ago, and offer some interesting prospects for RNA sequence design. They all in one way or the other rely at present on thermodynamic parameters corresponding to the nearestneighbor model and therefore structures that are known to be well predicted by energy minimization, for example the secondary structure of the guanine-binding riboswitch aptamer that is illustrated in the second test case example of previous section and in Figure 1, are the best to work with as inputs to these programs in order to achieve reliable results. Though exceptional cases exist, in general the upper range estimate Figure 2: Histogram comparison between the six selected programs available in Table 2 for the example test case that is designated as (1) in the Details of Use Section and for which the runtimes are reported in the first column of Table 2.

for the sequence length that these programs are useful for is around 150 nt [START_REF] Lange | Global or local? Predicting secondary structure and accessibility in mRNAs[END_REF]. It is expected that in future, having more experimental structures elucidated, the number of RNA sequences with a well-predicted secondary structure by energy minimization techniques will grow significantly and more biological systems involving RNAs will be designed by the aid of these programs.

Runtime can be a critical issue concerning the usage of these tools. A runtime comparison of six programs is provided in Table 2 for the two test case examples that were provided in the previous section in dot-bracket notation (the first is a toy-problem, the second is the structure of the guanine-binding riboswitch aptamer). The times reported are in minutes. Standard parameter values were used in the comparison. Because runtimes are measured in downloadable source code and cannot be measured in programs that require user interactive intervention such as with IncaRNAfbinv and RNAfbinv, we replaced IncaRNAfbinv that was discussed in the previous section by the simpler and less developed program RNAexinv. The justification is that RNAexinv is shape aware (preserving the same coarse-grain tree graph shape in the output as in the input) without the user's interactive selection of a fragment for preserving its secondary structure exactly like in IncaRNAfbinv, therefore RNAexinv contains the shape aware feature itself for inclusion in the comparison. RNAexinv is still much slower than the rest of the programs because it solves a more general inverse RNA folding problem that is shape aware. By our past experience, the programs RNAfbinv and IncaRNAfbinv are about 10% more computationally expensive than RNAexinv. Additionally, we inserted the program INFO-RNA because it is known to be the most computationally efficient among all programs, as is also observed in Table 2. Correspondingly, for the two test cases measured in Table 2 by 1000 runs, a histogram is plotted in Figure 2 for the first test case Figure 3: Histogram comparison between the five selected programs available in Table 2 for the example test case that is designated as (2) in the Details of Use Section and for which the runtimes are reported in the second column of Table 2. and in Figure 3 for the second test case in order to examine how far away the predicted structures of the designed sequences are from the input structure that is initially given. The distance between the two secondary structures was measured by the RNAdistance routine available in the Vienna RNA package that calculates by default the tree edit distance.

While it is impossible to draw conclusions from Table 2 and the associated Figures 2,3 as to which program is better for use, because a program such as RNAexinv belonging to the shape aware category (see Table 1) is expected to be much slower along with histograms that are more wide-spread compared to the rest of the programs in a notably beneficial way for its purpose, some trends can be observed. For example, INFO-RNA is the quickest as expected, but its histograms are more widespread compared to most of the other programs and this correlates with the known result that sequences designed with INFO-RNA tend to be more biased to high GC-contents [START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution[END_REF]. In contrast, the newly introduced antaRNA program from the same laboratory is both relatively fast and achieving histograms with results that are very close to the input structure. RNAiFold is also showing some fairly balanced outcomes between efficiency and proximity of the results to the input structure. Finally, RNAinverse shows impressively that although it was written more than 20 years ago and it features less constraints compared to the newer programs, it is still both fast and faithful to the input structure.

The above comparison is by no means exhaustive and can be supplemented by the additional references [START_REF] Garcia-Martin | RNAiFold 2.0: A web server and software to design custom and Rfam-based RNA molecules[END_REF][START_REF] Kleinkauf | antaRNA: Ant colony-based RNA sequence design[END_REF][START_REF] Kleinkauf | antaRNA-Multi-objective inverse folding of pseudoknot RNA using ant-colony optimization[END_REF][START_REF] Retwitzer | IncaRNAfbinv: A webserver for the fragment-based design of RNA sequences[END_REF]. These references could benefit a reader interested in the topic of run time comparisons, design capabilities and properties of the output sequences produced by each method. The RNAiFold webserver article [START_REF] Garcia-Martin | RNAiFold 2.0: A web server and software to design custom and Rfam-based RNA molecules[END_REF] contains a section on comparison with other software. Because this article does not include the Program Time (minutes) Figure 2 Time (minutes) Figure 3 antaRNA most recent methods antaRNA and IncaRNAfbinv, the interested reader can find more information about the performance of these algorithms in [START_REF] Kleinkauf | antaRNA: Ant colony-based RNA sequence design[END_REF][START_REF] Kleinkauf | antaRNA-Multi-objective inverse folding of pseudoknot RNA using ant-colony optimization[END_REF] and [START_REF] Retwitzer | IncaRNAfbinv: A webserver for the fragment-based design of RNA sequences[END_REF], respectively. It should also be noted that the use of tree edit distance to the target structure in Figure 3 as a performance measure may not consider that some of the methods included do not necessarily use the same energy model and dangle treatment as used herein.

The one used herein for computing the tree edit distance is the Turner 2004 model [START_REF] Mathews | Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure[END_REF] included in the Vienna RNA package 2 [START_REF] Lorenz | ViennaRNA Package 2.0[END_REF]. While NUPACK results could be moderately accurate in any case since ensemble defect optimization mitigates the slight differences between energy models, the results of RNAinverse, RNAiFold, INFO-RNA and the rest of the programs could be affected by the energy model of choice.

New Prospects: designed RNAs for structure-based search

As was mentioned in the Introduction, a major new application of inverse RNA folding programs is the discovery of novel, structured and functional RNAs in transcriptomic data. We briefly describe the concept and refer the interested reader to [START_REF] Dotu | Using RNA inverse folding to identify IRES-like structural subdomains[END_REF][START_REF] Retwitzer | An efficient minimum free energy structurebased search method for riboswitch identification based on inverse RNA folding[END_REF] for more information.

Sequence-based search tools like MPMATH [START_REF] Altschul | Gapped MPMATH and PSI-MPMATH: a new generation of protein database search programs[END_REF] have been used extensively for the detection of novel RNAs of interest, such as riboswitches, in newly sequenced data. They are easily available, highly efficient, and can partially address this task. However, when the search is restricted to only sequence-based considerations, it is rather limited. The idea to augment MPMATH search with inverse RNA folding for including structure-based considerations has been developed independently for identifying IRES-like structural subdomains [START_REF] Dotu | Using RNA inverse folding to identify IRES-like structural subdomains[END_REF] and riboswitch aptamer domains [START_REF] Retwitzer | An efficient minimum free energy structurebased search method for riboswitch identification based on inverse RNA folding[END_REF], where in the first reference the findings were also verified experimentally and in the second reference the experimental verifications are ongoing. In both of these works, this strategy has been shown to yield attractive candidates that are beyond the reach of well-established methods like Infernal [START_REF] Nawrocki | Infernal 1.1: 100-fold faster RNA homology search[END_REF]. Consequently, an idea was even suggested by the authors of Infernal to augment their own tool by the inverse RNA folding pre-processing step. Various combinations should be tried and in any case, it is expected that in the future, inverse RNA folding would become useful not only for the design of synthetic RNAs but also for the search of naturally occurring RNAs by the use of designed RNAs as a pre-processing step.

Concluding remarks

The various programs, especially the ones who are gaining experiences in biological meaningful problems and are being improved as a consequence by updated versions, should best be examined along with the constraints they allow and their orientation purposes. There are already several programs that were described in detail and offer both a webserver implementation and source-code availability, along with a proven experience in biological meaningful problems. Other programs should strive to achieve these goals. Practitioners should then select which program is more suitable for their needs according to the specific constraints and capabilities that are advertised in each one of the programs. department at Ben-Gurion University. His research interests include computational biology, RNA structure predictions and RNA structure probing by experimental methods.

Vladimir Reinharz is a computational scientist and a PhD graduate of the computer science department at McGill University. He is the recipient of an Azrieli and FQRNT postdoctoral fellowships to be pursued at Ben-Gurion University. His research interests include computational biology, RNA structure predictions, and scientific computing.

Yann Ponty is a computational scientist and a faculty member of the bioinformatics group in the computer science department (LIX) of École Polytechnique. His research interests include computational biology, RNA folding and visualization, combinatorics and discrete algorithms.

Jérôme Waldispühl is a computational scientist and a faculty member in the computer science department at McGill University. His research interests include computational biology, RNA and Protein structure predictions, as well as crowd-computing systems.

Danny Barash is a computational scientist and a faculty member in the computer science department at Ben-Gurion University. His research interests include computational biology, RNA structure predictions, and scientific computing.

incaRNAfbinv

incaRNAfbinv is a combination of two seperate programs, both have specific command line interfaces. It is recommended for most users to use the webserver as it already combines the two interfaces.

incaRNAtion generates the seed sequences later inserted into RNAfbinv. To run it, a Python distribution must be installed. The command line interface includes many fine tuning parameters while the main structure input must be inserted in a file. The input file must contain a target structure. In addition to the structure, a multiple sequence alignment (MSA) may be added to allow for sequence information. For the average user, the following are the ones that will be most used.

-d The path for an input file containing secondary structure and optional MSA.

-a A number between 0 and 1 used by the algorithm as a weight. 1 takes into account on the structure while 0 only considers the MSA.

-m Maximum penalty for an invalid pair.

-s gc This is followed by 2 numbers. The first, between 0 and 1, forces a given GC content while se second show the minimal number of output sequences required.

-gc max err A number between 0 and 1 with the maximal GC difference between the output sequences and the requested number. 0.1 by default. The output will be a list a seed sequences. Those seed sequences can later be insrted into RNAfbinv using the command line.

To use the incaRNAtion seeds, download the RNAfbinv extended version. The package includes a java GUI interface. The command line option allow for the same options as the GUI version. For the average user, the following are the ones that will be most used.

-i The number of simulated annealing iterations for a single sequence design.

-t Look ahead depth: The maximum number of consecutive mutations that generate a lower score sequence possible before a single simulated annealing iteration is over.

-c A starting sequence. This is where the incaRNAtion seed should be inserted.

To generate a single sequence, for the structure ((((...(((....)))...((((....))))...))))

allowing any sequence , starting from a given incaRNAtion seed, aiming at -23 dG (Kcal/mol) and target mutational robustness 0.8 given 50 iteration and a 4 nucleotide look ahead depth, the command should be a followed:

./RNAexinv -i 50 -t 4 -c <incaRNAtion seed>

Once the programs starts the following parameters should be inserted:

((((...(((....)))...((((....))))...)))) \ -23 \ 0.8

RNAiFold

The command line interface of RNAiFold has over 50 options allowing for an extremely fine tuning of the desired output. For the average user, the following are the ones that will be most used. There is two way to enable those options, or through a file, where the option name is on a line preceded by a # instead of a -, followed on the next line by the desired option. Usually, the option can be simply given as argument on the command line.

-RNAscdstr The target structure. Multiple target can be set, they must be on the same line separated by the pipe | symbol. The structures must have the same length.

-RNAseqcon The admissible sequences, in IUPAC format. It must be one string the same length as the structure.

-maxGCcont The maximal GC content admissible in the sequences.

-minGCcont The minimal GC content admissible in the sequences.

-TimeLimit The amount of time allowed to run (default 600 seconds).

-MAXsol The maximum number of solutions to be reached under the time limit Therefore, to find a maximum of 50 solutions, under an hour, for the structure ((((...(((....)))...((((....))))...))))

allowing any sequence with a mandatory GC base between the first and last nucleotide, and with a GC content between 60% and 70%, the command would be as follows:

./RNAiFold -RNAscdstr "((((...(((....)))...((((....))))...))))" \ -RNAseqcon GNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNC \ -maxGCcont 70 -minGCcont 60 -MAXsol 50 -TimeLimit 3600

antaRNA

The distribution of antaRNA as a Python2.7 executable and all options are given as arguments. A similar ensemble of constraints exists.

-Cstr The target structure in the dot bracket notation. A fuzzy notation can be used to define blocks allowed to base pair together using any lowercase and uppercase letter.

-Cseq The admissible sequences, in IUPAC format. It must be one string the same length as the structure.

-tGC Target GC content, in [0, 1], which also serves as a minimum.

-tGCmax Maximal GC content admissible in the sequences.

-tGCvar Variance (σ 2 ) in the case of normal distribution, -tGC serves as the expected value µ.

-t The amount of time allowed to run (default 600 seconds).

-n Number of solutions to be produced.

Therefore, to find a maximum of 50 solutions, under an hour, for the structure ((((...(((....)))...((((....))))...))))

allowing any sequence with a mandatory GC base between the first and last nucleotide, and with a GC content between 60% and 70%, the command would be as follows:

python antaRNA.py -Cstr "((((...(((....)))...((((....))))...))))" \ -Cseq GNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNC \ -tGC 0.6 -tGCmax 0.7 -n 50 -t 3600

Pseudoknotted structures can be considered with the -p parameter if pKiss mfe or HotKnots or IPKnot is installed,.

In addition, all parameters of the ant colony search algorithms can be directly modified through the command line, from the random seed to initiate the search -s, the number of ants exploring (-aps, default 10), the pheromone evaporation rate(-er, default 0.2), and a wealth of others.

NUPACK

The NUPACK program provides an ensemble of tools, design being the application for inverse folding. It has less options than the previous programs but with his focus for designing long sequences viable in vitro, it can extrapolate the energy parameters for a given concentration of sodium and magnesium.

The program loads the target structure and admissible sequences, in IUPAC format, from a file PREFIX.fold. The PREFIX can be any name chosen by the user but the extension .fold must be given. Additional parameters are:

-material which can be set as rna1995 to use Turner95 energy or rna1999 for Math-ews99 energy parameters.

-sodium The sodium concentration.

-magnesium The magnesium concentration.

-prevent The name of a file, which can contain one subsequence per line forbidden in the design.

-loadseed PREFIX.init A file containing one number, the random seed to be used. Each execution of the software will choose a different random seed, but the program is deterministic and will always return the same output for a given seed. Note that the name of the file must be the same as the one with the target sequence, followed by the extension .init.

To design any sequence with a mandatory GC base pair between the first and last nucleotide, for the structure ((((...(((....)))...((((....))))...)))) using Turner95 energy parameters should be done as follows.

First create a file PREFIX.fold, where PREFIX can be admissible file name, containing the two lines:

((((...(((....)))...((((....))))...)))) GNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNC Then, call the function design.

./design -material rna1995 PREFIX Note that the suffix .fold is not given. To generate a different sequence launch the program again.

1 .

 1 ((((...(((....)))...((((....))))...)))) 2. ((((((((...(.(((((.......))))).)........((((((.......))))))..)))))))) Programs W e b s e r v e r S o u r c e c o d e E x t e n s i o n t o p s e u d o k n o t s M u l t i -t a r g e t c a p a b i l i t y RemarksGeneral antaRNA[START_REF] Kleinkauf | antaRNA: Ant colony-based RNA sequence design[END_REF][START_REF] Kleinkauf | antaRNA-Multi-objective inverse folding of pseudoknot RNA using ant-colony optimization[END_REF]

-c

  Sequence constraint to output sequences To generate at least 50 seed sequences, for the structure ((((...(((....)))...((((....))))...)))) allowing any sequence with a mandatory GC base between the first and last nucleotide and a GC content of 70%, should be done as follows. First create a file INPUT FILE, where INPUT FILE can be admissible file name later given to the -d option, containing the line: ((((...(((....)))...((((....))))...)))) Then, call the incaRNAtion script: python IncaRNAtion -d INPUT_FILE -a 1 \ -m 20 -s_gc 0.7 59 \ -c GNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNC

Table 2 :

 2 Runtime comparison between six selected programs with availability of source code. The times are reported in minutes and for 1000 runs in each of the two input structure examples. The left column is for the toy-problem designated as[START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF] in the Details of Use Section, while the right column is for the guanine-binding riboswitch aptamer designated as (2) in the Details of Use Section.

		6.5	7.8
	RNAiFold	41.4	6.4
	INFO-RNA	0.5	0.8
	NUPACK	37.5	217
	RNAinverse	3.15	3.7
	RNAexinv	231	N/A

Acknowledgements

The authors would like to thank the project students Maya Bechler-Speicher and Sarah Damyahoo, co-guided by DB together with Chen Keasar from Ben-Gurion University who assisted with his expertise, for help with the program comparisons. This work was supported by the Lynn and William Frankel Center for Computer Sciences at Ben-Gurion University and ISF within the ISF-UGC joint research program framework [9/14]. VR is grateful for being supported by an Azrieli and FQRNT postdoctoral fellowships.

Key Points

• RNA design programs should be made user friendly and accessible to biologists as much as possible, both in terms of ease of use and simplification of the input and output such that it becomes understandable to the non-specialist.

• In most cases, a balanced tradeoff between efficiency and performance in terms of the quality of the designed sequences would be the best option for the design.

• From the algorithmic standpoint, the weighted sampling approach to sample the sequence space efficiently and the fragment based design approach are desired directions that can be further developed to yield more flexibility in the design procedure.

• Programs for RNA design should aim to accumulate practical experience in biological meaningful problems, be it experimental design or computational searches for novel non-coding RNAs.

Biographical Notes:

Alexander Churkin is a computational scientist, a research fellow in the physics department at Ben-Gurion University and an instructor at the Shamoon College of Engineeing. His research interests include computational biology, RNA structure predictions, and scientific computing.

Matan Drory Retwitzer is a computational scientist who is pursuing his Ph.D. studies in the computer science Supplementary material 5 Command Line Interfaces

RNAinverse

The command line interface of RNAinverse allows subtle optimizations while the main parameters are inserted upon software request. More advanced options exists for custom alphabet, energy parameters and base pairing. Those would not be discussed here as they are a very rare usecase. For the average user, the following are the ones that will be most used.

-T Rescale energy parameters for a given temperature.

-F Select the minimization algorithm. m for energy minimization or p for partition function.

-R The number of output sequence to output for the same structure. Negative number will force the software to continue until a perfect match.

--noGU Do not allow GU pairs.

--noClosingGU Do not allow GU pairs at the end of helices.

Therefore, to find a maximum of 50 solutions, using both partition and energy minimization algorithms, for the structure ((((...(((....)))...((((....))))...))))

allowing any sequence with a mandatory GC base between the first and last nucleotide, at 25 Celsius, the command would be as follows:

./RNAinverse -R50 -Fmp -T25

Once the software begins, it will request an input structure and starting sequence. Lowercase letters will be forced into the sequence while uppercase will be considered a starting sequence. If no sequence is inserted a random seed sequence will be used.

((((...(((....)))...((((....))))...)))) gNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNc