基于生物信息学及机器学习鉴定骨关节炎滑膜相关基因

杨延辑，廖文波

(1.遵义医科大学附属医院 脊柱外科，贵州 遵义 563099；2.遵义医科大学第二附属医院 骨外科，贵州 遵义 563000)

骨性关节炎(Osteoarthritis,OA)是由多种因素引起的以关节软骨退行性变,当骨关节炎发展到一定程度时形成不可逆疾病，世界上超过 10% 的人口患有OA，65岁以上人群患病率更高，OA 造成的经济和社会负担正在迅速增加，并严重影响患者的生活质量[1]。关节置换是关节炎目前终末期的主要治疗手段，其中膝骨关节炎(Knee osteoarthritis,KOA)患病率高、花费高是导致世界残疾的第四大原因[2]，目前临床诊断大多靠影像学、检验学、查体等；随着社会人口的老龄化 ,该病的发生率越来越高，目前的药物治疗大多是对症的，迄今为止还缺乏改善疾病的 OA 药物[3]；因此骨关节炎早期诊断、及时防治有其重要现实意义[4]。骨关节炎的发病机制目前尚不明确，以往骨关节炎的研究热点大多集中于软骨组织，现研究发现滑膜、成纤维样滑膜细胞和软骨下骨细胞也参与了骨关节炎的发病过程[5]。本研究拟通过生物信息学方法，对GEO公共数据库中有关OA滑膜的基因表达谱芯片数据进行差异表达基因分析、功能富集、尝试用机器学习的方式筛选骨关节炎关键基因，探讨OA发生发展中调节基因，为探索其诊断和治疗靶点提供生物学信息依据。

1 材料与方法

1.1 数据的下载及处理 本研究中从公共数据库GEO(GENE EXPRESSION OMNIBUS,GEO)下载骨关节炎滑膜相关芯片数据集。筛选条件为：①骨关节炎；②人类；③滑膜组织；④无药物及手术干预。基于筛选条件下载表达数据集GSE1919，GSE55235，GSE82017，GSE55457；GSE1919包含5个OA 滑膜组织样本和5个正常滑膜组织样本，GSE55235包含10个OA滑膜组织样本和10个正常滑膜组织样本，GSE82017包含10个OA滑膜组织样本和7个正常滑膜组织样本，GSE55457包含10个OA滑膜组织样本和10个正常滑膜组织样本，利用R软件读取下载相关基因的原始数据，对数据芯片进行预处理，将GSE1919，GSE55235两个数据集作为训练集进行数据合并，将GSE82017,GSE55457作为两个独立验证集，对样本中都不表达的探针、存在的缺失值或基因与多个探针存在对应关系等特殊情况，进行标准化、探针过滤、缺失值填充以及探针合并等，并对合并数据表达矩阵进行log2对数转换，对数据进行标准化处理。

1.2 差异基因及功能富集 利用R软件(4.1.2版本)Limma包进行差异表达分析，设定P<0.05,基因表达差异倍数(FoId change，Fc)绝对值≥1为筛选条件，筛选获得差异表达基因(Differentially expressed gene,DEG)。利用cluster-Profiler、org.Hs.eg.db及enrichplot包对获得的差异基因进行基因本体论(Gene Ontology ,GO)和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)富集分析，以P<0.05为筛选标准；并使用ggplot2包进行可视化处理。

1.3 机器学习筛选特征基因 机器学习(Machine learning,ML)是一门研究计算机如何模拟人类进行学习数据分析的科学[6]，在处理高维度、大批量数据等方面较传统方法有显著优势，目前已在临床获得广泛应用，对上述获得的差异基因通过R软件使用lasso回归算法[7]和SVM-RFE 支持向量机递归特征消除算法[8]两种机器学习方式筛选特征基因，LASSO 是一种回归算法，SVM-RFE 由Guyon等[9]提出，用于癌症分类中基因的选择，在临床已得到应用，如自闭症识别[10]、前列腺组织病理学分级等[11]。对两种方式获得的基因结果取交集，应用R软件“venneuler”包绘制韦恩图对结果进行可视化。

1.4 诊断标志物的验证 对上述取交集获得的特征基因，在GSE82017,GSE55457两个独立的数据集中来验证获得的基因作为诊断标志物的价值，通过绘制接受者操作特征曲线(Receiver operating characteristic，ROC)曲线评价其诊断价值，以P<0.05为阈值来确定。计算ROC曲线下面积(Area under the ROC，AUC)，AUC的取值在0～1，AUC越大，说明预测性能越好。

2 结果

2.1 差异基因筛选 对GSE1919及GSE55235两个数据进行合并，矫正、标准化处理后得到8 920个基因，进行差异分析共获得474个差异基因，其中226个表达上调基因，248个表达下调基因，绘制火山图并对表达排名前100的差异基因绘制热图(见图1)。

(红色高表达，蓝色低表达)

2.2 GO及KEGG富集结果 GO富集分析发现生物过程(Biological process,BP)主要富集在正向调节白细胞的激活，正向调节细胞激活、附着力，细胞趋化性，对脂多糖的反应，淋巴细胞活化的正调节等，分子功能(Molecular function,MF)中含胶原蛋白细胞外基质、等离子外侧膜、膜筏、膜微区、胞吞泡富集明显，细胞组分(Cellular component,CC)显示与受体配体活性、信号受体激活剂活性、糖胺聚糖结合、硫化合物结合等相关；KEGG富集分析显示其在脂质和动脉粥样硬化、1型人类 T 细胞白血病病毒感染、类风湿关节炎、MAPK信号通路、细胞-细胞因子受体相互作用、NF-kappa B 信号通路、破骨细胞分化、AGE-RAGE 信号通路、肿瘤坏死因子信号通路、IL-17信号通路、趋化因子信号通路等(见图2)。

图2 GO、KEGG富集分析结果

2.3 机器学习筛选OA标志物 对获得的474差异基因进一步通过机器学习的方式筛选特征基因，lasso算法选择lambda.min参数[12]获得14个特征基因分别为BCL6、DDIT4、KLF9、GADD45A、SIK1、HNRNPA1、MRC2、MTHFD2、KDELR3、CX3CR1、SCRG1、SLC2A3、ABCC3、ABL2，通过SVM-RFE算法得到19个特征基因分别为SELL、SCRG1、HLA-DMB、BCL6、UCP2、HLA-DMA、MIR6883、MXRA5、DDIT4、NFIL3、TNFAIP3、SORL1、MIR8071-2、STC1、NEDD9、KLF4、FAM107A、SNORD10、MAFF(见图3)，将两种方式获得的基因取交集，获得3个特征标志基因BCL6、DDIT4、SCRG1，绘制韦恩图(见图4)。

图3 机器学习Lasso结果及VSM-RFE结果

图4 两种机器学习结果取交集得到的3个基因

2.4 特征基因的验证价值 对获得的BCL6，DDIT4，SCRG1特征基因，在GSE82017,GSE55457两个独立的数据集中进行外部验证3个基因的诊断价值。ROC曲线结果显示BCL6在两个数据集的AUC值分别为0.9及0.83；SCRG1在两个数据集的AUC值分别为0.829及0.93；DDIT4的AUC值分别为0.757及0.7，3个诊断基因AUC值均>0.5，结果显示其均有较高诊断价值(见图5)。

A:GSE82107中验证结果;B:GSE55457中验证结果。

3 讨论

骨性关节炎是退行性关节疾病，多因素复杂作用所致；目前仍缺乏一种能够得到广泛认可的、临床上普遍使用的OA诊断的生物标记物。为了探索骨关节炎的诊断标志物，本研究通过从公共数据库GEO下载数据集GSE1919，GSE55235，GSE82017，GSE55457；共包含35个OA滑膜组织样本和，32个正常滑膜组织，将GSE1919，GSE55235作为训练集合并数据后筛选获得DEGs；对DEGs进行富集分析，其富集结果中现已有大量研究证实核转录因子kB(Nuclear factor-kapa B,NF-kB)信号通路的激活在OA的发生，发展中起着重要作用; NF-kB信号通路的激活增加炎症因子的释放，导致关节软骨的破坏，以及滑膜炎的发生[13]，MAPK 信号通路在早期 OA 中的调节作用，抑制MAPK信号通路可缓解OA的进展[14]，使用lasso回归算法和SVM-RFE 支持向量机递归特征消除算法得到的结果取交集获得3个特征基因，分别为BCL6，DDIT4，SCRG1。BCl6是位于3q27染色体的转录抑制因子[15]，其参与细胞的增殖、分化、凋亡及炎症过程等[16]。目前研究发现，BCL6 表达水平的升高提示患者滑膜炎症加重状态，在类风湿关节炎患者中BCL6的变化与患者病情程度及膜炎病理评分呈正相关[17]。90%的骨关节炎患者存在滑膜病变，且滑膜的病变程度与关节严重的疼痛和功能障碍有关[18]，滑膜炎症可促进促炎因子和疼痛神经递质的产生[19]。此外，滑膜炎的发生还可能促进软骨退变，滑膜炎通过释放促炎介质和软骨破坏因子引起软骨损伤，从而在OA中发挥引发作用[20]；Hou等[21]研究证实，BCL6高表达可以抑制NFATc1的活化来调节RANKL，而在OA中OPG/RANK/RANKL 调节系统的功能障碍与软骨下骨的组织学改变之间存在密切相关性，RANKL 不仅参与影响软骨下骨的通路，而且在骨关节炎软骨中表达，RANKL 的产生会导致软骨破坏[22]，并且XU等[23]通过实验证实，抑制RANKL 诱导的破骨细胞生成，能改善软骨细胞炎症减缓OA的进展，综上推测BCL6可能在OA起着一定作用，或可成为潜在治疗靶点，但目前BCL6大多集中在肿瘤，在其他疾病领域的研究近期才逐渐增多，仍需进一步临床实验研究来证实。目前已有研究发现，SCRG1在人关节软骨中特异性表达[24]，间充质干细胞(MSCs)是非造血基质细胞，具有自我更新和分化成间充质细胞的能力[25]，研究发现SCRG1正向调节hMSC自我更新、迁移和成骨分化和成软骨分化[26]，SCRG1是一种促进软骨基因表达的转录因子，ADAMTS9-AS2高表达可以逆转 miR-942-5p对SCRG1的抑制[27]。新近研究发现，SCRG1通过 Wnt5a 促进脐带间充质干细胞的成软骨分化[28]；骨关节炎可导致软骨撕裂和软骨细胞丢失，软骨变薄，磨损退化[4]，而SCRG1在人关节软骨中特异性表达，其可能在OA的发生发展中起到一定作用，但目前整体关于SCRG1在各个疾病领域的研究还相对较少，其在骨关节炎中的具体作用机制还需进一步研究。DDIT4是一种肿瘤相关蛋白，也是一种代谢和免疫相关蛋白，应激反应蛋白，可通过mTOR 途径调节蛋白质合成、自噬和细胞凋亡；在化疗、缺氧和DNA损伤等应激条件下高度表达[29-30],DDIT4也被称为发育和DNA损伤反应调控(REDD1)，关节软骨退化是OA发生的主要因素之一，涉及mTOR信号传导增加和自噬缺陷，REDD1 在正常人关节软骨中高度表达，并在衰老和 OA 期间减少，REDD1是mTOR的抑制剂；通过建立OA实验小鼠研究发现，REDD1缺乏会增加小鼠实验性OA的严重程度，并且在滑膜、半月板和软骨下骨中观察到类似的差异[31]，另有研究发现REDD1 在正常人关节软骨中高度表达，并在衰老和 OA 期间减少，通过实验观察发现REDD1在人和小鼠软骨以及小鼠半月板和滑膜中的表达随着年龄的增长而降低[32]。通过从铁死亡相关基因从FerrDb数据库(http://zhounan.org/ferrdb)[33]查询DDI4是铁死亡相关基因；研究发现REDD1 表达被发现是许多病理的早期生物标志物，包括炎症性疾病，REDD1 的表达与细胞凋亡、活性氧(ROS)的堆积和导致组织损伤的炎症激活有关，其可通过NF-κB、一氧化氮合酶和超氧化物歧化酶/谷胱甘肽过氧化物酶/NAPDH 氧化酶途径激活炎症和 ROS 产生[34]。并且新进研究发现，mTORC1作为mTOR的复合物，实验证明铁螯合通过多种途径抑制 mTORC1，而铁对 mTORC1 的激活至关重要，其中证实REDD1 通路部分参与铁螯合诱导的 mTORC1 抑制[35]；铁死亡作为是近年发现的一种新的程序性细胞死亡类型，这种新型死亡方式与活性氧(ROS)及细胞内铁相关，铁螯合剂可以抑制这一过程[36]，现在越来越多的研究发现铁死亡与OA的发生、发展相关，并且认为铁代谢异常是OA发展的影响因素之一[37]，综上,DDIT4可能成为临床在铁死亡方向作为OA治疗及诊断的一个新标志物。

综上所述，在本研究中使用综合生物信息学分析和机器学习方法筛选出BCL6，SCRG1，DDIT4鉴定为OA相关的特征基因，它们可能在OA的发生发展中充当重要的调节因子；并且在2个独立验证集ROC曲线中AUC值来验证3个相关基因，均获得较好的诊断价值，其中BCL6，SCRG1较DDIT4具有更高的诊断价值，DDIT4作为铁死亡相关基因，本研究通过生物信息学分析也证明铁死亡相关基因在OA发生发展中的作用，为OA的发病机制和治疗提供新的见解，但本研究是基于生物信息学及机器学习方法进行分析与解读，还需进一步的实验来进一步验证。