多芯片整合分析腰椎间盘突出症发生相关的分子通路及关键基因研究*

刘权祥，陈 钱，齐明宇，庄新明，郭建平△，李增新，程 维，代起志

(1.北华大学附属医院骨外二科，吉林吉林 132011；2.吉林省吉林市中心医院内分泌骨代谢科 132011；3.吉林大学第一医院脊柱外科，长春 132021)

腰椎间盘突出症(lumbar disc herniation，LDH)是一组以椎间盘(intervertebral disc，IVD)退行性变，间盘内纤维环破裂，髓核在切应力作用下向内突，进一步压迫神经根或马尾神经而引起的以神经症状为主的临床常见病和多发病[1]。新近的流行病学研究提示，在亚洲人群中，LDH的年总发病率为15 877/10万(其中男性为13 181/10万；女性为18 588/10万；男女比为1.00∶1.41)。此外，每5年患者的发病率和年治疗费用分别增加7.6%和14.7%。年龄矫正后的报道提示，LDH发病率随年龄增长而增加，75～79岁人群发病率高达42.6%。其中，年龄大于65岁患者占31.0%，且其所占医疗费用高达40.1%[2-3]。

正常的IVD内，由于髓核包含蛋白聚糖在内的亲水性物质，能够将髓核物质包裹在水分子中可加固其物理稳定性，同时可产生一定的黏弹性将IVD内的轴向载荷转换成环向应力，并由周围的致密结缔组织包裹，保证其结构完整性，并提供了独特的承受高压的能力，极大地缓冲了对锥体、纤维环外组织及椎管内神经的挤压作用力。反复长期磨损及其他营养和环境因素的作用下，IVD开始退化，失去髓核的保护作用，一方面可促使锥体压缩载荷和这些力纵向传递到椎体的软骨终板上，并径向地传递到周围的肌肉组织、神经及血管；此外，不均匀的压力施加在椎骨上，相关结构变化可导致不稳定脊柱节段的适应性重塑，并形成骨赘骨刺，加重关节平面不平衡及损伤[3]。另一方面，髓核慢性脱出时可进一步刺激脊柱韧带的肥大和骨化，伴有毛细血管增生和结缔组织的浸润，促使椎盘软组织的空间更加狭窄[3]。

分子机制方面，既往研究提示，慢性脊髓压迫会导致慢性脊髓实质内血流量减少，同时慢性缺血会激活独特的免疫反应， M1和M2巨噬细胞会在压缩部位聚集，进一步促进内皮细胞功能障碍和血脊髓屏障破坏，同时部分脊髓缺血后可增加髓内压力，导致神经元丢失、脱髓鞘，进一步促使小血管扁平化而导致的血供减少。新近的观点认为，LDH的疾病进展是由环境和遗传风险因素共同决定的，不同基因型和环境因素之间的相互作用机制也不尽不同。MARTIROSYAN等[4]综述LDH的发生、发展与基因的关联，认为遗传因素和环境诱导的基因改变占IVD病因的75%，特别是在炎性、降解性、稳态和结构系统多样性。ZHANG等[5]对128例LDH患者和132例年龄和性别匹配的健康人进行对照研究，利用基质辅助激光解吸/电离及快速质谱分析技术分析了3个基因中9个多态位点，发现FasL-844C/T(rs763110)和CASP9-1263A>G(rs4645978)两个位点的多态性与LDH密切相关，揭示其可能是LDH的潜在风险位点。PERERA等[6]评估慢性腰痛患者LDH候选基因的单核苷酸多态性(SNV)与LDH严重程度的关系发现，蛋白聚糖(aggrecan，ACAN)的核苷酸多态性如rs2272023、rs35430524、rs2882676、rs2351491、rs938609、rs3825994、rs1042630、rs698621和rs3817428与LDH严重程度密切相关。此外，KURZAWSKI等[7]认为电压门控钠通道Nav1.7相关调控亚基(sodium channel,voltage-gated,type Ⅸ,alpha subunit,SCN9A)rs676030多态性可能与症状性椎间盘突出症患者的疼痛强度有关[7]； SADOWSKA等[8]发现白细胞介素(IL)-15基因表达与LDH患者年龄相关，α1干扰素(IFNα1)、IL-6、IL-15、瞬时受体电位阳离子通道6(transient receptor potential cation channel-6，TrPC6)基因表达与IVD退变等级密切相关[8]。

近年来，随着全基因组关联分析(genome-wide association study，GWAS)的推进，IDH疾病相关风险基因在已知风险基因中被重新鉴定出来。这些研究不仅提高对IVDD各种病理生理学要素背后众多遗传变异的认识，同时有助于推进IDH患者的个性化治疗和药物治疗策略。然而，整体而言，IDH病因复杂且具体机制尚未探明，进行进一步的生物信息学分析有助于对该疾病有更加深入的认识，以及进一步地解决多维度的科学问题。

1 资料与方法

1.1数据获取与预处理

1.1.1数据获取 经检索，本课题组从基因表达综合数据库(Gene Expression Omnibus，GEO，网址为http://ncbi.nlm.nih.gov/geo/)中下载得到GSE23130、GSE17077和GSE15227 3组芯片的原始数据[9]，共得到57个样品，其中对照组35个(轻度退化)，试验组22个(严重退化)[10]。所有的芯片组织来自LDH患者的IVD组织，组织分类均采用Thompson分级方法(严重退化的为Thompson等级Ⅳ和Ⅴ级；轻度退化的为Ⅰ～Ⅲ级)。同时，通过匹配3组芯片对应的平台注释文件(GPL1352[U133_X3P] Affymetrix Human X3P Array；Affymetrix,Santa Clara,CA)，对探针号进行基因名转换，有利于鉴别出合适的基因。

1.1.2数据预处理 本研究利用Affy函数读取芯片的原始表达值、稳健多阵列平均值法(robust multi-array average，RMA)对原始数据进行归一化、K-邻近值法(k-nearest-neighbor，KNN)填补缺失值等处理[11]。其中，RMA是一种用于从Affymetrix数据创建表达式矩阵的算法，主要原理包括：(1)对原始强度值背景校正、log2变换；(2)分位数归一化；(3)对归一化数据进行线性拟合，以获得每个阵列上的每个探针集的表达式值[11]。KNN算法是机器学习中常用的经典算法之一，其处理缺失值原理是：当第m行、第n列数据缺失，则选取和基因m相似、第n个样品数据完好的k个基因的第n个实验进行替代[12]。最后，对于多个探针对应一个基因名的表达矩阵，课题组选取平均值作为基因表达值，而一个探针对应多个基因名的数据将进行删减。

1.2方法

1.2.1差异基因分析 通过Combat函数对3组芯片的表达值进行异质性检验，消除因为不同实验而产生的数据偏移和批次效应[13]。同时，利用主成分分析(principal component analysis，PCA)对Combat前后的表达值进行可视化分析，观察数据整体分布情况。最后，利用LIMMA差异表达分析对试验组和对照组芯片表达值进行差异评估[11]。其中，删除批效应可以减少对数据本身的依赖性、降低分析的错误率，并提高可重复性；而LIMMA函数是经典的差异分析方法，基于线性模型来分析设计的实验表达差异，主要应用于芯片微阵列、RNA seq、定量PCR和许多蛋白质技术产生的高维度数据。

接下来，利用Benjamini-Hochberg方法对表达数据进行数据错误发现率(false discovery rate,FDR)矫正，并进一步推算基因表达倍数(genes expressional fold-changes，FC)[12]。本次研究设定|log2 FC|>2且FDR矫正后P<0.05作为筛选表达差异基因的标准。

1.2.2差异基因的生物学功能分析 利用DAVID在线分析软件(DAVID database；http://david.abcc.ncifcrf.gov/)进一步分析了基因本体(gene ontology，GO)功能和KEGG(kyoto encyclopedia of genes and genomes，KEGG)通路信息[14]。DAVID数据库是国际上权威的集基因注释、可视化和集成探索的功能数据库，有利于整合探索基因通路和功能富集程度。一般认为富集P<0.05的GO功能和KEGG通路为有统计学意义的基因富集功能和通路。

1.2.3基于差异基因的蛋白互作网络分析 导出差异基因后，课题组基于STRING在线分析软件 (V10.5;http://string-db.org/)对基因进行了蛋白互作网络(protein-protein interaction network，PPI)分析[15]。STRING平台整合了现有的实验结论、数据库、文本挖掘库、基因表达库、基因连接数据、基因融合及基因共表达数据库等信息，是目前权威的蛋白功能预测工具。此外，本次研究还通过Cytoscape软件(V3.5.1;http://cytoscape.org/)对蛋白互作网络进行可视化及关键基因筛选分析。

2 结 果

2.1差异表达分析结果 基于前述方法，3组芯片57个样品(对照组35个，试验组22个)进行RMA背景矫正、归一化处理、KNN算法处理缺失值、Combat函数去除批次效应后共筛选出149个差异表达基因同时满足|log2 FC|>2(即表达值差异大于4倍)及FDR矫正后的P<0.05的选择标准。利用PCA分析可视化3组芯片整合过程中去除批次效应前后的数据分组情况，见图1。图1A示3组芯片数据分布较散乱，同一组不同芯片之间数据距离比较大；而图1B示在去除批次效应后芯片间数据相对集中，有利于后续分析。图2的热图展现了差异前50 的基因表达值的情况，颜色深浅差异代表了表达值的差异。

A：去除异质性前；B：去除异质性后

图1 PCA显示Combat函数去除批次效应前后的3组芯片的整体数据分布

上栏颜色：分组信息；绿色：对照组；黄色：试验组；热图每一小格的颜色：表达值；蓝色：低表达；红色：高表达

图2 差异基因的表达热图

2.2差异基因生物学功能分析 通过DAVID生物在线分析平台，此次研究进一步分析了差异基因的最显著的GO功能和KEGG通路。其中，GO功能主要包括生物学进程(Biological processes,BP)、细胞组成(cellular components，CC)、分子功能(molecular functions，MF)3个主要方面。BP功能方面，筛选的差异基因主要与GO:0006614～SRP依赖的膜翻译蛋白(富集数目13，P=5.84×10-12)、GO:0000184～核转录的mRNA分解过程(富集数目13，P=9.89×10-11)及GO:0006413～翻译起始(富集数目13，P=5.17×10-10)等功能密切相关。CC功能方面，差异基因主要与GO:0070062～细胞外泌体(富集数目59，P=2.36×10-14)、GO:0022625～胞浆大核糖体亚基(富集数目9，P=3.97×10-8)、GO:0031012～细胞外基质(富集数目15，P=4.72×10-8)等密切相关。而在MF方面则主要与GO:0003735～核糖体的结构组成(富集数目13，P=9.49×10-8)、GO:0008137～NADH脱氢酶(泛醌)活性(富集数目22，P=7.96×10-5)及GO:0046961～质子转运ATP酶活性、旋转机制(富集数目4，P=0.005 4)等相关。差异基因本体功能的可视化结果见图3。

图形形状：代表不同的功能类型(圆形：基因的生物学进程，三角形：细胞学组成，正方形：分子功能)；图形大小：代表富集基因的数目；图形颜色：代表富集的统计值，转换为负log10的P值进行绘图

图3 基因本体功能分析结果

A：KEGG通路富集分析结果；图形形状：代表不同的功能类型(圆形：基因的生物学进程，三角形：细胞学组成，正方形：分子功能)；图形大小：代表富集基因的数目；图形颜色：代表富集的统计值，转换为负log10的P值进行绘图；B：差异基因进行蛋白互作网络分析的结果；基因名字体大小：代表网络中关键度的大小

图4 KEGG通路富集分析及蛋白互作网络

KEGG通路方面，主要与HSA03010：核糖体转录翻译活动(富集数目13，P=8.34×10-9)、HSA000 190：氧化磷酸化(富集数目12，P=7.30×10-8)、HSA045 12:ECM受体相互作用(富集数目10，P=1.20×10-4)、HSA04060：细胞因子-细胞因子受体相互作用(富集数目8，P=0.001 5)、HSA04370:血管内皮生长因子(VEGF)信号通路(富集数目7，P=0.041)、 HSA01100：代谢途径(富集数目20，P=0.033)及HSA04140:自噬的调节(富集数目3，P=0.048)等明显相关。KEGG通路富集结果见图4A。

2.3差异基因蛋白互作网络分析 根据STRING数据库和Cytoscape软件分析的网络关系图，课题组基于网络连接度大小筛选出网络中的关键基因。其中，在此PPI网络中，泛素-A 52残基核糖体蛋白融合产物1(ubiquitin A-52 residue ribosomal protein fusion product 1，UBA52；核心紧密度：0.40；相关度：22)、大核糖体蛋白-P0(large ribosomal protein P0，RPLP0；核心紧密度：0.40；相关度：18)、核糖体蛋白L3(ribosomal protein L3，RPL3；核心紧密度：0.38；相关度：17)、大核糖体蛋白P2(large ribosomal protein P2，RPLP2；核心紧密度：0.37；相关度：16)、核糖体蛋白L27(ribosomal protein L27，RPL27；核心紧密度：0.40；相关度：16)等基因在网络中具有较高的连接度，被认为是LDH疾病差异基因网络中的枢纽基因。差异基因蛋白互作网络见图4B。

3 讨 论

LDH是一种复杂的脊柱疾病，具有高致残率，给社会和家庭带来沉重的医疗负担。本文通过多芯片整合分析，通过Thompson分级将IVD分为试验组(Ⅳ～Ⅴ级)和对照组(Ⅰ～Ⅲ级)，利用RMA、KNN算法对原始芯片数据进行背景化和归一化处理、Combat函数去除3组芯片间的异质性、LIMMA函数筛选差异基因等一系列严格的数据处理流程。多芯片整合分析结果提示，HSA03010：核糖体转录翻译、HSA000 190：氧化磷酸化、HSA045 12:ECM受体相互作用等通路及UBA52、RPLP0、RPL3、RPLP2、RPL27等核糖体相关枢纽基因被认为是LDH疾病发生、发展的重要通路及调控因子。

本研究提示，核糖体活动及相关基因可能是LDH疾病进展的重要调控介质。然而，其具体机制尚不明确。KOBAYASHI等[16]通过构建机械压迫神经根的体内模型观察神经根压迫引起的轴突血流紊乱在腰椎运动神经元中的作用。此研究中，在压迫开始1周后脊髓的运动神经元出现中央染色质溶解，3周后出现腹角神经元凋亡；同时，作者观察到在此过程中核糖体参与了明显的病理活动，可能与LDH或腰椎管狭窄症引起的神经根压迫时产生的轴突反应延伸到脊髓内的运动神经元，而导致神经元染色质溶解和细胞死亡有关。同样，ALAMIN等[17]利用实时PCR检测靶向16S核糖体RNA基因，然后用扩增子测序分析切除的IVD组织，发现在慢性疼痛的LDH组织中核糖体基因表达明珠增高，而这些改变与研究前期发现的低毒力微生物感染无关。那么，为何核糖体在LDH疾病占据如此重的作用呢？SLOMNICKI等[18]新近的研究提示，在神经元损伤初期，为了保护核仁在细胞内的调控地位不破坏核仁完整性、细胞存活和对神经营养素的信号应答，神经元内将通过降低核糖体含量、RNA应激及核糖体本身的募集而减少蛋白质合成。然而，在损伤末期，当核糖体蛋白从树突的神经元中耗尽时则需要强健的核糖体来支持树突生长和维持蛋白质合成。因此，SLOMNICKI等[18]也认为RNA应激颗粒形成及神经元树突状损伤是神经退行性疾病的早期表现，通过干扰核糖体稳态可能启动或放大神经变性相关联的回路。POIROT等[19]认为随着神经元的延伸，核糖体蛋白将形成复杂的互作网络，并且通过微界面进行进一步的通信。在此过程中，在类似于“分子突触”相互作用的神经元中，核糖体蛋白形成了一个与简单的脑类似的分子网络，部分核糖体表面蛋白被认为是支配功能性核糖体的调控位点，通过“蛋白间交流”连接成合适的循环网络来协调复杂的大分子运动和翻译过程，并为信息传递和处理开辟了新的视角[19]。在枢纽基因中，RPL3属于核糖体相关的功能大亚基核心蛋白，主要与RNA-蛋白质相互作用使rRNA形成其功能上正确的构象，从而参与核糖体不同功能区域之间及核糖体和细胞因子之间的通信等密切相关[20]。编码基因UBA52不仅是泛素库的泛素供体，而且也是核糖体蛋白复合物的调节器。其中，UBA52 N端含有泛素的融合蛋白和C末端的核糖体蛋白L40 (ribosomal protein L40，RPL40)结合参与了核糖体泛素化修饰[21]。然而，ARTERO-CASTRO等[22]认为RPLP0和RPLP2不仅与核糖体活动相关还可能参与活性氧(reactive oxygen species，ROS)积累及对应的丝裂原活化蛋白激酶1/细胞外信号调节激酶1/2(MAPK1/Erk2)信号通路激活，且与内质网应激和细胞的自噬等密切相关。核糖体蛋白基因RPL27进化过程中保存良好，相对保守且与维持核糖体蛋白的特异性3D结构相关[22]。

在氧化应激方面，XIAO等[23]认为氧化应激反应在LDH疾病进展中发挥重要作用，调控细胞外基质蛋白的再生，增加突出IVD的新生胶原和聚集蛋白聚糖含量，而抑制组织氧化磷酸化水平可明显改善LDH，从而促进IVD突出后IVD高度的恢复。CHANG等[24]认为术前术后降低体内氧化应激水平将明显改善LDH患者的手术效果和病情。在ECM相互作用方面，GUTERL等[25]认为ECM参与了LDH的纤维化基质形成，且参与椎体、IVD组织三维结构重构，改变相应的机械性能，可能是LDH疼痛的根本原因。同样，HIROSE等[26]的研究也证实IVD的ECM代谢及重塑在LDH的病因和发病机制中起着至关重要的作用。

基于多芯片整合分析，本研究发现无论在蛋白互作网络还是通路富集分析结果中均提示核糖体活动或相关枢纽基因(包括UBA52、RPLP0、RPL3、RPLP2、RPL27等核糖体相关基因)与LDH疾病的进展密切相关，机制方面拟参与疾病所致的神经损伤的病理过程。此外，氧化应激和ECM相互作用通路也同样显著富集于LDH疾病的差异表达基因中，可能与聚集蛋白聚糖代谢、胶原形成及纤维化重构密切相关。然而，本研究尚且缺乏实验进一步验证且LDH发病机制复杂，在具体的机制方面也将有待于后续的深入研究。