CN110021347B - 一种基于miRBase数据库的动物有参的miRNA数据分析方法 - Google Patents

Abstract

本发明公开了一种基于miRBase数据库的动物有参的miRNA数据分析方法，其特征在于，包括如下步骤：准备步骤；下机数据过滤步骤；基因组比对步骤；小RNA注释步骤；miRNA特征分析步骤；miRNA表达量分析步骤；miRNA功能和通路分析步骤；结果整理步骤。本发明的有益效果在于：分析结果全面，包含涉及到的miRNA分析内容以及其他测到的小RNA信息注释，实现自动化分析，只需提供相关数据库信息即可自动快速进行所有分析内容以及结果统计和整理。所有操作步骤可见，方便错误查询，在进行每一步分析时，都会记录所用到的命令行和参数，以及运行中产生的日志结果，一旦程序运行出错，可以快速检查错误。

Description

一种基于miRBase数据库的动物有参的miRNA数据分析方法

技术领域

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

背景技术

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

miRNA的作用机制在动物和植物之间存在明显差异，且有的物种有丰富的miRNA参考数据，但有的物种缺乏参考数据，甚至有些物种没有参考基因组信息，这些情况下的miRNA测序的数据分析方法十分不同。现有的miRNA数据分析工具功能单一且独立，无法实现自动化数据分析，不能进行批量的miRNA数据分析。并且在进行小RNA测序时，除了miRNA还会测到其他的各种小RNA，需要使用不同的数据库信息，对这些小RNA进行注释。目前也没有现成的工具进行相关分析；或者将这些步骤进行关联的工具。

发明内容

为了克服现有技术所存在的上述缺陷，本发明的目的在于提供一种自动化分析方法，通过整合恰当的分析软件，对于在miRBase数据库中有参考数据的动物miRNA测序数据，进行全套的自动化数据分析以及结果整理。

为了实现本发明的目的之一，所采用的技术方案是：一种基于miRBase数据库的动物有参的miRNA数据分析方法，包括如下步骤：

步骤一，准备步骤：

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

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

下机后的原始数据，去除接头，然后过滤低质量序列，过滤序列用于后续分析；

步骤三，基因组比对步骤：

调用miRDeep2中的mapper.pl脚本进行与基因组序列的比对分析，获得每条序列在基因组上的位置信息，并整理比对到每条染色体上的小RNA的序列数量；

步骤四，小RNA注释步骤：

首先调用quantifier.pl将这些序列与miRBase数据中该物种的miRNA序列进行比对，注释出miRNA序列，并得到其表达量，

然后对于piRBase数据库中有相关信息的物种，再使用Blast对其余序列与piRBase数据库比对，注释出piRNA信息，

再使用Blast将其余序列与Rfam数据库比对，注释其他如rRNA，tRNA,snRNA,snoRNA等非编码RNA信息，

再使用mireap对剩余的小RNA序列，进行新的miRNA结构预测。最后对所有的小RNA序列的注释结果进行汇总；

步骤五，miRNA特征分析步骤：

对上一步注释到的miRNA进行着重分析，首先数对miRNA的序列特征进行分析：包括miRNA碱基偏好性分析，即不同长度的miRNA的首位碱基的偏好性，以及所有miRNA每个位置上的碱基偏好性；

统计测到的miRNA序列中发生的碱基编辑事件；将该物种成熟miRNA序列与近缘物种进行blast比对，筛选出物种间保守的miRNA，并标记其相似度；

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

步骤六，miRNA表达量分析步骤：

对miRNA进行表达定量以及差异分析，

根据比对到该物种miRNA成熟体的序列数量，统计测到的miRNA的Reads Count值；

然后对这些miRNA表达量进行样本间矫正，通过miRNA表达密度分布图形，考察miRNA在样本间的表达模式；

然后采用DESeq对miRNA进行差异表达分析，并按照差异倍数(FoldChange>2)和显著性(Pvalue<0.05)筛选差异表达的miRNA，采用R语言的ggplot2软件包绘制差异表达miRNA的火山图和/或MA图；和/或采用Pheatmap包对差异表达miRNA的表达量绘制热图；

步骤七，miRNA功能和通路分析步骤：

以目标物种的mRNA的3’UTR序列为目标序列，使用miRanda软件对差异表达的miRNA序列，进行靶基因位点搜索；

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

步骤八，最终包含结果整理步骤：

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

在本发明的一个优选实施例中，所述步骤一中的读取的文件中包括：下机数据位置以及对应的样本名和分组名，用于差异分析的分组，分析结果保存路径，任务名称，物种简称，测序接头序列，该物种名miRNA的成熟体和前体序列，piRNA位点信息文件，基因组序列及其index文件的位置，用于画图的染色体列表，用于功能注释的基因注释文件，mRNA的3’UTR序列，GTF文件信息。

在本发明的一个优选实施例中，所述步骤二中的过滤低质量序列为：以5个碱基长度为窗口对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃。将过滤后的数据进行去重，获得无重复的序列，并标记所有序列数量。同时对原始数据和过滤数据量进行统计，并以柱状图展示不同长度的序列的数量分布特征。

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

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

实现自动化分析：只需提供相关数据库信息即可自动快速进行所有分析内容以及结果统计和整理。

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

附图说明

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

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

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

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

具体实施方式

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

参见图1，本发明的详细步骤包括：

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

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

其次是将高质量无冗余的测序数据与基因组进行比对。得到所有序列在基因组上的位置信息，并判断测到的数据是否是该物种序列信息。作为示例的，采用预设的miRDeep2中的mapper.pl脚本的默认参数进行基因组比对。统计每个样本中比对到每条染色体正负链上的小RNA序列的数量。

然后将测到的小RNA序列依次与miRNA、piRNA、ncRNA数据库进行比对，以尽可能注释到所有的小RNA信息。作为示例的，先用quantifier.pl将这些序列与miRBase数据中该物种的miRNA序列进行比对，注释出miRNA序列，并得到其表达量，然后对于piRBase数据库中有相关信息的物种，再使用Blast对剩余的序列与piRBase数据库比对，注释出piRNA信息，再使用Blast将剩余的序列与Rfam数据库比对，注释其他如rRNA，tRNA,snRNA,snoRNA等非编码RNA信息。再使用mireap对剩余的序列，进行新的miRNA结构预测。最后对所有的小RNA序列的注释结果进行汇总。

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

然后进行miRNA表达定量以及差异分析。作为示例的，根据比对到每个miRNA上的序列数量作为miRNA的表达量。并采用CPM标准化进行样本间的表达量矫正，使用R语言画出每个样本中miRNA的表达密度分布图，以考察各个样本的miRNA表达模式。然后使用DESeq软件进行miRNA表达差异分析。使用R语言筛选显著差异表达结果(默认筛选FoldChange>2且Pvalue<0.05的结果)，使用ggplot2包绘制差异表达miRNA的火山图(直观了解差异miRNA的分布情况)和MA图(评估文库标准化的好坏)。采用Pheatmap包对差异表达miRNA的表达量绘制热图，参见图2-4。

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

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

Claims (3)

1.一种基于miRBase数据库的动物有参的miRNA数据分析方法，其特征在于，包括如下步骤：

步骤一，准备步骤：

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

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

下机后的原始数据，去除接头，然后过滤低质量序列，过滤序列用于后续分析；

步骤三，基因组比对步骤：

调用miRDeep2中的mapper.pl脚本进行与基因组序列的比对分析，获得每条序列在基因组上的位置信息，并整理比对到每条染色体上的小RNA的序列数量；

步骤四，小RNA注释步骤：

首先调用quantifier.pl将这些序列与miRBase数据中对应物种的miRNA序列进行比对，注释出miRNA序列，并得到其表达量，

然后对于piRBase数据库中有相关信息的物种，再使用Blast对其余序列与piRBase数据库比对，注释出piRNA信息，

再使用Blast将其余序列与Rfam数据库比对，注释rRNA，tRNA,snRNA,snoRNA等非编码RNA信息，

再使用mireap对剩余的小RNA序列，进行新的miRNA结构预测；

最后对所有的小RNA序列的注释结果进行汇总；

步骤五，miRNA特征分析步骤：

对上一步注释到的miRNA进行着重分析，首先数对miRNA的序列特征进行分析：包括miRNA碱基偏好性分析，即不同长度的miRNA的首位碱基的偏好性，以及所有miRNA每个位置上的碱基偏好性；

统计测到的miRNA序列中发生的碱基编辑事件；将该物种成熟miRNA序列与近缘物种进行blast比对，筛选出物种间保守的miRNA，并标记其相似度；

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

步骤六，miRNA表达量分析步骤：

对miRNA进行表达定量以及差异分析，

根据比对到对应物种miRNA成熟体的序列数量，统计测到的miRNA的Reads Count值；

然后对这些miRNA表达量进行样本间矫正，通过miRNA表达密度分布图形，考察miRNA在样本间的表达模式；

然后采用DESeq对miRNA进行差异表达分析，并按照差异倍数：FoldChange>2和显著性：Pvalue<0.05筛选差异表达的miRNA，采用R语言的ggplot2软件包绘制差异表达miRNA的火山图和/或MA图；和/或采用Pheatmap包对差异表达miRNA的表达量绘制热图；

步骤七，miRNA功能和通路分析步骤：

以目标物种的mRNA的3’UTR序列为目标序列，使用miRanda软件对差异表达的miRNA序列，进行靶基因位点搜索；

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

步骤八，最终包含结果整理步骤：

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

2.如权利要求1所述的一种基于miRBase数据库的动物有参的miRNA数据分析方法，其特征在于，所述步骤1中的读取的文件中包括：下机数据位置以及对应的样本名和分组名，用于差异分析的分组，分析结果保存路径，任务名称，物种简称，测序接头序列，该物种名miRNA的成熟体和前体序列，piRNA位点信息文件，基因组序列及其index文件的位置，用于画图的染色体列表，用于功能注释的基因注释文件，mRNA的3’UTR序列，GTF文件信息。

3.如权利要求1所述的一种基于miRBase数据库的动物有参的miRNA数据分析方法，其特征在于，所述步骤2中的过滤低质量序列为：以5个碱基长度为窗口对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃；

将过滤后的数据进行去重，获得无重复的序列，并标记所有序列数量；同时对原始数据和过滤数据量进行统计，并以柱状图展示不同长度的序列的数量分布特征。