RNAFOLD(1)
π RNAFOLD(1) π
okay bestie, this is the User Commands page for RNAfold π β¨
β¨ NAME β¨
RNAfold - manual page for RNAfold 2.7.0 π
π§ SYNOPSIS π§
RNAfold [OPTIONS] [<input0.fa>] [<input1.fa>]...
π DESCRIPTION π
RNAfold 2.7.0 πΈ
Calculate minimum free energy secondary structures and partition function of RNAs π§¬
okay so RNAfold reads RNA sequences, works out their minimum free energy (mfe) structure, and prints that structure in bracket notation along with its free energy π. Unless you tell it otherwise with commandline arguments, it takes input from stdin or from an input file, and spits output to stdout π¦. If you give it the -p option, it also computes the partition function (pf) and base pairing probability matrix, and prints the free energy of the thermodynamic ensemble, the frequency of the mfe structure in the ensemble, and the ensemble diversity to stdout. Main character energy, literally π
It also makes PostScript files with plots of the resulting secondary structure graph and a "dot plot" of the base pairing matrix π. The dot plot is a matrix of squares. In the upper right half, the area of each square is proportional to the pairing probability, and the lower left half has one square for each pair in the minimum free energy structure. For each pair i-j with probability p>10E-6 there is a line of this form π«Ά:
i j sqrt(p) ubox
and it lives in the PostScript file, so you can easily pull the pair probabilities out π.
Sequences can be given in a simple text format where each sequence sits on a single line π. Output files are named "rna.ps" and "dot.ps" π.
Bestie, be careful: existing files of the same name will be overwritten! π
It's also possible to feed in sequence data in FASTA format π§. In this case, the first word of the FASTA header becomes the prefix for your output file names. The PostScript files "prefix_ss.ps" and "prefix_dp.ps" get made for the structure and dot plot, respectively.
Heads up: once you've given FASTA input, all following sequences must be FASTA too π.
To avoid problems with certain operating systems and/or file systems, the prefix will ALWAYS be sanitized! π This swaps any special character in the prefix for a filename delimiter. See the --filename-delim option if you want to change the delimiting character π.
The program keeps reading new sequences until it hits a line with just the single character @, or an end of file (EOF) condition π₯.
-h,--helpβ print help and exit π--detailed-helpβ print help, including all details and hidden options, and exit π--full-helpβ print help, including hidden options, and exit π«§-V,--versionβ print version and exit π·οΈ-v,--verboseβ be verbose. (default=off) Lowers the log level setting so even INFO messages get passed through π£οΈ.
π₯ I/O Options π€
Command line options for input and output (pre-)processing πΈ
-i,--infile=filenameβ read a file instead of reading from stdin. By default RNAfold reads input from stdin or from the file(s) that follow the RNAfold command. With this one you name the input file(s) you want data read from. Heads up: any additional files you supply to RNAfold still get processed too π .-o,--outfile[=filename]β print output to a file instead of stdout. This lets you write all output to output files instead of printing to stdout π. The default filename isRNAfold_output.foldif no FASTA header comes before the input sequences and the--auto-idfeature is inactive. Otherwise, output files use the schemeprefix.fold, where the prefix comes from the sequence id, like the FASTA header π§Ύ. You can give one output file name for all the data generated from the input by typing a filename right after this option. If a file with that same name already exists, any output gets appended to it π·.
Note: any special characters in the filename get replaced by the filename delimiter, so there is no way (yet) to pass a whole directory path through this option π
ββοΈ. (See also the --filename-delim parameter.)
-j,--jobs[=number]β split batch input into jobs and start processing in parallel using multiple threads π§΅. A value of0means use as many parallel threads as there are computation cores available. (default=`0') By default, input is processed serially, meaning one sequence at a time. With this switch, you can start computing many sequences from the input in parallel instead π. RNAfold makes as many parallel computation slots as you specify and hands out the input sequences from your input file(s) to the open slots.
Note: this increases memory consumption, since input alignments have to stay in memory until an empty compute slot opens up, and each running job needs its own dynamic programming matrices π.
--unorderedβ don't try to keep output in order with input while parallel processing is happening π. (default=off) When parallel input processing (the--jobsflag) is on, the order that input gets processed in depends on the host machine's job scheduler. So any output to stdout or files will most likely not follow the order of the matching input π€·ββοΈ. By default RNAfold uses a specialized data structure to keep results in input order, but that costs memory, since all output has to be held in memory until chunks of consecutive, ordered output are available. Set this flag and RNAfold won't buffer individual results. It prints them as soon as they've been computed β‘.--noconvβ don't automatically swap nucleotide "T" for "U" π. (default=off)--auto-idβ automatically generate an ID for each sequence. (default=off) The default mode is to determine an ID from the input sequence data whenever the input file format allows it, and sequence IDs usually come in the FASTA header πͺͺ. If this flag is active, RNAfold ignores any IDs pulled from the input and generates one for each sequence instead. Each generated ID is a prefix plus an increasing number π’. You can also use this flag to add a FASTA header to the output even if the input doesn't have one β¨.--id-prefix=STRINGβ prefix for automatically generated IDs (as used in output file names). (default=sequence') If you set this, every sequence gets prefixed with your string, so the output files follow this naming scheme:prefix_xxxx_ss.ps(secondary structure plot),prefix_xxxx_dp.ps(dot-plot),prefix_xxxx_dp2.ps` (stack probabilities), etc., where xxxx is the sequence number π¨.
Note: setting this parameter implies --auto-id π.
--id-delim=CHARβ change the delimiter between prefix and increasing number for automatically generated IDs (as used in output file names). (default=`_') This swaps out the default "_" between the prefix string and the increasing number in your auto-generated IDs π©·.--id-digits=INTβ specify the number of digits of the counter in automatically generated alignment IDs. (default=4') Auto-generated alignment IDs get an increasing number starting with 1. That number is always left-padded with leading zeros, so it takes up a fixed width π. Use this parameter to set the width to whatever you need. We allow numbers in the range[1:18]. This option implies--auto-id` π.--id-start=LONGβ specify the first number in automatically generated IDs. (default=1') Auto-generated sequence IDs get an increasing number, usually starting at 1. This lets you pick the first number you want π. Note: negative numbers are not allowed π«. Note: setting this parameter means ignoring any IDs retrieved from the input data, i.e. it activates the--auto-id` flag π«Ά.--filename-delim=CHARβ change the delimiting character used in sanitized filenames. (default=ID-delimiter') This changes the character used while sanitizing filenames, aka swapping out invalid characters π§Ή. Note: the default delimiter ALWAYS is the first character of the "ID delimiter" you supplied through the--id-delimoption. If the delimiter is a whitespace character or empty, invalid characters get simply removed rather than substituted π«§. Right now, these count as illegal for filenames: backslash'', slash'/', question mark'?', percent sign'%', asterisk'*', colon':', pipe symbol'|', double quote'"', and triangular brackets'<'and'>'` π¨.--filename-fullβ use full FASTA header to create filenames. (default=off) This turns off the default behavior of limiting output filenames to the first word of the sequence ID π§ . Consider this example: an input with the FASTA header>NM_0001 Homo Sapiens some geneusually produces output files with the prefix "NM_0001", leaving out the extra data in the header. So you get something like "NM_0001_ss.ps" for secondary structure plots. With this flag set, nothing gets truncated, and output filenames receive the full FASTA header data as their prefix π. Note, however, that invalid characters (such as whitespace) will be swapped for a delimiting character or simply removed (see also the--filename-delimoption) π.--log-level=levelβ set log level threshold. (default=`2')
By default, any log messages are filtered so only warnings (level 2) or errors (level 3) get printed π This setting lets you pick the log level threshold yourself. Higher values mean fewer messages, so less info but way more quiet β¨ Log-level 5 turns off everything, even errors and other critical stuff π
--log-file[=filename]β Sends log messages to a file instead of stderr π (default=RNAfold.log)--log-timeβ Adds a time stamp to log messages β° (default=off)--log-callβ Includes the file and line of the function making the log call π (default=off)
π§ͺ Algorithms π§ͺ
Pick any extra algorithms you want included in the calculations π The MFE and a structure representative are always calculated no matter what, so you're covered either way π«Ά
-p,--partfunc[=INT]β Calculates the partition function and base pairing probability matrix π (default=1)
Along with the MFE structure, it prints a coarse pseudo bracket notation of the pair probabilities, followed by the ensemble free energy. The notation uses the letters ., ,, |, {, }, (, and ) for bases that are essentially unpaired, weakly paired, strongly paired without preference, weakly upstream (downstream) paired, or strongly up- (down-)stream paired, respectively π§©
On the next line you get the centroid structure derived from the pair probabilities, along with its free energy and distance to the ensemble β¨ Last but not least, it prints the frequency of the mfe structure and the structural diversity (the mean distance between the structures in the ensemble) π For the deep dive, check the description of vrna_pf(), mean_bp_dist(), and vrna_centroid() in the RNAlib documentation π
Heads up! Unless you also specify-d2or-d0, the partition function and mfe calculations will use a slightly different energy model. See the discussion of dangling end options below β οΈ
An extra value passed to this option changes how the partition function gets calculated, so pay attention to these:
-p0β Calculates the partition function but not the pair probabilities, saving about 50% in runtime πββοΈ It prints the ensemble free energydG=-kT ln(Z).-p2β Computes stack probabilities, i.e. the probability that a pair(i,j)and the immediately enclosed pair(i+1,j-1)form at the same time, on top of the regular pair probabilities π A second postscript dot plot namedname_dp2.ps(ordot2.psif the sequence has no name) gets made. It shows pair probabilities in the upper right half and stack probabilities in the lower left π¦
--betaScale=DOUBLEβ Sets the scaling of the Boltzmann factors π‘οΈ (default=1.)
This value scales the thermodynamic temperature in the Boltzmann factors on its own, separate from the temperature of the individual loop energy contributions. The Boltzmann factors then become exp(- dG/(kT*betaScale)), where k is the Boltzmann constant, dG is the free energy contribution of the state, and T is the absolute temperature π§ͺ
-S,--pfScale=DOUBLEβ In the pf calculation, usesscale*mfeas an estimate for the ensemble free energy, which helps avoid overflows π (default=1.07)
The default is 1.07, and useful values are 1.0 to 1.2 π You'll occasionally need this for long sequences π§΅
--MEA[=gamma]β Computes the MEA (maximum expected accuracy) structure π― (default=1.)
The expected accuracy comes from the pair probabilities. Each base pair (i,j) gets a score of 2*gamma*p_ij, and the score of an unpaired base is the probability of not forming a pair π The parameter gamma is the vibe check between correctly predicted pairs and unpaired bases. So with small gamma values, the MEA structure only keeps pairs with super high probability π€ Using --MEA automatically implies -p so the pair probabilities get computed π«Ά
-c,--circβ Treats the RNA molecule as circular instead of linear π (default=off)--ImFeelingLuckyβ Returns exactly one stochastically backtracked structure π² (default=off)
This one computes the partition function and returns exactly one secondary structure, randomly sampled from the Boltzmann equilibrium according to how probable it is in the ensemble π
--bppmThreshold=cutoffβ Sets the threshold/cutoff for which base pair probabilities show up in the postscript output ποΈ (default=1e-5)
By tweaking the threshold, you control which base pair probabilities end up in the output π
By default, only the ones exceeding 1e-5 in probability show up as squares in the dot plot. Change the threshold to any other value and you can increase or decrease the amount of data β¨
-g,--gquadβ Incorporates G-Quadruplex formation into the structure prediction algorithm π (default=off)
π§΅ Structure Constraints π§΅
Command line options that let you interact with the structure constraints feature of this program π«§
--maxBPspan=INTβ Sets the maximum base pair span π (default=-1)-C,--constraint[=filename]β Calculates structures subject to constraints π (default='')
The program reads the sequence first, then a string containing constraints on the structure, encoded with these symbols:
.β no constraint for this base|β the corresponding base has to be pairedxβ the base is unpaired<β base i is paired with a base j>i>β base i is paired with a base j<i(and)β matching brackets, meaning base i pairs base j
With the exception of |, constraints disallow all pairs that conflict with the constraint π« This is usually enough to enforce the constraint, but every now and then a base may stay unpaired even with constraints in place π PF folding ignores constraints of type |.
--batchβ Uses constraints for multiple sequences π (default=off)
Normally, constraints from an input file only apply to a single input sequence. So RNAfold stops computing and quits after the first input sequence is processed. With this switch, RNAfold processes multiple input sequences and applies the same constraints to each one π
--canonicalBPonlyβ Removes non-canonical base pairs from the structure constraint π§Ή (default=off)--enforceConstraintβ Enforces base pairs given by round brackets()in the structure constraint πͺ (default=off)--shape=filenameβ Uses SHAPE reactivity data to guide structure predictions π§¬--shapeMethod=methodβ Selects the SHAPE reactivity data incorporation strategy π₯ (default=D)
These methods convert SHAPE reactivities into pseudo energy contributions:
Dβ Converts using the linear equation according to Deigan et al 2009. Derived pseudo energy terms get applied to every nucleotide involved in a stacked pair. This method is recognized by a capitalDin the provided parameter, so--shapeMethod="D"is the default setting π The slopemand the interceptbcan be set to non-default values if needed, otherwise m=1.8 and b=-0.6. To alter them, e.g. m=1.9 and b=-0.7, use a parameter string like--shapeMethod="Dm1.9b-0.7". You can also provide just one of the two, like--shapeMethod="Dm1.9"or--shapeMethod="Db-0.7".Zβ Converts SHAPE reactivities to pseudo energies according to Zarringhalam et al 2012 πΈ Reactivities get turned into pairing probabilities using linear mapping. Deviations from the observed pairing probabilities get penalized during the folding recursion. You can adjust the magnitude of the penalties by changing the factor beta, e.g.--shapeMethod="Zb0.8".Wβ Applies a given vector of perturbation energies to unpaired nucleotides according to Washietl et al 2012 π¦ Perturbation vectors can be calculated usingRNApvmin.
--shapeConversion=methodβ Selects the method for SHAPE reactivity conversion πͺ (default=O)
This parameter is useful when dealing with SHAPE incorporation according to Zarringhalam et al. These methods convert SHAPE reactivities into the probability that a certain nucleotide is unpaired:
Mβ Uses linear mapping according to Zarringhalam et al.Cβ Uses a cutoff approach to split nucleotides into paired and unpaired (e.g."C0.25")Sβ Skips the normalizing step, since the input data already represents probabilities of being unpaired rather than raw reactivity values πLβ Uses a linear model to convert the reactivity into a probability of being unpaired (e.g."Ls0.68i0.2"for a slope of 0.68 and an intercept of 0.2)Oβ Uses a linear model to convert the log of the reactivity into a probability of being unpaired (e.g."Os1.6i-2.29"for a slope of 1.6 and an intercept of -2.29) π§
--motif=SEQUENCE,STRUCTURE,ENERGYβ Specifies the stabilizing energy of a ligand binding to a particular sequence/structure motif π
Some ligands that bind RNAs need and/or induce particular sequence and structure motifs. For example, they might bind to an internal loop or a small hairpin loop. If a binding free energy is known for those cases, the binding (and therefore the stabilizing effect of the ligand) can be included in the folding recursions π§ Interior loop motifs are written as concatenations of 5' and 3' motifs, separated by an & character. Energy contributions must be given in kcal/mol. See the manpage for an example usage of this option π
--commands=filenameβ Reads additional commands from a file ποΈ
Commands include hard and soft constraints, plus structure motifs in hairpin and internal loops that need to be treated differently. You can also set commands for unstructured and structured domains π
β¨ Energy Parameters β¨
okay so the energy parameter sets can be tweaked or loaded from your own input files π
-T,--temp=DOUBLEβ rescale the energy parameters to a temperature of temp C. it defaults to 37C π‘οΈ ((default=`37.0'))-P,--paramFile=paramfileβ read energy parameters from paramfile instead of the default set π
bestie, different sets of energy parameters for RNA and DNA should come with your distribution. check the RNAlib documentation for the file format details π the placeholder file name 'DNA' lets you load DNA parameters without having to point at any input file π§¬
-4,--noTetraβ skip the special tabulated stabilizing energies for tri-, tetra- and hexaloop hairpins π« (default=off)
this one is mostly for testing π§ͺ
--salt=DOUBLEβ set the salt concentration in molar (M). the default is 1.021M π§-m,--modifications[=STRING]β allow modified bases inside the RNA sequence string π§¬β¨ ((default=`7I6P9D'))
so this is giving modified-base realness π€© the modified bases in your input have to be marked with their one-letter code. use the corrections and/or pairing partner rules when they're available. if you don't add any extra arguments, every available correction gets applied π or you can narrow it down to certain codes. for example, -mP6 will only correct for pseudouridine and m6A bases.
these are the currently supported one-letter codes and their energy corrections:
'7': 7-deaza-adenonsine (7DA)'I': Inosine'6': N6-methyladenosine (m6A)'P': Pseudouridine'9': Purine (a.k.a. nebularine)'D': Dihydrouridine
--mod-file=STRINGβ pull in additional modified base data from a JSON file π
β¨ Model Details β¨
okay this part is where you can tweak the energy model and the pairing rules with the parameters below π§π
-d,--dangles=INTβ decide how to treat "dangling end" energies for bases sitting next to helices in free ends and multi-loops π ((default=`2'))
so with -d1, only unpaired bases can take part in at most one dangling end. with -d2 that check gets skipped, and dangling energies get added for bases next to a helix on both sides no matter what. -d2 is the default for mfe and partition function folding (-p). -d0 ignores dangling ends altogether, which is mostly for debugging π and -d3 lets mfe folding allow coaxial stacking of adjacent helices in multi-loops. heads up, the implementation won't allow coaxial stacking of the two enclosed pairs in a loop of degree 3 right now, and it only works for mfe folding β οΈ
note: with-d1and-d3, only the MFE computations use this setting. partition function uses the-d2setting, so dangling ends get treated differently between the two π
--noLPβ produce structures without lonely pairs (helices of length 1) π ββοΈ (default=off)
for partition function folding, this only disallows pairs that can only show up isolated. other pairs might still pop up now and then as helices of length 1 π€·ββοΈ
--noGUβ do not allow GU pairs π« (default=off)--noClosingGUβ do not allow GU pairs at the end of helices π« (default=off)--nsp=STRINGβ allow other pairs on top of the usual AU, GC, and GU pairs π
the argument is a comma separated list of the extra pairs you want to allow. if the first character is a -, then AB implies that AB and BA are both allowed pairs. for example, --nsp="-GA" will allow GA and AG pairs. nonstandard pairs get a stacking energy of 0 π«
--energyModel=INTβ set the energy model β‘
this one is rarely used. it folds sequences from the artificial ABCD... alphabet, where A pairs with B, C pairs with D, and so on. use the energy parameters for GC (--energyModel 1) or AU (--energyModel 2) pairs π§ͺ
--helical-rise=FLOATβ set the helical rise of the helix in units of Angstrom π ((default=`2.8'))
use with caution! β οΈ this value gets reset to 3.4 automatically if you load DNA parameters with -P DNA and don't give any other value.
--backbone-length=FLOATβ set the average backbone length for looped regions in units of Angstrom 𦴠((default=`6.0'))
use with caution! β οΈ this value gets reset to 6.76 automatically if you load DNA parameters with -P DNA and don't give any other value.
β¨ Plotting β¨
these flags change the default behavior of structure layouts and pairing probability plots π¨π¦
--noPSβ skip the postscript drawing of the mfe structure πΌοΈ (default=off)--noDPβ don't produce the dot-plot postscript file with base pair or stack probabilities π (default=off)
when you combine this with -p, the flag turns off creation of the individual dot-plot files. that means the computed base pair probability output gets left out, but centroid and MEA structure prediction still happen β¨
-t,--layout-type=INTβ pick the layout algorithm π§© (default=`1')
this picks the algorithm that works out the nucleotide coordinates. these are the ones available right now:
'0': simple radial layout'1': Naview layout (Bruccoleri et al. 1988)'2': circular layout'3': RNAturtle (Wiegreffe et al. 2018)'4': RNApuzzler (Wiegreffe et al. 2018)
β¨ REFERENCES β¨
if you use this program in your work, you might wanna cite these ππ«Ά
- R. Lorenz, S.H. Bernhart, C. Hoener zu Siederdissen, H. Tafer, C. Flamm, P.F. Stadler and I.L. Hofacker (2011), "ViennaRNA Package 2.0", Algorithms for Molecular Biology: 6:26
- I.L. Hofacker, W. Fontana, P.F. Stadler, S. Bonhoeffer, M. Tacker, P. Schuster (1994), "Fast Folding and Comparison of RNA Secondary Structures", Monatshefte f. Chemie: 125, pp 167-188
- R. Lorenz, I.L. Hofacker, P.F. Stadler (2016), "RNA folding with hard and soft constraints", Algorithms for Molecular Biology 11:1 pp 1-13
- M. Zuker, P. Stiegler (1981), "Optimal computer folding of large RNA sequences using thermodynamic and auxiliary information", Nucl Acid Res: 9, pp 133-148
- J.S. McCaskill (1990), "The equilibrium partition function and base pair binding probabilities for RNA secondary structures", Biopolymers: 29, pp 1105-1119
- I.L. Hofacker & P.F. Stadler (2006), "Memory Efficient Folding Algorithms for Circular RNA Secondary Structures", Bioinformatics
- A.F. Bompfuenewerer, R. Backofen, S.H. Bernhart, J. Hertel, I.L. Hofacker, P.F. Stadler, S. Will (2007), "Variations on {RNA} Folding and Alignment: Lessons from Benasque", J. Math. Biol.
- D. Adams (1979), "The hitchhiker's guide to the galaxy", Pan Books, London
the mfe structure calculations are based on the dynamic programming algorithm originally developed by M. Zuker and P. Stiegler π§ and the partition function algorithm is based on work by J.S. McCaskill π‘
the energy parameters are taken from:
- D.H. Mathews, M.D. Disney, D. Matthew, J.L. Childs, S.J. Schroeder, J. Susan, M. Zuker, D.H. Turner (2004), "Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure", Proc. Natl. Acad. Sci. USA: 101, pp 7287-7292
- D.H Turner, D.H. Mathews (2009), "NNDB: The nearest neighbor parameter database for predicting stability of nucleic acid secondary structure", Nucleic Acids Research: 38, pp 280-282
β¨ EXAMPLES β¨
okay the examples are where it all clicks π₯π
single line sequence input, partition function and MEA structure
bestie, start by running this command, which gives you the partition function and the MEA structure for one sequence at a time π¦
$ RNAfold --MEA -d2 -p
the program will then ask you for a sequence. type in CGACGTAGATGCTAGCTGACTCGATGC and hit ENTER, and the output will look something like this π
CGACGUAGAUGCUAGCUGACUCGAUGC
(((.((((.......)).)))))....
minimum free energy = -1.90 kcal/mol
(((.((((.......))},})))....
free energy of ensemble = -2.86 kcal/mol
(((.(.((.......))..)))).... { 0.80 d=2.81}
(((.((((.......))).)))).... { -1.90 MEA=22.32}
frequency of mfe structure in ensemble 0.20997; ensemble diversity 4.19
so line one just repeats your sequence back to you π line two is the mfe structure in dot bracket notation, followed by its minimum free energy π next up are the pairing probabilities for each nucleotide, shown in a pseudo dot-bracket notation, and then the free energy of the ensemble β¨ the next two lines are the centroid structure with its free energy and its distance to the ensemble, plus the MEA structure with its free energy and its maximum expected accuracy π― and the last line gives you the frequency of the mfe representative in the whole ensemble of secondary structures, along with the ensemble diversity π for the full scoop on how the calculation works and how to read the output, check the RNAlib reference manual π
FASTA file sequence input
since version 2.0, you can also feed it a FASTA file π let's say you've got a file with two sequences in FASTA format, like this:
$ cat sequences.fa
>seq1
CGGCUCGCAACAGACCUAUUAGUUUUACGUAAUAUUUG
GAACGAUCUAUAACACGACUUCACUCUU
>seq2
GAAUGACCCGAUAACCCCGUAAUAUUUGGAACGAUCUA
UAACACGACUUCACUCUU
to compute the MFE for both sequences, use this command π§
$ RNAfold < sequences.fa
and you'll get output that looks like this π
>seq1
CGGCUCGCAACAGACCUAUUAGUUUUACGUAAUAUUUGGAACGAUCUAUAACACGACUUCACUCUU
.((.(((...((((..(((((........)))))))))...))).))................... ( -5.40)
>seq2
GAAUGACCCGAUAACCCCGUAAUAUUUGGAACGAUCUAUAACACGACUUCACUCUU
.......((((..............))))........................... ( -2.00)
β¨ CONSTRAINT EXAMPLES β¨
okay so you can also feed RNAfold secondary structure constraints alongside the sequence π Let's take the first sequence from the last example and make the outermost helix unpaired, meaning no base pairs (2,47) or (3,46). Your input file should look like this:
$ cat sequence_unpaired.fa
>seq1
CGGCUCGCAACAGACCUAUUAGUUUUACGUAAUAUUUG
GAACGAUCUAUAACACGACUUCACUCUU
.xx...................................
.......xx...................
Now run it with the structure constraint option -C and watch it serve β¨
$ RNAfold -C < sequence_unpaired.fa
>seq1
CGGCUCGCAACAGACCUAUUAGUUUUACGUAAUAUUUGGAACGAUCUAUAACACGACUUCACUCUU
....(((...((((..(((((........)))))))))...)))...................... ( -4.20)
That gives you the minimum free energy and a representative structure for the RNA, given that nucleotides 2,3,46 and 47 must not be involved in any base pair π«Ά For more on constrained folding, check the details of the -C option and the RNAlib reference manual.
Since version 2.2, the ViennaRNA Package distinguishes hard and soft constraints, which is literally why predictions can handle so many options: maximal base pair span, SHAPE reactivity data, and RNA-ligand binding to hairpin or interior loop motifs. Iconic behavior π
Restricting the maximal span of a base pair
There's a convenience command-line option that lets you cap the distance (j - i + 1) between two nucleotides i and j that form a base pair. For instance, a limit of 600nt can be set like this:
$ RNAfold --maxBPspan 600
Guide structure prediction with SHAPE reactivity data
You can use SHAPE reactivity data to steer secondary structure prediction:
$ RNAfold --shape=reactivities.dat < sequence.fa
Here, reactivities.dat is a two-column text file with sequence positions (1-based) and normalized reactivity values (usually between 0 and 2). Missing values can be left out, or given a negative score:
$ cat reactivities.dat
9 -999 # No reactivity information
10 -999
11 0.042816 # normalized SHAPE reactivity
12 0 # also a valid SHAPE reactivity
15 0.15027 # Missing data for pos. 13-14
...
42 0.16201
π¨ Important: when you feed it SHAPE reactivity data, RNAfold will only process the first sequence in the input file!
Complex structure constraints and grammar extensions
Structure constraints that go beyond what pseudo-dot-bracket notation can express can be provided in a so-called command file:
$ RNAfold --commands=constraints.txt < sequence.fa
The command file syntax is a generalization of the constraints used in UNAfold/mfold. Each line starts with a one or two letter command, followed by command parameters. For structure constraints, that's a single command character followed by three or four numbers. Optional auxiliary modifier characters can also limit a constraint to specific loop types. For base pair specific constraints, the loop types are:
- E β pairs in exterior loops
- H β closing pairs of hairpin loops
- I β closing pairs of interior loops
- i β enclosed pairs of interior loops
- M β closing pairs of multibranch loops
- m β enclosed pairs of multibranch loops
Nucleotide-wise constraints can be limited to their loop context using the matching uppercase characters. The default is to apply a constraint to all loop types, which is A. Pairing constraints for single nucleotides can also be limited to upstream (U) or downstream (D) orientation. Here's the full command file spec π
F i 0 k [TYPE] [ORIENTATION] # Force nucleotides i...i+k-1 to be paired
F i j k [TYPE] # Force helix of size k starting with (i,j) to be formed
P i 0 k [TYPE] # Prohibit nucleotides i...i+k-1 to be paired
P i j k [TYPE] # Prohibit pairs (i,j),...,(i+k-1,j-k+1)
P i-j k-l [TYPE] # Prohibit pairing between two ranges
C i 0 k [TYPE] # Nucleotides i,...,i+k-1 must appear in context TYPE
C i j k # Remove pairs conflicting with (i,j),...,(i+k-1,j-k+1)
E i 0 k e # Add pseudo-energy e to nucleotides i...i+k-1
E i j k e # Add pseudo-energy e to pairs (i,j),...,(i+k-1,j-k+1)
UD m e [LOOP] # Add ligand binding to unstructured domains with motif
# m and binding free energy e
# [LOOP] = { E, H, I, M, A }
# [TYPE] = [LOOP] + { i, m }
# [ORIENTATION] = { U, D }
Again, by default RNAfold only processes the first sequence in the input when constraints come from a command file. To apply the exact same constraints to every sequence in a multi-FASTA file, use the batch mode command-line option:
$ RNAfold --constraint=constraints.txt --batch < sequences.fa
Ligand binding contributions to specific hairpin/interior loop motifs
A convenience function lets you specify a hairpin/interior loop motif where a ligand binds with a particular binding free energy dG. Here's an example that adds a theophylline binding motif π§ The free energy contribution of this motif, dG=-9.22kcal/mol, is derived from k_d=0.32umol/l, taken from Jenison et al. 1994. The structure motif consists of a symmetric interior loop of size 6, followed by a small helix of 3 base pairs, and a bulge of 3 nucleotides, but the entire structure can still be represented by one interior loop.
See the motif description below, where the & character splits the motif into a 5' and a 3' part. The first line gives the sequence motif, the second line shows the actual structure motif of the aptamer pocket, and the third line is the interior loop motif that fully encapsulates the theophylline aptamer:
GAUACCAG&CCCUUGGCAGC
(...((((&)...)))...)
(......(&).........)
To use this in the folding recursions of RNAfold, you only need to provide the motif itself and the binding free energy:
$ RNAfold --motif="GAUACCAG&CCCUUGGCAGC,(...((((&)...)))...),-9.22" < sequences.fa
Adding the --verbose option to the call above also prints the sequence position of each motif found in the MFE structure. If interior-loop-like motifs are provided, two intervals get printed, denoting the 5' and 3' part respectively. Ok, that's a vibe πΈ
Ligand binding contributions to unpaired segments of the RNA structure
The extension of the RNA folding grammar with unstructured domains makes it easy to incorporate ligands that bind to unpaired stretches of an RNA structure π¦ To model these interactions, you only need two parameters: (i) a sequence motif in IUPAC notation that specifies where the ligand binds, and (ii) a binding free energy that can be derived from the association/dissociation constant of the ligand.
With those two in hand, you can model the competition between regular intramolecular base pairing and ligand interaction just by writing a simple command file of this form:
UD m e [LOOP]
Here, m is the motif string in uppercase IUPAC notation, e is the binding free energy in kcal/mol, and [LOOP] is an optional loop type restriction. See also the command file specification above.
For instance, say a protein has a 4-nucleotide footprint binding AAAA, with a binding free energy e = -5.0 kcal/mol, and you want to restrict binding to exterior and multibranch loops. Your command file would look like this:
$ cat commands.txt
UD AAAA -5.0 ME
And the corresponding call to RNAfold to compute the MFE and equilibrium probabilities becomes:
$ RNAfold --commands=commands.txt -p < sequence.fa
The resulting MFE plot will be annotated to display the binding site(s) of the ligand, and the base pair probability dot-plot is extended to include the probability that a particular nucleotide is bound by the ligand. main character energy, fr fr π
β¨ POST-TRANSCRIPTIONAL MODIFICATION EXAMPLES β¨
Many RNA molecules harbor (post-transcriptional) modifications 𧬠These modified bases often change the pairing behavior or energy contribution of the loops they're part of. To account for that effect (to a certain degree), you can use additional correcting energy parameters for loops containing the modified bases. The literature has a few stacking and some terminal mismatch energies, and some of them are already provided within the ViennaRNA Package.
The --modification and --mod-file command-line parameters can be used to apply these parameters in predictions. The former lets you select a subset of the implemented modified base corrections, while the latter lets the prediction programs read energy parameters for modified bases from a user-provided JSON file π
Consider, for instance, the following tRNA sequence, with dihydrouridines and pseudouridines annotated by their one-letter codes D and P:
$ cat tRNAphe.fa
>tRNAphe
GCCGAAAUAGCUCAGDDGGGAGAGCGPPAGACUGAAGAPCUAAAGGDCCCUGGUPCGAUCCCGGGUUUCGGCACCA
Now, a prediction that supports the destabilizing effect of D and the stabilizing effects of P within base pair stacks can be done like this:
$ RNAfold --modifications=DP < tRNAphe.fa
>tRNAphe
GCCGAAAUAGCUCAGDDGGGAGAGCGPPAGACUGAAGAPCUAAAGGDCCCUGGUPCGAUCCCGGGUUUCGGCACCA
(((((((..((((........)))).(((((.......)))))....(((.((......)).)))))))))).... (-23.37)
β¨ AUTHOR β¨
Ivo L Hofacker, Walter Fontana, Sebastian Bonhoeffer, Peter F Stadler, Ronny Lorenz π«Ά
β¨ REPORTING BUGS β¨
If in doubt, our program is right and nature is at fault π Comments should be sent to rna@tbi.univie.ac.at.
That's the whole page, bestie! π§π This is RNAfold 2.7.0, last updated February 2026, so go be iconic with your next fold πβ¨
π¦ where this page comes from
- ubuntu resolute: vienna-rna 2.7.0+dfsg-2 Β· non-free
- ubuntu stonking: vienna-rna 2.7.0+dfsg-2 Β· non-free
- π girlified by claude-haiku-5-5 on : a language model's rewrite of the original, under its terms