环状RNA(circRNA)是一类特殊的非编码RNA分子，呈封闭环状结构，不受RNA外切酶影响，相对于线性RNA分子更加稳定不容易降解^[1]。近年在功能上的研究表明，circRNA分子含有大量的microRNA(miRNA)结合位点，在细胞中起到miRNA海绵的作用从而消除或减弱miRNA对其靶基因的抑制，提高靶基因的表达水平，这种调控模式被称为竞争性内源RNA(ceRNA)^[2]。ceRNA是一种全新的基因表达调控模式，涉及许多RNA分子，包括mRNA、circRNA和miRNA等。有研究表明，ceRNA在肿瘤的发生发展中起着重要作用。

椎间盘变性是腰背痛和坐骨神经痛最常见的原因，而椎间盘退化的主要原因是年龄的增长，然而近年来椎间盘退化的发病群体却有年轻化的趋势^[3]。目前治疗椎间盘退化的方法包括药物治疗、物理治疗和外科手术，但这些都不是最理想的。因此，我们在基因分子的层面上做出假设：在椎间盘变性疾病中，circRNA和mRNA之间存在潜在的ceRNA关系，我们可以通过调节RNA之间的竞争性关系来准确高效地预诊治疗椎间盘变性。

通过miRNA桥梁构建了一个包含mRNA与circRNA两类节点的ceRNA网络，利用聚类分析、模块挖掘等生物信息学方法进行网络分析并且从网络中挖掘出与疾病相关的重要模块，利用Gene Ontology(GO)基因功能注释、KEGG通路富集等方法对circRNA的功能进行了预测。为了探索重要模块中的蛋白编码基因与circRNA的分子功能，在NCBI PubMed中搜集大量文献来证明重要模块中基因的功能及作用。在ceRNA表达调控关系的基础上推测出circRNA潜在的生物学意义和功能，从根本上预测椎间盘变性中基因表达水平和分子调控机制，为椎间盘变性相关药物的开发和临床预诊治疗提供了重要的分子标志物。

1 实验材料与方法 1.1 数据来源及处理

在NCBI GEO数据库中下载了人类椎间盘变性疾病的芯片表达谱数据(GSE67567)，其中包含同样本的mRNA和circRNA的表达谱。总共10个样本，包括五个正常样本和五个髓核退化疾病样本。其中mRNA表达谱来自GPL15314 Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version)平台，circRNA表达谱来自GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1平台。在GENCODE数据库中下载人类编码蛋白转录序列Fasta格式文件，并且在Circbase数据库中下载人类circRNA的序列文件。运用blast软件对mRNA序列与芯片探针序列进行了比对，得到探针-mRNA的对应关系后进行平台转换后得到基因的表达谱。与此同时，对circRNA的探针表达谱进行平台转换得到circRNA表达谱。

1.2 构造竞争性内源RNA(ceRNA)网络

在mirbase数据库上下载人类成熟的miRNA序列，利用miRanda^[4]软件预测miRNA-mRNA, 与circRNA-miRNA的的靶向关系。Score选取默认值>140，Energy选取更为严谨的域值< 0。

1.2.1 ceRNA关系预测

这两类关系对中的miRNA进行超几何检验。具体计算公式如下：

$ P - value = 1 - \sum\limits_{i = 0}^{r - 1} {\frac{{\left( {\frac{t}{i}} \right)\left( {\frac{{m - t}}{{n - i}}} \right)}}{{\left( {\frac{m}{n}} \right)}}} $

其中m代表人类miRNA总数，t表示与mRNA相互作用的miRNA的数目，n表示与circRNA相互作用的miRNA的数目，r表示mRNA与circRNA之间共享的miRNA数目。通过对mRNA与miRNA的靶向关系对和circRNA与miRNA的关系对做超几何检验，我们得到ceRNA对(p < 0.05)。

研究发现circRNA的过量表达会降低游离miRNA分子的数目，进而导致miRNA下游靶基因的高表达^[5]。共表达是一种发现调控关系的常用方法，我们在正常样本和疾病样本中分别计算mRNA与circRNA的pearson相关性，得到2 539 352个边对关系和对应的Pearson相关系数cor值。考虑到算法的统计学意义和网络的规模等因素，我们认为mRNA与circRNA在疾病与正常样本中变化程度大于0.5的关系对符合共表达关系。数学公式如下:

$ \left| {{\rm{cor}}\left( {{\rm{case}}} \right){\rm{ - cor}}\left( {{\rm{control}}} \right)} \right| > 0.5 $

根据最新研究结果，circRNA和mRNA之间共表达的增加或缺失可能是导致疾病的原因，所以我们同样假设这种竞争关系可以导致椎间盘变性疾病的发生^[2]。

为了使网络更加趋向于疾病特异性和相关性，既考虑到了通过Pearson相关性构成的共表达网络，又兼顾到了miRanda预测的miRNA潜在靶点构成的ceRNA网络，选取他们关系对的交集构成mRNA与circRNA的ceRNA二步网络进行网络功能与拓扑分析，为后续的研究做基础。

1.3 模块挖掘与功能注释通路富集

运用cytoscape 3.2.0插件jActiveModules进行模块挖掘，参数选取默认值(Number of Modules:5, Overlap Threshold:0.8)。并且使用了cytoscape3.2.0插件ClueGO对模块中包含的基因GO功能注释和KEGG通路富集。

方法流程图见图 1.

2 结果分析 2.1 表达谱处理和分析

运用blast软件对mRNA序列和芯片探针序列进行了比对，进行平台转换后我们得到含有11 502个编码基因的表达谱。同时对circRNA的探针表达谱进行平台转换得到含有2 793个circRNA的表达谱。我们通过聚类分别对mRNA的表达谱和circRNA的表达谱进行进一步分析，结果表明两个表达谱在疾病和正常样本中的差异性明显，样本数据较好，见图 2。

2.2 网络分析结果

通过以上算法生成了一个由mRNA与circRNA两类节点组成的二步网络。网络包含2 349个mRNA，264个circRNA，5 705个ceRNA关系对。其中有两类节点，黄色圆形为mRNA，蓝绿色菱形为circRNA。我们构建的网络服从幂率分布(见图 3(b))。

对网络进行拓扑属性分析，其中包括度分布、网络介数、最短路径和网络亲密度。由图 3(a)可见网络服从幂率分布，说明网络中存在一小部分的基因或者circRNA连接着大部分节点，网络连通性较强，同时也说明这小部分节点在网络调控中起重要作用。

2.3 识别ceRNA网络重要模块

使用cytoscape3.2.0的jActiveModules插件进行模块挖掘，发现了五个circRNA与mRNA连接紧密的模块。模块的规模都很相似，模块节点数为32-70，边为37-132。利用Cytoscape的ClueGO插件分别对每一个模块中的所有基因进行功能注释和通路富集，选取了P值< 0.05的GOTerm和Pathway。其中主要包含了几类生物学过程：相关炎症反应、骨发育和损伤、蛋白质结构和稳定性等等。

2.4 模块展示与分析 2.4.1 模块一

模块一(见图 4(a))包括36个节点、45条边，包含了19个mRNA和7个circRNA。其中编码蛋白的基因富集到蛋白质稳定期(见图 4(c))。通过对椎间盘蛋白质组学差异分析，我们发现正常人和腰椎间盘退变患者的蛋白质种类存在显著差异，腰椎间盘退变患者的蛋白质含量也随着病变程度的加深而逐渐减少^[6]。同时我们也发现了内源性凋亡信号通路的调控，研究表明人椎间盘髓核细胞的凋亡具有内在性和外源性^[7]。同时椎间盘变性与骨的发育相关，此模块也富集到了骨发育功能^[8-9]。

因PAPSS2是AR brachyolmia的疾病基因，并且PAPSS2突变可以产生从brachyolmia到脊椎干骺发育不良的表型分级。PAPS是体内活性硫酸根供体，它的合成需要在PAPS合成酶(PAPss)的催化下进行，PAPS合成酶2(PAPSS2)的突变与大骨节病、骨关节炎等多种骨骼疾病的发生发展有关^[10-16]。在此模块中，GARS与PAPSS2同时富集到骨发育的生物学过程。基于ceRNA理论，推断GARS基因与PAPSS2及他们靶向的circRNAs(has_circ_0013204，has_circ_0000077，has_circ_0073027，has_circ_0001234，has_circ_0007311)，可能通过ceRNA调控关系诱发椎间盘变性等骨疾病^[17]。发现模块二中的PAPSS2与has_circ_0000077的竞争性调控机制改变了共表达的强度从而影响了基因的表达和相应蛋白的功能。在正常样本中他们的共表达相关值为-0.633 56，在疾病状态下他们的表达相关值为0.642 102，在正常和疾病中体现了明显的共表达差异，并且由于竞争性关系从负相关变为正相关也暗示了竞争性RNA调控机制的表达趋势。并且他们的超几何P值统计学显著为1.02×10^-6，FDR矫正值为7.913×10^-3，证明了通过hsa-let-7a-2-3p^[18]介导的PAPSS2与has_circ_0000077之间紧密的ceRNA关系(见图 4(b))。其中得到PAPSS2基因的Foldchange值为1.899 515 547在疾病中显著上调。

2.4.2 模块二

模块二(见图 5(a))包括41个节点、44条边，包含了22个mRNA和19个circRNA。我们发现了”转化生长因子β受体信号通路、淋巴细胞分化正调节、白细胞分化正调节、T细胞选择阳性、跨质膜导入、白细胞活化负调节、细胞酰胺代谢过程负调节等生物学过程(P值< 0.05)。转化生长因子β1可以改善退变的椎间盘细胞活性，从而刺激胶原基因表达来延缓椎间盘组织的退变，因此对退变早期的椎间盘尤其是髓核具有修复功能，可逆转椎间盘的退变^[19]。也有实验证实腺病毒载体介导人转化生长因子β1基因转染退行性变的椎间盘细胞能够有效提高髓核细胞内转化生长因子β1含量，阻逆椎间盘退行性变^[20-22]。

同时椎间盘变性还伴随着相关的炎症和免疫反应，许多文献都指出在腰椎椎间盘突出症中，髓核边缘有明显的炎症和自身免疫反应，在破裂型间盘突出中有比较明显的免疫炎症反应，这种反应的主要特征是巨噬细胞浸润^[23]；髓核组织暴露于自身免疫系统中可以激活T、B细胞，引起自身免疫反应，血管内皮生长因子参与了退变椎间盘新生血管的形成，并可能在椎间盘突出与腰腿痛机制中与T、B细胞协同作用引起自身免疫反应；髓核边缘有散在的炎性细胞和较多巨噬细胞，但小血管增生和T淋巴细胞浸润不明显，提示有炎性反应，但自身免疫反应不如椎间盘突出症明显^[24-25]。这两种免疫调控机制分别与淋巴细胞分化的正向调节、白细胞分化的正向调节、阳性T细胞选择有密切联系^[26-30]。

其中参与这些与免疫炎症相关的生物学过程的基因CD74已经被文献证实参与了椎间盘内部调控过程。文献[31-34]中解释说明MIF / CD74途径可以代表治疗椎间盘退变的关键靶标，MIF可以通过CD74促进炎症反应，因为用其拮抗剂ISO-1抑制MIF的功能可以减轻MIF诱导的炎症反应并且发挥有效的治疗效果，以此来减轻椎间盘疾病患者的疼痛，缩短疾病的治愈周期^[34-40]。与CD74同时富集到几个重要的免疫炎症相关通路的PRKCZ基因和他们靶向的circRNAs(has_circ_0000077，has_circ_0081342，has_circ_0000389，has_circ_0000815)很可能同时作用，遵循ceRNA的竞争调控关系来引发椎间盘疾病(见图 5(c))。

药物方面有文献证实烟酰胺有助于减轻压力对椎间盘的损伤，能够促进压力损伤后的椎间盘恢复。烟酰胺可以抑制IL-1β诱导的椎间盘细胞凋亡，改善IL-1β导致的椎间盘能量代谢障碍，烟酰胺可以促进髓核细胞的增殖并抑制白细胞介素1β诱导的髓核细胞凋亡，通过对烟酰胺的负调控来影响椎间盘的变性^[41-43]。

模块二中的CD74与has_circ_0001320也同样可能通过ceRNA的关系改变了共表达关系的强度。在正常样本中他们的共表达相关值为0.868 441，在疾病状态下他们的表达相关值为-0.049 57。在正常和疾病中同样有着明显的共表达差异，并且由于竞争性关系从负相关变为正相关同样暗示了竞争性RNA调控机制的表达趋势。他们的超几何P值统计学显著为2.62×10^-6，FDR矫正值为0.015 134，表明通过hsa-miR-17-3p^[44]介导的CD74与has_circ_0001320有着潜在的ceRNA关系(图 5(b))。CD74基因的Foldchange值为0.389 194 149在疾病中显著下调。

图 5(a)为模块二网络图，我们把文献证明已知疾病相关基因用红色圆形节点表示；把富集到相同生物学过程的基因用黄色圆形节点表示；把黄色和红色节点的一步邻居circRNAs用紫色菱形节点表示并且把它们共同的一步邻居用绿色菱形节点表示；其他的基因用蓝色圆形节点表示，circRNA用橙色菱形节点表示。图 5(b)举例一个ceRNA关系对的MiRanda靶点展示。图 5(c)模块中所有编码基因Goterm功能和Pathway通路富集结果，颜色代表P值大小。

3 讨论及结论

1) 没有约束差异基因是为了排除常规差异分析方法的不确定性，以免影响实验结果。取而代之的是分析circRNA和mRNA之前差异的共表达相关性，从而更准确的分析两者的竞争性关系。

2) 通过cytoscape的一个网络模块挖掘插件识别出了五个与椎间盘变性高度相关的功能模块并且在功能模块中通过网络拓扑关系发现疾病潜在的编码基因和circRNA。

3) 改进了单一的共表达方法，更关注于基因在疾病和差异样本中差异的共表达调控关系。

4) 为椎间盘疾病的预测、药物靶点治疗和circRNA功能的探索提供了新思路。