Barrick Lab :: ToolsRNAStructureMutualInformation

---+!! RNA Structure Mutual Information

%TOC%

---++ Overview: What do these programs do?

[[http://en.wikipedia.org/wiki/Mutual_information][Mutual information]] between columns in a sequence alignment of structured RNAs can provide evidence for base pairs in a [[http://en.wikipedia.org/wiki/RNA_structure][secondary structure model]]. For example, if there is a [[http://en.wikipedia.org/wiki/Base_pair][Watson-Crick base pair]] between two columns, then one may have a perfectly balanced distribution of bases A, U, C, and G at each position, but only observe GC, CG, AU, UA pairs. Since you can predict the identity of the base in the second column from the first, there is more overall information (less [[http://en.wikipedia.org/wiki/Shannon_entropy][Shannon entropy]]) in the sequence string than would be predicted from each column considered alone. This correlation is also commonly referred to as "base covariation".

Raw mutual information scores between columns in RNA sequence alignments are typically poor indicators of base interactions, primarily because: (1) Real RNA alignments contain a limited number of sequences. There is a chance of observing elevated mutual information between columns by chance due to this limited sampling. (2) If the base present in one column rarely varies, then there are few opportunities to observe covariation with other columns. This means raw mutual information scores including the conserved column are depressed. (3) the sequences in an RNA alignment usually have shared evolutionary histories such that correlations between positions occur due to common descent rather than structural constraints.

These MI Perl scripts run an analysis pipeline that estimates the _statistical significance_ of mutual information scores between columns in an alignment, corrected for these considerations. In addition to providing evidence for base-paired stems in RNA structures, this procedure has been used to predict Non-Watson-Crick pairs based on unusual covariation patterns and support the identification of new RNA structure motif examples. It is important to note that this entire analysis relies on the quality of the input alignment, i.e. it can only find correlations between columns _as they are aligned_ and requires a fair number of divergent sequences in the input alignment.

More information on this approach can be found in the following publication where this approach was applied to ten riboswitch families. If you use predictions from these programs, please cite:

<div style="border-color: grey; border-style: solid; border-width: 1px; padding:1px;">
Barrick, J.E. and Breaker R.R. 2007. The structures, distributions, and mechanisms of metabolite-binding riboswitches. _Genome Biology_ *8*:R239. [[http://www.ncbi.nlm.nih.gov/sites/entrez?Db=pubmed&Cmd=ShowDetailView&TermToSearch=17997835][&laquo;PubMed&raquo;]]
</div>

---++ Installation

The mutual information Perl scripts have been developed and tested on <nop>MacOS X. It should be straightforward to install the required programs on any Unix-style system. If you encounter problems, <a href='mailto:admin@barricklab.org?subject=TWiki&#32;Feedback&#32;on&#32;Lab.ToolsRNAStructureMutualInformation'>let me know</a>.

---+++ Install a modified version of rate4site

The program <b>[[http://www.tau.ac.il/~itaymay/cp/rate4site.html][rate4site]]</b> is used to infer a phylogenetic tree with per-site substitution rates from the observed sequence alignment. In order to properly function on RNA alignments, <b>rate4site</b> requires a minor source code modification to deal with gap characters as a separate state. For convenience, I have included a modified version of the complete source release for download here.

| [[%ATTACHURL%/rate4site.tgz][%ICON{download}% Download <b>rate4site</b>]] | version 2.01 (Nov06) modified | 18 December 2008 |

The included Makefile is for compiling under Windows. Instead, compile using this command:

<code><div style="border-color: grey; border-style: solid; border-width: 1px; padding:1px;">
g++ -Dunix -DDOUBLEREP -o rate4site -O3 *.cpp
</div></code>

Finally, copy the new <b>rate4site</b> executable to a bin directory (e.g. =/usr/local/bin=) or add its location to your $PATH so that the mutual information scripts can invoke it.

---+++ Install the esl-weight program from Infernal

Chances are that you already have Infernal installed if you routinely work with RNA alignments. If not, download Infernal, compile, and install according to the included instructions.

| [[http://infernal.janelia.org/][%ICON{external}% Official Infernal Site]] |

These scripts only need the <b>esl-weight</b> utility program that is included in the easel subpackage. This program will be compiled by default, but may NOT be installed. For infernal-1.0rc5, you can find the binary at =infernal-1.0rc5/easel/miniapps/esl-weight= in your Infernal package. Manually install this program into your path (e.g. copy to =/usr/local/bin=) or add its location to your $PATH so that the mutual information scripts can employ it.

---+++ Install <nop>BioPerl

These scripts use modules for reading phylogenetic trees from [[http://www.bioperl.org/][BioPerl]]. Download and install BioPerl according to the official instructions.

| [[http://www.bioperl.org/][%ICON{external}% Official BioPerl Site]] |

---+++ Install mutual information scripts

Finally, the download and extract the mutual information Perl scripts themselves. You will (naturally) need [[http://www.perl.org][Perl]] to be installed on your system to use these.

| [[%ATTACHURL%/mi.tgz][%ICON{download}% Download <b>MI Scripts</b>]] | version 1 | 27 December 2008 |

You should be able to run these Perl scripts from their current location or add them to your $PATH.

---++ Usage

The analysis pipeline can be run using the single script =mutual_information_significance.pl=.

*Usage*
<code><div style="border-color: grey; border-style: solid; border-width: 1px; padding:1px;">
mutual_information_significance.pl -i stockholm.stk -o stockholm.mi.stk [-r 200 -n 300]
</div></code>

The input RNA alignment file in Stockholm format (=stockholm.stk=) *MUST* have an RF line. See the example alignments that are included and the Infernal documentation for a description of this line. An RF line will be generated by default if the Infernal program =cmalign= was used to generate the Stockholm file. You may want to alter or construct this line yourself. Only columns that contain non-gap characters in the RF line will be considered when removing redundant sequences.

The input alignment is processed before MI is calculated. First, identical sequences in the alignment are removed. Second, sequences that share the most identity with other sequences are removed until fewer than a specified number of sequences remain in the alignment. Third, sequence weights are calculated according to the GSC algorithm to de-emphasize closely related sequences. Fourth, all columns that are >50% gaps (taking into account sequence weights) are removed. These steps reduce the number of columns and sequences that must be considered in further calculations and typically do not affect the calculated MI significance scores.

The two parameters that you may want to adjust are =-n= to control the maximum number of sequences allowed after pruning to the most diverse and =-r= to specify how many different random alignments to generate to estimate the _p_-value significance of the actual MI score between each pair of columns (default = 200). The more random alignments used, the better the precision of the estimated _p_-value. For production runs, a value between 1000 to 10000 should be used. Note that a randomization procedure is used to estimate _p_-values, so they may differ slightly if you run the same procedure twice on the same input file.

The output file (=stockholm.mi.stk=) contains the pruned alignment with additional per-column annotation showing column pairs sorted by the significance of the mutual information between them.

Be forewarned that depending on (1) the number of sequences, (2) the number of columns in your alignment, and (3) the number of resamplings requested for estimating <i>p</i>-values that this procedure can be *extremely slow* and the intermediate resampling of alignments can require *large amounts of free disk space*. If you must interrupt operation of the script, it can usually be called from within the same working directory later and execution will pick up where it left off rather than restarting, if possible. If this procedure is too intensive for your input alignment, it is possible to parallelize the calculation of mutual information from each resampled alignment (Each alignment in the =resampled-tree= directory must be used to generate a MI file in the =resampled-mi= directory).

---++ Example

Several example Stockholm alignments are provided in the Examples folder. To "quickly" test this procedure and generate example output, use the smallest one: =SAM-II.stk=. From within this path, run:

<code><div style="border-color: grey; border-style: solid; border-width: 1px; padding:1px;">
mutual_information_significance.pl -i SAM-II.stk -o SAM-II.mi.stk
</div></code>

This command should take less than 30 minutes to complete. Several intermediate files and folders will be created. Open the resulting output file =SAM-I.mi.stk= in a text editor. Other files can be deleted after the command completes.
Attachments
Topic attachments
I	Attachment	History	Action	Size	Date	Who	Comment
tgz	mi.tgz	r1	manage	111.9 K	2008-12-27 - 16:24	JeffreyBarrick
tgz	rate4site.tgz	r1	manage	218.9 K	2008-12-18 - 20:39	JeffreyBarrick
This topic: Lab > ToolList > ToolsRNAStructureMutualInformation
Topic revision: r4 - 2008-12-27 - JeffreyBarrick