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

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

背景技术：

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

miRNA的作用机制在动物和植物之间存在明显差异，且有的物种有丰富的miRNA参考数据，但有的物种缺乏参考数据，甚至有些物种没有参考基因组信息，这些情况下的miRNA测序的数据分析方法十分不同。由于不同物种中的miRNA有一定的保守性，因此对于没有本物种参考miRNA数据的测序结果，也可以进行分析。但是目前还没有针对无参考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的3’UTR序列或mRNA序列为目标序列，使用miRanda软件或psRobot软件对差异表达的miRNA序列，进行靶基因位点搜索；

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

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

对miRNA碱基偏好性进行分析；

步骤七，结果整理步骤：

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

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

在本发明的一个优选实施例中，所述sRNA分类注释步骤当中，所述miRNA的命名方式为采用物种简称-miRNA家族名称的命名方式。

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

在本发明的一个优选实施例中，所述miRNA差异分析步骤中，所述绘制图像包括采用R语言的ggplot2软件包绘制差异表达miRNA的火山图、MA图；采用Pheatmap包对差异表达miRNA的表达量绘制热图。

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

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

针对无参考miRNA数据的miRNA测序数的分析方法。

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

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

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

附图说明

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

图2为准备文件示意图。

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

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

图5为本发明的热图示意图。

图6为每个序列首位碱基的分布情况示意图。

图7为所有序列每一位碱基的分布情况示意图。

具体实施方式

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

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

首先是对下机数据进行过滤和数量统计。本发明实施例中，对下机数据进行去除接头和低质量序列的过滤处理，得到高质量的测序结果。作为示例地，去除接头序列，并通过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预测的步骤。

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

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

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

参见图6、7，图6为每个序列首位碱基的分布情况；图7为所有序列每一位碱基的分布情况；

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