基于生物信息学的骨关节炎差异基因谱预测

祝乃强 侯静怡 刘金欣 马桂云

(1承德医学院附属医院脊柱外二科，河北 承德 067000;2承德医学院)

骨关节炎(OA)多见于中老年人，65岁以上老年人患病率高达50%〔1〕。近年来，研究发现OA发病机制并不是一组靶点或一些基因间的相互作用关系出现了明显改变，而是一组靶点或一些基因间的相互作用关系出现了明显改变从而导致其相关功能模块或生物网络发生故障，其中包括软骨细胞外基质的酶降解，新细胞外基质合成缺乏，细胞死亡与凋亡，软骨细胞的异常活化与分化等〔2，3〕。因此，寻找新的更有意义的生物标志物及其参加的信号通路，在预防与治疗OA疾病发生和发展过程中起到重要作用。随着高通量技术的快速发展，生物数据呈现指数级的增长，以基因芯片为代表的高通量生物检测技术得到广泛应用，其可在短时间内获得大量的疾病相关基因信息，使得疾病在基因层面进行全面而细致的分析成为可能〔4〕。本研究采用生物信息学方法分析OA基因芯片数据，目的在于寻找OA病理过程中的关键基因及其参与信号通路，为OA早期诊断、早期治疗干预、药物研发提供新的思路。

1 材料与方法

1.1数据收集与预处理 以“Osteoarthritis”作为检索词，在NCBI的基因表达汇编(GEO)数据库检索相应的基因芯片，下载符合实验要求的4套数据集，其系列号为GSE12021、GSE55235、GSE55457和GSE55584，均来自GPL96平台，芯片类型为Affymetrix Human Genome U133a芯片。GSE12021芯片数据集由Huber等〔5〕于2008年提交，该数据集包括13例正常滑膜组织样本，20例OA组织样本，20例类风湿性关节炎滑膜组织样本；GSE55235、GSE55457和GSE55584数据集是由Woetzel等〔6〕于2014年提交的，共包含20例正常滑膜组织样本，26例OA滑膜组织样本，33例类风湿性关节炎滑膜组织样本。本研究共选取4个数据集中的正常滑膜组织样本和OA滑膜组织样本作为后续分析的样本。

1.2差异表达基因的筛选分析 利用R语言中的 “sva”分析包〔7〕对GSE12021、GSE55235、GSE55457和GSE55584数据集中的数据进行合并及批次校正，利用“limma”分析包〔8〕对上述4个数据集进行分析，以P<0.05，|log2FC|>1为筛选条件选择数据集的差异表达基因。

1.3差异基因表达的基因本体论(GO)功能富集和京都基因与基因组百科全书(KEGG)通路富集分析 GO分析用于高通量基因组或转录组数据鉴定的生物学机制，包括分子功能(MF)、生物学过程(BP)、细胞组分(CC)〔9〕；KEGG 数据库用于识别候选目标的功能及生物相关性〔10〕。本研究利用R语言中的“clusterProfiler”富集分析包〔11〕对差异基因进行GO功能富集分析和KEGG通路富集分析。

1.4蛋白-蛋白互作网络(PPI)构建及关键基因分析(PPI) 将所有的差异基因导入到STRING线上数据库(https：//string-db.org/)〔12〕，进行蛋白相互作用分析，并利用Cytoscape软件对PPI进行可视化。应用Cytoscape软件中的Cytohubba插件对网络进行拓扑学分析，以度值(degree)前10位的差异基因作为PPI中的关键基因。DisGeNET数据库(http：//www.disgenet.org/web/DisGeNET/me，Version 5.0)〔13〕是包含与人类疾病相关的基因和变体的平台之一。在 DisGeNET 数据库中选用基因进行检索，将上述靶点依次导入到数据库中，获取靶点类型信息(Protein Class)。

1.5关键模块分析及功能富集分析 利用Cytoscape软件中的插件复杂分子检测(MCODE)对PPI进行模块分析，标准设置：节点得分截止值=0.2，K核心值(K-Core)=2，度值截止=2，计算MCODE的得分值；以MCODE分数>6分作为显著性模块的筛选标准，并且利用DAVID(https：//david.ncifcrf.gov/)在线数据库对筛选出的模块中的基因进行通路富集分析〔14〕。

2 结 果

2.1差异表达基因的筛选结果 GSE12021、GSE55235、GSE55457和GSE55584的数据集共包含33例正常滑膜组织样本和40例OA滑膜组织样本，筛选出符合条件的差异基因共320个，其中上调基因172个，下调基因148个。

2.2差异基因的GO功能富集和KEGG通路富集分析 将差异基因进行GO富集分析，其中上调的差异基因共有151个条目富集，主要包括加双氧酶活动、铁离子结合、细胞黏合力、泛素蛋白连接酶集合等，；下调的差异基因共有216个条目富集，主要包括抑制转运因子结合、核苷酸糖跨膜转运活动，激素受体结合等，详细内容见图1、图2。

图1 上调差异基因的GO功能富集分析

图2 下调差异基因的GO功能富集分析

KEGG通路富集结果表明，上调的差异基因共有159个条目富集，主要包括丝裂原活化蛋白激酶(MAPK)信号通路、脂肪酸代谢通路、凋亡通路等；下调的基因共有133个条目富集，主要包括Janus激酶/信号转导与转录激活子(JAK-STAT)信号通路、原发性免疫通路、造血细胞谱系通路等，详细内容见图3、4。

图3 上调差异基因的KEGG通路富集分析

2.3PPI及关键基因分析 将筛选出的差异基因导入到STRING在线数据库中，得到差异基因的PPI，该网络包含313个节点，2 160条边。如图5所示，该网络包含167个上调基因和144个下调基因。应用Cytoscape软件的Cytohubba插件，分析网络的拓扑学参数，且以degree值的大小表示该基因在此网络中的重要性，如若degree越大，表明与其关联的基因越多，则其基因越重要。

图4 下调差异基因的KEGG通路富集分析

(白色节点代表上调基因，灰色节点代表下调基因)图5 差异基因的蛋白-蛋白互相作用网络

以degree前10位(degree > 40)的基因作为网络的关键基因，所筛选的关键基因为：连环蛋白(CTNNB)1，糖原合成酶激酶(GSK)-3β，信号转导和转录激活因子(STAT)1，Ras同系物集组成员(RHO)C，组蛋白去乙酰化酶(HDAC)9，早老素(PSEN)1，赖氨酸脱甲基酶(KDM)5C，β-分泌酶(BACE)1，JAK 3和卡林蛋白(CUL1)；所涉及的OA过程中有信号通路、转录因子、蛋白(钙结合蛋白)、酶(蛋白酶、激酶、酶调节剂、连接酶)。

2.4功能模块分析 应用Cytoscape插件MCODE寻找差异基因PPI的关键模块，其中符合MCODE得分大于4分的模块共有3个，模块1的MCODE得分为7.867，16个节点和59条边。见图6。模块1中的基因进行功能富集分析，模块1中基因主要参与6条富集通路，其中BP 3个条目，包括转录/DNA模板化、神经元分化负调剂、RNA聚合酶Ⅱ启动子的转录正调节；CC 1个条目-质膜外侧和MF 2个条目，包括RNA聚合酶Ⅱ核心启动子序列、转录因子活性-DNA结合；KEGG通路富集分析显示，模块1中的基因主要参与5条信号通路，如JAK-STAT信号通路、趋化因子、磷脂酰肌醇3激酶(PI3K)-丝氨酸-苏氨酸激酶(Akt)信号通路等。见表1。

(浅色节点代表上调基因，深色节点代表下调基因)图6 功能模块的PPI

表1 功能模块的功能富集和通路富集分析

FDR：错误发现率

3 讨 论

OA是最常见的一种全关节疾病，肥胖、年龄、过度运动、炎症、遗传及创伤等多种因素与OA的发生密切联系〔15〕。目前，OA的病因及发病机制尚不明确。分析OA的潜在发病机制对于诊断预后和鉴定药物治疗的靶点至关重要。由于高通量测序和微阵列技术可同时提供有关人类基因组中数千种基因表达水平的信息，因此该方法已广泛用于预测OA的潜在诊断和治疗靶点。

本研究进行差异基因分析表明OA发病机制的复杂性，其并不是简单的单个基因/靶点单独其作用，而是多基因/靶点相互作用的复杂结果。KEGG通路富集结果显示，这些差异基因主要涉及MAPK信号通路，脂肪酸代谢通路，凋亡通路、JAK-STAT信号通路，原发性免疫通路，造血细胞谱系通路等。MAPK信号通路能够调控软骨细胞的增殖、凋亡、胞外基质代谢、炎症因子分泌等过程，在OA病理过程中发挥重要作用〔16，17〕。Ge等〔16〕研究表明，3-磷酸肌醇依赖性蛋白激酶(PDK)1可通过p38 MAPK信号通路促进OA的软骨细胞凋亡。而Peng等〔18〕研究进一步发现OA中的双特异性磷酸酶(DUSP)1过表达抑制了MAPK信号通路的激活和OA相关基质的表达。研究表明，JAK/STAT信号通路受到促炎细胞因子基因表达而显著上调的调控。Lim等〔19〕研究表明，MMP-13在IL-1β处理的人软骨肉瘤细胞系中的表达诱导了JAK2和STAT1/STAT2激活及MMP-13基因表达，其被泛酪氨酸激酶抑制剂AG490阻断。然而，也有研究发现〔20〕，用IL-1β处理这种软骨细胞系也激活了p38-MAPK/c-Fos/AP-1。

利用STRING在线数据库对这些差异基因进行PPI网络分析得到位于中心的关键靶点，其中，STAT1作为调节巨噬细胞炎症反应中不可或缺的组成部分，如若激活的STAT1被磷酸并易位到细胞核中以上调炎症的基因〔21〕。JAK3作为一种非受体酪氨酸激酶，可与细胞因子受体复合物的常见γ链(IL2RG)相互作用〔22〕。

PI3K-Akt信号通路对于生理和病理条件下细胞生长和存活的许多方面都至关重要〔23〕。研究表明，细胞内信号转导途径MAPK和PI3K-Akt已被证实在调节炎症介质和骨关节炎发病机制中的细胞外基质改变中都是至关重要的〔24，25〕。其中IL-1β能够激活MAPK的三种途径(ERK，p38，JNK)及PI3K/Akt信号通路，随着诱导分解代谢的表达，可进一步促进细胞外基质破坏〔26〕，表明MAPK和PI3K/Akt的调节在IL-1β处理的软骨细胞中分解代谢因子和软骨基质改变的调节中起关键作用〔27〕。趋化因子包含免疫细胞分泌的一小部分细胞因子或蛋白，且在免疫反应过程中调节和控制免疫系统细胞〔28〕。鉴于趋化因子在炎症期间白细胞运输和活化过程中起关键性作用，趋化因子很可能在OA的发育发挥重要作用〔29〕。研究表明，在OA液中检测到许多趋化因子可以促进细胞类型的转化，包括CCL1848和CX3CL1〔30，31〕。