CN107516021A - 一种基于高通量测序的数据分析方法 - Google Patents

Abstract

本发明公开一种基于高通量测序的数据分析方法。其中，所述方法包括：获取参考测序样品和对照测序样品，获得三种的RNA的差异DNA甲基化区域相关的表达水平出现差异的RNA，获得甲基化水平和表达水平高度相关的RNA，构成候选竞争性内源RNA调控关系对，生成候选竞争性内源RNA网络，筛选出hub节点，并将所述hub节点进行癌症的体细胞突变收录数据库注释。本发明提供的基于高通量测序的数据分析方法，将hub节点在癌症的体细胞突变收录数据库进行注释，发现癌症相关的竞争性内源RNA调控关系，预测癌症发生发展过程总的调控机制，提高了对癌症致病机理预测的准确性。

Description

一种基于高通量测序的数据分析方法

技术领域

本发明涉及生物医学技术领域，具体涉及一种基于高通量测序的数据分析方法。

背景技术

癌症是起源于上皮组织的恶性肿瘤，是全球发病和死亡的主要原因之一。癌症不仅极大的降低患者的生活质量，并且对患者造成身体和经济上的巨大伤害。早期分子诊断和综合治疗水平的提高使得其死亡率虽有下降，但是其发病机制目前尚不清楚。

随着高通量测序技术的迅速发展，高通量RNA测序成本的快速下降，作为一项颠覆性的技术二代测序已经产生了种类繁多的癌症预测的分子标记。在生物体的生命过程中，基因表达致力于将遗传信息转变为具有生物活性和功能的蛋白质。而这种表达的改变可能导致蛋白质功能的改变甚至疾病的发生。长链非编码RNA((long noncoding RNA,以下简称lncRNA)无论是在生理还是病理状态下都能调节许多重要的生物学进程。而环状RNA(circular RNA，以下简称circRNA)作为一类特殊的非编码RNA也已经成为最新的研究热点。大量研究表明，lncRNA和circRNA分子都富含微RNA(micro RNA，以下简称miRNA)结合位点，能够在细胞中充当miRNA海绵进而消除miRNA对其靶基因的抑制作用。DNA甲基化作为重要的表观遗传修饰之一，对基因表达具有重要的调控作用，在癌症的发生发展中可以作为一类关键分子标记。但仅仅对某种RNA或者DNA甲基化进行研究，并不足够来解释疾病中的生物学过程和调控机制。

因此，如何提出一种方法，分析癌症发生发展过程中DNA甲基化与多种RNA的关联来全局的阐述癌症的致病机制，从而提高分析结果的准确性成为业界亟待解决的重要课题。

发明内容

针对现有技术中的缺陷，本发明提供一种基于高通量测序的数据分析方法。

本发明提出一种基于高通量测序的数据分析方法，包括：

获取参考测序样品和对照测序样品，其中，所述参考测序样品来源于正常样本，所述对照测序样品来源于异常样本；

根据对所述参考测序样品和所述对照测序样品分别进行高通量测序获取的三种RNA中每个RNA的第一表达水平、每个胞嘧啶C位点的第一DNA甲基化水平和所述每个RNA的第二表达水平、所述每个胞嘧啶C位点的第二DNA甲基化水平，以及第一预设规则获得所述三种RNA中与差异DNA甲基化区域相关的所述表达水平出现差异的RNA；

根据对所述参考测序样品和所述对照测序样品进行高通量测序获得所述三种RNA中每个RNA的甲基化水平、所述第一表达水平和所述第二表达水平，以及第二预设规则获得甲基化水平与表达水平高度相关的RNA；

根据对所述参考测序样品和所述对照测序样品进行高通量测序获取到的每个miRNA的序列信息和所述三种RNA中每个RNA的序列信息，以及第三预设规则获得构成候选竞争性内源RNA调控关系对；

根据所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA、所述甲基化水平与表达水平高度相关的RNA、所述构成候选竞争性内源RNA调控关系对，以及第四预设规则生成候选竞争性内源RNA网络；

根据所述候选竞争性内源RNA网络，以及第五预设规则获得hub节点，并将所述hub节点进行癌症的体细胞突变收录数据库注释。

本发明提供的基于高通量测序的数据分析方法，通过对差异DNA甲基化区域相关的三种的RNA的表达水平的差异分析，甲基化水平和表达水平的相关性分析，挑选出构成候选竞争性内源RNA调控关系对，并形成候选竞争性内源RNA网络，从而筛选出hub节点，将hub节点在癌症的体细胞突变收录数据库进行注释，发现癌症相关的竞争性内源RNA调控关系，预测癌症发生发展过程总的调控机制，提高了对癌症致病机理预测的准确性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一实施例基于高通量测序的数据分析方法的流程示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。下列实施例中未注明具体条件的实验方法，通常按照常规条件或按照厂商所建议的条件实施检测。

图1为本发明一实施例基于高通量测序的数据分析方法的流程示意图，如图1所示，本发明实施例提供的基于高通量测序的数据分析方法，包括：

S101、获取参考测序样品和对照测序样品，其中，所述参考测序样品来源于正常样本，所述对照测序样品来源于异常样本；

具体地，可以通过癌症患者的癌旁正常组织获得正常样本，同时通过所述癌症患者的癌组织获得异常样本，在对所述正常样品和异常样本进行处理分别获得参考测序样品和对照测序样品，以供后续进行高通量测序。

S102、根据对所述参考测序样品进行高通量测序获取的三种RNA中每个RNA的第一表达水平和每个胞嘧啶C位点的第一DNA甲基化水平、对所述对照测序样品获得的所述每个RNA的第二表达水平和所述每个胞嘧啶C位点的第二DNA甲基化水平，以及第一预设规则获得所述三种RNA中与差异DNA甲基化区域相关的所述表达水平出现差异的RNA；

具体地，对所述参考测序样品处理后进行高通量测序获取三种RNA中每个RNA的第一表达水平和每个胞嘧啶C位点的第一DNA甲基化水平，对所述对照测序样品进行高通量测序获取所述每个RNA的第二表达水平和每个胞嘧啶C位点第二DNA甲基化水平，并根据所述每个RNA的所述第一表达水平和所述第二表达水平，所述每个胞嘧啶C位点的所述第一DNA甲基化水平和所述第二DNA甲基化水平，以及第一预设规则获得所述三种RNA中与差异DNA甲基化区域相关的所述表达水平出现差异的RNA。其中，所述三种的RNA可以是长链非编码RNA(long noncoding RNA,以下简称lncRNA)、信使RNA(messenger RNA，以下简称mRNA)和环状RNA(circular RNA，以下简称circRNA)。

所述第一预设规则可以包括：首先，通过DESeq或者EBSeq软件对所述每个RNA的所述第一表达水平和所述第二表达水平进行差异性分析，获得所述RNA的表达差异倍数值和FDR值，如果某一所述RNA的所述表达差异倍数值大于等于第一阈值并且所述FDR值小于第二阈值，那么所述RNA为表达水平出现差异的RNA；其中，所述第一阈值和所述第二阈值可以根据实际经验进行设定，本发明实施例不做限制，例如所述第一阈值为2，所述第二阈值为0.01。其次，采用MOABS软件，对所述参考参考测序样品和所述对照测序样品中相同区域的所述第一DNA甲基化水平和所述第二DNA甲基化水平进行差异化分析，该软件基于二项层次模型，可以获得所述区域的单位点覆盖度和平均差异，如果所述单位点覆盖度大于第三阈值且所述平均差异大于第四阈值，那么所述区域构成差异DNA甲基化区域(Differentialmethylation region，以下简称DMR)；其中，所述第三阈值和所述第四阈值可以根据实际经验进行设定，本发明实施例不做限制，例如所述第三阈值为10，所述第四阈值为0.2；最后，如果所述DMR的起始位置和终止位置与所述表达水平出现差异的RNA的起始位置和终止位置上下游2000bp存在交集，那么所述RNA为所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA。

S103、根据对所述参考测序样品和所述对照测序样品进行高通量测序获得所述三种RNA中每个RNA的甲基化水平、所述第一表达水平和所述第二表达水平，以及第二预设规则获得甲基化水平与表达水平高度相关的RNA；

具体地，对所述参考测序样品和所述对照测序样品处理后进行高通量测序，可以获得所述三种RNA中每个RNA的甲基化水平，根据所述每个RNA的所述甲基化水平、所述第一表达水平和所述第二表达水平以及第二预设规则可以获得甲基化水平与表达水平高度相关的RNA。

所述第二预设规则可以包括：对所述每个RNA的所述甲基化水平、所述第一表达水平和所述第二表达水平进行皮尔森相关性分析，例如采用R语言的cor.test函数来进行皮尔森相关性分析，可以获得每个RNA其表达水平和甲基化水平的相关性值和显著性值；如果某一所述RNA的所述相关性值的绝对值大于第五阈值并且显著性值小于第六阈值，那么该RNA的甲基化水平与表达水平高度相关，所述RNA即为所述甲基化水平与表达水平高度相关的RNA。所述第五阈值和所述第六阈值可以根据实际经验进行设定，本发明实施例不做限制，例如所述第五阈值为0.7，所述第六阈值为0.01。

S104、根据对所述参考测序样品和所述对照测序样品进行高通量测序获取到的每个miRNA的序列信息和所述三种RNA中每个RNA的序列信息，以及第三预设规则获得构成候选竞争性内源RNA调控关系对；

具体地，利用miRDeep2软件处理经过高通量测序的所述参考测序样品和所述对照测序样品的数据，再经过参考基因比对，获得每个miRNA的序列信息，对所述三种RNA进行参考基因比对，获得每个RNA的序列信息。根据所述每个miRNA的序列信息、所述每个RNA的序列信息，以及第三预设规则获得构成候选竞争性内源RNA调控关系对。

所述第三预设规则包括：将每个miRNA的序列信息和所述每个RNA的序列信息提供给miRNA的第一预设数量的靶点预测软件对miRNA的靶RNA进行预测，例如采用miranda和RNAhybrid，分别对miRNA的靶lncRNA,mRNA,circRNA进行预测。可理解的是，采用的所述靶点预测软件的数量越多，获得的预测结果也会越准确。每个所述靶点预测软件可以预测出所述每个miRNA的靶RNA，如果某一个RNA被所有的所述靶点预测软件预测为某个miRNA的所述靶RNA，那么所述RNA成为所述miRNA的候选靶点。如果任意两个成为候选靶点的RNA之间共享miRNA的数量大于第七阈值，并且所述两个成为候选靶点的RNA之间的校正后的超几何分布P值小于第八阈值，那么所述两个成为候选靶点的RNA构成候选竞争性内源RNA调控关系对。其中，所述第七阈值和所述第八阈值根据实际经验进行设定，本发明实施例不做限定。

例如，所述第七阈值等于5，所述第八阈值等于0.01，有两个成为候选靶点的RNA，a和b，a与30个miRNA相关，b与20个miRNA相关，a和b共享10个miRNA，所述miRNA的数量共有2000个，根据上述数据30、20、10和2000，经过编程可以计算得出a和b之间的超几何分布P值。对a和b的超几何分布P值进行错误发现率校正得到校正后的超几何分布p值。如果所述校正后的超几何分布p值小于0.01，由于a和b共享miRNA的数量为10，因此a和b构成候选竞争性内源RNA调控关系对。

S105、根据所述与差异甲基化区域相关的所述表达水平出现差异的RNA、所述甲基化水平与表达水平高度相关的RNA、所述构成候选竞争性内源RNA调控关系对，以及第四预设规则生成候选竞争性内源RNA网络；

具体地，根据步骤102获得的所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA、步骤103获得的所述甲基化水平与表达水平高度相关的RNA、步骤104获得的所述构成候选竞争性内源RNA调控关系对，以及第四预设规则生成候选竞争性内源RNA网络。

所述第四规则包括：所有的所述构成候选竞争性内源RNA调控关系对可以形成一个集合，从这个集合中挑选出部分RNA调控关系对，挑选出的所述RNA调控关系对中的至少一个RNA既是所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA，又是所述甲基化水平与表达水平高度相关的RNA，所述挑选出的部分RNA调控关系对可以形成一个网络，在这个网络中每个RNA与周围至少一个RNA相连，这个网络被称为候选竞争性内源RNA网络。

S106、根据所述候选竞争性内源RNA网络，以及第五预设规则获得hub节点，并将所述hub节点进行癌症的体细胞突变收录数据库注释。

具体地，在所述候选竞争性内源RNA网络中，每个RNA都是一个节点，根据第五预设规则在所述候选竞争性内源RNA网络中筛选出hub节点。所述第五预设规则包括在所述候选竞争性内源RNA网络中，与临近的第一预设数量RNA相连的所述RNA为所述hub节点，所述第一预设数量根据实际经验进行设定，本发明实施例不做限制，例如所述第一预设数量为10个。或者按照与临近的RNA相连的所述RNA的数量的多少由高到低对所述RNA进行排序，排名前第九阈值的所述RNA为所述hub节点，所述第九阈值根据实际经验进行设定，本发明实施例不做限制，例如所述第九阈值为10％。

获得的所述hub节点，由于其高连通度这一特性，它的改变通常被认为容易导致疾病的发生和发展，将所述hub节点进行癌症的体细胞突变收录数据库(Catalogue ofsomatic mutations in cancer，以下简称COSMIC)注释，即对于所述hub节点，统计每个hub节点与所述COSMIC中的基因的关联情况，所述hub节点如果与所述COSMIC中的基因的关联度非常高或者所述hub节点本身就是一个在所述COSMIC中已经注释的基因，那么说明所述hub节点与癌症的发生发展有着密不可分的作用。

本发明提供的基于高通量测序的数据分析方法，通过对三种的RNA的差异DNA甲基化区域相关的表达水平出现差异的分析，甲基化水平和表达水平的高相关分析，挑选出构成候选竞争性内源RNA调控关系对，并形成候选竞争性内源RNA网络，从而筛选出hub节点，将hub节点在癌症的体细胞突变收录数据库进行注释，发现癌症相关的竞争性内源RNA调控关系，预测癌症发生发展过程总的调控机制，提高了对癌症致病机理预测的准确性。

在上述实施例的基础上，进一步地，所述第一预设规则包括：

对所述参考测序样品和所述对照测序样品中相同的RNA的所述第一表达水平与所述第二表达水平进行差异性分析，获得所述RNA的表达差异倍数值和FDR值，若判断获知所述表达差异倍数值大于等于第一阈值且所述FDR值小于第二阈值，则所述RNA为表达水平出现差异的RNA；

对所述参考测序样品和所述对照测序样品中相同区域的所述第一DNA甲基化水平和所述第二DNA甲基化水平进行差异化分析，获得所述区域的单位点覆盖度和平均差异，若判断获知所述单位点覆盖度大于第三阈值且所述平均差异大于第四阈值，则所述区域构成差异DNA甲基化区域；

根据所述表达水平出现差异的RNA与所述差异DNA甲基化区域获得所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA。

具体地，首先，通过DESeq或者EBSeq软件对所述每个RNA的所述第一表达水平和所述第二表达水平进行差异性分析，获得所述RNA的表达差异倍数值和FDR值，如果某一所述RNA的所述表达差异倍数值大于等于第一阈值并且所述FDR值小于第二阈值，那么所述RNA为表达水平出现差异的RNA；其中，所述第一阈值和所述第二阈值可以根据实际经验进行设定，本发明实施例不做限制，例如所述第一阈值为2，所述第二阈值为0.01。

其次，采用MOABS软件，对所述参考测序样品和所述对照测序样品中相同区域的所述第一DNA甲基化水平和所述第二DNA甲基化水平进行差异化分析，该软件基于二项层次模型，可以获得所述区域的单位点覆盖度和平均差异，如果所述单位点覆盖度大于第三阈值且所述平均差异大于第四阈值，那么所述区域构成差异DNA甲基化区域(Differentialmethylation region，以下简称DMR)；其中，所述第三阈值和所述第四阈值可以根据实际经验进行设定，本发明实施例不做限制，例如所述第三阈值为10，所述第四阈值为0.2；最后，如果所述DMR区域的起始位置和终止位置与所述表达水平出现差异的RNA的起始位置和终止位置上下游2000bp存在交集，那么所述RNA为所述与差异甲基化区域相关的所述表达水平出现差异的RNA。

在上述各实施例的基础上，进一步地，所述第二预设规则包括：

将所述每个RNA的所述甲基化水平、所述第一表达水平和所述第二表达水平进行皮尔森相关性分析获得所述每个RNA的表达水平与甲基化水平之间相关性值和显著性值，若判断获知所述相关性值的绝对值大于第五阈值且所述显著性值小于第六阈值，则所述RNA为甲基化水平与表达水平高度相关的RNA。

具体地，对所述每个RNA的所述甲基化水平、所述第一表达水平和所述第二表达水平进行皮尔森相关性分析，例如采用R语言的cor.test函数来进行皮尔森相关性分析，可以获得每个RNA的表达水平与甲基化水平之间相关性值和显著性值；如果某一所述RNA的所述相关性值的绝对值大于第五阈值并且所述显著性值小于第六阈值，那么该RNA的甲基化水平与表达水平高度相关，所述RNA即为所述甲基化水平与表达水平高度相关的RNA。所述第五阈值和所述第六阈值可以根据实际经验进行设定，本发明实施例不做限制，例如所述第五阈值为0.7，所述第六阈值为0.01。

在上述各实施例的基础上，进一步地，所述第三预设规则包括：

根据所述每个miRNA的序列信息和所述每个RNA的序列信息，通过第一预设数量的靶点预测软件，分别获得每个靶点预测软件预测出的所述每个miRNA的靶RNA，若判断获知所述RNA被所有所述靶点预测软件预测为某个miRNA的靶RNA，那么所述RNA成为所述miRNA的候选靶点；

若判断获知任意两个成为候选靶点的RNA之间共享miRNA的数量大于第七阈值，并且所述两个成为候选靶点的RNA之间的校正后的超几何分布P值小于第八阈值，那么所述两个成为候选靶点的RNA构成候选竞争性内源RNA调控关系对。

具体地，将每个miRNA的序列信息和所述每个RNA的序列信息提供给miRNA的第一预设数量的靶点预测软件对miRNA的靶RNA进行预测，例如采用miranda和RNAhybrid，分别预测对miRNA的靶lncRNA,mRNA,circRNA进行预测。可理解的是，采用的所述靶点预测软件的数量越多，获得的预测结果也会越准确。每个所述靶点预测软件可以预测出所述每个miRNA的靶RNA，如果某一个RNA被所有的所述靶点预测软件预测为所述靶RNA，那么所述RNA成为所述miRNA的候选靶点。如果任意两个成为候选靶点的RNA之间共享miRNA的数量大于第七阈值，并且所述两个成为候选靶点的RNA之间的校正后的超几何分布P值小于第八阈值，那么所述两个成为候选靶点的RNA构成候选竞争性内源RNA调控关系对。其中，所述第七阈值和所述第八阈值根据实际经验进行设定，本发明实施例不做限定。

例如，所述第七阈值等于5，所述第八阈值等于0.01，有两个成为候选靶点的RNA，a和b，a与30个miRNA相关，b与20个miRNA相关，a和b共享10个miRNA，所述miRNA的数量共有2000个，根据上述数据30、20、10和2000，经过编程可以计算得出a和b之间的超几何分布P值。对a和b的超几何分布P值进行错误发现率校正得到校正后的超几何分布的p值。如果所述校正后超几何分布p值小于0.01，由于a和b共享miRNA的数量为10，因此a和b构成候选竞争性内源RNA调控关系对。

在上述各实施例的基础上，进一步地，所述第四预设规则包括：

从所有的所述构成候选竞争性内源RNA调控关系对中挑选出部分RNA调控关系对，生成候选竞争性内源RNA网络；其中，挑选出的所述RNA调控关系对满足所述RNA调控关系对中的至少一个RNA为所述与差异甲基化区域相关的所述表达水平出现差异的RNA和所述甲基化水平与表达水平高度相关的RNA。

具体地，所有的所述构成候选竞争性内源RNA调控关系对可以形成一个集合，从这个集合中挑选出RNA调控关系对，挑选出的所述RNA调控关系对中的至少一个RNA既是所述与差异甲基化区域相关的所述表达水平出现差异的RNA，又是所述甲基化水平与表达水平高度相关的RNA，所述挑选出的部分RNA调控关系对可以形成一个网络，在这个网络中每个RNA与周围至少一个RNA相连，这个网络被称为候选竞争性内源RNA网络。

在上述实施例的基础上，进一步地，所述第五预设规则包括：

在所述候选竞争性内源RNA网络中，与临近的第一预设数量RNA相连的所述RNA为所述hub节点或者按照与临近的RNA相连的所述RNA的数量的多少由高到低对所述RNA进行排序，排名前第九阈值的所述RNA为所述hub节点。

具体地，在所述候选竞争性内源RNA网络中，与临近的第一预设数量RNA相连的所述RNA为所述hub节点，所述第一预设数量根据实际经验进行设定，本发明实施例不做限制，例如所述第一预设数量为10个。或者按照与临近的RNA相连的所述RNA的数量的多少由高到低对所述RNA进行排序，排名前第九阈值的所述RNA为所述hub节点，所述第九阈值根据实际经验进行设定，本发明实施例不做限制，例如所述第九阈值为10％。

在上述各实施例的基础上，进一步地，所述RNA的甲基化水平通过对RNA区域上多个甲基化位点进行加权算法计算获得。

具体地，获取每个所述RNA在CG、CHG和CHH位点的甲基化水平，再采用加权算法计算每个所述RNA的甲基化水平，即所述甲基化水平等于每个所述RNA所有位点的甲基化覆盖度除以每个所述RNA所有位点的总覆盖度。

在上述各实施例的基础上，进一步地，还包括：

在基于Starbase的竞争性内源RNA关系中查找出所述注释过的所述hub节点对应的竞争性内源RNA关系。

具体地，Starbase是一个存储了高通量实验数据支持的microRNA靶标并整合多个靶标预测软件的交集和调控关系的综合型数据库，在基于Starbase的竞争性内源RNA调控关系对中，提取注释过的所述hub节点，以便对这些与癌症发生发展高度相关的所述hub节点与哪些相连的RNA起到相互的调控作用进行研究。

在上述各实施例的基础上，进一步地，所述方法还包括：

在候选竞争性内源RNA网络中查找出对应的所述注释过的所述hub节点对应的竞争性内源RNA关系。

具体地，在获得的候选竞争性内源RNA网络中，提取注释过的所述hub节点，以便对这些与癌症发生发展高度相关的所述hub节点与哪些相连的RNA起到相互的调控作用进行研究。

在上述各实施例的基础上，进一步地，所述三种RNA包括lncRNA、mRNA和circRNA。

本发明提供的基于高通量测序的数据分析方法，通过对三种的RNA的差异DNA甲基化区域相关的表达水平出现差异的分析，甲基化水平和表达水平的高相关分析，挑选出构成候选竞争性内源RNA调控关系对，并形成候选竞争性内源RNA网络，从而筛选出hub节点，将hub节点在癌症的体细胞突变收录数据库进行注释，发现癌症相关的竞争性内源RNA调控关系，预测癌症发生发展过程总的调控机制，提高了对癌症致病机理预测的准确性。

下面通过一具体的实施例对本发明提出的基于高通量测序的数据分析方法进行说明。

优选地，本发明实施例采用6组RRBS测序数据进行分析，包括每组数据均包括lncRNA测序数据，mRNA测序数据、circRNA测序和miRNA测序数据。其中，mRNA,lncRNA的表达水平为RPKM值，circRNA的表达水平为read计数。上述6组数据，分别通过三个不同个体的正常组织和疾病组织获得。

分别将来自于相同个体的2组数据进行对比，获得3组表达水平的差异RNA和差异DNA甲基化区域，合并3组结果后，共获得3249个表达水平出现差异的mRNA，1098个表达水平出现差异的lncRNA，165个表达水平出现差异的circRNA，以及7858个DMR。通过判断每个所述DMR区域的起始位置和终止位置与所述表达水平出现差异的RNA的起始位置和终止位置上下游2000bp是否存在交集，获得与DMR相关的表达水平出现差异RNA集合，包含1996个mRNA、599个lncRNA和63个circRNA。

通过计算每组数据的所有RNA的甲基化水平，可以发现lncRNA和mRNA在基因体区域的甲基化水平相对偏高，同时转录起始位点附近呈现出低甲基化，而circRNA的甲基化水平并未呈现出特定的趋势。将每组数据中的每个RNA的所述RNA的甲基化水平与RNA表达水平进行皮尔森相关性分析，获得相关性值的绝对值大于0.7，显著性值小于0.01的RNA。将六组数据中满足上述条件的RNA合并，获得甲基化水平与表达水平高度相关的RNA共包括1321个mRNA，515个lncRNA，13个circRNA。

利用miRDeep2软件处理miRNA测序数据，在经过参考基因比对，获得miRNA的序列信息。通过其他测序数据经过各自的比对软件处理后获得lncRNA、mRNA和circRNA的序列信息。根据miRNA的序列信息，lncRNA、mRNA和circRNA的序列信息和采用miranda和RNAhybrid预测miRNA靶点，结果如表1所示，可见circRNA具有最多的miRNA，暗示了circRNA更易充当miRNA海绵来行使作用。

表1 有靶点的miRNA数目预测统计

类型	miRNA	circRNA	lncRNA	mRNA
					已知的	1098	1096	595	788
预测的	173	33	83	101
					总计	1271	1129	678	889

对于所有mRNA，lncRNA，circRNA，根据所述第三预设规则，满足任意两个RNA之间存在5个共享miRNA且校正后的超几何分布P值小于0.01的、能够构成候选竞争性内源RNA(Competing endogenous RNA，以下简称ceRNA)调控关系对，包含141415个候选ceRNA调控关系对。

对所有mRNA，lncRNA，circRNA进行筛选，同时满足既是所述与DMR相关的表达水平出现差异的RNA，又是甲基化水平与表达水平高度相关的RNA，共包括323个mRNA和19个lncRNA。从所有的ceRNA关系对中提取323个mRNA和19个lncRNA以及与这些RNA相关的RNA，构成候选竞争性内源RNA网络。

在上述候选竞争性内源RNA网络中共包含17291个癌症相关ceRNA关系对，包含1862个RNA，其中1107个mRNA，125个lncRNA，630个circRNA。按照与临近的RNA相连的所述RNA的数量的多少由高到低对上述1862个RNA进行排序，排名前10％的RNA被提取出来，共计186个RNA，包括90个mRNA，3个lncRNA和93个circRNA，上述186个RNA通常被认为是候选的致病位点。

进一步发现有3个RNA注释在COSMIC癌症数据库中，其中有一个基因与癌症的关系是已经证实的，有2个基因虽然注释在数据库中但并未与该癌症类型相关。通过癌症相关ceRNA网络找到显著影响COSMIC中的基因的RNA(与COSMIC中相关基因的连通度>5)。他们可能是更有潜力影响癌症发生发展的RNA。

最后，分别从候选竞争性内源RNA网络以及starbase数据库的ceRNA网络中提取COSMIC中被注释的这3个RNA和显著影响COSMIC基因的RNA的调控关系子网。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (10)

1.一种基于高通量测序的数据分析方法，其特征在于，包括：

获取参考测序样品和对照测序样品，其中，所述参考测序样品来源于正常样本，所述对照测序样品来源于异常样本；

根据对所述参考测序样品进行高通量测序获取的三种RNA中每个RNA的第一表达水平和每个胞嘧啶C位点的第一DNA甲基化水平、对所述对照测序样品获得的所述每个RNA的第二表达水平和所述每个胞嘧啶C位点的第二DNA甲基化水平，以及第一预设规则获得所述三种RNA中与差异DNA甲基化区域相关的所述表达水平出现差异的RNA；

根据对所述参考测序样品和所述对照测序样品进行高通量测序获得所述三种RNA中每个RNA的甲基化水平、所述第一表达水平和所述第二表达水平，以及第二预设规则获得甲基化水平与表达水平高度相关的RNA；

根据对所述参考测序样品和所述对照测序样品进行高通量测序获取到的每个miRNA的序列信息和所述三种RNA中每个RNA的序列信息，以及第三预设规则获得构成候选竞争性内源RNA调控关系对；

根据所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA、所述甲基化水平与表达水平高度相关的RNA、所述构成候选竞争性内源RNA调控关系对，以及第四预设规则生成候选竞争性内源RNA网络；

根据所述候选竞争性内源RNA网络，以及第五预设规则获得hub节点，并将所述hub节点进行癌症的体细胞突变收录数据库注释。

2.根据权利要求1所述的方法，其特征在于，所述第一预设规则包括：

对所述参考测序样品和所述对照测序样品中相同的RNA的所述第一表达水平与所述第二表达水平进行差异性分析，获得所述RNA的表达差异倍数值和FDR值，若判断获知所述表达差异倍数值大于等于第一阈值且所述FDR值小于第二阈值，则所述RNA为表达水平出现差异的RNA；

对所述参考测序样品和所述对照测序样品中相同区域的所述第一DNA甲基化水平和所述第二DNA甲基化水平进行差异化分析，获得所述区域的单位点覆盖度和平均差异，若判断获知所述单位点覆盖度大于第三阈值且所述平均差异大于第四阈值，则所述区域构成差异DNA甲基化区域；

根据所述表达水平出现差异的RNA与所述差异DNA甲基化区域获得所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA。

3.根据权利要求1所述的方法，其特征在于，所述第二预设规则包括：

将所述每个RNA的所述甲基化水平、所述第一表达水平和所述第二表达水平进行皮尔森相关性分析获得所述每个RNA的表达水平与甲基化水平之间的相关性值和显著性值，若判断获知所述相关性值的绝对值大于第五阈值且所述显著性值小于第六阈值，则所述RNA为所述甲基化水平与表达水平高度相关的RNA。

4.根据权利要求1所述的方法，其特征在于，所述第三预设规则包括：

根据所述每个miRNA的序列信息和所述每个RNA的序列信息，通过第一预设数量的靶点预测软件，分别获得每个靶点预测软件预测出的所述每个miRNA的靶RNA，若判断获知所述RNA被所有所述靶点预测软件预测为某个miRNA的靶RNA，那么所述RNA成为所述miRNA的候选靶点；

若判断获知任意两个成为候选靶点的RNA之间共享miRNA的数量大于第七阈值，并且所述两个成为候选靶点的RNA之间的校正后的超几何分布P值小于第八阈值，那么所述两个成为候选靶点的RNA构成所述候选竞争性内源RNA调控关系对。

5.根据权利要求1所述的方法，其特征在于，所述第四预设规则包括：

从所有的所述构成候选竞争性内源RNA调控关系对中挑选出部分RNA调控关系对，生成候选竞争性内源RNA网络；其中，挑选出的所述RNA调控关系对满足所述RNA调控关系对中的至少一个RNA为所述与差异DNA甲基化区域相关的所述表达水平出现差异的RNA和所述甲基化水平与表达水平高度相关的RNA。

6.根据权利要求1所述的方法，其特征在于，所述第五预设规则包括：

在所述候选竞争性内源RNA网络中，与临近的第一预设数量RNA相连的所述RNA为所述hub节点或者按照与临近的RNA相连的所述RNA的数量的多少由高到低对所述RNA进行排序，排名前第九阈值的所述RNA为所述hub节点。

7.根据权利要求1所述的方法，其特征在于，所述RNA的甲基化水平通过对RNA区域上多个甲基化位点进行加权算法计算获得。

8.根据权利要求1所述的方法，其特征在于，还包括：

在基于Starbase的竞争性内源RNA关系中查找出所述注释过的所述hub节点对应的竞争性内源RNA关系。

9.根据权利要求1所述的方法，其特征在于，还包括：

在所述候选竞争性内源RNA网络中查找出对应的所述注释过的所述hub节点对应的竞争性内源RNA关系。

10.根据权利要求1至9任一项所述的方法，其特征在于，所述三种RNA包括lncRNA、mRNA和circRNA。