基于生物信息学探讨骨关节炎患者LncRNA表达及CeRNA互作网络

陈锋 章晓云,2 陈跃平 廖建钊 周毅 李嘉骏 吴炜锋

(1广西中医药大学附属瑞康医院，广西 南宁 530011；2江西中医药大学)

骨关节炎(OA)是一种以关节软骨退变钙化、软骨下骨重塑、骨赘形成、滑膜炎症和血管侵袭为特征的慢性退行性疾病〔1～3〕。伴随全球人口老龄化的不断发展，OA发病率逐年递增，已成为全球范围内的公共健康问题〔4〕。研究表明，OA的发生与创伤、遗传、衰老、肥胖等多种因素有关，但具体发病机制尚不明确。目前临床上非手术治疗方法包括口服硫酸氨基葡萄糖、皮质类固醇、透明质酸钠、多西环素和基质金属蛋白酶抑制剂等，虽然这些药物均有不同程度的疗效，但均无法完全阻断OA的进展〔5,6〕；而手术治疗虽能矫正关节畸形、改善关节功能，但存在术后感染、血栓和二次手术等潜在风险〔7,8〕。因此，全面挖掘OA的潜在发病机制，寻找新的治疗靶点，具有非常重要的意义。

长链非编码RNA(LncRNA)是一类长度大于200nt且不能编码蛋白的RNA。首次发现时认为其没有编码功能，被定义成遗传暗物质。但近年研究证实LncRNA在生命活动过程中发挥着重要的调控作用，与疾病的发生发展密切相关〔9〕。研究表明LncRNA在OA发展过程中可作为竞争性内源RNA(CeRNA)调节微小RNA(miRNA)的表达模式和生物学特征〔10～12〕。即LncRNA可以海绵吸附miRNA来影响信使RNA(mRNA)的表达，从而参与OA的免疫反应、炎症反应、细胞代谢等生物学过程。因此本研究通过GEO数据库挖掘与OA相关的芯片,借助R语言分析差异基因,并构建LncRNA-miRNA-mRNA网络来深入研究LncRNA在CeRNA调控网络中的作用机制。

1 资料和方法

1.1筛选OA的相关芯片 以“Osteoarthritis”、“Genome-wide Expression Profiles”为关键词，在GEO数据库中检索与OA相关的芯片,获得编号为GSE51588的芯片矩阵文件和编号为GPL13497的基因注释文件,该芯片中包含50个样本,包括40例OA患者的关节组织和10例健康对照。

1.2芯片数据分析 利用Perl对数据进行基因重注释，并利用R语言校正、分类，limma包对编码蛋白质的mRNA及LncRNA进行差异分析，设置差异倍数(FC)>0.5，P<0.05，最后运用pheatmap包对存在差异的LncRNA绘制差异热图。

1.3LncRNA、miRNA、mRNA互作预测 利用mircode数据库筛选与差异LncRNA互作的miRNA。再从公共数据库(miRTarBase、TargetScan、miRDB)预测与miRNA互作的靶基因，取3个数据库均能预测到的靶基因，与步骤1.2所得的差异mRNA取交集。整合得到LncRNA-miRNA-mRNA互作网络，并导入Cytoscape绘制CeRNA网络图。

1.4蛋白互作网络构建 为进一步探究OA的Hub基因，将差异基因导入String数据库中，限定研究物种为“Homo Sapiens”获得蛋白互作关系，设置连接评分>0.97，再导出蛋白互作关系的数据文件，选取文件中的node1、node2和连接评分三列导入到Cytoscape中，最后利用Cytoscape软件的“Network Analyzer”工具进行可视化处理，构建PPI网络图，并根据度值筛选出Hub基因。

1.5Hub基因与核心LncRNA的相关性分析 选取CeRNA调控网络中度值最高的LncRNA与其互作的Hub基因，保留各自对应的疾病组数据，利用R语言对数据进行校正后，设置相关系数的绝对值>0.4，P<0.01，最后运用limma包对存在统计学意义的基因绘图。

1.6GO与KEGG富集分析 利用DAVID数据库对差异基因进行GO功能富集分析，研究OA的主要生物功能；对差异基因进行KEGG通路富集分析，研究OA的主要信号通路，P<0.05代表富集结果显著。最后利用R语言绘制GO和KEGG富集分析气泡图。

2 结 果

2.1差异mRNA及LncRNA 利用R语言对芯片进行基因重注释和数据校正后，再对其进行差异分析。结果显示，与健康对照者相比，OA患者的膝关节组织中共存在3 292个明显改变的mRNA，其中上调1 476个，下调1 816个；共存在67个明显改变的LncRNA，其中上调45个，下调22个。在上调和下调的LncRNA中分别选取最显著的5个来绘制差异热图，见图1。

图1 差异LncRNA在不同样本中的表达

2.2CeRNA网络 利用mircode等数据库共预测到208个miRNA和1 936个mRNA，将预测到的mRNA与芯片中的3 292个差异mRNA取交集，再删除不在该互作关系中的miRNA、LncRNA，整合后得到LncRNA-miRNA-mRNA互作网络。将该网络导入Cytoscape中绘制CeRNA网络，见图2。网络中共包括310个节点(54个LncRNA节点、38个miRNA节点、218个mRNA节点)和1 085条边，每条边表示LncRNA与miRNA、miRNA与mRNA之间的关系。该网络中度值排名前5位的LncRNA为钾电压门控通道亚家族Q成员1反义链(KCNQ1OT)1、X 染色体失活特异性转录本基因(XIST)(度值36)、反义Opa反应蛋白5RNA1 (OIP5-AS)1(度值33)、淋巴细胞白血病缺失基因(DLEU)1(度值27)、慢性淋巴细胞白血病上调(CLLU)1(度值27)。

图2 CeRNA调控网络

2.3PPI网络 借助String数据库和Cytoscape软件构建PPI网络图，见图3。图中共包含30个节点、41条边，其中节点代表蛋白基因，边则代表蛋白基因间的互作关系，节点大小及颜色深浅表示其度值的大小，节点越大相应的度值越大，颜色由黄变橙相应的度值则由小变大，边的粗细反映连接评分，边越粗，评分越高。度值排名前5的蛋白基因为信号转导和转录激活剂(STAT)3(度值6)、生长因子受体结合蛋白(GRB)2(度值6)、丝裂原激活的蛋白激酶(MAPK)1(度值6)、叉头家族(FOX)O3(度值5)、原癌基因MYC(度值5)。这些度值较大的蛋白基因在整个网络中起着关键的作用，为该网络的Hub基因。

图3 蛋白互作网络

2.4共表达结果 选取CeRNA中度值最高的LncRNA(KCNQ1OT1、XIST)与Hub基因(STAT3、GRB2、MAPK1、FOXO3、MYC)分析其表达的相关性。除XIST与MAPK1之间呈正相关，其余均呈负相关，选取负相关部分存在统计学意义的数据进行展示，见图4。

图4 共表达结果

2.5GO富集分析和KEGG富集分析 为探究更具体的功能模式,利用DAIVD数据库对关系网络中的218个差异基因进行GO和KEGG富集分析。GO富集分析共得到237个条目,KEGG富集分析共得到45条信号通路,根据P值排序,选取P值最显著的前20个条目以气泡图的形式展现,见图5。图中横坐标为该通路富集的基因占人体总基因的比例，气泡越大代表该通路富集的基因数越多，颜色越红则代表富集程度越显著。分析结果显示,骨关节炎的生物过程主要包括对DNA结合转录因子活性的调控、对营养水平的反应、肌细胞的增殖等；信号通路主要涉及PI3K-AKT信号通路、MAPK信号通路、mTOR信号通路等。

图5 GO和KEGG信号通路可视化结果

3 讨 论

在OA的研究过程中，针对骨代谢、基质降解、氧化应激和免疫炎症反应等方面的治疗均已应用于临床，但其疗效并不理想〔13〕。而LncRNA在介导OA发病机制中起着重要调控作用，通过基因测序技术发现在软骨分化过程中表达量发生改变的LncRNA高达数千个〔14,15〕。因此，深入研究LncRNA在OA发生发展中的作用及其调控机制意义重大。

本研究基于GEO数据库联合DAVID功能富集分析软件，通过生物信息技术筛选调控差异表达miRNA的上游LncRNA及miRNA作用的下游靶基因，旨在探讨OA发展中LncRNA表现出的调控网络模式。CeRNA网络分析得出有54个LncRNA、38个miRNA和218个mRNA与OA发生发展相关。对差异表达的mRNA进行GO和KEGG富集分析，得出这些基因所富集的功能主要包括DNA结合转录因子活性的调控、营养水平的反应、肌细胞的增殖等生物学过程，涉及PI3K-AKT、MAPK、mTOR等信号通路，这些通路在OA发生发展中均起着重要的作用〔16〕。

通过分析共得到5个关键LncRNA(KCNQ1OT1、XIST、OIP5-AS1、DLEU1、CLLU1)，与之互作的mRNA在GO和KEGG通路分析中均显示与OA高度相关。这些结果表明，在OA的发生发展过程中，上述LncRNA可能发挥着重要作用。研究表明，KCNQ1OT1与miR-200a存在潜在结合位点，而在CeRNA网络中Hub基因FOXO3又是miR-200a的潜在靶基因。KCNQ1OT1可通过作用miR-200a影响FOXO3的表达，进而调控细胞自噬、衰老、凋亡和炎症反应等〔17〕；XIST在细胞增殖、分化和基因组维持中起着关键作用，研究表明降低XIST的表达，可显著降低骨髓间充质干细胞(BMSC)的增殖和迁移能力，从而导致炎症反应发生〔18～20〕。研究表明沉默OIP5-AS1能显著抑制细胞增殖活性，通过上调OIP5-AS1可调节miR-410的表达，并进一步调控靶基因KLF10，从而介导PI3K/AKT信号通路，促进细胞增殖〔21,22〕；DLEU1和CLLU1可在转录中、转录后、表观遗传等方面调控基因的表达，参与胶质瘤、骨肉瘤、宫颈癌等疾病的发展进程，但其在OA中的潜在作用尚无相关研究，这就为今后通过作用DLEU1和CLLU1治疗OA提供新的方向和线索〔23,24〕。

本研究结果说明OA的作用机制无法依靠某个特定靶点来调控，且各基因间存在错综复杂的互作关系：LncRNA与miRNA结合，通过隔离mRNA而下调miRNA的结合，最终促进mRNA的表达。STAT3是一类参与成骨细胞增殖、分化、炎症和免疫反应的特异性基因，通过激活STAT3信号通路不仅能刺激MSC成骨分化，还参与Cdc42等生物因子介导的软骨损伤过程。这表明STAT3在OA的软骨下骨硬化和软骨损伤过程中发挥了重要作用，因此抑制STAT3表达是阻止关节损伤进入不可逆阶段的重要干预措施,对治疗OA意义重大〔25,26〕；GRB2作为一种衔接蛋白，经研究证实与PI3K/AKT/NF-κB信号通路的转导存在密切关联，沉默GRB2能够抑制上述通路的活化，进而减少组织的炎性浸润发挥保护作用，GRB2有望成为OA治疗的重要靶点〔27〕；MAPK1参与了细胞增殖分化、应激反应等活动，在细胞外信号转导与细胞反应中发挥着关键作用〔28〕；MYC作为一种转录因子，可调控CyclinD1、CyclinE及P27等细胞周期相关基因的表达，与软骨细胞的增殖、分化和凋亡密切相关〔29〕。本研究表明通过调节MAPK1、FOXO3及MYC水平减少细胞凋亡，并通过调节STAT3及GRB2表达，降低炎症因子等物质的水平，从而改善OA的炎症症状及软骨退变情况。这些特异性指标蛋白的筛选可为进一步研究OA的临床或动物实验提供可靠的检测指标，不仅能有效缩短实验周期，还可降低实验成本。但阐明LncRNA-miRNA-mRNA间的互作关系及LncRNA在OA中的作用机制是未来研究的重点。