基于WGCNA筛选骨关节炎的生物标志物

龙通华，张红参，许冬梅，曾志华，傅子原，关皓邺

(1.右江民族医学院研究生学院，2.右江民族医学院附属医院康复医学科，3.广西高校重点实验室右江流域特色民族药研究重点实验室，广西百色 533000)

骨关节炎(osteoarthritis，OA)是最常见的慢性退行性骨骼系统疾病之一，多见于中老年人。骨关节炎的病理特征是软骨退化、滑膜发炎、软骨下改变及骨赘形成，最终导致关节功能丧失[1]，因此骨关节炎是导致老年人关节疼痛、身体残疾和生活质量下降的主要原因之一[2]。目前，骨关节炎的发病机制仍不清楚，缺少可用于诊断的生物标志物，因此，利用生物信息学方法探索骨关节炎的生物标志物及其潜在的作用机制对于骨关节炎的诊断、治疗和预后至关重要。

有研究发现，骨关节炎的发病机制与基因表达调控有关[3-4]。基因芯片技术、基因测序技术和生物信息学分析在基因组水平的数据分析中被广泛应用，可识别与骨关节炎发生发展相关的差异表达基因(differentially expressed genes，DEGs)和功能通路。加权基因共表达网络分析(weighted gene co-expression network analysis，WGCNA)是一种基于高通量基因表达谱的算法[5]。它被广泛用于识别各种疾病的基因共表达网络，以揭示基因之间的相关性，并找到疾病显著相关的基因模块[6-8]。基因表达综合数据库(gene expression omnibus，GEO, http: //www.ncbi.nlm.nih.gov/geo)是一个高通量微阵列和下一代测序序列功能基因组数据库[9]。基因本体论(gene ontology，GO)是注释基因和分析基因生物学过程的主要生物信息学工具[10]。京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)是一个整合了基因组、化学和系统功能信息的数据库，并对基因的功能进行系统化分析[11]。

本研究基于生物信息学分析方法，通过GEO数据库下载含有20个骨关节炎和18个正常对照样本的高通量测序数据集GSE114007进行差异分析，利用WGCNA将差异基因分为不同的基因模块，并筛选出与骨关节炎显著相关的基因模块进行GO和KEGG富集分析，寻找其潜在的作用机制。蛋白互作(protein-protein interaction，PPI)网络筛选出DDIT3和BCL6作为关键基因。最后从GEO数据库下载GSE57218和GSE169077表达谱数据对训练集GSE114007分析结果进行验证。

1.1 数据下载和预处理本研究从GEO数据库搜索并下载与骨关节炎相关的转录组测序数据和表达谱数据，得到GSE114007、GSE57218和GSE169077三个数据集，其中GSE114007作为训练集，GSE57218和GSE169077作为验证集。GSE114007是转录组测序的counts数据，有20个骨关节炎样本和18个正常对照样本，GSE57218有33个骨关节炎样本和7个正常对照样本， GSE169077有6个骨关节炎样本和5个正常对照样本。利用R语言软件对表达谱数据进行log2转换和归一化处理，根据平台注释信息将探针ID转换基因名，数据中多个基因探针ID对应单个基因名，将数据进行求和取最大值处理。

1.2 筛选差异基因使用R语言软件的“DESeq2”包[12-13]对GSE114007数据的骨关节炎样本与正常对照样本进行差异分析，筛选DEGs。DEGs筛选标准为矫正P(adj.P.Val)<0.05和|log2FC|>1。

1.3 构建DEGs共表达基因模块利用“WGCNA”包[5]对DEGs构建共表达网络，首先pickSoftThreshold函数计算出相关系数为0.86时对应的软阈值为10。以软阈值为10将差异表达基因构建成邻接矩阵，之后转变为拓扑重叠矩阵[14]，以识别基因模块。基因模块的最小基因数为30，合并相似基因模块的切割高度的阈值设为0.25[5]。

1.4 GO和KEGG富集分析为了进一步分析与OA最相关的基因模块的基因功能，我们使用“clusterProfiler”包[15]分别对与OA最相关的基因模块的基因进行GO富集分析和KEGG通路富集分析，筛选条件为P<0.05。

1.5 PPI分析和关键基因的筛选通过STRING在线数据11.5(https://string-db.org/cgi/input.pl)构建与骨关节炎最相关的基因模块基因的蛋白质-蛋白质相互作用网络关系，以识别蛋白质之间的相互作用，最低互动分数为0.9。并通过Cytoscape软件v3.8.2版本进行可视化PPI网络[16]，利用cytoHubba插件的Degree算法筛选出前十个基因作为关键基因[17]。

1.6 关键基因的验证在10个关键基因中有两个关键基因(DDIT3和BCL6)在骨关节炎中尚未见相关研究报道。利用t检验计算出关键基因在GSE57218和GSE169077表达谱数据中OA组和正常组的表达水平，验证关键基因(DDIT3和BCL6)表达水平的差异。

1.7 统计学方法根据数据分布特点，采用非参数检验或t检验分析两组间差异的统计学意义。所有分析均使用R软件4.0.5版本进行分析，adj.P.Val或P<0.05为有统计学意义。

2.1 差异表达基因筛选在骨关节炎和正常对照两组中共筛选出1752个DEGs，其中上调基因927个，下调基因825个，DEGs的结果展示如图1。

2.2 DEGs共表达网络分析我们选取相关系数为0.86时对应的软阈值为10(如图2A所示)构建无尺度网络进行WGCNA分析。通过WGCNA分析将DEGs构建了6个不同的基因模块，分别是棕色基因模块、蓝色基因模块、蓝绿色基因模块、绿色基因模块、黄色基因模块及灰色基因模块，如图2B所示。根据图3A所示的基因模块与骨关节炎相关性分析表明，这些基因模块都与骨关节炎的发生发展相关。图3B显示了每个基因模块与OA的相关显著性，其中棕色基因模块与OA的相关关系最为显著。图3C显示了棕色模块基因与OA相关的显著性，因此将棕色基因模块作为关键基因模块。

2.3 关键基因模块的功能富集分析GO富集分析结果表明，棕色基因模块的基因主要富集在脂类代谢调控、T细胞活化、细胞凋亡、炎症反应调控、组织缺氧反应、细胞饥饿反应等生物学过程(图4)。KEGG通路富集分析结果表明，棕色模块基因主要富集在MAPK信号通路、破骨细胞通路、肿瘤信号通路、脂肪细胞因子信号通路、mTOR信号通路、昼夜节律等信号通路(图5)。

2.4 关键基因模块的PPI网络分析采用STRING在线数据库以及Cytoscape软件构建棕色模块基因的PPI(图6A)，利用cytoHubba插件的degree算法筛选出PPI前10个基因(图6B)，分别是JUN原癌基因(jun proto-oncogene, JUN)、CCAAT增强子结合蛋白β(CCAAT enhancer binding protein beta，CEBPB)、Fos原癌基因(Fos proto-oncogene，FOS)、RELA原癌基因(RELA proto-oncogene，RELA)、泛素C(Ubiquitin C，UBC)、CCAAT增强子结合蛋白α(CCAAT enhancer binding protein alpha，CEBPA)、激活转录因子3(activating transcription factor 3，ATF3)、DNA 损伤诱导的转录因子3(DNA damage inducible transcript 3，DDIT3)、BCL6转录抑制因子(BCL6 transcription repressor，BCL6)和干扰素调节因子4(interferon regulatory factor 4，IRF4)。

2.5 关键基因的表达水平在cytoHubba插件的degree算法筛选出PPI前10个基因中，DDIT3和BCL6基因目前未见其与骨关节炎相关的报道研究，因此最终选取DDIT3和BCL6作为关键基因。图7显示了DDIT3和BCL6在骨关节炎与正常对照两组间的表达量水平情况，这两个关键基因在骨关节炎组样本中相对于正常对照组表达量显著下调。

2.6 关键基因的验证在验证集GSE57218中，DDIT3和BCL6这两个关键基因的表达量如图8A、B所示，DDIT3和BCL6的表达量在OA组中相对于正常对照组显著下调。在另一个验证集GSE169077中，DDIT3和BCL6的表达量如图8C、D所示，其表达量在OA组中相对于正常对照组也显著下调。

骨关节炎是一种退行性骨关节疾病，以软骨退化和破坏、关节滑膜炎、软骨下骨改变及骨赘形成为病理特征[1]，导致日常生活能力严重下降，给患者和医疗保健系统带来沉重的经济负担[2]。目前，骨关节炎的发病机制仍不清楚，了解其发病机制对骨关节炎的诊断、治疗和预后具有重要的意义[18-19]，本研究从GEO数据库下载GSE114007的转录组数据，利用“DESeq2”包进行差异分析，筛选出1752个DEGs，其中上调基因927个，下调基因825个。WGCNA分析方法将DEGs构建基因共表达网络，得到6个不同表达模式的基因模块，其中棕色基因模块(共281个基因)与OA的相关系数最高，且呈负相关，表明了棕色基因模块是参与骨关节炎发生发展的关键基因模块。GO富集分析表明，棕色模块基因主要富集在脂类代谢调控、T细胞活化、细胞凋亡、炎症反应和组织缺氧反应等生物调控过程，在这些生物学过程中，脂类代谢、细胞凋亡和炎症反应已经被证明与OA的发展有关[20-23]。KEGG通路富集分析表明，MAPK信号通路、破骨细胞通路、肿瘤信号通路、脂肪细胞因子信号通路和mTOR信号通路与OA的发病机制有关，而这些信号通路已经被证明在OA发病机制中起着重要的作用[24-25]。正常软骨组织中的代谢受严格的调节，以维持合成和降解之间的平衡，有研究认为骨关节炎发生时，这种平衡被打破，在骨关节炎的软骨退化后其代谢和基因表达谱都会发生变化。

PPI网络分析筛选出互作关系最强的10个基因JUN、CEBPB、FOS、RELA、UBC、CEBPA、ATF3、DDIT3、BCL6、IRF4。有研究发现JUN通过调节软骨细胞合成代谢和分解代谢基因的表达在OA软骨破坏中发挥重要作用[26]。LU等人[27]进一步研究发现JUN参与软骨细胞的凋亡，且与OA的发生有紧密关系。FOS是一种具有调控细胞的增殖、分化和转化生理功能的原癌基因[28]，其可上调MMP1、抑制TIMP1启动子的转录，从而降解软骨组织，是其致骨关节炎发生的机制之一[29]。ATF3与多种疾病的发病机制有关，包括代谢性、免疫性和炎症性疾病，研究发现ATF3通过核因子-kB(NF-kB)途径软骨细胞炎性细胞因子的表达参与OA的发病过程[30]。此外，CEBPB、RELA、UBC、CEBPA和IRF4这些基因也陆续被发现与OA的发病机制有关[31-33]。

DDIT3是一种内质网应激通路中的关键转录因子，研究发现DDIT3通过内质网应激通路调控细胞凋亡和炎症反应过程[34]，XIONG等人[35]进一步研究发现在软骨发育过程中DDIT3通过内质网应激通路参与细胞凋亡，以维持正常的代谢。YU等人[36]研究发现DDIT3通过诱导SOX9的表达来调控软骨细胞分化，YANG等人[37]进一步研究阐明了DDIT3通过SIRT1-AKT途径促进软骨细胞的自噬，这些研究表明DDIT3对软骨的代谢、分化和凋亡具有极其重要的作用。GO和KEGG富集分析显示DDIT3参与脂类代谢、细胞凋亡和MAPK信号通路，因此DDIT3有可能通过脂类代谢、细胞凋亡和MAPK信号通路参与OA的发生发展。此外，PPI网络分析表明，DDTI3与JUN、FOS、ATF3存在着相互作用网络关系，因此DDIT3与这些基因存在相互调节的作用，又因JUN、FOS、ATF3基因与OA的发病机制有关，所以DDIT3可能通过这些基因进而参与OA的发病机制。有研究发现FOS和DDIT3存在相互调控关系共参与软骨形成和分化[38]，LI等[39]研究发现DDIT3和JUN通过形成新的复合物参与细胞凋亡的调控。此外，DDIT3通过IL-1β和MMP3介导软骨细胞分解代谢和细胞凋亡反应[40]。

BCL6是一种转录抑制因子，在细胞增殖、炎症反应和凋亡中起到重要的作用。GO富集分析显示BCL6参与炎症反应、免疫反应、衰老和细胞凋亡等生物学过程。有研究发现BCL6对成骨细胞分化有重要的作用[41]，另外，有研究发现BCL6对破骨细胞分化也起到重要的作用[42]。许多研究发现BCL6在类风湿性关节炎发病机制中扮演着重要的角色[43]。PPI网络分析显示BCL6与JUN、FOS及IRF4存在蛋白互作关系，因此BCL6有可能通过JUN、FOS和IRF4之间相互作用参与OA的发病机制。CHAWEEWANNAKORN等[44]研究发现BCL6通过下调FOS影响破骨细胞调节因子的生成。这些研究表明BCL6与骨关节炎的发生有着密切关联。

本研究又在GSE57218和GSE169077两个数据集中进行验证，DDIT3和BCL6在两个验证集的骨关节炎组中的表达量也是显著低于正常对照组，进一步验证了本研究的结果，因此，我们可以到得DDIT3和BCL6与骨关节炎的发生发展有密切关联。虽然本研究结果在另外两个表达谱数据得到了验证，但我们研究还有不足之处。其一DDIT3和BCL6这两个关键基因没有进一步实验验证，其二是对DDIT3和BCL6关键基因参与OA发病的具体机制没有进一步地研究分析，只是通过GO、KEGG富集分析和现有文献研究结果进行关联探寻其参与OA的潜在发病机制。

综上所述，本研究结果表明DDIT3和BCL6与骨关节炎的发生发展有密切关联，DDIT3和BCL6可能通过MAPK信号通路、mTOR信号通路、脂类代谢参与骨关节炎的发生发展，因此DDIT3和BCL6可作为新的骨关节炎候选生物标志物。