基于生物信息学分析鉴定激素性股骨头坏死的关键生物标志物及相关免疫细胞浸润

葛彦志，黄培洁，瞿杭波，陈佳烨，陈祖祥，陈文娴，夏清馨，童培建,3

股骨头坏死（osteonecrosis of the femoral head,ONFH）是导致髋关节疼痛和患者残疾的一个重要原因。目前将ONFH分为两类：创伤性和非创伤性。创伤性ONFH 主要是由股骨头血液供应的急性机械破坏引起的；非创伤性ONFH 中激素性ONFH（steroidinduced ONFH,SIONFH）是最常见的形式[1]。然而，SIONFH 发生的确切机制仍不清楚。目前，关于SIONFH的发病机制主要有以下五种理论：①脂质代谢紊乱理论；②骨髓间充质干细胞降低成骨潜能的理论；③血液供应不足理论；④炎症与细胞凋亡的理论；⑤基因多态性和非编码RNA 理论。SIONFH可能是由以上多种机制共同作用的结果[2]。由于骨坏死组织难以恢复、高残疾率，以及部分患者仍需长期使用糖皮质激素等原因，SIONFH在医学界受到广泛关注。SIONFH 的发生发展会对股骨头造成不可逆的损伤，从而严重影响患者的生活质量。因此，研究SIONFH的发病机制对其预防、诊断和治疗等方面都具有重要意义。

生物信息学是生物学与信息学的交叉科学，其研究对象主要集中在基因和蛋白质两个方面，并且在生命科学的研究中发挥着至关重要的作用[3]。在研究基因表达谱与生物学功能之间可能存在的联系时主要运用基因芯片技术，该技术可一次性对上万个基因的表达谱进行检测，极大推动了生物信息学技术的进步[4-6]。

本研究利用生物信息学方法对来自基因表达综合数据库（gene expression omnibus,GEO）中SIONFH患者外周血清中的基因芯片数据进行差异表达基因（differentially expressed genes,DEGs）分析，并计算SIONFH患者外周血清中免疫细胞之间的相关性，通过挖掘到的DEGs分析其分子作用关系，构建疾病相关DEGs 的分子调控网络，从而进一步探究SIONFH发生的分子机制、SIONFH患者外周血清中浸润性免疫细胞的组成特点及其在疾病发生发展过程中扮演的角色。

1 资料与方法

1.1 SIONFH样本数据的获取

由GEO 数据库在线检索人类基因芯片样本，以“Femoral head necrosis”为检索策略，获得GPL15207平台上提供的芯片数据GSE123568，所用实验平台为美国昂飞公司（Affymetrix,Inc）的人类全基因芯片数据。共纳入40 例样本，其中SIONFH 组30 例，非SIONFH组10例。

1.2 基因数据统计方法

使用Perl语言（版本：5.34.0）对探针ID数据进行注释，并转换为基因名称。使用R 语言（版本：3.4.3）编辑软件进行芯片数据预处理和分析，使不同样本之间归一化并具有可比性，并将基因表达原始数据进行标准化。通过BiocManager 的“Affy”包及RMA函数进行背景校正和归一化处理后，对数据进行汇总从而获取校正后表达水平的标准化数据[4]。通过安装R语言中的“limma”包对基因表达矩阵进行差异计算，并用贝叶斯方法进行多重检验校正。通过两组间表达量的比值（fold change,FC），并对其取以2为底的对数即log2FC 和P值筛选DEGs（筛选条件：|log2FC|≥1 和校正P＜0.05，负值代表基因下调，正值代表基因上调），并R语言中的“ggplot2”和“heatmap”包分别绘制DEGs 的火山图和热图，从而直观地对DEGs进行可视化[7,8]。

1.3 DEGs的基因本体论（gene ontology,GO）富集分析

GO数据库是一种对基因及其蛋白质产物功能进行系统描述的数据库，在分析基因间的功能相似性和基于高通量生物学数据分析疾病相关的生物学功能及通路上被广泛应用[9]。使用R语言的“clusterProfiler”包，对上调及下调DEGs分别进行GO富集分析。GO富集分析分为三部分：细胞成分（cellular component,CC）、生物过程（biological process,BP）和分子功能（molecular function,MF）。错误发现率（false discovery rate,FDR）＜0.25 和P＜0.05 被认为其功能具有显著差异[10]。

1.4 SIONFH相关靶基因的筛选

本研究中，利用关键词“Femoral head necrosis”筛选SIONFH 相关的靶基因信息并来源于以下5 个数据库：①GeneCards（https://www.genecards.org），获得SIONFH 相关靶标基因2 082个；②DrugBank[11]，该数据库未鉴定出与SIONFH 相关的基因；③TTD（therapeutic target database,http://db.idrblab.net/ttd），未发现SIONFH 相关基因；④OMIM（online mendelian inheritance in man®,http://www.omim.org）[12]，从中发现102 个与SIONFH 相关的基因；⑤PharmGKB（https://www.pharmgkb.org），确定了126 个预测基因。所有获取的靶基因删除重复数据后与DEGs 进行交集，绘制韦恩图（http://bioinformatics.psb.ugent.be/webtools/Venn），共收集62个与SIONFH相关的靶基因，并进行进一步分析。

1.5 蛋白质网络的构建（protein-protein interaction networks,PPI）

STRING在线数据库（http://string-db.org）是一个对已知和预测的蛋白进行互作分析的数据库和网络资源。Cytoscape 3.4.0 软件（http://www.cytoscape.org）可以将生物分子交互网络与高通量基因表达数据和其他的分子状态信息整合在一起，以进行分子可视化分析[13]。本研究利用STRING数据库在线检索，预测靶蛋白质之间直接和间接的相互作用。将筛选出的DEGs 输入到STRING 网站数据库中，选取交互作用评分＞0.4 分（中等置信度）的相互作用关系，以构建SIONFH组和非SIONFH组之间相关DEGs的互作用网络。最后，将STRING数据库中得到的蛋白质相互作用结果导入Cytoscape 软件中，进行网络分析及可视化操作，建构可视化的分子交互作用网络，并且对蛋白质之间、蛋白质和DNA 之间等交互作用进行分析。利用软件中的CytoNCA[14]插件同时计算各蛋白之间相互关联紧密程度，并进一步筛选出与SIONFH关联紧密的关键靶蛋白[15]。

1.6 获取免疫细胞矩阵及免疫细胞浸润的相关性分析

使用Perl软件将基因表达数据整理成行名为基因名、列名为样品名的基因矩阵，采用R语言BioManager中的“limma”包对SIONFH 患者和非SIONFH 组外周血清中的mRNA 表达谱矩阵进行校正。用CIBERSORT 法对人类免疫细胞亚型的表达矩阵进行去卷积分析，计算22种免疫细胞的相对比例，且每个样本获得一个P值，根据P＜0.05筛选样本，获得免疫细胞组成矩阵。利用R语言“graphics”包的barplot函数绘制两组样本中每种免疫细胞构成比的柱状图；利用R语言中的“corrplot”包对SIONFH患者外周血清中的免疫细胞进行相关性分析；利用R语言中的“vioplot”包对SIONFH 患者和非SIONFH 组外周血清中免疫细胞占比进行比较分析并绘制小提琴图[16]。

2 结果

2.1 数据标准化

由于基因原始数据表达值的中位数值显示不一，需要采用R 语言软件中的“Affy”包及RMA 函数对基因表达原始数据进行背景校正和归一化处理，最后对数据进行汇总从而获取校正后基因表达水平的数据，校正前后的数据如图1所示。

图1 GSE123568样本原始数据（A）和校正后数据（B）的箱线图

2.2 DEGs

选取|log2FC|≥1 及P＜0.05 作为阈值，从而得到SIONFH 的384 个DEGs。图2A、B 分 别为所筛选DEGs的火山图及热图，表1为排名前10的上调和下调DEGs列表。

表1 排名前10的上下调DEGs

图2 SIONFH相关的DEGs表达结果

2.3 GO功能注释

使用DAVID在线工具对SIONFH组与非SIONFH组患者外周血清中的DEGs进行GO 富集分析（包括CC、BP 和MF）。上调的DEGs 在生物学过程中主要富集在免疫反应、细胞防御反应、模式识别受体信号通路、神经炎症反应等方面；在细胞组分中主要富集在分泌颗粒膜、质膜外侧和内溶酶体等方面；在分子功能中主要富集在模式识别受体活性、趋化因子受体活性、免疫受体活性等方面（图3A）。下调的DEGs 在生物学过程中主要富集在组蛋白修饰、共价染色质修饰、肽基赖氨酸修饰和RNA 剪接等方面；在细胞组分中主要富集在核斑点、中心粒和前核糖体等方面；在分子功能中主要富集在组蛋白结合、丝氨酸/苏氨酸激酶活性和转录共调节活性等方面（图3B）。

图3 上调DEGs（A）和下调DEGs（B）GO功能富集分析条形图

2.4 靶基因的获取及DEGs蛋白质相互作用分析

通过对5个疾病相关数据库的筛选，最终筛选出2310 个SIONFH 相关基因，与384 个DEGs 取交集后共发现62 个疾病相关靶基因（图4A）。随后，应用STRING数据库在线检索、预测蛋白质之间的相互作用，采用Cytoscape 软件及其CytoNCA 插件，最终筛选出与疾病相关性最高的前3 个蛋白靶点（图4B），分别为精氨酸酶1（Arginase-1,ARG1）、Toll 样受体2（toll-like receptors-2,TLR2）和前列腺素内过氧化物合成酶2（prostaglandin-endoperoxide synthase-2,PTGS2）。

图4 SIONFH靶基因的获取和核心基因的鉴定

2.5 免疫细胞浸润结果

免疫细胞占比柱状图中展示了22种免疫细胞的浸润比例（图5A），与非SIONFH 组比较，SIONFH 组患者外周血清中的中性粒细胞（39.5%±11.2%）、单核细胞（21.2%±8.0%）和CD8+T 细胞（14.2%±6.7%）浸润水平较高；静息的CD4+记忆T 细胞（＜0.01%）、树突细胞（dendritic cell,DC）（＜0.01%）和滤泡辅助性T 细胞（T cells follicular helper,Tfh）（＜0.01%）浸润水平较低。相关性热图（图5B）结果显示活化的DC 与Tfh 呈正相关（r=0.670,P＜0.05），即当活化的DC 减少时，Tfh 也相对减少。CD8+T 细胞与中性粒细胞呈负相关（r=-0.500,P＜0.05），即当CD8+T 细胞减少时，中性粒细胞相对增多。通过小提琴图对SIONFH 组与非SIONFH 组间免疫细胞浸润进行差异分析，相较于非SIONFH 组，SIONFH 组活化的DC（P=0.002）浸润程度显著降低，即SIONFH 组活化的DC 较 非SIONFH 组明显减少；而SIONFH 组的记忆B 细胞（P=0.014）浸润程度显著提高，即SIONFH 组的记忆B 细胞较非SIONFH 组明显增加（图5C）。

图5 SIONFH组及非SIONFH组外周血清中22种免疫细胞占比的柱状图（A）、热图（B）和小提琴图（C）

3 讨论

3.1 生物信息学的发展

生物信息学作为一门通过利用数理和信息科学的观点、原理和方法来研究生命现象的学科，在生命科学的研究中发挥着极其重要的作用[3,7,17]。为了从分子层面了解疾病的发生发展过程，通过基因芯片表达谱数据进行分析，从而将获得的信息数据和生物学过程联系起来，以解释基因的功能。

3.2 SIONFH的基因表达谱研究

SIONFH的发生是由于使用激素类药物，从导致股骨头血液供应受损，骨细胞和骨髓成分死亡，进而导致结构改变、股骨头塌陷和关节功能障碍。患者最终需要通过THA 进行治疗，但严重损害了患者的生活质量[2]。既往研究表明，可以通过基因表达来分析SIONFH患者的发病原因并进行探讨，可以从不同的全基因组表达谱中获得互补的结果。然而，尽管大量独立的全基因组方法已被用于SIONFH 组与非SIONFH组独立样本中的DEGs分析，但该疾病发生的潜在分子机制尚不清楚。本研究获得了SIONFH组与非SIONFH组外周血清中的基因表达分析，筛选出与SIONFH相关的靶基因，并通过生物信息学分析探讨了SIONFH 存在的潜在关键基因，通过CIBERSORT法对SIONFH 患者的外周血清进行免疫细胞浸润分析，初步探讨免疫细胞在SIONFH中的作用。因此，可以从SIONFH患者与非SIONFH患者之间的基因表达分析中发现一些可能与SIONFH发病机制或治疗存在密切联系的DEGs。在以往的生物信息学分析中，通过加权基因共表达网络分析（weighted correlation network analysis,WGCNA）方法对芯片GSE123568分析后发现94 个靶基因，对其进行一系列的分析，如PPI互作、基因富集分析、转录因子富集分析及竞争性内源RNA。本研究不仅对该芯片进行DEGs 分析及功能富集分析，进一步从SIONFH相关的5个数据库进行发掘。同时，通过免疫细胞浸润分析，进一步对SIONFH的生物信息学研究进行探究。

3.3 SIONFH相关数据讨论

本研究中在GO富集分析（包括CC、MF及BP）中发现中性粒细胞脱粒显著上调，并且参与免疫反应的活化中性粒细胞也显著上调。以往研究表明，在先天免疫中起关键效应的细胞是中性粒细胞，其可被迅速招募从而保护宿主免遭入侵病原体的侵害。中性粒细胞杀死病原体的方式主要有两种：①在吞噬作用后再在细胞内杀死病原体；②在细胞外通过脱颗粒和释放中性粒细胞胞外陷阱来杀死病原体。所有杀死病原体的策略都需要细胞毒性蛋白质和蛋白酶的参与，而这些蛋白质和蛋白酶在中性粒细胞发育过程中都被包装在细胞质颗粒中[18]。而Toll 样受体（toll-like receptors,TLR）作为中性粒细胞表面的7种先天免疫受体之一，有研究表明，TLR信号传导由病原体相关的分子模式触发，这些模式介导完善的细胞因子驱动途径，与IRF3/IRF7 一起激活核因子κB（nuclear factor kappa-B,NF-κB）。在TLR2、TLR3等刺激过程中抑制或消除半胱天冬酶-8，随后导致受体相互作用蛋白激酶3（receptor-interacting protein kinase 3,RIP3）程序性坏死，这种坏死是通过含有Toll样/白细胞介素受体1（Toll/interleukin-1 receptor,TIR）结构域的接头诱导干扰素-β（TIR-domain-containing adaptor inducing interferon-β,TRIF）或髓样分化因子88（myeloid differentiation factor 88,MyD88）信号转导的，最终导致细胞程序性坏死[19]。本研究构建的PPI网络中，TLR2 显著上调，表明中性粒细胞表面受体TLR2可能通过细胞程序性坏死，从而引起SIONFH的发生。此外，在生物过程中上调的GO 还有IL-8[20]。IL-8可刺激破骨细胞生成或抑制成骨细胞形成，是骨破坏的介质[21]。因此，结合先前的研究结果推测，SIONFH 的发生可能与IL-8 对骨细胞的影响有关。在生物过程中下调的GO有组蛋白修饰。研究表明，组蛋白具有从核小体突出的尾部，这些尾部中的许多残基可以进行翻译后修饰，从而影响DNA 的翻译过程（包括染色质压缩、核小体动力学和转录），而组蛋白核心修饰不仅可以直接调节转录，还可以影响DNA 修复、复制和细胞状态变化等过程[22]。虽然增强干细胞的分化能力是一个难题，但目前可以通过组蛋白修饰来增强干细胞的分化能力，从而为损伤修复和再生提供更坚实的基础[23]。在GO 富集分析中，组蛋白修饰下调与股骨头坏死有密切联系，可能是由于组蛋白修饰下调，导致干细胞分化能力减弱，损伤修复和再生减弱，从而对SIONFH造成影响。

3.4 本研究的局限性

首先，本研究仅从数据库及血清层面进行生物信息学分析，与临床实际情况存在一定差别，有待进一步实验验证；其次，样本量相对较少，对结果造成一定影响，今后将通过增大样本量进一步验证。

4 结论

多种途径参与SIONFH发生发展的过程，炎症因子和细胞因子可能是SIONFH发病中的重要环节。基于PPI网络，ARG1、TLR2和PTGS2被确定为SIONFH患者外周血清中的潜在核心靶点，与SIONFH的发生和发生有着密切的联系。免疫细胞浸润结果表明SIONFH 患者的外周血清的免疫细胞浸润模式与非SIONFH患者相比具有一定的特异性。其中，活化的中性粒细胞、单核细胞和DC 对SIONFH 的不同作用可以对SIONFH 的研究方向起到提示作用。从生物信息学分析SIONFH可能的发病原因与机制，为未来的科学研究和临床治疗提供了思路和参考。

【利益冲突】所有作者均声明不存在利益冲突