


1: INTRODUCTION
2: INSTALLATION
3: PREPROCESSING THE READS
4: ALIGNING THE READS TO THE GENOME
5: FILTERING UBIQUITOUS ALIGNMENTS
6: FILTERING BY ANNOTATION
7: PREDICTING MICRORNAS
8: MAKING PERMUTED CONTROLS
9: THE MIRDEEP SOFTWARE PACKAGE
10: THE BLASTPARSED FORMAT
11: MIRDEEP OUTPUT FORMAT



1: INTRODUCTION


This document is an introduction to the miRDeep package.

Section 2 describes the installation of dependencies.

Section 3-8 is a tutorial that reports and explains the exact command lines that were used to 
produce the set of nematode microRNA predictions presented in the miRDeep article. The files used
(except for very large files) are included in the 'demo_full' and 'demo_limited' directories, so that
the user can reproduce every step and gain familarity with miRDeep (we strongly recommend this).

Section 9 describes the scripts that constitute the miRDeep software package.

Section 10-11 describes the special formats used.




2: INSTALLATION

Several dependencies are required to run miRDeep.


First, the NCBI package can be downloaded from the site:

http://www.ncbi.nlm.nih.gov/Ftp/

The package, when unzipped, contains all executables and documentation.
Notice that numerous alternative mapping tools, specialized for short query sequences, now exist.
Users are of course very welcome to use these in the miRDeep pipeline, and will likely find that
it requires a minimum of format parsing.



Second, the Vienna package must be downloaded from the site:

http://www.tbi.univie.ac.at/~ivo/RNA/

This site contains detailed installation information for both UNIX and Windows users. If you are a
UNIX user, you can download the source from this site and type:

	tar -zxvf ViennaRNA-1.6.tar.gz
	cd ViennaRNA-1.6
	./configure
	make
and (as root)
	make install



Third, if you plan to use miRDeep's Randfold option (completely optional), you should download
Sean Eddy's SQUID library from the site:

http://selab.janelia.org/software.html


The source distribution can be unpacked on a UNIX system with the following commands:

	gunzip squid.tar.gz   Uncompresses the archive.
	tar xf squid.tar      Unpacks the archive. 
	                      (makes a new directory, squid)
	cd squid              Moves into the distribution toplevel directory.
	./configure               Configures the software for your system.
	make                      Builds the binaries. 
	make check                (optional) Runs the testsuite.	
	make install              (optional) Installs the software. 
	                          (You may need to be root to install.)
	make clean                (optional) Cleans up.



Four, if you plan to use miRDeep's Randfold option, you should download Randfold from the site:

http://bioinformatics.psb.ugent.be/supplementary_data/erbon/nov2003/

(make sure to download version 2 (C version))

The source distribution can be unpacked on a UNIX system with the following commands:

	tar xzvf randfold.tar.gz (create randfold-2.0 dirextory)
	cd randfold-2.0
	(optional) edit the Makefile
	make
	cp randfold /usr/local/bin (or whatever destination)


Fifth, you should set paths to include the location of the downloaded miRDeep scripts, as well
as the directories where you have the blastall, megablast, formatdb and RNAfold executables.
If you plan to use the Randfold option, you should also set paths to include the randfold
executable.

This is done by adding the lines:

export PATH=~/location:$PATH

to the .bashrc or equivalent file ('location' designates the desired location in the
file system).



3: PREPROCESSING THE READS

Before reads are mapped to the genome, they must be preprocessed.

This has already been done for the 454 reads used in the nematode tutorial. The preprocessing step
vary widely depending on the source of the deep sequencing reads (raw output from sequencing machine,
GEO database etc.). Because of this, miRDeep does not include a catch-all tool for this step.
However, we describe exactly what format should be attained and give an reference file, 454_total.fa.

First, the deep sequencing reads should have the adapters removed from 5' and 3' ends (if present).

Second, deep sequencing reads shorter than 15 nt should be discarded, since they will otherwise flood
the mapping output.

Third, the deep sequencing reads must be parsed into fasta format. Also, redundancy should be removed
such that reads with identical sequence are represented with a single fasta entry. Therefore, each deep
sequence id must end with a '_x' and a integer, with the integer indicating the number of times the exact
sequence was detected in the deep sequencing dataset. Last, all of the fasta ids should be unique. One way
to ensure this is to include a running number in the id.

For reference, see the 454_total.fa file.



4: ALIGNING THE READS TO THE GENOME

The next step is to align the deep sequencing reads to the genome. This step is computationally
demanding, and may be omitted by the user for the purpose of this tutorial. Note that not all files
mentioned are included in the 'demo_limited' directory, since some are very big.


The nematode genome is formatted for a blast search:

formatdb -i genome_cel.fa -p F -o T


The deep sequencing reads are aligned to the genome using megablast. -W is alignment seed length, -D is
output format and -p is minimum local identity. Note that the 'blastout' file is not included in the demos
since it is in the gigabyte range:

megablast -d genome_cel.fa -i 454_total.fa -o blastout -W 12 -D 2 -p 100


The blast output is formatted to a custom tabular separated format:

blastoutparse.pl blastout > blastparsed




5: FILTERING UBIQUITOUS ALIGNMENTS

The reads are now filtered such that only perfect alignments (full length, 100% identity) are
retained. The script allows for a designated number of unaligned nucleotides in the 3'end using
the -a option. The -c 5 options specifies that alignments from reads that align to more than
five positions in the genome should be discarded:

filter_alignments.pl blastparsed -c 5 > blastparsed_excision





6: FILTERING BY ANNOTATION

Note that this step can normally be omitted. Previous predictions indicate that the accuracy of
miRDeep increases only slightly when reads are filtered by annotation. However, since both the
C. elegans genome and the dataset used has been mined for small RNAs, this step is used in the
nematode predictions.

The ids of reads that overlap with UCSC annotations of rRNA, scRNA, snRNA, snoRNA and tRNA are
found:

overlap.pl blastparsed_excision ucsc_ncRNA -b > ids_overlap_ucsc_ncRNA


In a similar way, reads overlapping with UCSC annotations of coding region is found. Reads that
align perfectly to rfam sequences rRNA, scRNA, snRNA, snoRNA and tRNA are found using blast. Reads
that are identical to 21Us are found by sequence matching. All of these ids are combined:

cat ids_overlap_* | sort -u > ids_overlap_total


Only alignments where the read id is not included in the ids_overlap_total file are retained. -g
designates that lines where the query read id is included in the 'ids_overlap_total' file are
discarded:

blastparselect.pl blastparsed_excision -g ids_overlap_total > blastparsed_excision_filtered

Reads can also (optionally) be filtered, such that only reads that have one or more remaining
alignments are kept. -b designates that the output should be fasta entries and not alignments:

filter_alignments.pl blastparsed_excision_filtered -b 454_total.fa > 454_aligned.fa





7: PREDICTING MICRORNAS


Using the remaining alignments as guidelines, the potential precursor sequences are excised from the
genome:

excise_candidate.pl genome_cel.fa blastparsed_excision_filtered > precursors.fa


The secondary structures of the potential precursors are predicted using RNAfold. -noPS means that no
graphical output is produced:

cat precursors.fa | RNAfold -noPS > structures


The signatures are generated by aligning the remaining reads to the potential precursors. -b means that
only perfect alignments are retained, and output is sorted according to the potential precursors to
generate signature format (note that it is necessary for the prediction to do this sorting):

auto_blast.pl 454_aligned.fa precursors.fa -b > signatures


Predictions are made. Note that the miRDeep.pl does both the filtering of unplausible Dicer products
and the scoring of plausible Dicer products. Option -s designates the file that contains the microRNA
that should be considered for nucleus conservation purposes. mature_metazoan_no_cel_v10.fa contains
all miRBase version 10.0 metazoan mature microRNA sequences, except for the C. elegans microRNAs.
Option -y designates that potential precursors are tested for their relative stability using Randfold:

miRDeep.pl signatures structures -s mature_metazoan_no_cel_v10.fa -y > predictions 


Note that 117 sequences are predicted to be microRNA precursors. Of these, 103 are known ones, while
14 are new candidates.



The default option of miRDeep should be sufficient for prediction of miRNAs in genomes not previously
mined for small RNAs. However, using the advanced options will often improve sensitivity. Further,
some options give transparancy to the prediction procedure:


-s (for instance -s mature.fa) designates a fasta file with mature miRNA sequences to be considered for
conservation purposes. Since miRNA seeds are often conserved over vast phylogenetic distances, this
option will often improve sensitivity, even when predicting in model systems with few phylogenetic neighbours.

-t prints the sequences that are discarded by miRDeep rather than the ones that are predicted. The option
also prints out basic information on why the sequences was discarded. If miRDeep does not predict known
miRNAs in a sample, this option can reveal why the sensitivity is not satisfying. This can for instance
be because of a too strict cut-off.

-u gives limited output in the form of ids of the predicted potential precursors.

-v (for instance, -v 5) varies the cut-off. This is the single most important and difficult option. In
general, if you are missing a large fraction of known miRNAs that you know are in the dataset, you should  
lower the cut-off. If you make the permuted controls (see below) and get many false positives, you should
raise the cut-off. As a rule of thumb, the cut-off should probably not be lowered below 0 or 1 unless the
dataset is _highly_ enriched in miRNA sequences.

-x is a sensitive option for use with Sanger sequences obtained through conventional cloning. It uses a
lowered cut-off (-5) for added sensitivity, and emphazises read counts more in the absence of a star sequence.
It does this by allowing the read count score contribution of potential precursors without a star sequence
to exceed 0.

-y uses the Randfold application to resolve borderline candidates. Using this option increases the computing
time significantly, but the option can also increase accuracy signicantly.

-z scores the stems below the potential precursor sequence for potential Drosha recognition. In general,
more base pairings between the 10 nt flanking the potential precursor suggest a tighter stem and gives a
higher score contribution. For each number of base pairings b in the stem, the probability p(b) of a
miRNA precursor hairpin or a background hairpin having p(b) was estimated, and log-odd score contributions
were calculated. When this option is used, we recommend that a slightly higher cut-off is used (as default
1.5).




8: MAKING PERMUTED CONTROLS


This script permutes the ids and structures of an RNAfold output file.

permute_structure.pl structures > structures_permuted


The number of false positives can now be estimated by doing a regular mirdeep.pl run with the
'structures_permuted' file. Of course, a large number of permutations should be done for an accurate
estimation:

miRDeep.pl signatures structures_permuted -s mature_metazoan_no_cel_v10.fa -y > controls 







9: THE MIRDEEP SOFTWARE PACKAGE


The miRDeep software package consists of seven documented Perl scripts that should be run
sequentially by the user. miRDeep can be run on Linux or Windows platforms or any other
system that supports Perl.

a. blastoutparse.pl is used to parse standard NCBI blast output format into a custom tabular
separated format ('blastparsed', see below).

b. blastparselect.pl cleans the output from blastoutparse.pl.

c. filter_alignments.pl filters the alignments of deep sequencing reads to a genome. It
filters when only a limited part of a read is aligned. It can also filter reads that are
aligning multiple times (user-specified) to the genome. The basic input is a file in
blastparsed format.

d. overlap.pl can be used (user specified) to remove reads that align to the genome in
positions that overlap with selected annotation tracks provided by the user (for example
known rRNAs, tRNAs, etc). The basic input is a file in blastparsed format and an annotation
file in standard gff format.

e. excise_candidate.pl cuts out potential precursor sequences from a genome using aligned
reads as guidelines. The basic input is a file in blastparsed format and a genome fasta
file. The basic output is also fasta format.

f. miRDeep.pl is the core algorithm. Several files are given as input. The first is a file
in blastparsed format giving information on reads aligning to the potential precursors.
The second is a RNAfold output file giving information on the sequence, structure and absolute
stability of the potential precursors. The -s option designates a fasta file containing
the known mature miRNA sequences, whose nuclei should be considered for conservation purposes.
Command line options allow the cut-off to be varied and different verbose functions to be
enabled. Basic output is the total information on the predicted miRNA precursors. 

g. permute_structure.pl permutes the id and sequence/structure combinations of an RNAfold
output file. This is used to do the permutation controls.






10: THE BLASTPARSED FORMAT



The blastparsed format contains most of the information from ncbi blast output in tabulated fields.
One line corresponds to one local alignment:

The first field is the query id
The second field is the query length
The third field is the range of nucleotides in the query covered by the alignment
The fourth field is the subject id
The fifth field is the subject length
The sixth field is the range of nucleotides in the subject covered by the alignment
The seventh field is the blast e-value
The eigth field is the percentage of identical nucleotides in the alignment
The ninth field is the bit score
The tenth field (can contain whitespaces) are the strands of the query/subject in the alignment 

Example:

fire_9_x1       22      1..22   chrII   15279308        11537795..11537816      4e-05   1.00    44.1    Plus / Plus

Note that the blast_parsed format always uses nucleotide positions that are relative to the 5' of
the given strand. This means that for a sequence of length n, the first nucleotide on the minus
strand base pairs with the n't nucleotide on the plus strand. All nucleotides are 1-indexed.






11: MIRDEEP OUTPUT FORMAT


With the default option, miRDeep outputs the potential precursors that pass the cut-off.
With the -t options, miRDeep outputs the potential precursors that do not pass the cut-off.
For each potential precursor, miRDeep outputs one entry. Entries are separated by three newlines.


Each entry begins with the scoring fields. If the potential precursor was discarded by the initial
filtering before being scored, "signature problem" or "structure problem" will be reported. Otherwise,
the score contributions are reported:


A) 'score_nucleus' is the score contribution of the presence or absence of a conserved seed.
The following line will indicate which known miRNAs share the seed.

B) 'score_star' is the score contribution of the presence or absence of a star strand.

C) 'score_randfold' is the score contribution of the presence or absence of a relatively
energetically stable structure. This is only calculated for borderline potential precursors,
and only if the Randfold option is used.

D) 'score_mfe' is the score contribution of the minimum free energy of the potential
precursor.

E) 'score_freq' is the score contribution of the frequency of reads that are in correspondence
with Dicer processing

F) 'score' is the total score.



After the scoring fields, general sequence and structure information is reported:


G) all 'beg' and 'end' designates the begin and end positions of a feature in the 'pri'
sequence (the sequence excised from the genome). All positions are 1-indexed.

H) all 'seq' designates the sequence of a given feature.

I) all 'struct' designates the structure of a given feature in dot-bracket nomenclature.

J) 'freq' is the frequency of reads that are in correspondence with Dicer processing

K) 'flank_first' is the sequence to the 5' side of 'pre' in 'pri'.

L) 'flank_second' is the sequence to the 3' side of 'pre' in 'pri'.

M) 'loop' is the loop sequence, positionen between 'mature' and 'star'.

O) 'mature' is the mature sequence.

P) 'pre' is the potential precursor sequence, consisting of the the mature, the loop and the
star sequence.

R) 'pri' is the sequence that was excised from the genome, using the reads as guidelines. It
consists of the 'pre' sequence and the two flanks.

Q) 'star' is the star sequence.


At the very end of each entry is the signature, in blast_parsed format.
