miRNA预测工具miRDeep-P2

之前讲过预测植物miRNA的一款软件miR-PREFER, 今天在介绍一款软件miRDeep-p2, 也叫miRDP2

 

  • 安装

　　在此之前，应安装一下软件

　　Bowite, Bowtie2, Vienna (RNA二级结构预测软件大礼包)

安装以上软件以后，在mirdp2下载最新版的miRDP2，以及ncRNA_rfam.tar.g

 1 tar -xf miRDP2-v1.1.4.tar

 2mv 1.1.4 miRDP2-v1.1.4 

在TestData下载测试数据集--TestData.tar.gz

 

  • miRNA数据处理

（1）去接头，长度选择在18-30 bp，选用cutadapt

  (2) 去低质量reads， 可以用fastp

（3）将fastq 文件转成fasta文件，并去除冗余序列，每个reads的编号：read0_x29909，x后面表示相同的序列数，最后要保证FASTA中的每个序列都唯一。

可以选用以下脚本（将.fq 放在一个文件夹）：

 1 #!/usr/bin/env python
 2 
 3 import os,re
 4 from collections import defaultdict
 5 
 6 li = os.listdir(os.getcwd())
 7 oli = filter(lambda x: x.endswith(".fa"),li)
 8 oli.sort()
 9 
10
11 for fil in oli:
12     info = defaultdict(int)
13     with open(fil) as f,\
14         open("%s.fa" %fil,"w") as o:
15         while 1:
16             name = f.readline()
17             seq = f.readline()
18             plus = f.readline()
19             qual = f.readline()
20             if name == '':
21                 break
22             info[seq.strip()] +=1
23         count = -1
24         for k,v in info.items():
25             count +=1
26             o.write(">read%s_x%s\n%s\n" %(count,v,k))

 

  • 运行 

再次之前，修改一下miRDP2-v1.1.4_pipeline.bash中的一个参数，因为我的RNAfold跑不通，所以修改

 RNAfold --noPS  中的 --noPS参数。为-noPS  

 

新建文件夹，用于存放测试数据

 1 mkdir miRDP2_Test 

将下载的测试数据以及Rfam文件上传到改文件夹，并压缩

 1 tar xf ncRNA_rfam.tar.gz 2 tar xf TestData.tar.gz 

建立索引

1 bowtie-build -f ./TestData/TAIR10_genome.fa ./TestData/TAIR10_genome.fa
3 ##为Rfam建立索引，一定得在流程中script/index 下目录下 
5 bowtie-build -f ./ncRNA_rfam.fa miRDP2-v1.1.4/scripts/ram_index

其中ncRNA_rfam.fa 为Rfam中非编码RNA (包括rRNA, tRNA,snRNA, and snoRNA), 也可以从Rfam上自行下载所有RNA.fa序列，并根据RNA类型进行分类合并。

运行流程

1 miRDP2-v1.1.4_pipeline.bash -g ./TestData/TAIR10_genome.fa -x ./TestData/TAIR10_genome -f -i ./TestData/GSM2094927.fa -o ./ 

2
3 #-g 基因组序列 
4 #-x 索引 
5 #-f sRNA-seq 为fasta格式 
6 #-i 输入RNA文件，多个文件用逗号隔开 
7 #可选 
8 #-L：reads匹配到最少的位置，默认15， 以防有重复序列 
9 #-M：bowtie 的错配，默认为0

结果：

  • miRNA预测结果: GSM2094927-15-0-10_filter_P_prediction, 每列的内容分别为，“染色体编号”，“所在链”，“代表性的短读编号”，“前体编号”，“成熟miRNA位置”，“前体位置”，“成熟序列”，“前体序列 ”
  • 日志文件: script_log和 script_err, 在运行出错时用于排除

 

 

  • 软件大概步骤

（1）将reads 比对到ncRNA seq，和known miR mature seq得到 rfam_reads.aln, known_miR.aln

利用脚本 preprocess_reads.pl 对上述 rfam_reads.aln, known_miR.aln 过滤reads，得到 *.fa 以及 *-precessed.fa，*.total_reads

（2）mapping filtered reads

将 *-precessed.fa 比对参考基因组， 得到 *_processed.aln

用 convert_bowtie_to_blast.pl 将 *_processed.aln --》 *-processed.bst （

用 filter_alignments.pl 过滤掉比对到一定次数以上（默认15）的reads， *-processed.bst ---》 *-processed_filter${len}.bst

（3）根据比对位置，提取上下游一定长度序列作为前提序列，并预测二级结构

利用 excise_candidate.pl ，将 *-processed_filter${len}.bst --》 *_precursors.fa

利用 RNAfold 软件 预测2级结构， *_precursors.fa --》 _structures

（4）提取不是ncRNA的reads 作为signature preparation

将 *.fa 比对到参考基因组， 得到 *.aln

利用convert_bowtie_to_blast.pl 将 *.aln --》*.bst

用 filter_alignments.pl 过滤掉比对到一定次数以上（默认15）的reads， 将 *.bst ---》 *_filter${len}.bst

用 filter_alignments.pl 将  *_filter${len}.bst --》 *_filtered.fa

准备 reads signature file

对 *_precursors.fa 进行bowtie-build 建库

将  *_filtered.fa 比对到 *_precursors.fa， 得到  *_precursors.aln

利用convert_bowtie_to_blast.pl 将 **_precursors.aln --》*_precursors.bst 

 将 *_precursors.bst  --〉*_signatures

 （5）miRDP core algorithm

利用  mod-miRDP.pl 将 *_signatures， *_structures --》_predictions

 

------END------

关注下方公众号可获得更多精彩

posted @ 2020-02-17 16:28  斩毛毛  阅读(3509)  评论(2)  编辑  收藏  举报