之前讲过预测植物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
(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
结果:
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------
关注下方公众号可获得更多精彩
手机扫一扫
移动阅读更方便
你可能感兴趣的文章