基于基因芯片技术结合CIBERSORT反卷积算法研究骨关节炎的免疫机制及潜在中药靶向预测①

杨雷 伍搏宇 熊辉齐 新宇 简功辉 杨少锋 聂颖 马露

（湖南中医药大学，长沙 410208）

骨关节炎（osteoarthritis，OA）是一种以关节破坏、骨质增生和关节软骨退行性病变为特征的慢性疾病。目前机制尚不清楚，多种因素均可导致OA［1-3］。既往认为OA是一种“磨损”疾病，而目前却认为OA是一种免疫相关疾病，免疫复合物介导的补体激活、破骨细胞激活以及细胞因子平衡失调在OA病情进展中处于重要地位，备受关注［4-5］。早期OA患者影像学改变并不明显，而滑膜炎却是所有关节炎最早出现的共有特征［6］。滑膜组织作为OA免疫系统的损伤对象，是持续性关节退行性变的主要驱动因素，尤其免疫细胞浸润，包括浆细胞、M2型巨噬细胞、活化树突状细胞等多种免疫细胞滑膜浸润参与并介导了OA的自身免疫应答，诱导促炎细胞因子、趋化因子及蛋白酶分泌，扰乱免疫平衡进而加速软骨和骨侵蚀，在OA发生发展、恶化过程中起重要作用［7-10］。因此，从免疫浸润角度研究OA发病机制是针对OA早期治疗、防止恶化的先决条件，可为OA精准靶向免疫治疗提供新途径。

生物信息学是分子生物学和信息技术的新兴交叉学科，可揭示疾病分子机制。基因芯片作为一种新兴技术，已被用于生物信息的高效、大规模获取，可广泛收集疾病表达谱数据，尤其对多种免疫细胞间的动态调节研究具有指导作用。本文采用生物信息学工具对通用基因芯片数据库（gene expression omnibus，GEO）中OA患者和正常人群滑膜数据进行分析，通过CIBERSORT反卷积算法分析OA滑膜中的免疫浸润情况，并根据差异的免疫浸润基因预测相关免疫调节中药，为OA免疫机制研究和免疫调节相关中药筛选提供参考。

1 材料与方法

1.1 芯片数据采用GEO（www.ncbi.nlm.nih.gov/geo/）数据库搜索OA滑膜基因芯片数据，筛选GSE55235数据，由Thomas Häupl提供（包含10个正常人群滑膜组织、10个OA患者滑膜组织），设正常滑膜为对照组，OA患者滑膜数据为实验组，根据GPL96-57554（［HG-U133A］Affymetrix Human Genome U133A Array）平台将探针转换为基因名称。

1.2 差异表达基因获取通过R语言程序“limma”安装包对基因探针进行两样本t检验，并根据筛选条件|log2FC|≥1和校正P<0.05进行差异基因分析，其中log2FC<0代表基因表达下调，log2FC>0代表基因表达上调。

1.3 核心基因获取和注释将两组差异基因通过在线平台STRING（www.string-db.org/cgi/input.pl）获取蛋白互作（protein-protein interaction，PPI）网络模型，并将结果导入Cytoscape软件，采用Cytoscape插件工具cyto Hubba中的MCC方法识别重要基因，将重要基因进行可视化，并列出排名前10的作为核心基因。采用R语言cluster Profiler对差异基因进行基因本体论（genetic ontology，GO）分析与京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析，并探究其潜在生物学过程（biological process，BP）、细胞成分（cellular component，CC）、分子功能（molecular function，MF）及重要信号通路。

1.4 免疫细胞浸润矩阵获取采用R语言程序并链接CIBERSORT反卷积算法进行活化记忆性CD4 T细胞、浆细胞、活化NK细胞、调节性T细胞等22种免疫细胞转录特征的模拟计算［11］。为了结果准确，设定模拟次数为1 000次，采用Kruskal-Wallis秩和检验对P<0.05的数据进行后续分析。

1.5 免疫细胞相关性和差异分析为分析各类免疫细胞相关性，对CIBERSORT反卷积算法以P<0.05基础上筛选的可信样本数据中进行不同免疫细胞间Pearson相关系数计算，并采用秩和检验比较OA组与对照组差异。

1.6 OA免疫浸润过程潜在中药预测将两组差异的核心靶基因和GO富集与免疫相关BP导入在线平台（www.coremine.com/medical/），筛选具有潜在免疫调节作用的药物，采用秩和检验，以P<0.01表示差异具有统计学意义，并对中药进行功效、归经统计。

1.7 主要观察指标通过GEO数据库获取正常人群滑膜和OA患者滑膜组织差异基因及免疫细胞浸润分布，并进行潜在中药预测，主要观察指标：①正常滑膜组织和OA滑膜组织免疫细胞浸润的分布差异及相关性；②正常滑膜组织和OA滑膜组织差异的靶基因；③OA滑膜组织免疫浸润相关BP和可能信号通路；④挖掘调节OA滑膜组织局部免疫微环境的潜在治疗中药。

2 结果

2.1 OA患者差异表达基因结果对来自GSE55235的微阵列进行筛选，共得到10个健康志愿者和10个OA患者滑膜组织芯片，正常滑膜组织和OA滑膜组织差异基因共716个，其中上调基因382个，下调基因334个（图1）。

图1 OA滑膜组织基因芯片差异表达基因热图和火山图Fig.1 Heat map and volcano map of OA synovium-related differential expression genes

2.2 核心靶基因获取结果将716个差异表达基因导入STRING绘制PPI，并将分析结果的TSV文件导入Cytoscape软件进行可视化展示（图2），经cyto Hubba计算列取前10位作为核心基因，分别为IL-6、CXCL8、JUN、VEGFA、IL-1β、MMP9、ITGB2、FOS、APOB、CXCL12（表1）。

表1 PPI模型中OA组和正常组差异核心靶基因Tab.1 Core differential target genes in PPI networks of OA and normal groups

图2 OA滑膜组织差异基因PPI及核心蛋白基因模块Fig.2 OA synovium tissue differential gene PPI networks map and Hub genes

2.3 差异基因GO和KEGG富集分析采用R软件cluster Profiler对OA差异基因进行BP和KEGG分析，上调和下调的差异基因分别为382、334个，GO分析得到BP 374、879条，上调基因共富集到34条通路，下调基因共富集到39条通路，对比上调基因和下调基因GO和KEGG富集分析结果，发现上调基因与免疫密切相关，GO主要涉及白细胞正向调节、吞噬作用、T细胞呈递激活等BP（图3A）。KEGG富集的相关通路主要涉及吞噬小体、TNF信号通路、IL-17信号通路等相关信号通路（图3B），而下调基因GO主要涉及对氧水平的反应、对氧水平下降的反应、节律过程等非免疫BP（图3C），相关信号通路主要涉及毒性心肌炎相关通路、利什曼病相关信号通路等转导机制（图3D）。

图3 OA差异基因GO及KEGG富集分析Fig.3 GO and KEGG enrichment analysis of differential genes related to OA

2.4 免疫浸润细胞分布采用CIBERSORT反卷积法以P<0.05为筛选条件对芯片进行筛选，共得到20个可信样本，热图左侧10个为正常滑膜组织，右侧10个为OA滑膜组，未活化的记忆CD4 T细胞和活化的肥大细胞含量在OA滑膜组织中不同程度减少，而M0型巨噬细胞、未活化的肥大细胞含量在OA滑膜中均出现不同程度增加（图4A），图4B进一步展示了22种免疫细胞在各样本中的含量分布。

图4 正常人群和OA患者免疫细胞浸润分布Fig.4 Infiltration of immune-associated cells in normal and OA samples

2.5 免疫细胞相关性及差异分析未活化的CD4记忆T细胞和活化的肥大细胞（r=0.66）、未活化的树突状细胞和M1型巨噬细胞（r=0.54）、M1型巨噬细胞和活化的CD4记忆T细胞（r=0.49）呈正相关；而未活化的肥大细胞和活化的肥大细胞（r=-0.78）、记忆B细胞和活化的B细胞（r=-0.67）、未活化的CD4记忆T细胞和M2型巨噬细胞（r=-0.62）呈负相关（图5A）。通过小提琴图较为直观地展现正常人群和OA患者滑膜组织免疫浸润细胞差异，以P<0.05表示差异有统计学意义，浆细胞（P=0.019）、M0型巨噬细胞（P=0.038）和未活化的肥大细胞（P<0.001）在OA滑膜组织中含量较高，而未活化的CD4记忆T细胞（P<0.001）、活化的NK细胞（P=0.026）、活化的肥大细胞（P<0.001）在OA滑膜组织中均出现不同程度降低（图5B）。

图5 OA组织中免疫浸润细胞相关性及差异小提琴图Fig.5 Correlation analysis and violin plot of differences of immune cells in OA group

2.6 干预OA滑膜组织免疫浸润机制的中药挖掘通过Coremine Medical预测和挖掘潜在干预OA免疫浸润的中药，以P<0.01作为筛选条件，共得到123味中药，与白细胞黏附、白细胞增殖、吞噬作用、抗原呈递与表达、T细胞呈递等BP及IL-6、CXCL8、JUN、VEGFA、IL-1β、MMP9等核心靶基因密切相关，将结果导入Cytoscape进行可视化分析，共得到节点137个，边226条，青风藤、姜黄、雷公藤、藤黄等度值较高，提示与OA免疫浸润机制密切相关，将这些中药进行归经，功效分析显示OA免疫浸润相关药物均归属肝、肾、肺经，药物功效以清热类、补虚类、祛风湿类、活血化瘀类、解表类等为主，提示以上药物可能作为OA免疫治疗的潜在候选中药（图6）。

图6 OA免疫浸润机制潜在中药功效及雷达图Fig.6 Potential traditional Chinese medicine efficacy map and radar map of osteoarthritis immune infiltra⁃tion mechanism

3 讨论

滑膜组织是OA免疫系统的主要入侵对象，本研究将正常人群和OA患者滑膜组织差异基因进行PPI分析，筛选出10个与OA免疫相关的核心靶基因IL-6、CXCL8、JUN、VEGFA、IL-1β、MMP9、ITGβ2、FOS、APOB、CXCL12，除ITGβ2和APOB基因未见报道以外，其他均发现参与OA滑膜免疫浸润过程。VEGFA可介导早期OA滑膜中性粒细胞、淋巴细胞、肥大细胞浸润诱导血管炎症［12］；JUN和FOS在OA滑膜中高表达，参与炎症细胞弥散浸润［13］；MMP9参与CD4+T细胞调控巨噬细胞浸滑膜病程；IL在OA患者滑膜组织中主要由白细胞、巨噬细胞、滑膜纤维细胞等产生［14］，OA患者滑液、滑膜、软骨下骨和软骨中均不同程度增高［15］；此外，过多IL产生可抑制软骨中Ⅱ型胶原生成，诱导基质金属蛋白酶类产生，加速软骨侵蚀［16］。本研究预测出的IL-6、IL-1β在OA进程中备受关注，已展开大量相关研究［17-18］。而趋化因子（CXCL）是一类具有化学趋化作用的细胞因子，主要由滑膜巨噬细胞产生，可对中性粒细胞、淋巴细胞、单核细胞等细胞产生趋化作用，且与血管炎发展密切相关，介导早期OA出现滑膜炎症［19］。CXCL8可影响JAK-STAT、NF-κB和JNK-MAPK信号通路并增强其他促炎细胞因子表达促进细胞凋亡并抑制软骨细胞增殖［20］；CXCL12可驱动树突状细胞从血液募集至发生炎症反应的滑膜诱导OA慢性炎症［21］。本研究筛选出的大部分核心靶基因与目前OA免疫相关研究结果相符，进一步证实研究结果的可靠性。对OA滑膜组织差异基因进行GO和KEGG富集分析，结果显示上调基因与OA免疫机制密切相关，主要涉及多种免疫细胞，进一步说明OA是由多种炎症细胞及因子共同参与的免疫过程，存在错综复杂的转导过程，而KEGG富集预测出大部分信号通路均介导了OA免疫过程，这些通路在诱导OA滑膜炎症因子渗出、免疫细胞浸润、补体系统激活、体液免疫和细胞免疫过程中均发挥重要作用［22］。

为深入研究免疫细胞在OA滑膜组织的作用，采用CIBERSORT反卷积法对两组免疫相关差异基因进行筛选，CIBERSORT反卷积算法是基于线性支持向量回归（SVR）的可靠机器学习方法，广泛用于评估22种免疫细胞相对含量及动态调节过程，在噪声、未知混合物含量和密切相关的细胞类型方面优于其他方法。最终结果显示浆细胞、M0型巨噬细胞和未活化的肥大细胞在OA滑膜组织中含量较高，而未活化的CD4记忆T细胞、活化的NK细胞、活化的肥大细胞在OA滑膜组织不同程度降低，以上研究结果部分得到相关基础研究证实，如OA滑膜组织中肥大细胞浸润程度较高，与OA患者结构损伤有关，其机制可能与肥大细胞释放大量TNF-α及IL-6等细胞因子诱导炎症反应，刺激滑膜增生和软骨侵蚀［23］；记忆CD4+T细胞积累是关节局部炎症反应过程中的常见现象，可能与慢性OA形成有关［24］；调节性T细胞在OA滑膜中含量丰富，其水平与炎症因子水平相关［25］；而滑液中过量的CD19+CD24hiC7+B细胞可加速软骨破坏［26］；M1型巨噬细胞产生炎症因子与产生抗炎因子的M2型巨噬细胞比例失衡是促进OA进展的关键机制之一，调节M2型巨噬细胞有望对延缓OA起积极治疗作用［27-28］。上述文献证据结合课题组分析结果表明巨噬细胞、浆细胞、肥大细胞等在OA中发挥重要作用，应成为进一步研究重点。

OA属于中医学“痹症”范畴，中医对其认识较早，长时间临床实践显示中药具有良好有效性及安全性，且中药具有多成分、多靶点、多途径作用机制，从中药中寻找针对OA免疫机制整体调节的药物应作为研究重点［29-31］。OA中医病机属于年老体衰，肝肾亏虚，筋骨失养，外感风寒湿邪气，痹阻脉络，流注关节［32］，证属“本虚标实”，以肝肾亏虚为主，风寒湿邪为标实。本次预测药物的归经和功效分析，大部分药物归经多属于肺、肝肾经且药物多为补虚类、祛风湿类、活血化瘀类、解表类等，进一步说明本研究药物预测结果符合OA中医病机，且本研究预测药物以藤类药物居多，中医学受“取类比象”思维影响，认为藤类药物形条达，缠绕蔓延，如同人体经络纵横交错走行，具有通畅经络、舒通筋脉作用，善治痹证久治不愈，甚则变形之顽痹，完全符合OA病症特点。本研究从免疫浸润角度将OA差异基因中免疫相关BP及作用于PPI核心靶基因的中药进行反向预测，筛选出青风藤、姜黄、雷公藤、藤黄等中药在免疫调节干预OA过程中处于核心地位，现代药理学也证实以上药物可通过调节免疫细胞发挥OA治疗作用，如青风藤可抑制TNF-α、IL-β、IL-6等炎症因子产生，减少炎症反应对关节软骨的损伤［33］。此外，青风藤还可抑制滑膜组织血管新生和血管炎症，减少免疫细胞浸润对滑膜的刺激；姜黄中有效成分姜黄素具有强大的抗氧化和抗炎作用，可抑制IL-1β诱导的NF-κB介导的炎症和细胞凋亡，保护软骨［34］；雷公藤除抗炎、镇痛、抗血管新生外，还对滑膜有免疫调节功能，抑制T细胞增殖、诱导T细胞凋亡、恢复Th1/Th2细胞平衡［35-36］。以上文献支持本研究结果的可信度，也为中药干预骨关节炎免疫机制候选中药提供参考。

目前美国风湿病学会（ACR）、欧洲抗风湿联盟（EULAR）和中国2018版《骨关节炎诊疗指南》均推荐针对OA治疗的药物，包括非甾体类抗炎药（NSAIDs）、阿片类药、糖皮质激素、抗焦虑药、软骨保护药等，促使消化道出血、脑、肾和心血管疾病、关节感染、药物依赖性等风险性增加，长期使用饱受争议，治疗目的主要是减轻或消除疼痛，纠正畸形，改善或恢复关节功能，但却未对免疫系统进行调节，远期效果仍不理想［37-40］。对于晚期OA患者关节软骨缺损，组织工程学方法的自体软骨移植修复是目前研究热点［41］。但关节周围组织免疫微环境引发的自身免疫排斥反应也是组织工程材料关节植入需考虑的关键问题。中医药是中国人民几千年来的智慧结晶，其独特魅力受到广泛关注［42-43］。因此，从中草药中寻找干预OA免疫失衡的潜在中药，对早期OA患者减缓病程具有重要意义，对于晚期OA患者关节组织工程材料的特性选择也具有指导意义。

本研究依赖于GEO数据库，是既往发布数据集的二次挖掘和分析，实验结果可能与既往实验结论不同，原因可能为样本量过少导致数据分析偏倚。此外，CIBERSORT反卷积算法分析是基于有限的遗传数据，这些数据可能会偏离细胞异型相互作用、疾病的诱发因素以及疾病表型的可塑性，导致结果不准确。虽然本研究结果得到了多数文献研究支持，但研究结果的可靠性仍需进一步实验验证。