用Matlab筛选mirbase,一种基于miRBase数据库的植物有参的miRNA数据分析方法与流程…

本发明涉及转录组测序领域，具体涉及一种在miRBase数据库中有参考数据的植物miRNA测序的数据分析方法。

背景技术：

miRNA是一类由内源基因编码非编码单链RNA分子，在动植物中参与转录后基因表达调控。多数miRNA以单拷贝、多拷贝或基因簇的形式存在于基因组中。miRNA在很多物种中被广泛发现，且在进化进程中高度保守，因此研究miRNA的确切功能、目的靶基因、以及其作用机制，是转录组学数据分析中的重要一环，对于了解生物体内基因的表达调控机制有重要意义。

miRNA的作用机制在动物和植物之间存在明显差异，且有的物种有丰富的miRNA参考数据，但有的物种缺乏参考数据，甚至有些物种没有参考基因组信息，这些情况下的miRNA测序的数据分析方法十分不同。目前还没有针对植物小RNA分析的工具。植物miRNA和动物miRNA在生物体内的作用机制不同，保守性程度也不同。目前还没有现成的流程分析植物小RNA测序数据；尤其是没有自动化的分析平台实现植物小RNA测序结果的流程化分析工具，包括后续的sRNA注释，miRNA序列的特征分析，表达量分析和差异分析，靶基因位点分析，等各个步骤的自动化整合。

技术实现要素：

为了克服现有技术所存在的上述缺陷，本发明的目的在于提供一种基于miRBase数据库的植物有参的miRNA数据分析方法。

为了实现本发明的目的之一，所采用的技术方案是：

一种基于miRBase数据库的植物有参的miRNA数据分析方法，包括如下步骤：

步骤一、文件准备步骤：准备并读取config文件，软件读取相关信息后，会生成进行以下列出的所有分析步骤对应的shell脚本，按顺序运行即可，在运行同时每一步都会有运行日志，方便结果检查；

步骤二、下机数据过滤步骤：

下机后的原始数据，去除接头，然后过滤低质量序列，即：以5个碱基长度为窗口对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃。将过滤后的数据进行去重，获得无重复的序列，并标记所有序列数量并统计，过滤序列用于后续分析；

步骤三：sRNA分类注释步骤：

将去重后的序列与Rfam数据库进行blast比对，筛选出碱基错配数小于2的结果，注释出其中的非编码RNA序列，

将其余的小RNA序列与miRBase数据库中该物种的miRNA成熟体序列进行比对，筛选出碱基错配数小于2的结果，注释为已知的miRNA序列，同时计算测到的miRNA表达量，进行表达模式分析；

步骤四、miRNA差异分析步骤：

根据上一步注释到的miRNA信息以及表达量结果，使用DESeq进行差异表达分析，并按照差异倍数(FoldChange>2)和显著性(Pvalue<0.05)筛选差异表达的miRNA并绘制图像；

步骤五、miRNA功能和通路分析步骤：

以目标物种的mRNA序列为目标序列，使用psRNATarget或者psRobot软件对差异表达的miRNA序列，进行靶基因位点搜索；

对上一步预测到的miRNA靶基因进行GO功能和KEGG通路的富集分析，获得差异miRNA可能参与的功能和代谢通路；

步骤六、miRNA序列特征分析步骤：

miRNA碱基偏好性分析；

将该物种成熟miRNA序列与近缘物种进行blast比对，筛选出物种间保守的miRNA，并标记其相似度；

对检测到的已知miRNA进行家族归类，并查找相应miRNA家族在其他物种中的存在情况；

步骤七、结果整理步骤：

将所有用于生成miRNA结题报告的统计分析结果进行整理。

在本发明的一个优选实施例中，所述文件准备步骤中，所述的文件包括：下机数据位置以及对应的样本名和分组名、用于差异分析的分组、分析结果保存路径、任务名称、物种简称、测序接头序列、该物种名miRNA的成熟体序列、基因组序列及其index文件的位置、用于功能注释的基因注释文件，mRNA序列、GTF文件中的任意一种或多种。

在本发明的一个优选实施例中，所述下机数据过滤步骤当中，所述的统计为同时对原始数据和过滤数据量进行统计，并以柱状图展示不同长度的序列的数量分布特征。

在本发明的一个优选实施例中，所述sRNA分类注释步骤当中，还包括新的miRNA预测：使用mapper.pl将剩余的序列与基因组进行比对，并使用mireap.pl对比对上的序列进行新的miRNA预测，并使用RNAfold获得结构信息，最后对所有的小RNA序列的注释结果进行统计。

在本发明的一个优选实施例中，所述miRNA差异分析步骤当中，所述绘制图像包括：采用R语言的ggplot2软件包绘制差异表达miRNA的火山图(直观了解差异miRNA的分布情况)和、或MA图(评估文库标准化的好坏)、采用Pheatmap包对差异表达miRNA的表达量绘制热图。

在本发明的一个优选实施例中，所述miRNA碱基偏好性分析为不同长度的miRNA的首位碱基的偏好性，以及所有miRNA每个位置上的碱基偏好性。

本发明的主要创新点在于：

针对植物miRNA的特点采用合适的分析方法。

结果全面，包含涉及到的miRNA分析内容以及其他测到的小RNA信息注释。

自动整理所有分析结果，完成各个部分分析之后，自动对结果进行统计，可视化，以及归类整理，使结果排布一目了然，直接用于报告生成。

所有操作步骤可见，方便错误查询，在进行每一步分析时，都会记录所用到的命令行和参数，以及运行中产生的日志结果，一旦程序运行出错，可以快速检查错误。

附图说明

图1为本发明的流程示意图。

图2为运行日志示意图。

图3为本发明的MA示意图。

图4为本发明的火山图示意图。

图5为本发明的psRobot结果示意图。

图6为本发明的blast结果示意图。

图7为结果目录示意图。

图8为原始数据以及相关统计图表示意图。

图9为去冗余序列以及相关统计图表示意图。

图10为各种小RNA注释结果示意图。

图11为miRNA特征分析结果示意图。

图12为miRNA表达量相关分析结果示意图。

图13为miRNA靶基因预测相关分析结果示意图。

具体实施方式

以下通过实施例对本发明作进一步的说明，但这些实施例不得用于解释对本发明的限制。

在步骤S1)中接受用户的小RNA测序数据，以及相关的数据库信息，然后对所有的数据进行相关的分析，得到每个样本中所有小RNA的注释信息，并对miRNA进行序列特征分析和表达量分析，以及样本间差异表达分析，功能和通路富集分析。

首先是对下机数据进行过滤和数量统计。本发明实施例中，对下机数据进行去除接头和低质量序列的过滤处理，得到高质量的测序结果。作为示例地，采用perl语言脚本去除接头序列(filter_data.py脚本)，并通过5bp的滑动窗口，对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃过滤低质量序列。然后过滤掉长度小于18或者大于36bp的序列。然后对高质量数据的重复序列进行归纳，得到所有的无冗余序列。并对原始数据和高质量进行数量统计。

接下来先通过比对注释出ncRNA序列。作为示例的，使用Blast将这些序列与Rfam数据库比对，注释其他如rRNA，tRNA,snRNA,snoRNA等非编码RNA信息。然后使用perl脚本对结果筛选出碱基错配数小于2的结果，注释出其中的非编码RNA序列。

然后注释出miRNA序列。作为示例的，将其余的小RNA序列与miRBase数据库中该物种的miRNA成熟体序列进行Blast比对，筛选出碱基错配数小于2的结果，注释为已知的miRNA序列，同时计算测到的miRNA表达量，进行表达模式分析。

然后从剩余的序列预测新的miRNA信息。作为示例的，使用mapper.pl将剩余的序列与基因组进行比对，并使用mireap.pl对比对上的序列进行新的miRNA预测，并使用RNAfold获得结构信息。最后对所有的小RNA序列的注释结果进行统计。

对于之前检测到的保守miRNA序列根据其表达量，进行差异表达分析。作为示例的，使用DESeq进行差异表达分析，并按照差异倍数(FoldChange>2)和显著性(Pvalue<0.05)筛选差异表达的miRNA。同时采用R语言的ggplot2软件包绘制差异表达miRNA的火山图(直观了解差异miRNA的分布情况)和MA图(评估文库标准化的好坏)。采用Pheatmap包对差异表达miRNA的表达量绘制热图。

根据序列相似性，对筛选到的显著差异表达的miRNA进行靶基因预测。作为示例的，以本物种的mRNA的3’UTR序列为目标序列,使用psRNATarget或者psRob软件对差异表达的miRNA序列，进行靶基因位点搜索。然后使用R语言通过超几何检验计算靶基因富集到哪些GO功能和KEGG代谢通路上，从而了解这些差异miRNA所发挥的功能。

psRobot结果参见图5，为miRNA：

vvi-miR156a与mRNA：VIT_03s0097g00240.t01的关系预测结果以及位点详情。

blast结果示例参见图6，第一列为本物种miRNA，第二列为与之最匹配的其他物种的miRNA。

本发明还对预测到的保守的miRNA序列进行序列特征分析，包括碱基偏好性分析，保守性分析和家族分析。

作为示例的，采用perl脚本，先对不同长度的miRNA序列，分别统计第一位碱基的种类分布数量；以及所有miRNA每个位置上的碱基种类分布数量，并使用R语言画图展示结果。然后将该物种的miRNA序列与近缘物种进行比对，找出物种间存在的保守性miRNA，并标记之间的相似度。根据每个miRNA的家族信息，找出在近缘物种中是否包含对应家族的miRNA信息。

最终整理所有的分析结果，所所有分析内容按类别排放在不同的目录下。作为示例的，将原始数据单独存放；将数据过滤的统计结果，序列长度分布图形单独存放；将所有小RNA的注释结果及注释结果统计都单独存放；将miRNA序列特征分析结果单独存放；将miRNA表达量以及差异表达相关的分析内容单独存放；将差异表达的miRNA对应的靶基因预测结果，以及功能和通路富集分析结果单独存放，结果目录排布参见7-13。