使用mirDeep2进行miRNA-seq数据分析

软件安装

首先从GitHub上下载最新的miRDeep2

git clone https://github.com/rajewsky-lab/mirdeep2.git mirdeep2.0.1.2
cd mirdeep2.0.1.2/

使用install.pl脚本进行安装

perl install.pl

会有如下的提示信息

提示信息

可以按照他的要求，直接使用source ~/.bashrc加载环境变量，然后再次运行perl install.pl就会帮你解决依赖关系，依赖工具如下

1. bowtie short read aligner
2. Vienna package with RNAfold
3. SQUID library
4. randfold
5. Perl package PDF::API2

数据分析流程

miRDeep2处理数据主要用到了三个脚本: miRDeep2.pl, mapper.pl和quantifier.pl， 需要提供如下的数据集:

1. 参考基因组的FASTA文件
2. miRBase中该物种的成熟miRNA
3. mRBase中该物种的前体miRNA
4. 高通量测序结果的FASTA文件

假如你已经有了如下文件

文件名	描述信息
cel_cluster.fa	参考基因组的FASTA文件
mature_ref_this_species.fa	miRBase中该物种的成熟miRNA
mature_ref_other_species.fa	miRBase中该物种邻近物种的成熟miRNA
precursors_ref_this_species.fa	mRBase中该物种的前体miRNA
reads.fa	高通量测序结果的FASTA文件

第一步: 建立索引

bowtie-build cel_cluster.fa cel_cluster

第二步: 将read回帖到参考基因组

mapper.pl reads.fa -c -j -k TCGTATGCCGTCTTCTGCTTGT  -l 18 -m -p cel_cluster \
  -s reads_collapsed.fa -t reads_collapsed_vs_genome.arf -v

各个参数的含义如下:

• -c: 表示输入文件是fasta，
• -e fastq: 表示输入文件是fastq
• -h 如果不是fasta，用该参数处理成fasta
• -j 移除ATCGUNatcgun以外的字符
• -k: 表示去除接头序列
• -l 18 剔除长度在18 bp以下的序列
• -m 合并相同的reads
• -p bowite索引
• -s 处理后的read
• -t 处理后比对文件
• -d 如果要处理多个样本，则指定配置文件

第三步（可选）: 快速进行定量。如果不需要预测新的miRNA， 可以用直接用miRBase数据库进行定量

quantifier.pl -p precursors_ref_this_species.fa -m mature_ref_this_species.fa \
  -r reads_collapsed.fa -t cel -y 16_19

输出结果为miRNA_expressed.csv, 记录每个样本的每个miRNA的count数，结果同样可以用网页打开expression_16_19.html查看

第四步: 鉴定新的miRNA，并进行定量

miRDeep2.pl reads_collapsed.fa cel_cluster.fa reads_collapsed_vs_genome.arf \
  mature_ref_this_species.fa mature_ref_other_species.fa \
  precursors_ref_this_species.fa -t C.elegans 2> report.log

这一步要求的参考基因组的序列不能有'ATCGN'以外的字符，没遇到报错就万事大吉，遇到报错就用tr解决吧

第五步: 浏览结果

最后可以打开results.html查看结果。

参考资料

https://github.com/rajewsky-lab/mirdeep2/blob/master/TUTORIAL.md