來源:https://drmirdeep.github.io/mirdeep2_tutorial.html
Extended miRDeep2 tutorial with step by step instructions. It will cover the mapper.pl for preprocessing and mapping, the miRDeep2.pl for de-novo prediction and the quantifier.pl for expression profiling.
Disclaimer
This tutorial comes with no warranty and demands common sense of the reader. I am not responsible for any damage that happens to your computer by using this tutorial. For comments or questions just create an 'issue' here
Apparently you will need a github account for that.
Preface
This is a step by step guide for a full small RNA sequencing data analysis using the miRDeep2 package and its patched files. The first part will describe the general workflow to do de-novo miRNA predictions based on a small RNAseq data seq and the second part will focus on expression analysis with the quantification module.
Installation instructions
If you haven't installed them yet you can obtain the main package frommiRDeep2and the patched files frompatchby clicking on 'Clone or download' and then on 'Download Zip'. Extract the zipped files and then open a command line window. If you have git installed you can obtain the packages also directly from the command line by typing
·git clone https://github.com/rajewsky-lab/mirdeep2.git·
and
git clone https://github.com/Drmirdeep/mirdeep2_patch.git
To install the miRDeep2 package enter the directory to which the package was extracted to. If you extracted the folder on the Desktop then typing
cd ~/Desktop/mirdeep2
will bring you to the mirdeep2 folder. Then typing
perl install.pl
will start the installer and download and install third party software. In particular. bowtie version 1, RNAfold, randfold and the perl packages PDF:API and TTF will be installed. Please follow the instructions on the screen. When mirdeep2 was installed successfully please open a new terminal window and just type
miRDeep2.pl
If you see the miRDeep2 usage instructions on the screen you can continue to install the patch. Otherwise something went wrong during the installation. In case everything worked fine you can now enter the directory of the mirdeep2_patch by typing
cd ~/Desktop/mirdeep2_patch
if you extracted the patched file to the Desktop. Typing
bash patchme.sh
will add the patched files to your mirdeep2 installation. After that we can start with some miRNA data analysis. Download the miRBase reference files for version 21 by typing
mirbase.pl 21
This will download the hairpin.fa.gz and mature.fa.gz for version 21 to directory ~/mirbase/21/
(The ~ sign will be expanded by your computer to your home directory).
If you want the gff files as well then you need to type
mirbase.pl 21 1
The second argument can be anything but 0 which will tell the script to also get the gff files from mirbase. For miRNA quantification we next extract the miRNAs for our species of interest. For that you need to know the 3-letter code of miRBase for your species. For humans this will be 'hsa' and mouse would be mmu. To extract the mature sequences from the mirbase file we downloaded before you just need to type
extract_miRNAs.pl ~/mirbase/21/mature.fa.gz hsa > mature_ref.fa
and to get the hairpin sequences by typing
extract_miRNAs.pl ~/mirbase/21/hairpin.fa.gz hsa > hairpin_ref.fa
For running miRDeep2 for de-novo miRNA prediction it is beneficial to supply also mature miRNAs from related species. You could for example use mouse 'mmu' and chimp 'ptr' as related species. For extracting those miRNAs you can type
extract_miRNAs.pl ~/mirbase/21/mature.fa.gz mmu,ptr > mature_other.fa
Tutorial files
In case you want to follow the tutorial with example data you can get the files from heretutorial_filesand type
unzip drmirdeep.github.io-master.zip
change to the unzipped directory with
cd drmirdeep.github.io-master
and continue with the tutorial.
The data and what else you will need
Usually you will have gotten a small RNA sequencing data file from a collaborator that wants you to analyze the data file. Before you can start with any kind of analysis you should either already know the small RNA sequencing adapter that was used for the sequencing of the sample or ask your collaborator to sent it to you. If you don't clip the adapter then the majority of the reads having an adapter are likely not to be aligned to anywhere.
Once you know the adapter sequence you should do a simple check to see how many of your sequences contain the adapter. This you can do by typing
grep -c TGGAATTC example_small_rna_file.fastq
where 'TGGAATTC' are the first 8 nucleotides of the adapter that has been used for this sample. Replace it with your own sequencing adapter. MicroRNAs have a mean length of 22 nucleotides in animals so if you have sequenced one of those it will likely have the sequencing adapter attached to it. If the resulting number of sequences with an adapter is around 70% of the number of your input sequences the data set can be considered as reasonably good. Note: In case that only adapters have been sequenced predominantly you will also get a high number which is obviously not good. If you only get 10% of sequences with an adpater then likely something went wrong during the sequnecing library preparation or your sample doesn't contain too many small RNAs. For novel miRNA prediction we need to map the reads against a reference database which has to be indexed by bowtie 1. For this we take a reference database file, lets call it refdb.fa (This can be a genome file or simply a file with scaffolds) and build a bowtie index by typing
bowtie-build refdb.fa refdb.fa
The first argument is here the actual file to index, the second argument is the prefix for the bowtie index files. You can name it differently but for ease of use I use the same name as my reference database file. Depending on the input file size it can take several hours (for the human genome for example) to be indexed. However, the bowtie website has already some prebuild index files for download. If you decide to download index files you will also need to download the fasta file with which the index was build. Otherwise the results in the miRDeep2 prediction will be not reliable.
Data preprocessing for novel miRNA prediction
Since the miRDeep2 package was designed as a complete solution for miRNA prediction and quantification it also contains data preprocessing routines that will also clip the sequencing adapter. The main function of the mapping module is the mapping of the preprocessed reads file to reference database. The reference database is typically an annotated genome sequence but can also be simply a scaffold assembly if no genome is available. The scaffolds itself however should be at least 200 nucleotides long so that a sane miRNA precursor plus some flanking region fits into it. Apart from clipping adapters the module does sanity checks on your sequencing reads and also collapse read sequences to reduce the file size which will save computing time.
mapper.pl example_small_rna_file.fastq -e -h -i -j -k TGGAATTC -l 18 -m -p refdb.fa -s reads_collapsed.fa -t reads_vs_refdb.arf -v -o 4
What does this command do? The first argument needs to be your sequencing file. Typically, this will be a fastq file. The format of the fastq file is designated by specifying option '-e'. If your file is in fasta format already you specify option '-c' instead. If your reads file is not in fasta format you need to specify option '-h' which advises the mapper module to parse your file to fasta format. Option -i will convert RNA to DNA and option '-j' will remove sequences that contain characters other than ACGTN.
Now comes the actual adapter clipping which is only done if a adapter sequence is given by option -k. Only the first 6 nucleotides of this sequence will be used to search for an exact match in the sequencing reads. Option '-m' will collapse the reads to remove redundancy and decrease the file size. A sequnecing read seen 10 times in your raw file will occur only once in the collapsed file and have a _x10 in its identifier.
After that the reads will be mapped to the given refence genome which index file was specified by option '-p'. Option -s indicates the preprocessed read file name which is output by the mapper module and option -t is the file name of the read mappings to the reference database ('refdb.fa') in miRDeep2's arf format. A mapping file in arf format can be easily obtained from a standard bowtie 1 output file (This is NOT in 'sam' format but a proprietary bowtie text file format) by typing
convert_bowtie_output.pl reads_vs_refdb.bwt > reads_vs_refdb.arf
However, if you used the mapping module then the mapped output file is already in arf format.
Identification of known and novel miRNAs
For predicting novel miRNAs the miRDeep2 module from the package is called with a collapsed reads file and a reference genome file in fasta format. For better prediction results reference files of miRNAs and related miRNAs should be given since miRDeep2 considers predicted miRNAs with conserved seeds in other species more reliable that miRNAs with non-conserved seeds.
miRDeep2.pl reads_collapsed.fa refdb.fa reads_vs_refdb.arf mature_ref.fa mature_other.fa hairpin_ref.fa -t hsa 2>report.log