通過構建競爭性內源RNA網絡識別潛在的椎間盤疾病相關circRNA

陳佳新,宋 超,張玥欣,尹佳琦,李姍珊,李春權,張 建

(哈爾濱醫(yī)科大學 大慶校區(qū)信息學院，黑龍江 大慶 163319)

環(huán)狀RNA(circRNA)是一類特殊的非編碼RNA分子，呈封閉環(huán)狀結構，不受RNA外切酶影響，相對于線性RNA分子更加穩(wěn)定不容易降解[1]。近年在功能上的研究表明，circRNA分子含有大量的microRNA(miRNA)結合位點，在細胞中起到miRNA海綿的作用從而消除或減弱miRNA對其靶基因的抑制，提高靶基因的表達水平，這種調控模式被稱為競爭性內源RNA(ceRNA)[2]。ceRNA是一種全新的基因表達調控模式，涉及許多RNA分子，包括mRNA、circRNA和miRNA等。有研究表明，ceRNA在腫瘤的發(fā)生發(fā)展中起著重要作用。

椎間盤變性是腰背痛和坐骨神經痛最常見的原因，而椎間盤退化的主要原因是年齡的增長，然而近年來椎間盤退化的發(fā)病群體卻有年輕化的趨勢[3]。目前治療椎間盤退化的方法包括藥物治療、物理治療和外科手術，但這些都不是最理想的。因此，我們在基因分子的層面上做出假設：在椎間盤變性疾病中，circRNA和mRNA之間存在潛在的ceRNA關系，我們可以通過調節(jié)RNA之間的競爭性關系來準確高效地預診治療椎間盤變性。

通過miRNA橋梁構建了一個包含mRNA與circRNA兩類節(jié)點的ceRNA網絡，利用聚類分析、模塊挖掘等生物信息學方法進行網絡分析并且從網絡中挖掘出與疾病相關的重要模塊，利用Gene Ontology(GO)基因功能注釋、KEGG通路富集等方法對circRNA的功能進行了預測。為了探索重要模塊中的蛋白編碼基因與circRNA的分子功能，在NCBI PubMed中搜集大量文獻來證明重要模塊中基因的功能及作用。在ceRNA表達調控關系的基礎上推測出circRNA潛在的生物學意義和功能，從根本上預測椎間盤變性中基因表達水平和分子調控機制，為椎間盤變性相關藥物的開發(fā)和臨床預診治療提供了重要的分子標志物。

1 實驗材料與方法

1.1 數(shù)據來源及處理

在NCBI GEO數(shù)據庫中下載了人類椎間盤變性疾病的芯片表達譜數(shù)據(GSE67567)，其中包含同樣本的mRNA和circRNA的表達譜。總共10個樣本，包括五個正常樣本和五個髓核退化疾病樣本。其中mRNA表達譜來自GPL15314 Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version)平臺，circRNA表達譜來自GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1平臺。在GENCODE數(shù)據庫中下載人類編碼蛋白轉錄序列Fasta格式文件，并且在Circbase數(shù)據庫中下載人類circRNA的序列文件。運用blast軟件對mRNA序列與芯片探針序列進行了比對，得到探針-mRNA的對應關系后進行平臺轉換后得到基因的表達譜。與此同時，對circRNA的探針表達譜進行平臺轉換得到circRNA表達譜。

1.2 構造競爭性內源RNA(ceRNA)網絡

在mirbase數(shù)據庫上下載人類成熟的miRNA序列，利用miRanda[4]軟件預測miRNA-mRNA,與circRNA-miRNA的的靶向關系。Score選取默認值>140，Energy選取更為嚴謹?shù)挠蛑?0。

1.2.1 ceRNA關系預測

這兩類關系對中的miRNA進行超幾何檢驗。具體計算公式如下：

其中m代表人類miRNA總數(shù)，t表示與mRNA相互作用的miRNA的數(shù)目，n表示與circRNA相互作用的miRNA的數(shù)目，r表示mRNA與circRNA之間共享的miRNA數(shù)目。通過對mRNA與miRNA的靶向關系對和circRNA與miRNA的關系對做超幾何檢驗，我們得到ceRNA對(p <0.05)。

研究發(fā)現(xiàn) circRNA的過量表達會降低游離miRNA分子的數(shù)目，進而導致miRNA下游靶基因的高表達[5]。共表達是一種發(fā)現(xiàn)調控關系的常用方法，我們在正常樣本和疾病樣本中分別計算mRNA與circRNA的pearson相關性，得到2 539 352個邊對關系和對應的Pearson相關系數(shù)cor值?？紤]到算法的統(tǒng)計學意義和網絡的規(guī)模等因素，我們認為mRNA與circRNA在疾病與正常樣本中變化程度大于0.5的關系對符合共表達關系。數(shù)學公式如下:

|cor(case)-cor(control)|>0.5

根據最新研究結果， circRNA和mRNA之間共表達的增加或缺失可能是導致疾病的原因，所以我們同樣假設這種競爭關系可以導致椎間盤變性疾病的發(fā)生[2]。

為了使網絡更加趨向于疾病特異性和相關性，既考慮到了通過Pearson相關性構成的共表達網絡，又兼顧到了miRanda預測的miRNA潛在靶點構成的ceRNA網絡，選取他們關系對的交集構成mRNA與circRNA的ceRNA二步網絡進行網絡功能與拓撲分析，為后續(xù)的研究做基礎。

1.3 模塊挖掘與功能注釋通路富集

運用cytoscape 3.2.0插件jActiveModules 進行模塊挖掘，參數(shù)選取默認值(Number of Modules:5, Overlap Threshold:0.8)。并且使用了cytoscape3.2.0插件 ClueGO對模塊中包含的基因GO功能注釋和KEGG通路富集。

方法流程圖見圖1.

圖1 方法流程圖Fig.1 Method flowchart

2 結果分析

2.1 表達譜處理和分析

運用blast軟件對mRNA序列和芯片探針序列進行了比對，進行平臺轉換后我們得到含有11 502個編碼基因的表達譜。同時對circRNA的探針表達譜進行平臺轉換得到含有2 793個circRNA的表達譜。我們通過聚類分別對mRNA的表達譜和circRNA的表達譜進行進一步分析，結果表明兩個表達譜在疾病和正常樣本中的差異性明顯，樣本數(shù)據較好，見圖2。

圖2 表達譜聚類圖Fig.2 Clustering map of expression profile

注：最下方淺灰色條代表正常樣本，深灰色條代表髓核退化樣本。顏色由淺入深代表表達量的標準化數(shù)值大小。

2.2 網絡分析結果

通過以上算法生成了一個由mRNA與circRNA兩類節(jié)點組成的二步網絡。網絡包含2 349個mRNA，264個circRNA，5 705個ceRNA關系對。其中有兩類節(jié)點，黃色圓形為mRNA，藍綠色菱形為circRNA。我們構建的網絡服從冪率分布(見圖3(b))。

對網絡進行拓撲屬性分析，其中包括度分布、網絡介數(shù)、最短路徑和網絡親密度。由圖3(a)可見網絡服從冪率分布，說明網絡中存在一小部分的基因或者circRNA連接著大部分節(jié)點，網絡連通性較強，同時也說明這小部分節(jié)點在網絡調控中起重要作用。

圖3 網絡分析圖Fig.3 Network analysis chart

2.3 識別ceRNA網絡重要模塊

使用cytoscape3.2.0的jActiveModules插件進行模塊挖掘，發(fā)現(xiàn)了五個circRNA與mRNA連接緊密的模塊。模塊的規(guī)模都很相似，模塊節(jié)點數(shù)為32-70，邊為37-132。利用Cytoscape的ClueGO插件分別對每一個模塊中的所有基因進行功能注釋和通路富集，選取了P值<0.05的GOTerm和Pathway。其中主要包含了幾類生物學過程：相關炎癥反應、骨發(fā)育和損傷、蛋白質結構和穩(wěn)定性等等。

2.4 模塊展示與分析

2.4.1 模塊一

模塊一(見圖4(a))包括36個節(jié)點、45條邊，包含了19個mRNA和7個circRNA。其中編碼蛋白的基因富集到蛋白質穩(wěn)定期(見圖4(c))。通過對椎間盤蛋白質組學差異分析，我們發(fā)現(xiàn)正常人和腰椎間盤退變患者的蛋白質種類存在顯著差異，腰椎間盤退變患者的蛋白質含量也隨著病變程度的加深而逐漸減少[6]。同時我們也發(fā)現(xiàn)了內源性凋亡信號通路的調控，研究表明人椎間盤髓核細胞的凋亡具有內在性和外源性[7]。同時椎間盤變性與骨的發(fā)育相關，此模塊也富集到了骨發(fā)育功能[8-9]。

因PAPSS2是AR brachyolmia的疾病基因，并且PAPSS2突變可以產生從brachyolmia到脊椎干骺發(fā)育不良的表型分級。PAPS是體內活性硫酸根供體，它的合成需要在PAPS合成酶(PAPss)的催化下進行，PAPS合成酶2(PAPSS2)的突變與大骨節(jié)病、骨關節(jié)炎等多種骨骼疾病的發(fā)生發(fā)展有關[10-16]。在此模塊中，GARS與PAPSS2同時富集到骨發(fā)育的生物學過程?；赾eRNA理論，推斷GARS基因與PAPSS2及他們靶向的circRNAs(has_circ_0013204，has_circ_0000077，has_circ_0073027，has_circ_0001234，has_circ_0007311)，可能通過ceRNA調控關系誘發(fā)椎間盤變性等骨疾病[17]。發(fā)現(xiàn)模塊二中的PAPSS2與has_circ_0000077的競爭性調控機制改變了共表達的強度從而影響了基因的表達和相應蛋白的功能。在正常樣本中他們的共表達相關值為-0.633 56，在疾病狀態(tài)下他們的表達相關值為0.642 102，在正常和疾病中體現(xiàn)了明顯的共表達差異，并且由于競爭性關系從負相關變?yōu)檎嚓P也暗示了競爭性RNA調控機制的表達趨勢。并且他們的超幾何P值統(tǒng)計學顯著為1.02×10-6，F(xiàn)DR矯正值為7.913×10-3，證明了通過hsa-let-7a-2-3p[18]介導的PAPSS2與has_circ_0000077之間緊密的ceRNA關系(見圖4(b))。其中得到PAPSS2 基因的Foldchange值為1.899 515 547在疾病中顯著上調。

注：圖(a)中文獻證明已知疾病相關基因用紅色圓形節(jié)點表示；把富集到相同生物學過程的基因用黃色圓形節(jié)點表示；把黃色和紅色節(jié)點的一步鄰居circRNAs用紫色菱形節(jié)點表示并且把它們共同的一步鄰居用綠色菱形節(jié)點表示；其他的基因用藍色圓形節(jié)點表示，circRNA用橙色菱形節(jié)點表示。圖4(b)舉例一個ceRNA關系對的MiRanda靶點展示。圖4(c)模塊中顏色代表P值大小。彩圖見電子版(http://swxxx.alljournals.cn/ch/index.aspx.2018年第4期)。

2.4.2 模塊二

模塊二(見圖5(a))包括41個節(jié)點、44條邊，包含了 22個mRNA和19個circRNA。我們發(fā)現(xiàn)了” 轉化生長因子β受體信號通路、淋巴細胞分化正調節(jié)、白細胞分化正調節(jié)、T細胞選擇陽性、跨質膜導入、白細胞活化負調節(jié)、細胞酰胺代謝過程負調節(jié)等生物學過程(P值<0.05)。轉化生長因子β1可以改善退變的椎間盤細胞活性，從而刺激膠原基因表達來延緩椎間盤組織的退變，因此對退變早期的椎間盤尤其是髓核具有修復功能，可逆轉椎間盤的退變[19]。也有實驗證實腺病毒載體介導人轉化生長因子β1基因轉染退行性變的椎間盤細胞能夠有效提高髓核細胞內轉化生長因子β1含量，阻逆椎間盤退行性變[20-22]。

同時椎間盤變性還伴隨著相關的炎癥和免疫反應，許多文獻都指出在腰椎椎間盤突出癥中，髓核邊緣有明顯的炎癥和自身免疫反應，在破裂型間盤突出中有比較明顯的免疫炎癥反應，這種反應的主要特征是巨噬細胞浸潤[23]；髓核組織暴露于自身免疫系統(tǒng)中可以激活T、B 細胞，引起自身免疫反應，血管內皮生長因子參與了退變椎間盤新生血管的形成，并可能在椎間盤突出與腰腿痛機制中與 T、B 細胞協(xié)同作用引起自身免疫反應；髓核邊緣有散在的炎性細胞和較多巨噬細胞，但小血管增生和T淋巴細胞浸潤不明顯，提示有炎性反應，但自身免疫反應不如椎間盤突出癥明顯[24-25]。這兩種免疫調控機制分別與淋巴細胞分化的正向調節(jié)、白細胞分化的正向調節(jié)、陽性T細胞選擇有密切聯(lián)系[26-30]。

其中參與這些與免疫炎癥相關的生物學過程的基因CD74已經被文獻證實參與了椎間盤內部調控過程。文獻[31-34]中解釋說明MIF / CD74途徑可以代表治療椎間盤退變的關鍵靶標，MIF可以通過CD74促進炎癥反應，因為用其拮抗劑ISO-1抑制MIF的功能可以減輕MIF誘導的炎癥反應并且發(fā)揮有效的治療效果，以此來減輕椎間盤疾病患者的疼痛，縮短疾病的治愈周期[34-40]。與CD74同時富集到幾個重要的免疫炎癥相關通路的PRKCZ基因和他們靶向的circRNAs(has_circ_0000077，has_circ_0081342，has_circ_0000389，has_circ_0000815)很可能同時作用，遵循ceRNA的競爭調控關系來引發(fā)椎間盤疾病(見圖5(c))。

藥物方面有文獻證實煙酰胺有助于減輕壓力對椎間盤的損傷，能夠促進壓力損傷后的椎間盤恢復。煙酰胺可以抑制IL-1β誘導的椎間盤細胞凋亡，改善IL-1β導致的椎間盤能量代謝障礙，煙酰胺可以促進髓核細胞的增殖并抑制白細胞介素1β誘導的髓核細胞凋亡，通過對煙酰胺的負調控來影響椎間盤的變性[41-43]。

模塊二中的CD74與has_circ_0001320也同樣可能通過ceRNA的關系改變了共表達關系的強度。在正常樣本中他們的共表達相關值為0.868 441，在疾病狀態(tài)下他們的表達相關值為-0.049 57。在正常和疾病中同樣有著明顯的共表達差異，并且由于競爭性關系從負相關變?yōu)檎嚓P同樣暗示了競爭性RNA調控機制的表達趨勢。他們的超幾何P值統(tǒng)計學顯著為2.62×10-6，F(xiàn)DR矯正值為0.015 134，表明通過hsa-miR-17-3p[44]介導的CD74與has_circ_0001320有著潛在的ceRNA關系(圖5(b))。CD74基因的Foldchange值為0.389 194 149在疾病中顯著下調。

圖5 模塊二功能分析圖Fig.5 Functional analysis chart for Module 2

圖5(a)為模塊二網絡圖，我們把文獻證明已知疾病相關基因用紅色圓形節(jié)點表示；把富集到相同生物學過程的基因用黃色圓形節(jié)點表示；把黃色和紅色節(jié)點的一步鄰居circRNAs用紫色菱形節(jié)點表示并且把它們共同的一步鄰居用綠色菱形節(jié)點表示；其他的基因用藍色圓形節(jié)點表示，circRNA用橙色菱形節(jié)點表示。圖5(b)舉例一個ceRNA關系對的MiRanda靶點展示。圖5(c)模塊中所有編碼基因Goterm功能和Pathway通路富集結果，顏色代表P值大小。

3 討論及結論

1)沒有約束差異基因是為了排除常規(guī)差異分析方法的不確定性，以免影響實驗結果。取而代之的是分析circRNA和mRNA之前差異的共表達相關性，從而更準確的分析兩者的競爭性關系。

2)通過cytoscape的一個網絡模塊挖掘插件識別出了五個與椎間盤變性高度相關的功能模塊并且在功能模塊中通過網絡拓撲關系發(fā)現(xiàn)疾病潛在的編碼基因和circRNA。

3)改進了單一的共表達方法，更關注于基因在疾病和差異樣本中差異的共表達調控關系。

4)為椎間盤疾病的預測、藥物靶點治療和circRNA功能的探索提供了新思路。