骨性关节炎滑膜核心基因筛选和生物信息学分析

张美龄 徐亚东 江小华

骨性关节炎(osteoarthritis, OA)是一种累及全关节的致残性疾病，随着全球人口老龄化和肥胖增加等因素的影响，发生率还在逐年上升,给个人健康及社会经济造成很大负担[1]。OA主要病理特征包括关节软骨退化、滑膜炎及关节结构的改变[2]。目前OA的发病机制还不是完全清晰，关节置换是疾病终末期的有效治疗方法，然而预后依然较差,因此疾病预防和早期OA诊断更为重要。此外，滑膜的病理改变可能发生在软骨变化之前，滑膜组织在基因层面的变化，可能是疾病发生、发展的关键所在，随着高通量微阵列技术的发展为深化基因-疾病关系的研究带来了希望[3，4]。本研究通过对公共数据库GEO来源的OA滑膜组织微阵列数据分析，包括筛选差异表达基因、功能富集分析、ROC曲线验证，以期探索OA潜在诊断和治疗靶点，同时为疾病分子机制研究提供生物信息学依据。

材料与方法

1.数据来源：从国际公开数据库GEO (https:∥www.ncbi.nlm.nih.gov/geo/)下载骨性关节炎相关数据集GSE55235、GSE55457、GSE55584，物种选择智人，3个微阵列数据集都基于GPL96平台[5，6]。GSE55235包含10个OA滑膜组织样本和10个正常滑膜组织样本，GSE55457包含10个OA滑膜组织样本和10个正常滑膜组织样本，GSE55584包含6个OA滑膜组织样本。

2.差异表达基因筛选：整合GSE55235、GSE55457、GSE55584微阵列数据集共26个OA滑膜组织样本，20个正常滑膜组织样本，并使用在线分析工具Network Analyst (https:∥www.networkanalyst.ca/)进行标准化处理，并进行主成分分析，再以P<0.05，|Log2(fold change, FC)|≥1为条件鉴定差异表达基因[7]。

3.差异表达基因的GO和KEGG通路分析：基因本体论(gene ontology, GO)术语包括生物过程(biological process, BP)、细胞组分(cellular components, CC)、分子功能(molecular function, MF)。使用R软件(3.6.3版本)Org.Hs.eg.db包(ID转换)、Ggplot2包(用于可视化)、Cluster-Profiler包(数据分析)进行GO和KEGG信号通路富集分析，调整后P<0.1为筛选标准。

4.蛋白-蛋白相互作用网络构建及核心蛋白筛选：使用在线分析工具STRING (https:∥string-db.org/)构建基因编码蛋白质的相互作用(protein-protein interaction, PPI)网络(可信度=0.4)并下载其TSV格式文件，标注基因上下调关系，利用细胞全景观软件Cyto-scape(3.8.3版本)展示DEGs的调控网络，并使用软件插件“Cytohubba”筛选最大相关性算法(maximal clique centrality，MCC)中位列前5个基因作为核心基因[8,9]。

5.核心基因的表达差异：使用在线分析工具Network Analyst进行核心基因的差异表达分析，以箱式图形式展示。

6.核心基因诊断效能评价:使用R软件(3.6.3版本)的PROC包(用于分析)、Ggplot2包(用于可视化)绘制5个核心基因的受试者工作特征(receiver operating characteristic curve, ROC)曲线。

结 果

1.数据预处理和差异表达基因筛选：整合3个芯片数据集共得到12545个基因的mRNA表达矩阵，并进行主成分分析显示,OA组和对照组间样本自然分开(图1)。根据P<0.05，|Log2FC|≥1标准共鉴定出327个DEGs，包括219个上调基因和108个下调基因(图2)。

图1 微阵列数据集的主成分分析括号内百分数表示OA组与对照组间差异中可以解释全面分析结果的百分数，样本与样本之间的距离越小越相似，反之差异越大

图2 差异表达基因火山图每个点代表1个基因，红色219个上调基因，蓝色108个下调基因，灰色表示表达量差异无统计学意义的基因

2.差异表达基因的功能富集分析：通过Cluster-Profiler包富集分析，219个上调DEGs主要参与含胶原的细胞外基质、质膜外侧、内质网腔的组成；参与细胞外基质组织、白细胞迁移等生物过程；发挥信号蛋白受体激活、糖胺聚糖结合、胞外基质结构组成等分子功能，参与破骨细胞分化、轴突引导等KEGG通路(图3A)。108个下调DEGs主要参与转录因子复合物的组成；参与糖皮质激素反应、皮质类固醇反应、脂多糖反应等生物过程；发挥DNA结合转录激活子激活、细胞因子激活、趋化因子激活等分子功能;富集于白细胞介素-17(interleukin-17, IL-17)、肿瘤坏死因子(tumor necrosis factor, TNF)、卡波西肉瘤相关疱疹病毒感染等KEGG通路(图3B)。

图3 差异表达基因GO和KEGG富集柱状图A.上调基因功能富集分析结果；B.下调基因功能富集分析结果

3.PPI网络构建及关键基因鉴定：基于在线分析工具STRING构建(303个节点、873条边组成)蛋白质调控网络，由Cyto-scape软件进行PPI可视化展示(图4A),利用Cytohubba插件MCC算法鉴定前5个核心基因分别是白介素-6(interleukin-6, IL-6)、C-X-C基序趋化因子配体8(c-x-c motif chemokine ligand8, CXCL8)、血管内皮生长因子A(vascular endothelial growth factor A, VEGFA)、toll样受体3(toll-like receptor, TLR3)、toll样受体7(toll-like receptor, TLR7,图4B)。

图4 蛋白质-蛋白质相互作用网络A.红色表示上调的基因编码蛋白质，绿色表示下调的基因编码蛋白质；B.Cytohubba筛选5个核心基因编码蛋白质

4.核心基因的表达：使用Network Analyst数据库分析发现，与对照组比较，TLR3和TLR7 滑膜mRNA水平在OA组表达升高(图5A、B)；CXCL8、IL-6、VEGFA 滑膜mRNA水平则在OA组表达下降(图5C～E)，以P<0.05为差异有统计学意义。

图5 核心基因表达水平箱式图A.TLR7；B.TLR3；C.CXCL8；D.IL-6；E.VEGFA

5.诊断效能评估：基于R语言绘制核心基因的ROC曲线其中IL-6(AUC=0.810)、CXCL8(AUC=0.846)、VEGFA(AUC=0.921)、TLR3基因(AUC=0.785)、TLR7基因(AUC=0.875，图6)，AUC均>0.7具有较好的OA诊断效能。

图6 核心基因诊断OA的ROC曲线A.表达上调基因TLR7、TLR3；B.表达下调基因CXCL8、IL-6、VEGFA

讨 论

骨性关节炎是一种慢性退行性关节疾病，由机体体质和机械等多因素复杂作用所致，常需要体格检查、实验结果、影像学来支持OA诊断[10，11]。基质金属蛋白酶(matrix metalloproteinases, MMPs)、C反应蛋白(C-reactive protein, CRP)、ⅡA型胶原氨基末端前肽(procollagen ⅡA amino-terminal pro-peptide, PⅡANP)等为目前报道有限数量的生物学标志物[12]。尚需探索更多高敏感度和特异性的诊断靶点来协助OA预防和治疗，而且有证据表明滑膜炎发生在影像学变化之前，可能是OA发病和结构进展的独立驱动因素[3]。

本研究首先对3个GEO来源的骨性关节炎滑膜微阵列数据集进行主成分分析，发现骨性关节炎组和正常对照组间样本自然分开，提示样本来源可靠，接着构建PPI网络筛选得到327个差异表达基因，功能聚类分析发现DEGs主要集中于含胶原的细胞外基质、白细胞迁移、细胞因子激活、趋化因子激活等基因本体术语；主要参与破骨细胞分化，轴突引导、IL-17、TNF等KEGG信号通络。有实验研究表明，在骨性关节炎初期，软骨失去完整性，胶原外基质降解，炎性介质产生，并作用于临近滑膜释放IL-7、TNF-α等细胞因子[13,14]。这些结果与本研究发现具有较高吻合度，提示使用生物信息学方法可以从海量数据信息中，挖掘并寻找规律，为疾病研究提供一种较为高效的分析手段[15]。接着筛选关键调控基因以及ROC曲线验证发现5个核心基因，其中TLR7、TLR3表达量上调，VEGFA、IL-6、CXCL8表达量下调，且AUC均>0.7(0.765～0.921)被认为具有较高诊断效能。

VEGFA为血管内皮生长因子A,属于VEGF家族中的一员，编码肝素结合蛋白，对生理性和病理性血管生成均至关重要，VEGF以多种亚型形式存在，依赖于mRNA选择性剪切。有研究显示，VEGFA参与了OA病理过程，在有症状OA滑膜中基因和蛋白水平表达降低，具有促进血管生成和疼痛的作用[16,17]。IL-6作为广泛表达于非特异性炎性介质，在OA免疫、炎症、软骨破坏、软骨重建中发挥重要作用[18]。有研究表明，IL-6在OA中似乎起着双向调控作用，既能促进白细胞介素-1(interleukin-1, IL-1)、TNF受体、MMP阻滞剂下调炎性反应，又能提高免疫细胞功能和炎性反应，IL-6差异性表达至少说明OA存在慢性改变状态[19]。CXCL8为C-X-C基序趋化因子配体8，该基因编码的蛋白质属于CXC趋化家族的一员，是炎性反应的主要介质，一项Meta分析显示，CXCL8作为滑膜液生物学标志物具有较高的诊断价值[20]。本研究发现CXCL8 mRNA水平在OA组中表达下降，提示CXCL8表达量可能与OA严重程度相关，在疾病中可能潜在具有炎症和稳态双重功能，其具体机制有待于进一步研究。

TLR3、TLR7编码的蛋白是Toll样受体家族成员，TLRs在病原体识别和固有免疫激活中发挥基本作用，其中TLR3可直接绑定β干扰素TIR结构域衔接蛋白激活干扰素调节因子(interferon regulatory factor，IRFs)，进而产生炎性细胞因子、趋化因子和蛋白酶以促进OA病程进展。此外有研究表明OA与2个TLR3单核苷酸多态性(single nucleotide polymorphism, SNP)位点显著相关，并且与TLR7的SNP位点也具有一定相关性，这些与OA风险遗传相关的证据表明TLRs在OA发生、发展中具有重要作用。

综上所述，本研究通过整合(GSE55235、GSE55584、GSE554573)微阵列数据集进行全基因组学分析，发现TLR7、TLR3、CXCL8、IL-6、VEGFA可能是骨性关节炎滑膜损害的关键基因，并初步探讨了基因在疾病中的表达规律及诊断效能，这些都有助于阐明OA发生的分子机制，进一步基因功能聚类分析提示破骨细胞分化、细胞因子激活、TNF信号通路等可能参与了OA的发展过程。