多芯片整合分析腰椎間盤突出癥發(fā)生相關的分子通路及關鍵基因研究*

劉權祥，陳 錢，齊明宇，莊新明，郭建平△，李增新，程 維，代起志

(1.北華大學附屬醫(yī)院骨外二科，吉林吉林 132011；2.吉林省吉林市中心醫(yī)院內(nèi)分泌骨代謝科 132011；3.吉林大學第一醫(yī)院脊柱外科，長春 132021)

腰椎間盤突出癥(lumbar disc herniation，LDH)是一組以椎間盤(intervertebral disc，IVD)退行性變，間盤內(nèi)纖維環(huán)破裂，髓核在切應力作用下向內(nèi)突，進一步壓迫神經(jīng)根或馬尾神經(jīng)而引起的以神經(jīng)癥狀為主的臨床常見病和多發(fā)病[1]。新近的流行病學研究提示，在亞洲人群中，LDH的年總發(fā)病率為15 877/10萬(其中男性為13 181/10萬；女性為18 588/10萬；男女比為1.00∶1.41)。此外，每5年患者的發(fā)病率和年治療費用分別增加7.6%和14.7%。年齡矯正后的報道提示，LDH發(fā)病率隨年齡增長而增加，75～79歲人群發(fā)病率高達42.6%。其中，年齡大于65歲患者占31.0%，且其所占醫(yī)療費用高達40.1%[2-3]。

正常的IVD內(nèi)，由于髓核包含蛋白聚糖在內(nèi)的親水性物質(zhì)，能夠?qū)⑺韬宋镔|(zhì)包裹在水分子中可加固其物理穩(wěn)定性，同時可產(chǎn)生一定的黏彈性將IVD內(nèi)的軸向載荷轉(zhuǎn)換成環(huán)向應力，并由周圍的致密結(jié)締組織包裹，保證其結(jié)構完整性，并提供了獨特的承受高壓的能力，極大地緩沖了對錐體、纖維環(huán)外組織及椎管內(nèi)神經(jīng)的擠壓作用力。反復長期磨損及其他營養(yǎng)和環(huán)境因素的作用下，IVD開始退化，失去髓核的保護作用，一方面可促使錐體壓縮載荷和這些力縱向傳遞到椎體的軟骨終板上，并徑向地傳遞到周圍的肌肉組織、神經(jīng)及血管；此外，不均勻的壓力施加在椎骨上，相關結(jié)構變化可導致不穩(wěn)定脊柱節(jié)段的適應性重塑，并形成骨贅骨刺，加重關節(jié)平面不平衡及損傷[3]。另一方面，髓核慢性脫出時可進一步刺激脊柱韌帶的肥大和骨化，伴有毛細血管增生和結(jié)締組織的浸潤，促使椎盤軟組織的空間更加狹窄[3]。

分子機制方面，既往研究提示，慢性脊髓壓迫會導致慢性脊髓實質(zhì)內(nèi)血流量減少，同時慢性缺血會激活獨特的免疫反應， M1和M2巨噬細胞會在壓縮部位聚集，進一步促進內(nèi)皮細胞功能障礙和血脊髓屏障破壞，同時部分脊髓缺血后可增加髓內(nèi)壓力，導致神經(jīng)元丟失、脫髓鞘，進一步促使小血管扁平化而導致的血供減少。新近的觀點認為，LDH的疾病進展是由環(huán)境和遺傳風險因素共同決定的，不同基因型和環(huán)境因素之間的相互作用機制也不盡不同。MARTIROSYAN等[4]綜述LDH的發(fā)生、發(fā)展與基因的關聯(lián)，認為遺傳因素和環(huán)境誘導的基因改變占IVD病因的75%，特別是在炎性、降解性、穩(wěn)態(tài)和結(jié)構系統(tǒng)多樣性。ZHANG等[5]對128例LDH患者和132例年齡和性別匹配的健康人進行對照研究，利用基質(zhì)輔助激光解吸/電離及快速質(zhì)譜分析技術分析了3個基因中9個多態(tài)位點，發(fā)現(xiàn)FasL-844C/T(rs763110)和CASP9-1263A>G(rs4645978)兩個位點的多態(tài)性與LDH密切相關，揭示其可能是LDH的潛在風險位點。PERERA等[6]評估慢性腰痛患者LDH候選基因的單核苷酸多態(tài)性(SNV)與LDH嚴重程度的關系發(fā)現(xiàn)，蛋白聚糖(aggrecan，ACAN)的核苷酸多態(tài)性如rs2272023、rs35430524、rs2882676、rs2351491、rs938609、rs3825994、rs1042630、rs698621和rs3817428與LDH嚴重程度密切相關。此外，KURZAWSKI等[7]認為電壓門控鈉通道Nav1.7相關調(diào)控亞基(sodium channel,voltage-gated,type Ⅸ,alpha subunit,SCN9A)rs676030多態(tài)性可能與癥狀性椎間盤突出癥患者的疼痛強度有關[7]； SADOWSKA等[8]發(fā)現(xiàn)白細胞介素(IL)-15基因表達與LDH患者年齡相關，α1干擾素(IFNα1)、IL-6、IL-15、瞬時受體電位陽離子通道6(transient receptor potential cation channel-6，TrPC6)基因表達與IVD退變等級密切相關[8]。

近年來，隨著全基因組關聯(lián)分析(genome-wide association study，GWAS)的推進，IDH疾病相關風險基因在已知風險基因中被重新鑒定出來。這些研究不僅提高對IVDD各種病理生理學要素背后眾多遺傳變異的認識，同時有助于推進IDH患者的個性化治療和藥物治療策略。然而，整體而言，IDH病因復雜且具體機制尚未探明，進行進一步的生物信息學分析有助于對該疾病有更加深入的認識，以及進一步地解決多維度的科學問題。

1 資料與方法

1.1數(shù)據(jù)獲取與預處理

1.1.1數(shù)據(jù)獲取 經(jīng)檢索，本課題組從基因表達綜合數(shù)據(jù)庫(Gene Expression Omnibus，GEO，網(wǎng)址為http://ncbi.nlm.nih.gov/geo/)中下載得到GSE23130、GSE17077和GSE15227 3組芯片的原始數(shù)據(jù)[9]，共得到57個樣品，其中對照組35個(輕度退化)，試驗組22個(嚴重退化)[10]。所有的芯片組織來自LDH患者的IVD組織，組織分類均采用Thompson分級方法(嚴重退化的為Thompson等級Ⅳ和Ⅴ級；輕度退化的為Ⅰ～Ⅲ級)。同時，通過匹配3組芯片對應的平臺注釋文件(GPL1352[U133_X3P] Affymetrix Human X3P Array；Affymetrix,Santa Clara,CA)，對探針號進行基因名轉(zhuǎn)換，有利于鑒別出合適的基因。

1.1.2數(shù)據(jù)預處理 本研究利用Affy函數(shù)讀取芯片的原始表達值、穩(wěn)健多陣列平均值法(robust multi-array average，RMA)對原始數(shù)據(jù)進行歸一化、K-鄰近值法(k-nearest-neighbor，KNN)填補缺失值等處理[11]。其中，RMA是一種用于從Affymetrix數(shù)據(jù)創(chuàng)建表達式矩陣的算法，主要原理包括：(1)對原始強度值背景校正、log2變換；(2)分位數(shù)歸一化；(3)對歸一化數(shù)據(jù)進行線性擬合，以獲得每個陣列上的每個探針集的表達式值[11]。KNN算法是機器學習中常用的經(jīng)典算法之一，其處理缺失值原理是：當?shù)趍行、第n列數(shù)據(jù)缺失，則選取和基因m相似、第n個樣品數(shù)據(jù)完好的k個基因的第n個實驗進行替代[12]。最后，對于多個探針對應一個基因名的表達矩陣，課題組選取平均值作為基因表達值，而一個探針對應多個基因名的數(shù)據(jù)將進行刪減。

1.2方法

1.2.1差異基因分析 通過Combat函數(shù)對3組芯片的表達值進行異質(zhì)性檢驗，消除因為不同實驗而產(chǎn)生的數(shù)據(jù)偏移和批次效應[13]。同時，利用主成分分析(principal component analysis，PCA)對Combat前后的表達值進行可視化分析，觀察數(shù)據(jù)整體分布情況。最后，利用LIMMA差異表達分析對試驗組和對照組芯片表達值進行差異評估[11]。其中，刪除批效應可以減少對數(shù)據(jù)本身的依賴性、降低分析的錯誤率，并提高可重復性；而LIMMA函數(shù)是經(jīng)典的差異分析方法，基于線性模型來分析設計的實驗表達差異，主要應用于芯片微陣列、RNA seq、定量PCR和許多蛋白質(zhì)技術產(chǎn)生的高維度數(shù)據(jù)。

接下來，利用Benjamini-Hochberg方法對表達數(shù)據(jù)進行數(shù)據(jù)錯誤發(fā)現(xiàn)率(false discovery rate,FDR)矯正，并進一步推算基因表達倍數(shù)(genes expressional fold-changes，F(xiàn)C)[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數(shù)據(jù)庫是國際上權威的集基因注釋、可視化和集成探索的功能數(shù)據(jù)庫，有利于整合探索基因通路和功能富集程度。一般認為富集P<0.05的GO功能和KEGG通路為有統(tǒng)計學意義的基因富集功能和通路。

1.2.3基于差異基因的蛋白互作網(wǎng)絡分析 導出差異基因后，課題組基于STRING在線分析軟件 (V10.5;http://string-db.org/)對基因進行了蛋白互作網(wǎng)絡(protein-protein interaction network，PPI)分析[15]。STRING平臺整合了現(xiàn)有的實驗結(jié)論、數(shù)據(jù)庫、文本挖掘庫、基因表達庫、基因連接數(shù)據(jù)、基因融合及基因共表達數(shù)據(jù)庫等信息，是目前權威的蛋白功能預測工具。此外，本次研究還通過Cytoscape軟件(V3.5.1;http://cytoscape.org/)對蛋白互作網(wǎng)絡進行可視化及關鍵基因篩選分析。

2 結(jié) 果

2.1差異表達分析結(jié)果 基于前述方法，3組芯片57個樣品(對照組35個，試驗組22個)進行RMA背景矯正、歸一化處理、KNN算法處理缺失值、Combat函數(shù)去除批次效應后共篩選出149個差異表達基因同時滿足|log2 FC|>2(即表達值差異大于4倍)及FDR矯正后的P<0.05的選擇標準。利用PCA分析可視化3組芯片整合過程中去除批次效應前后的數(shù)據(jù)分組情況，見圖1。圖1A示3組芯片數(shù)據(jù)分布較散亂，同一組不同芯片之間數(shù)據(jù)距離比較大；而圖1B示在去除批次效應后芯片間數(shù)據(jù)相對集中，有利于后續(xù)分析。圖2的熱圖展現(xiàn)了差異前50 的基因表達值的情況，顏色深淺差異代表了表達值的差異。

A：去除異質(zhì)性前；B：去除異質(zhì)性后

圖1 PCA顯示Combat函數(shù)去除批次效應前后的3組芯片的整體數(shù)據(jù)分布

上欄顏色：分組信息；綠色：對照組；黃色：試驗組；熱圖每一小格的顏色：表達值；藍色：低表達；紅色：高表達

圖2 差異基因的表達熱圖

2.2差異基因生物學功能分析 通過DAVID生物在線分析平臺，此次研究進一步分析了差異基因的最顯著的GO功能和KEGG通路。其中，GO功能主要包括生物學進程(Biological processes,BP)、細胞組成(cellular components，CC)、分子功能(molecular functions，MF)3個主要方面。BP功能方面，篩選的差異基因主要與GO:0006614～SRP依賴的膜翻譯蛋白(富集數(shù)目13，P=5.84×10-12)、GO:0000184～核轉(zhuǎn)錄的mRNA分解過程(富集數(shù)目13，P=9.89×10-11)及GO:0006413～翻譯起始(富集數(shù)目13，P=5.17×10-10)等功能密切相關。CC功能方面，差異基因主要與GO:0070062～細胞外泌體(富集數(shù)目59，P=2.36×10-14)、GO:0022625～胞漿大核糖體亞基(富集數(shù)目9，P=3.97×10-8)、GO:0031012～細胞外基質(zhì)(富集數(shù)目15，P=4.72×10-8)等密切相關。而在MF方面則主要與GO:0003735～核糖體的結(jié)構組成(富集數(shù)目13，P=9.49×10-8)、GO:0008137～NADH脫氫酶(泛醌)活性(富集數(shù)目22，P=7.96×10-5)及GO:0046961～質(zhì)子轉(zhuǎn)運ATP酶活性、旋轉(zhuǎn)機制(富集數(shù)目4，P=0.005 4)等相關。差異基因本體功能的可視化結(jié)果見圖3。

圖形形狀：代表不同的功能類型(圓形：基因的生物學進程，三角形：細胞學組成，正方形：分子功能)；圖形大小：代表富集基因的數(shù)目；圖形顏色：代表富集的統(tǒng)計值，轉(zhuǎn)換為負log10的P值進行繪圖

圖3 基因本體功能分析結(jié)果

A：KEGG通路富集分析結(jié)果；圖形形狀：代表不同的功能類型(圓形：基因的生物學進程，三角形：細胞學組成，正方形：分子功能)；圖形大小：代表富集基因的數(shù)目；圖形顏色：代表富集的統(tǒng)計值，轉(zhuǎn)換為負log10的P值進行繪圖；B：差異基因進行蛋白互作網(wǎng)絡分析的結(jié)果；基因名字體大小：代表網(wǎng)絡中關鍵度的大小

圖4 KEGG通路富集分析及蛋白互作網(wǎng)絡

KEGG通路方面，主要與HSA03010：核糖體轉(zhuǎn)錄翻譯活動(富集數(shù)目13，P=8.34×10-9)、HSA000 190：氧化磷酸化(富集數(shù)目12，P=7.30×10-8)、HSA045 12:ECM受體相互作用(富集數(shù)目10，P=1.20×10-4)、HSA04060：細胞因子-細胞因子受體相互作用(富集數(shù)目8，P=0.001 5)、HSA04370:血管內(nèi)皮生長因子(VEGF)信號通路(富集數(shù)目7，P=0.041)、 HSA01100：代謝途徑(富集數(shù)目20，P=0.033)及HSA04140:自噬的調(diào)節(jié)(富集數(shù)目3，P=0.048)等明顯相關。KEGG通路富集結(jié)果見圖4A。

2.3差異基因蛋白互作網(wǎng)絡分析 根據(jù)STRING數(shù)據(jù)庫和Cytoscape軟件分析的網(wǎng)絡關系圖，課題組基于網(wǎng)絡連接度大小篩選出網(wǎng)絡中的關鍵基因。其中，在此PPI網(wǎng)絡中，泛素-A 52殘基核糖體蛋白融合產(chǎn)物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)等基因在網(wǎng)絡中具有較高的連接度，被認為是LDH疾病差異基因網(wǎng)絡中的樞紐基因。差異基因蛋白互作網(wǎng)絡見圖4B。

3 討 論

LDH是一種復雜的脊柱疾病，具有高致殘率，給社會和家庭帶來沉重的醫(yī)療負擔。本文通過多芯片整合分析，通過Thompson分級將IVD分為試驗組(Ⅳ～Ⅴ級)和對照組(Ⅰ～Ⅲ級)，利用RMA、KNN算法對原始芯片數(shù)據(jù)進行背景化和歸一化處理、Combat函數(shù)去除3組芯片間的異質(zhì)性、LIMMA函數(shù)篩選差異基因等一系列嚴格的數(shù)據(jù)處理流程。多芯片整合分析結(jié)果提示，HSA03010：核糖體轉(zhuǎn)錄翻譯、HSA000 190：氧化磷酸化、HSA045 12:ECM受體相互作用等通路及UBA52、RPLP0、RPL3、RPLP2、RPL27等核糖體相關樞紐基因被認為是LDH疾病發(fā)生、發(fā)展的重要通路及調(diào)控因子。

本研究提示，核糖體活動及相關基因可能是LDH疾病進展的重要調(diào)控介質(zhì)。然而，其具體機制尚不明確。KOBAYASHI等[16]通過構建機械壓迫神經(jīng)根的體內(nèi)模型觀察神經(jīng)根壓迫引起的軸突血流紊亂在腰椎運動神經(jīng)元中的作用。此研究中，在壓迫開始1周后脊髓的運動神經(jīng)元出現(xiàn)中央染色質(zhì)溶解，3周后出現(xiàn)腹角神經(jīng)元凋亡；同時，作者觀察到在此過程中核糖體參與了明顯的病理活動，可能與LDH或腰椎管狹窄癥引起的神經(jīng)根壓迫時產(chǎn)生的軸突反應延伸到脊髓內(nèi)的運動神經(jīng)元，而導致神經(jīng)元染色質(zhì)溶解和細胞死亡有關。同樣，ALAMIN等[17]利用實時PCR檢測靶向16S核糖體RNA基因，然后用擴增子測序分析切除的IVD組織，發(fā)現(xiàn)在慢性疼痛的LDH組織中核糖體基因表達明珠增高，而這些改變與研究前期發(fā)現(xiàn)的低毒力微生物感染無關。那么，為何核糖體在LDH疾病占據(jù)如此重的作用呢？SLOMNICKI等[18]新近的研究提示，在神經(jīng)元損傷初期，為了保護核仁在細胞內(nèi)的調(diào)控地位不破壞核仁完整性、細胞存活和對神經(jīng)營養(yǎng)素的信號應答，神經(jīng)元內(nèi)將通過降低核糖體含量、RNA應激及核糖體本身的募集而減少蛋白質(zhì)合成。然而，在損傷末期，當核糖體蛋白從樹突的神經(jīng)元中耗盡時則需要強健的核糖體來支持樹突生長和維持蛋白質(zhì)合成。因此，SLOMNICKI等[18]也認為RNA應激顆粒形成及神經(jīng)元樹突狀損傷是神經(jīng)退行性疾病的早期表現(xiàn)，通過干擾核糖體穩(wěn)態(tài)可能啟動或放大神經(jīng)變性相關聯(lián)的回路。POIROT等[19]認為隨著神經(jīng)元的延伸，核糖體蛋白將形成復雜的互作網(wǎng)絡，并且通過微界面進行進一步的通信。在此過程中，在類似于“分子突觸”相互作用的神經(jīng)元中，核糖體蛋白形成了一個與簡單的腦類似的分子網(wǎng)絡，部分核糖體表面蛋白被認為是支配功能性核糖體的調(diào)控位點，通過“蛋白間交流”連接成合適的循環(huán)網(wǎng)絡來協(xié)調(diào)復雜的大分子運動和翻譯過程，并為信息傳遞和處理開辟了新的視角[19]。在樞紐基因中，RPL3屬于核糖體相關的功能大亞基核心蛋白，主要與RNA-蛋白質(zhì)相互作用使rRNA形成其功能上正確的構象，從而參與核糖體不同功能區(qū)域之間及核糖體和細胞因子之間的通信等密切相關[20]。編碼基因UBA52不僅是泛素庫的泛素供體，而且也是核糖體蛋白復合物的調(diào)節(jié)器。其中，UBA52 N端含有泛素的融合蛋白和C末端的核糖體蛋白L40 (ribosomal protein L40，RPL40)結(jié)合參與了核糖體泛素化修飾[21]。然而，ARTERO-CASTRO等[22]認為RPLP0和RPLP2不僅與核糖體活動相關還可能參與活性氧(reactive oxygen species，ROS)積累及對應的絲裂原活化蛋白激酶1/細胞外信號調(diào)節(jié)激酶1/2(MAPK1/Erk2)信號通路激活，且與內(nèi)質(zhì)網(wǎng)應激和細胞的自噬等密切相關。核糖體蛋白基因RPL27進化過程中保存良好，相對保守且與維持核糖體蛋白的特異性3D結(jié)構相關[22]。

在氧化應激方面，XIAO等[23]認為氧化應激反應在LDH疾病進展中發(fā)揮重要作用，調(diào)控細胞外基質(zhì)蛋白的再生，增加突出IVD的新生膠原和聚集蛋白聚糖含量，而抑制組織氧化磷酸化水平可明顯改善LDH，從而促進IVD突出后IVD高度的恢復。CHANG等[24]認為術前術后降低體內(nèi)氧化應激水平將明顯改善LDH患者的手術效果和病情。在ECM相互作用方面，GUTERL等[25]認為ECM參與了LDH的纖維化基質(zhì)形成，且參與椎體、IVD組織三維結(jié)構重構，改變相應的機械性能，可能是LDH疼痛的根本原因。同樣，HIROSE等[26]的研究也證實IVD的ECM代謝及重塑在LDH的病因和發(fā)病機制中起著至關重要的作用。

基于多芯片整合分析，本研究發(fā)現(xiàn)無論在蛋白互作網(wǎng)絡還是通路富集分析結(jié)果中均提示核糖體活動或相關樞紐基因(包括UBA52、RPLP0、RPL3、RPLP2、RPL27等核糖體相關基因)與LDH疾病的進展密切相關，機制方面擬參與疾病所致的神經(jīng)損傷的病理過程。此外，氧化應激和ECM相互作用通路也同樣顯著富集于LDH疾病的差異表達基因中，可能與聚集蛋白聚糖代謝、膠原形成及纖維化重構密切相關。然而，本研究尚且缺乏實驗進一步驗證且LDH發(fā)病機制復雜，在具體的機制方面也將有待于后續(xù)的深入研究。