引用本文： 王海明, 张驰, 魏向阳, 魏全, 何成奇. 膝骨关节炎软骨细胞中差异基因表达的生物信息学分析. 华西医学, 2021, 36(5): 623-631. doi: 10.7507/1002-0179.202103320 复制
骨关节炎（osteoarthritis，OA）的发病过程主要为生物力学和生物学等综合因素共同作用下，软骨细胞、细胞外基质和软骨下骨三者合成-降解失衡所致[1-2]。其发病机制不明，现有的研究主要集中于软骨细胞，以及细胞外基质和软骨下骨的变化。OA 康复治疗目标以缓解疼痛、改善日常活动功能为主，其中物理治疗起着举足轻重的作用[1, 3]。随着细胞生物学、分子生物学等相关学科的交叉渗透，涉及 OA 疾病信号通路中生物标志物的研究成为关注的重点、热点[4-5]，现有研究提示 OA 病变机制复杂，单一信号通路可能无法明确阐述其发病机制[6]，因此，急需对 OA 复杂的信号转导网络关系进行筛选及深入探讨，这对于解释疾病背后的发病机制、物理治疗有效性的机制，以及未来疾病预防、诊治方面都重要的意义。生物信息学是生物学与信息学的交叉科学，其研究对象主要集中在基因和蛋白质 2 个方面，在生命科学的研究中发挥着至关重要的作用[7-8]。基因芯片技术用于研究基因表达谱与生物学功能之间可能存在的联系[9-11]，它的出现使我们可一次性对上万个基因的表达谱进行检测，极大地推动了生物信息学技术的进步[12]。本研究利用生物信息学方法对来自基因表达综合数据库（Gene Expression Omnibus，GEO）中 OA 患者软骨细胞基因芯片数据进行差异表达分析，随后进行京都基因和基因组数据库（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析[9-10]，采用文本挖掘分析其分子作用关系，构建疾病相关差异基因的分子调控网络，旨在进一步探究 OA 发生的分子机制。
1 资料与方法
1.1 软骨细胞样本数据的获取
在 GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）在线检索人类基因芯片样本，以“osteoarthritis OR cartilage OR chondrocyte*”为检索策略，满足以下纳入标准：① OA 疾病诊断符合美国风湿病学会诊断标准[13]；② 样本包含 OA 软骨细胞和对照组正常软骨细胞的检测数据；③ 软骨细胞样本与对照组样本数均需≥3 个，标本具有可重复性。获得 GPL570 平台上由 Dehne 等[14]提供的芯片数据系列 GSE16464，所用的实验平台为美国昂飞公司（Affymetrix, Inc）的 Affymetrix Human Genome U133 Plus 2.0 Array 人类全基因芯片数据。共纳入样本 6 个，其中 OA 软骨细胞组样本 3 个，对照组为正常膝关节软骨细胞样本 3 个（表1）。
表1
基因芯片数据 GSE16464 样本信息
1.2 基因分析
1.2.1 数据处理
对基因数据平台中存储的原始数据样本进行数据预处理，减少原始数据误差，增强进一步数据挖掘分析的信度[12, 14-15]。分析过程中我们借助基因本体论（Gene Ontology，GO）、KEGG 通路分析和基因/蛋白质相互作用关系检索工具（Search Tool for the Retrival of Interacting Genes/Proteins，STRING）等信息数据库和 R 语言、Perl 语言、Cytoscape 分析软件及 DAVID（Database for Annotation, Visualization and Integrated Discovery）（https://david.ncifcrf.gov/home.jsp）等分析工具。通过已知的信号传导通路及生化代谢反应通路，与在实验中得到的具体数据结合对其网络进行分析。
1.2.2 基因数据统计方法
使用 Perl 5.22.4 语言编辑软件将探针 ID 数据进行注释，并转换为基因名称（gene symbol）。使用 R 3.4.3 语言编辑软件进行芯片数据预处理和分析，使样本之间归一化具有可比性，基因表达原始数据进行标准化。通过 Affy 包中的 RMA（robust multi-arry avery）背景矫正和归一化处理后，对数据进行汇总从而获取矫正后表达水平的标准化数据[9]。通过 R 语言中线性回归模型软件包 limma 包对不同组的芯片进行差异计算，并用贝叶斯方法进行多重检验校正，通过倍比法（fold change，FC）和P值筛选获得差异基因[16]。差异基因的获得需同时满足以下条件：① |log2FC|>2；② P<0.05。然后，通过 GO、KEGG、蛋白相互作用网络分析采用超几何算法和 Benjamini 法对数据进行矫正分析。
1.2.3 基因 GO 数据库功能富集分析
GO 是一组预先定义好的、用来描述基因及其产物功能和行为标准术语，通过分析蛋白质术语之间的语义关系可以估计蛋白质之间的功能相似性。GO 数据库作为对基因及其蛋白质产物的功能进行系统描述的数据库，已经被广泛应用于分析基因（及其产物）间的功能相似性、基于高通量生物学数据分析疾病相关的生物学功能通路上，是目前最为成功的对生物学进行系统描述的工具[17]。我们通过 GO 分类号和 GO 数据相关分析工具将分类与具体基因联系起来，从而对该基因的功能分别在生物学过程、分子功能和细胞成分 3 个细胞生物学领域对基因及其产物的功能进行定义。
本研究所得到的差异基因通过 DAVID 数据库进行基因功能分化，应用 EASE（expression analysis systemic explore）方法选取 EASE<0.1 注释基因条目[18]。
1.2.4 蛋白质相互作用
我们利用 STRING 数据库[19]（https://string-db.org/）在线检索、预测蛋白质之间直接的物理相互作用和间接功能的相关性。将筛选出的差异基因输入到 STRING 10.5 数据库中，选取交互作用最小评分大于 0.4（中等置信度）的相互作用关系构建 OA 软骨细胞和正常组相关差异基因的蛋白相互作用网络。
最后，将 STRING 数据库中得到的蛋白质相互作用结果导入 Cytoscape 软件[20-21]中，进行网络分析及可视化操作，建构可视化的分子交互作用网络，并且对大规模蛋白质和蛋白质之间交互作用、蛋白质和 DNA 之间等交互作用的关联性进行分析。利用软件中 cytoHubba 插件同时计算各蛋白之间相互关联紧密程度的等级（degree）进一步筛选出 OA 软骨细胞的作用关键基因（hub gene）。
1.2.5 KEGG 通路富集分析
将本研究所得到的差异基因进行 KEGG 通路富集分析，通过对细胞内已知生物学过程的计算机化和将现有的基因功能信息解释标准化，对基因的功能进行注释和分析[22]，筛选出 OA 软骨细胞代谢的相关通路。
2 结果
2.1 数据标准化
基因原始数据表达值的中位数值呈现不均一状态（图 1a）。
图1
GSE16464 样本数据箱线图
a. 原始数据；b. 标准化后数据
采用 R 语言软件中的 Afyy 包 RMA 法对基因表达原始数据背景矫正和归一化处理后，对数据进行汇总从而获取矫正后表达水平的数据见图 1b。
2.2 差异表达基因
选取|log2FC|>2 及P<0.05 作为阈值，从而得到疾病组的差异表达基因 422 个。其中下调基因 191 个，上调基因 231 个。图 2 为差异表达基因的热图，图 3 为差异表达基因的火山图，表2 为排名前 30 差异表达基因列表。
图2
差异表达基因的热图
横轴为样品名称，纵轴为基因名称，其中亮蓝色 type 为正常组软骨细胞，红色 type 为 OA 软骨细胞，热图谱中蓝色代表下调基因，红色代表上调基因
图3
差异表达基因的火山图
红色圆点代表上调基因，绿色圆点代表下调基因
表2
排名前 30 的差异表达基因列表
2.3 差异表达基因 GO 注释
使用 DAVID 在线工具对正常软骨细胞与 OA 软骨细胞差异基因进行 GO 富集分析（包括生物学过程、细胞组分、分子功能 3 个方面）。生物学过程差异基因主要富集在细胞分裂、有丝分裂、姐妹染色单体内聚力、细胞增殖、RNA 聚合酶Ⅱ启动子转录调控、炎症反应、染色体分离、增强调控细胞外信号调节激酶（extracellular regulated protein kinases，ERK）1 和 ERK2 的级联反应等方面（图 4、表3）。细胞组分差异基因主要富集在细胞核、细胞外区、染色体着丝粒浓缩等方面（图 5）。分子功能差异基因主要富集在蛋白激酶绑定、丝氨酸/苏氨酸蛋白激酶活化、微管结合、蛋白激酶激活、激酶激活以及细胞因子活性等方面（图 6、表4）。
图4
差异表达基因 GO 富集分析-生物学过程
表3
GO 富集分析-生物学过程中的相关基因
图5
差异表达基因 GO 富集分析-细胞组分
图6
差异表达基因 GO 富集分析-分子功能
表4
GO 富集分析-分子功能中的相关基因
2.4 差异表达基因蛋白质相互作用分析
应用 STRING 数据库在线检索、预测蛋白质之间直接的物理相互作用和间接功能的相关性得到 267 个节点、972 条边（图 7）。使用 Cytoscape 的 cytoHbba 软件选取 degree 前 30 的节点（图 8）。其中CDK1、CCNB1、CDC20、NDC80、CDCA8、MAD2L1、TTK、NCAPG、PBK、KIF11、DLGAP5等基因是蛋白相互作用的重要联系节点。
图7
差异基因相关蛋白相互作用网络
圆点即节点为差异表达基因相关蛋白，边表示蛋白相互作用关系
图8
Cytoscape 的 cytoHbba 软件构建蛋白相互作用网络
方框即节点，其中蓝色代表前 30 位以外的差异表达基因相关蛋白，蓝色以外色彩代表前 30 位差异表达基因，从黄色到红色代表节点关系越来越密切，边表示蛋白相互作用关系
2.5 差异表达基因 KEGG 分析
使用 DAVID 在线工具对差异基因进行 KEGG 富集分析（图9）。卵母细胞减数分裂、细胞周期、丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信号通路、黄体酮调节卵母细胞成熟、生理周期、安非他命成瘾为最主要富集的信号通路。表5 为上述主要富集相关通路中的基因。
图9
差异基因进行 KEGG 富集分析结果
表5
主要富集通路中的基因
3 讨论
生物信息学从广义上讲是通过利用数理和信息科学的观点、原理和方法去研究生命现象的一门学科，在生命科学的研究中发挥着至关重要的作用[7, 23]。生物信息学通过以计算机为主要工具，开发各种生物信息学软件，对基因和蛋白质的序列、结构等相关信息进行分析和研究，建立相应的理论模型，以了解这些生物大分子所隐藏信息的生物学意义。如何充分利用这些数据，并通过对海量数据的分析、处理，揭示这些数据的内涵，得到对人类有价值的信息，将会成为有意义的研究方向。
基因芯片表达谱数据分析将获得的信息数据和生物学过程联系起来，解释基因的功能，从分子层面上了解疾病的发生发展。基于信息学的研究对于基因组信息进行挖掘，极大减轻了大规模筛选 OA 疾病后涉及的差异表达基因、蛋白、信号通路相关信息的工作量，对全面、系统、综合了解疾病的发生发展提供了新的方向。生物信息学的发展，为我们利用基因芯片技术和生物信息学分析方法，从基因水平或蛋白层面全面了解 OA 疾病的变化过程提供了可能。
目前尚缺乏对 OA 软骨细胞生物信息学分析的相关研究。我们对 OA 软骨细胞与正常组比较差异表达基因进行 GO 富集分析，结果显示这些基因大多数位于细胞核、细胞外区域，涉及染色体着丝粒浓缩。而就这些差异基因参与生物学过程而言，它们通过蛋白激酶活性的调节等方面参与细胞分裂、有丝分裂、细胞增殖、增强调控 ERK1 和 ERK2 的级联反应和炎症反应等过程，概括来说生物学过程主要涉及到遗传信息传递、蛋白代谢和炎症反应。进一步分析发现，膝 OA 软骨细胞中这些差异基因的分子功能主要富集在蛋白代谢过程和细胞因子活性等方面。在蛋白相互作用分析中的关键节点基因（CDK1、CCNB1、CDC20、NDC80、CDCA8、MAD2L1、TTK、NCAPG、PBK、KIF11、DLGAP5）也主要参与细胞周期、蛋白代谢。KEGG 富集分析显示，差异基因主要富集于卵母细胞减数分裂、细胞周期等细胞增殖相关通路和 MAPK 信号通路。上调、下调的差异基因生物信息学分析主要涉及细胞蛋白代谢途径，这与 OA 发生中合成、分解代谢耦联失衡的机制相一致[1-2]。在 KEEG 分析中 MAPK 通路[24]是主要富集通路之一，该通路是多种软骨细胞外刺激信号（包括应激、炎症性细胞因子、生长因子等）从细胞膜转导到细胞核内，与其他信号通路协同作用，使细胞对外界环境刺激作出实时的反应，参与软骨细胞生长、增殖、分化及凋亡等生理、病理的过程[24-25]。这一通路的分析与 GO 分析中炎症反应、细胞因子活性是 OA 软骨细胞差异基因富集的重要方面一致。许展仪等[26]对 OA 软骨与正常患者软骨进行生物信息学分析，发现他们差异基因也主要与慢性炎症反应、先天免疫反应有关。Liu 等[27]在分析 OA 患者与正常人群软骨中的长链非编码 RNA 生物富集分析、华芳等[28]对 OA 患者滑膜样本数据富集分析也发现细胞因子、趋化因子是其中重要的通路。在 OA 发生中，涉及软骨细胞的因子参与炎症反应和免疫反应过程，是当前关于发病机制的研究热点之一[6, 29-30]。细胞因子是细胞经刺激而合成、分泌的一类具有广泛生物学活性的小分子蛋白质，通过结合相应受体调节细胞生长、分化和效应，调控免疫应答。根据其在炎症反应中的作用又可分为促炎细胞因子和抗炎细胞因子。目前，OA 疾病被越来越多的学者认为是低度的慢性炎症性疾病[31-32]，发生机制中关于炎症因子研究逐步成为热点，其中较多研究的促炎细胞因子有肿瘤坏死因子 α、白细胞介素（interleukin，IL）-1、IL-6、IL-34 等[33-35]，而抗炎因子主要有 IL-10、IL-13、IL-4、转化生长因子等[6, 36-37]。这与我们在生物信息学分析中发现的炎症反应、细胞因子参与 OA 过程具有一致性，但现在 OA 研究较多的炎症因子、细胞因子与我们生物信息学分析中的差异基因并不完全一致，这些差异基因是否就是 OA 发生、发展的关键环节，这为我们后续进一步深入研究 OA 发生机制提供了理论上的新着眼点。
本研究存在一定的不足。本研究仅从 OA 软骨细胞层面进行生物信息学分析，可能与临床实际情况存在出入，还有待进一步的验证；研究中软骨细胞的样本虽然均一性较好，但是样本量不大（6 个样本），在一定程度上限制了结论的稳定性。
综上所述，多种途径参与了 OA 软骨细胞变化的过程，主要涉及到细胞周期、蛋白代谢基因/通路，炎症因子和细胞因子可能是 OA 发病中的重要环节。
骨关节炎（osteoarthritis，OA）的发病过程主要为生物力学和生物学等综合因素共同作用下，软骨细胞、细胞外基质和软骨下骨三者合成-降解失衡所致[1-2]。其发病机制不明，现有的研究主要集中于软骨细胞，以及细胞外基质和软骨下骨的变化。OA 康复治疗目标以缓解疼痛、改善日常活动功能为主，其中物理治疗起着举足轻重的作用[1, 3]。随着细胞生物学、分子生物学等相关学科的交叉渗透，涉及 OA 疾病信号通路中生物标志物的研究成为关注的重点、热点[4-5]，现有研究提示 OA 病变机制复杂，单一信号通路可能无法明确阐述其发病机制[6]，因此，急需对 OA 复杂的信号转导网络关系进行筛选及深入探讨，这对于解释疾病背后的发病机制、物理治疗有效性的机制，以及未来疾病预防、诊治方面都重要的意义。生物信息学是生物学与信息学的交叉科学，其研究对象主要集中在基因和蛋白质 2 个方面，在生命科学的研究中发挥着至关重要的作用[7-8]。基因芯片技术用于研究基因表达谱与生物学功能之间可能存在的联系[9-11]，它的出现使我们可一次性对上万个基因的表达谱进行检测，极大地推动了生物信息学技术的进步[12]。本研究利用生物信息学方法对来自基因表达综合数据库（Gene Expression Omnibus，GEO）中 OA 患者软骨细胞基因芯片数据进行差异表达分析，随后进行京都基因和基因组数据库（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析[9-10]，采用文本挖掘分析其分子作用关系，构建疾病相关差异基因的分子调控网络，旨在进一步探究 OA 发生的分子机制。
1 资料与方法
1.1 软骨细胞样本数据的获取
在 GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）在线检索人类基因芯片样本，以“osteoarthritis OR cartilage OR chondrocyte*”为检索策略，满足以下纳入标准：① OA 疾病诊断符合美国风湿病学会诊断标准[13]；② 样本包含 OA 软骨细胞和对照组正常软骨细胞的检测数据；③ 软骨细胞样本与对照组样本数均需≥3 个，标本具有可重复性。获得 GPL570 平台上由 Dehne 等[14]提供的芯片数据系列 GSE16464，所用的实验平台为美国昂飞公司（Affymetrix, Inc）的 Affymetrix Human Genome U133 Plus 2.0 Array 人类全基因芯片数据。共纳入样本 6 个，其中 OA 软骨细胞组样本 3 个，对照组为正常膝关节软骨细胞样本 3 个（表1）。
表1
基因芯片数据 GSE16464 样本信息
1.2 基因分析
1.2.1 数据处理
对基因数据平台中存储的原始数据样本进行数据预处理，减少原始数据误差，增强进一步数据挖掘分析的信度[12, 14-15]。分析过程中我们借助基因本体论（Gene Ontology，GO）、KEGG 通路分析和基因/蛋白质相互作用关系检索工具（Search Tool for the Retrival of Interacting Genes/Proteins，STRING）等信息数据库和 R 语言、Perl 语言、Cytoscape 分析软件及 DAVID（Database for Annotation, Visualization and Integrated Discovery）（https://david.ncifcrf.gov/home.jsp）等分析工具。通过已知的信号传导通路及生化代谢反应通路，与在实验中得到的具体数据结合对其网络进行分析。
1.2.2 基因数据统计方法
使用 Perl 5.22.4 语言编辑软件将探针 ID 数据进行注释，并转换为基因名称（gene symbol）。使用 R 3.4.3 语言编辑软件进行芯片数据预处理和分析，使样本之间归一化具有可比性，基因表达原始数据进行标准化。通过 Affy 包中的 RMA（robust multi-arry avery）背景矫正和归一化处理后，对数据进行汇总从而获取矫正后表达水平的标准化数据[9]。通过 R 语言中线性回归模型软件包 limma 包对不同组的芯片进行差异计算，并用贝叶斯方法进行多重检验校正，通过倍比法（fold change，FC）和P值筛选获得差异基因[16]。差异基因的获得需同时满足以下条件：① |log2FC|>2；② P<0.05。然后，通过 GO、KEGG、蛋白相互作用网络分析采用超几何算法和 Benjamini 法对数据进行矫正分析。
1.2.3 基因 GO 数据库功能富集分析
GO 是一组预先定义好的、用来描述基因及其产物功能和行为标准术语，通过分析蛋白质术语之间的语义关系可以估计蛋白质之间的功能相似性。GO 数据库作为对基因及其蛋白质产物的功能进行系统描述的数据库，已经被广泛应用于分析基因（及其产物）间的功能相似性、基于高通量生物学数据分析疾病相关的生物学功能通路上，是目前最为成功的对生物学进行系统描述的工具[17]。我们通过 GO 分类号和 GO 数据相关分析工具将分类与具体基因联系起来，从而对该基因的功能分别在生物学过程、分子功能和细胞成分 3 个细胞生物学领域对基因及其产物的功能进行定义。
本研究所得到的差异基因通过 DAVID 数据库进行基因功能分化，应用 EASE（expression analysis systemic explore）方法选取 EASE<0.1 注释基因条目[18]。
1.2.4 蛋白质相互作用
我们利用 STRING 数据库[19]（https://string-db.org/）在线检索、预测蛋白质之间直接的物理相互作用和间接功能的相关性。将筛选出的差异基因输入到 STRING 10.5 数据库中，选取交互作用最小评分大于 0.4（中等置信度）的相互作用关系构建 OA 软骨细胞和正常组相关差异基因的蛋白相互作用网络。
最后，将 STRING 数据库中得到的蛋白质相互作用结果导入 Cytoscape 软件[20-21]中，进行网络分析及可视化操作，建构可视化的分子交互作用网络，并且对大规模蛋白质和蛋白质之间交互作用、蛋白质和 DNA 之间等交互作用的关联性进行分析。利用软件中 cytoHubba 插件同时计算各蛋白之间相互关联紧密程度的等级（degree）进一步筛选出 OA 软骨细胞的作用关键基因（hub gene）。
1.2.5 KEGG 通路富集分析
将本研究所得到的差异基因进行 KEGG 通路富集分析，通过对细胞内已知生物学过程的计算机化和将现有的基因功能信息解释标准化，对基因的功能进行注释和分析[22]，筛选出 OA 软骨细胞代谢的相关通路。
2 结果
2.1 数据标准化
基因原始数据表达值的中位数值呈现不均一状态（图 1a）。
图1
GSE16464 样本数据箱线图
a. 原始数据；b. 标准化后数据
采用 R 语言软件中的 Afyy 包 RMA 法对基因表达原始数据背景矫正和归一化处理后，对数据进行汇总从而获取矫正后表达水平的数据见图 1b。
2.2 差异表达基因
选取|log2FC|>2 及P<0.05 作为阈值，从而得到疾病组的差异表达基因 422 个。其中下调基因 191 个，上调基因 231 个。图 2 为差异表达基因的热图，图 3 为差异表达基因的火山图，表2 为排名前 30 差异表达基因列表。
图2
差异表达基因的热图
横轴为样品名称，纵轴为基因名称，其中亮蓝色 type 为正常组软骨细胞，红色 type 为 OA 软骨细胞，热图谱中蓝色代表下调基因，红色代表上调基因
图3
差异表达基因的火山图
红色圆点代表上调基因，绿色圆点代表下调基因
表2
排名前 30 的差异表达基因列表
2.3 差异表达基因 GO 注释
使用 DAVID 在线工具对正常软骨细胞与 OA 软骨细胞差异基因进行 GO 富集分析（包括生物学过程、细胞组分、分子功能 3 个方面）。生物学过程差异基因主要富集在细胞分裂、有丝分裂、姐妹染色单体内聚力、细胞增殖、RNA 聚合酶Ⅱ启动子转录调控、炎症反应、染色体分离、增强调控细胞外信号调节激酶（extracellular regulated protein kinases，ERK）1 和 ERK2 的级联反应等方面（图 4、表3）。细胞组分差异基因主要富集在细胞核、细胞外区、染色体着丝粒浓缩等方面（图 5）。分子功能差异基因主要富集在蛋白激酶绑定、丝氨酸/苏氨酸蛋白激酶活化、微管结合、蛋白激酶激活、激酶激活以及细胞因子活性等方面（图 6、表4）。
图4
差异表达基因 GO 富集分析-生物学过程
表3
GO 富集分析-生物学过程中的相关基因
图5
差异表达基因 GO 富集分析-细胞组分
图6
差异表达基因 GO 富集分析-分子功能
表4
GO 富集分析-分子功能中的相关基因
2.4 差异表达基因蛋白质相互作用分析
应用 STRING 数据库在线检索、预测蛋白质之间直接的物理相互作用和间接功能的相关性得到 267 个节点、972 条边（图 7）。使用 Cytoscape 的 cytoHbba 软件选取 degree 前 30 的节点（图 8）。其中CDK1、CCNB1、CDC20、NDC80、CDCA8、MAD2L1、TTK、NCAPG、PBK、KIF11、DLGAP5等基因是蛋白相互作用的重要联系节点。
图7
差异基因相关蛋白相互作用网络
圆点即节点为差异表达基因相关蛋白，边表示蛋白相互作用关系
图8
Cytoscape 的 cytoHbba 软件构建蛋白相互作用网络
方框即节点，其中蓝色代表前 30 位以外的差异表达基因相关蛋白，蓝色以外色彩代表前 30 位差异表达基因，从黄色到红色代表节点关系越来越密切，边表示蛋白相互作用关系
2.5 差异表达基因 KEGG 分析
使用 DAVID 在线工具对差异基因进行 KEGG 富集分析（图9）。卵母细胞减数分裂、细胞周期、丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信号通路、黄体酮调节卵母细胞成熟、生理周期、安非他命成瘾为最主要富集的信号通路。表5 为上述主要富集相关通路中的基因。
图9
差异基因进行 KEGG 富集分析结果
表5
主要富集通路中的基因
3 讨论
生物信息学从广义上讲是通过利用数理和信息科学的观点、原理和方法去研究生命现象的一门学科，在生命科学的研究中发挥着至关重要的作用[7, 23]。生物信息学通过以计算机为主要工具，开发各种生物信息学软件，对基因和蛋白质的序列、结构等相关信息进行分析和研究，建立相应的理论模型，以了解这些生物大分子所隐藏信息的生物学意义。如何充分利用这些数据，并通过对海量数据的分析、处理，揭示这些数据的内涵，得到对人类有价值的信息，将会成为有意义的研究方向。
基因芯片表达谱数据分析将获得的信息数据和生物学过程联系起来，解释基因的功能，从分子层面上了解疾病的发生发展。基于信息学的研究对于基因组信息进行挖掘，极大减轻了大规模筛选 OA 疾病后涉及的差异表达基因、蛋白、信号通路相关信息的工作量，对全面、系统、综合了解疾病的发生发展提供了新的方向。生物信息学的发展，为我们利用基因芯片技术和生物信息学分析方法，从基因水平或蛋白层面全面了解 OA 疾病的变化过程提供了可能。
目前尚缺乏对 OA 软骨细胞生物信息学分析的相关研究。我们对 OA 软骨细胞与正常组比较差异表达基因进行 GO 富集分析，结果显示这些基因大多数位于细胞核、细胞外区域，涉及染色体着丝粒浓缩。而就这些差异基因参与生物学过程而言，它们通过蛋白激酶活性的调节等方面参与细胞分裂、有丝分裂、细胞增殖、增强调控 ERK1 和 ERK2 的级联反应和炎症反应等过程，概括来说生物学过程主要涉及到遗传信息传递、蛋白代谢和炎症反应。进一步分析发现，膝 OA 软骨细胞中这些差异基因的分子功能主要富集在蛋白代谢过程和细胞因子活性等方面。在蛋白相互作用分析中的关键节点基因（CDK1、CCNB1、CDC20、NDC80、CDCA8、MAD2L1、TTK、NCAPG、PBK、KIF11、DLGAP5）也主要参与细胞周期、蛋白代谢。KEGG 富集分析显示，差异基因主要富集于卵母细胞减数分裂、细胞周期等细胞增殖相关通路和 MAPK 信号通路。上调、下调的差异基因生物信息学分析主要涉及细胞蛋白代谢途径，这与 OA 发生中合成、分解代谢耦联失衡的机制相一致[1-2]。在 KEEG 分析中 MAPK 通路[24]是主要富集通路之一，该通路是多种软骨细胞外刺激信号（包括应激、炎症性细胞因子、生长因子等）从细胞膜转导到细胞核内，与其他信号通路协同作用，使细胞对外界环境刺激作出实时的反应，参与软骨细胞生长、增殖、分化及凋亡等生理、病理的过程[24-25]。这一通路的分析与 GO 分析中炎症反应、细胞因子活性是 OA 软骨细胞差异基因富集的重要方面一致。许展仪等[26]对 OA 软骨与正常患者软骨进行生物信息学分析，发现他们差异基因也主要与慢性炎症反应、先天免疫反应有关。Liu 等[27]在分析 OA 患者与正常人群软骨中的长链非编码 RNA 生物富集分析、华芳等[28]对 OA 患者滑膜样本数据富集分析也发现细胞因子、趋化因子是其中重要的通路。在 OA 发生中，涉及软骨细胞的因子参与炎症反应和免疫反应过程，是当前关于发病机制的研究热点之一[6, 29-30]。细胞因子是细胞经刺激而合成、分泌的一类具有广泛生物学活性的小分子蛋白质，通过结合相应受体调节细胞生长、分化和效应，调控免疫应答。根据其在炎症反应中的作用又可分为促炎细胞因子和抗炎细胞因子。目前，OA 疾病被越来越多的学者认为是低度的慢性炎症性疾病[31-32]，发生机制中关于炎症因子研究逐步成为热点，其中较多研究的促炎细胞因子有肿瘤坏死因子 α、白细胞介素（interleukin，IL）-1、IL-6、IL-34 等[33-35]，而抗炎因子主要有 IL-10、IL-13、IL-4、转化生长因子等[6, 36-37]。这与我们在生物信息学分析中发现的炎症反应、细胞因子参与 OA 过程具有一致性，但现在 OA 研究较多的炎症因子、细胞因子与我们生物信息学分析中的差异基因并不完全一致，这些差异基因是否就是 OA 发生、发展的关键环节，这为我们后续进一步深入研究 OA 发生机制提供了理论上的新着眼点。
本研究存在一定的不足。本研究仅从 OA 软骨细胞层面进行生物信息学分析，可能与临床实际情况存在出入，还有待进一步的验证；研究中软骨细胞的样本虽然均一性较好，但是样本量不大（6 个样本），在一定程度上限制了结论的稳定性。
综上所述，多种途径参与了 OA 软骨细胞变化的过程，主要涉及到细胞周期、蛋白代谢基因/通路，炎症因子和细胞因子可能是 OA 发病中的重要环节。
