激素性股骨头坏死的关键基因筛选及发病机制探索：基于生物信息学分析

肖立康，董永辉，郑稼

激素性股骨头坏死（steroid-induced osteonecrosis of femoral head, SONFH）是一种非创伤性股骨头坏死（osteonecrosis of the femoral head, ONFH），一般有明确的激素应用史，早期可进行药物治疗及髓心减压术、带血管蒂骨移植等保髋手术，如未进行有效的治疗，80%的患者会在1~4年内发生股骨头塌陷[1]。股骨头一旦塌陷，保守治疗疗效不佳，只能进行全髋关节置换术[2]。随着糖皮质激素广泛应用于血液病、严重急性呼吸综合征、系统性红斑狼疮等免疫系统损伤相关疾病，SONFH的发病率逐年上升，且ONFH患者逐渐年轻化。SONFH的发病机制包括脂肪栓塞、骨髓内压力升高、骨髓间充质干细胞（bone marrow mesenchymal stem cells, BMSC）成骨和成脂分化不平衡及血管生成障碍等[3]。近年来随着基因研究的深入，SONFN的分子机制研究有了较大进展，但SONFH确切的分子发病机制至今尚无定论。因此，有必要继续寻找SONFN相关的潜在生物标志物，并探讨其分子发病机制，从而寻找更有效、更完善的ONFN诊疗方法。本研究拟采用生物信息学分析筛选SONFH新的潜在生物标志物，并探索其相关发病机制。

1 资料与方法

1.1 研究样本

本研究从高通量基因表达（gene expression omnibus, GEO）数据库（http://www.ncbi.nlm.nih.gov/geo）提取了来源于GPL15207平台的数据集GSE123568，其中包括30例SONFH患者和10例非SONFH患者（类固醇给药后）。

1.2 筛选差异表达基因（differential expressed gene,DEG）

使用R语言程序（4.1.1版）的LIMMA（微阵列数据的线性模型）软件包对样本进行分析，筛选样本中阈值标准|log2FC|＞1和P＜0.05的DEG。然后，在R语言程序中使用Pheatmap软件包（1.0.12版）绘制数据集DEG的热图。R语言程序中使用的参数包括分位数、Median polish和PM矫正。

1.3 DEG的富集分析

基因本体论（gene ontology, GO）为功能基因组学研究提供全面的资源，包括各个物种的生物过程（biological process, BP）、细胞成分（cell component,CC）和分子功能（molecular function, MF）。京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes, KEGG）整合了所有完整测序的基因组和部分基因组，同时附带最新的基因功能注释。使用David数据库（6.8版，https://david.ncifcrf.gov/）进行富集分析，从分子水平和功能水平研究DEG。

1.4 基因富集分析（gene set enrichment analysis,GSEA）

GSEA是一种不同于GO分析和KEGG分析的富集分析方法，可以协助找出表达差异不显著却同样值得关注的基因。运行R语言程序的ClusterProfiler软件包（3.10.1版）对所有基因进行GSEA，为筛选DEG的生物学功能提供另一种选择。用R语言程序的Enrichplot软件包（1.2.0版）绘制ES图、热图和Rank图。显著富集标准：P＜5%，FDR＜25%。

1.5 构建蛋白相互作用（protein-protein interaction,PPI）网络

通过搜索String数据库（http://www.string-db.org/）来预测PPI网络，构建了以综合得分＞0.4为截断点的DEG的PPI网络。应用Cytoscape软件的MCODE插件选出最显著模块，设置标准：Degree Cutoff=2、Max depth=100、Node Score Cutoff=0.2、K-core=2。再次使用David数据库对该模块中的基因进行GO分析和KEGG分析。

1.6 筛选关键基因

使用Cytoscape软件的Cytohubba插件运行12种方法筛选PPI网络中的关键基因，Degree、BottleNeck、Radiality、Stress、Closeness较适合人类。筛选出关键基因后再次构建PPI网络。

1.7 关键基因在BMSC中的表达

1.7.1 样本分组

本研究选择人BMSC［赛业（苏州）生物科技有限公司］进行验证，使用不同浓度的地塞米松（dexamethasone, DEX）在37℃、95% O2、5% CO2下处理BSMC 72 h，以模拟细胞水平SONFH的发生过程。将BMSC分为5组：①空白对照组，BMSC不进行任何处理；②二甲基亚砜（dimethyl sulfoxide, DMSO）组，用DEX的溶媒DMSO处理BMSC，排除因溶媒引起的mRNA表达变化；③10-8mol/L组，用10-8mol/L的DEX处理BMSC；④10-7mol/L组，用10-7mol/L的DEX处理BMSC；⑤10-6mol/L组，用10-6mol/L的DEX处理BMSC。

1.7.2 RNA提取

选择超纯总RNA快速提取试剂盒（北京博迈德基因技术有限公司）进行RNA提取，向培养板中加入裂解液RL后迅速摇匀，使RL充分与细胞接触，裂解细胞并灭活RNA酶，使用移液枪吹打混匀。

将得到的混匀样品室温下放置孵育5 min。每1 mL RL加0.2 mL氯仿，盖紧样品管盖，剧烈振荡15 s，并将其在室温下放置孵育3 min。

在4℃环境下设定转速为12 000 r/min（r=5 cm，下同），将样品离心10 min。取上清液转移至新管中，加入等量的70%乙醇溶液，颠倒混匀后转移至吸附柱RA中。然后在转速12 000 r/min下离心45 s，弃去废液，将吸附柱RA再次套回收集管中。加400 μL RE，12 000 r/min离心45 s，弃去废液。

取新的RNase-free离心管，加入10 μL DNase Ⅰ储存液中和70 μL RDD溶液，轻柔混匀后加入吸附柱RA中央，室温放置15 min后加入400 μL RE，在转速12 000 r/min下离心45 s，弃去废液。然后再加入500 μL漂洗液RW，在转速12 000 r/min下离心45 s，弃去废液。将吸附柱RA放回空的收集管中，再次12 000 r/min离心2 min后去除漂洗液。

取出吸附柱RA，放入RNase-free离心管，根据预期RNA产量向吸附柱RA中央加入适量RNase-free H2O，室温放置2 min后在转速12 000 r/min下离心1 min。使用超微量核酸分析仪NANO-400A（杭州奥盛仪器有限公司）测定RNA浓度，并记录。

1.7.3 cDNA合成

使用HiFiScript cDNA Synthesis Kit试剂盒（江苏康为世纪生物科技股份有限公司）进行cDNA合成。将RNA模板、Primer Mix、dNTP Mix、DTT、RT Buffer、HiFiScript和RNase-free H2O溶解并置于冰上备用。配制成20 μL反应体系进行基因组DNA去除和cDNA合成（表1）。将混合液涡旋震荡混匀，短暂离心使溶液收集到管底。42℃下孵育15 min，然后85℃下孵育5 min。待反应结束后短暂离心，然后置于冰上等待下一步操作。逆转录产物进行后续实时荧光定量PCR。

表1 去除基因组DNA和反转录反应体系

1.7.4 实时荧光定量PCR

使用Universal SYBR Green qPCR Super Mix试剂盒（美国EVERBRIGHT公司）进行实时荧光定量PCR，取2×Universal SYBR Green qPCR Super Mix、引物、cDNA、RNase-free H2O，恢复至室温，按照表2制备反应体系，轻轻涡旋混匀后转移至PCR管进行扩增。

表2 荧光定量体系

使用PCR检测系统对mRNA转录本进行定量。扩增条件：95℃作用5 min，然后分别以95℃作用5 s和60℃作用30 s，进行40个循环，建立熔解曲线，获得定量数据。以GAPDH为参照基因，进行RT-qPCR反应，重复3次。采用标准比较法（ΔΔCt）评价候选基因的表达稳定性。用2-ΔΔCt法计算相对靶基因表达水平，多样本间比较采用单因素方差分析，应用GraphPad Prism 9软件进行统计学分析及图像处理。本研究中使用的引物序列如表3所示。

表3 实时荧光定量PCR基因引物序列

2 结果

2.1 确定SONFH的DEG

数据集GSE123568中包含402个DEG，其中189个基因表达上调，213个基因表达下调。图1和图2分别显示了所有DEG的火山图和差异明显的前50个表达上调基因和前50个表达下调基因的热图。

图1 DEG的火山图

图2 差异基因表达热图

2.2 GO和KEGG富集分析

考虑到某些表达上调基因和表达下调基因可能参与同一个通路，所以本研究将所有DEG一起进行GO和KEGG富集分析，运行R语言程序ggplot包进行分析和可视化（图3）。GO富集分析结果表明，BP的DEG显著富集了炎症反应、免疫反应、血小板脱颗粒。对于CC，DEG在质膜、胞质、外泌体中显著富集。对于MF，DEG在蛋白结合、受体活性、趋化因子受体活性和铵盐跨膜转运蛋白活性中显著富集。KEGG富集分析结果表明，DEG主要富集于破骨细胞分化、细胞因子-细胞因子受体相互作用、趋化因子信号通路。

图3 GO与KEGG富集分析气泡图

2.3 GSEA

使用GSEA富集ONFH所有可能的通路（图4）。其中最显著富集的基因集与ONFH的发生呈正相关，包括细胞因子-细胞因子受体相互作用效应、趋化因子信号通路、自然杀伤细胞介导的细胞毒作用、JAK-STAT信号通路、Toll样受体信号通路、溶酶体。

图4 GSEA富集最明显的通路

图5 关键基因的筛选过程

图6 DEX处理后mRNA表达变化

2.4 构建PPI网络和确定关键基因

基于String数据库搜索到402个DEGs互相作用信息，利用Cytoscape软件聚集成一个由326个节点和1 193条边组成的网络（图5A），使用MCODE插件确定了最显著模块（图5B）。对最显著模块中的基因再次进行了富集分析。结果显示，最显著模块中的基因主要在GO中的免疫反应、Toll样受体（Toll-like receptors, TLR）1∶TLR2信号通路和三酰化细菌脂肽对细胞的反应中显著富集，在KEGG中主要在吞噬体、破骨细胞分化和Toll样受体信号通路等通路中富集（表4）。

表4 显著模块基因进行富集分析的结果

选择Cytohubba插件中5种最合适分类方法筛选中心基因（表5）。最后从中取交集，得到了8个关键基因，分别是TLR2、TYROBP、IRF8、SNCA、PECAM1、SLC4A1、HCK、CD14（图5C），对关键基因再次构建PPI网络（图5D）。

表5 Cytohubba插件5种方法筛选的关键基因

2.5 关键基因在BMSC中的表达

为了进一步验证上述关键基因，本研究选择了TLR2、IRF8、HCK、CD14、SNCA和SLC4A1进行了验证。

SONFH组CD14、SNCA、HCK、IRF8和SLC4A1在DEX处理后表达显著增加（P＜0.05），而TLR2的表达并无明显增加，在使用10-7mol/L DEX处理后其表达反而有所减少（图6）。CD14在使用10-8mol/L DEX处理后较DMSO处理后表达显著增加（P＜0.05），在使用10-7mol/L和10-6mol/L DEX分别处理后较10-8mol/L DEX处理后其表达均显著增加（P＜0.001）。SNCA在使用10-6mol/L DEX处理后较DMSO、10-8mol/L和10-7mol/L DEX分别处理后其表达均显著增加（P＜0.001）。SLC4A1在使用10-6mol/L DEX处理后较DMSO、10-8mol/L、10-7mol/L DEX分别处理后其表达均显著增加（P＜0.001）。IRF8在使用10-8mol/L DEX处理后较DMSO处理后表达显著增加（P＜0.05），在使用10-6mol/L DEX处理后较DMSO、10-8mol/L和10-7mol/L DEX分别处理后表达均显著增加（P＜0.001）。HCK在使用10-8mol/L DEX处理后较DMSO处理后表达显著增加（P＜0.001），在使用10-7mol/L DEX处理后较10-8mol/L DEX处理后表达却显著减少（P＜0.01），而在使用10-6mol/L DEX处理后较10-8mol/L和10-7mol/L DEX分别处理后表达均显著增加（P＜0.001）。

3 讨论

ONFH是股骨头缺血引起骨细胞和骨髓成分死亡及之后的修复，继而导致股骨头结构改变、股骨头塌陷，引起髋关节疼痛、功能障碍的疾病，按病因可分为创伤性和非创伤性两大类。随着糖皮质激素在临床上的广泛应用，以及一部分患者因病情需要长期服用激素，使得SONFH的发病率逐年上升。SONFH一旦发生，随着病程进展会造成股骨头不可逆损伤，致残率极高，严重影响患者的生活。本研究采用生物信息学分析筛选SONFH的DEG，并对其发病机制进行进一步探讨。

3.1 富集通路的具体作用

在富集分析过程中，本研究将表达上调基因和表达下调基因一同分析，避免遗漏表达上调基因与表达下调基因共同的通路。GO富集分析结果显示，DEG主要富集于免疫反应和炎症反应过程。免疫反应与SONFH发生、发展密切相关，肿瘤坏死因子-α激活巨噬细胞打破成骨与破骨的平衡最终可导致骨坏死[4]。同时，在BP中，血小板脱颗粒也明显富集，血清-糖皮质激素诱导激酶（serum and glucocorticoids-inducible kinase-1, SGK1）在DEG中表达显著上调，推测可能是遗传性或者获得性因素导致，继而出现个体对糖皮质激素敏感性不同的情况[5]。SGK1可调节血小板Orai1蛋白丰度，随后活化血小板，出现血小板脱颗粒，使局部处于高凝状态[6]。激素可直接损伤骨微血管内皮细胞，产生过多的活性氧自由基，致使血管损伤或栓塞，最终发生局部血供中断而坏死[7]。推测血小板脱颗粒可能与这种微血管损伤学说呈协同作用，共同导致SONFH。KEGG富集分析结果表明，DEG主要参与破骨细胞的分化。近年研究显示，破骨细胞在SONFH进展和股骨头塌陷的发生中起重要作用，股骨头坏死区的破骨细胞明显增多[8-9]。股骨头是糖皮质激素过量而产生不良反应的敏感解剖部位，前人利用小鼠模型证明，糖皮质激素可延长破骨细胞寿命，缩短成骨细胞寿命，骨质破坏与硬化骨的形成导致了股骨头塌陷这一结局[10]。研究破骨细胞分化过程中的分子机制对SONFN的早期诊断与治疗可以起到指导作用。

对于GSEA，本研究将表达矩阵中包括DEG在内的所有基因放在一起进行分析。结果表明，样本中患者的大部分基因主要富集于细胞因子-细胞因子受体相互作用、趋化因子信号通路、JAK-STAT信号通路、自然杀伤细胞的细胞毒性作用、Toll样受体信号通路和溶酶体。这些都被证实是糖皮质激素促进炎症反应的重要机制，从而参与了SONFH的发病机制。糖皮质激素促进Toll样受体信号通路激活是促进炎症的重要因素，根据GO富集分析结果，考虑是通过激活MyD88信号通路来实现的。此外，还涉及溶酶体途径，推测自噬参与了SONFH的发病机制。自噬是一把双刃剑，可以通过清除损伤细胞器来保护细胞，过度自噬则会导致细胞死亡[11]。糖皮质激素可通过诱导转录调控来保护破骨细胞破坏成骨细胞，进一步导致ONFH，然而这仍需要进一步的实验来验证。JAK-STAT信号通路在调节炎症反应过程中起着至关重要的作用，多种促炎细胞因子会通过此途径介导软骨的炎症与降解[12]。对PPI网络的显著模块基因进行GO和KEGG富集分析，结果与之前的分析一致。这些机制可能解释了SONFH发生的具体过程。

3.2 关键基因的具体作用

糖皮质激素可能是通过改变关键基因的表观遗传状态从而调控相关基因表达，最终发挥在人体内抗炎、允许、免疫抑制等作用[13]。为了筛选关键基因，本研究借助Cytohubba插件，选择了5种适用于人类的测序方法，结果取交集，得到了8个关键基因：TLR2、TYROBP、IRF8、SNCA、PECAM1、SLC4A1、HCK、CD14，其中SNCA和SLC4A1在SONFH组表达下调，其余6个基因皆表达上调。TLR2是最为重要的关键基因，已有文献报道同样位于细胞膜上的TLR4在小鼠SONFH模型发挥重要作用。根据上述结果推测，TLR2与TLR1协同介导的先天免疫反应，通过MyD88起作用，然后启动了股骨头的破坏程序。通常糖皮质激素能够抑制单核细胞表达和释放CD14（脂多糖），但是本研究结果显示CD14在SONFH组表达上调，这种现象可能与糖皮质激素剂量相关而产生不同的效应[14-15]。高表达的CD14可参与到TLR2介导的通路中导致核因子κB（nuclear factor κB, NF-κB）激活、细胞因子分泌和炎症反应[16]。SNCA（α-突触核蛋白）能参与多巴胺释放和运输的调节，可导致caspase-3激活减少并降低神经元对各种凋亡刺激的反应。TLR2可介导胞外SNCA聚集体的神经毒性和促炎效应，同时本研究筛选出表达下调的SNCA，这将间接导致激活增加，之前已有报道提出caspase-3参与了与ONFH相关的细胞凋亡，然而目前这种机制在骨组织中需要更多的实验来进行验证[17]。HCK（酪氨酸蛋白激酶）是非受体酪氨酸蛋白激酶家族成员，主要表达于髓系细胞，在调节先天免疫反应中起重要作用，包括中性粒细胞、单核细胞、巨噬细胞和肥大细胞功能、吞噬作用、细胞存活和增殖、细胞黏附和迁移。目前已经证实HCK抑制剂可改善与炎症相关的骨破坏，这可能是通过介导NLRP3（结节样受体家族蛋白3）炎性小体的激活来实现的[18-19]。TYROBP（酪氨酸蛋白激酶结合蛋白）在关键基因中仅次于TLR2，可与HCK结合进行接下来的一系列反应[20]。IRF8（干扰素调节因子8）可起到转录激活剂或阻遏物的作用，对破骨细胞的形成至关重要[21]。PECAM1（血小板内皮细胞黏附分子）是白细胞跨内皮迁移所需的细胞黏附分子，参与大多数炎症反应过程[22]。同时PECAM1对先天免疫和获得性免疫都是必不可少的，使组织能够对炎症刺激做出适当的反应，属于一种保护性机制，本课题组推测这是在ONFH发生过程中的代偿性高表达，在股骨头缺血时代偿性地促进新血管的生成。SLC4A1（阴离子转运蛋白）具有维持正常红细胞形状的功能，也是维持人类肾脏远端碳酸氢盐正常酸碱平衡所必需的，疾病组中SLC4A1下调导致的肾小管性酸中毒将需要骨骼释放钙盐来缓冲，间接影响股骨头的骨质状况。

3.3 验证结果的差异

本研究使用不同浓度DEX处理BMSC，在细胞水平模拟SONFH发生过程。然而DEX处理后BMSC中TLR2mRNA的表达结果并不令人满意，不排除技术原因导致；其次可能与TLR2在不同样本类型表达不同有关，本课题组之前得出的TLR2表达上调的结果来源于血清样本，本研究的验证过程则来源于细胞样本。结果中SNCAmRNA表达与之前分析相反，TLR2可介导胞外SNCA聚集体的神经毒性和促炎效应。Raffaele等[23]通过动物实验发现小鼠敲除SNCA后，成熟红细胞中的活性氧和一氧化氮合酶减少，因此推测SNCA可能在非神经组织中具有促进氧化应激的作用。SLC4A1基因编码的蛋白又被称为阴离子交换蛋白1（anion exchanger 1, AE1），在维持细胞稳态方面起重要作用。氧化应激会引起AE1结构的变化进而影响其转运功能，最终引起细胞损伤，在BMSC中检测到表达上调的SLC4A1可能是代偿性反应[24]。CD14、HCK和IRF8mRNA表达明显上调，与本研究分析结果一致。

3.4 优点与不足

上述标志物可能在SONFH发生、发展过程中产生重要作用。本研究通过GEO数据库进行深入的生物信息学分析，并初步进行了细胞水平表达的验证。本研究在String数据库中应用了高可信互相作用网络，对各个基因之间的相互作用进行了可靠的逻辑学分析。此外，本研究利用GSEA防止遗漏重要的通路，完善研究中SONFH的发病机制。但是，本研究也存在一定的局限性和不足。首先本研究验证这项研究的样本量偏少，相对缺乏说服力，今后将进行大样本及设计更严格的细胞生物实验来进一步验证。其次本研究进行的二次数据挖掘没有实验组和对照组研究对象的基本特征资料，缺乏糖皮质激素剂量、疗程等信息。因数据库的更新等原因使得本研究结果和结论存在一定的局限性，需进行进一步的研究。

4 结论

TLR2、TYROBP、IRF8、SNCA、PECAM1、SLC4A1、HCK、CD14共8个关键基因是SONFH发病过程中的潜在标志物。同时还发现脂多糖介导的免疫反应在SONFH发生发展的过程中发挥重要作用，CD14、HCK和IRF8可能是介导这一过程的主要标志物。本研究揭示了SONFH的可能途径，对SONFH的早期诊断、潜在治疗靶点的发现提供了新的参考。

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