基于基因共表達(dá)權(quán)重網(wǎng)絡(luò)分析骨質(zhì)疏松相關(guān)的關(guān)鍵LncRNA

黎永華 譚仁霆 黃世福 邱金龍

海南省儋州市中醫(yī)醫(yī)院骨傷科，海南 儋州571700

骨質(zhì)疏松癥是常見的骨代謝疾病。骨質(zhì)疏松癥患者常伴隨著脊柱、髖部和腕部骨折的發(fā)生及其他骨科圍手術(shù)期的并發(fā)癥[1]。目前的研究重點(diǎn)聚焦在骨質(zhì)疏松的發(fā)病機(jī)制，但卻忽略了骨質(zhì)疏松的治療手段[2]。而循環(huán)單核細(xì)胞作為外周細(xì)胞的一種，可分化為巨噬細(xì)胞，樹突狀細(xì)胞和破骨細(xì)胞等，并作為破骨細(xì)胞的前體細(xì)胞用于破骨細(xì)胞分化[3]。而破骨分化可以激活和凋亡某些細(xì)胞因子(如IL-1，IL-6)，這些相關(guān)因子對(duì)骨平衡和骨代謝至關(guān)重要[4]。因此，循環(huán)單核細(xì)胞的功能與骨質(zhì)疏松癥的發(fā)病機(jī)制密切相關(guān)。絕經(jīng)后卵巢功能下降，會(huì)導(dǎo)致雌激素分泌減少，隨后抑制RANKL蛋白的表達(dá)，從而導(dǎo)致單核細(xì)胞的破骨分化能力下調(diào)，阻礙了骨重塑的進(jìn)程[5-6]。LncRNA是一種長(zhǎng)度大于200個(gè)單位的非編碼RNA，參與多種疾病的發(fā)生發(fā)展進(jìn)程。有研究發(fā)現(xiàn)，LncRNA可以通過負(fù)反饋的方式反向調(diào)節(jié)骨髓間充質(zhì)干細(xì)胞(BMCs)的成骨分化，最終導(dǎo)致骨質(zhì)疏松的形成[7]。而機(jī)械應(yīng)力和金屬顆粒作為常見的外源性因素，可以誘導(dǎo)LncRNA的表達(dá)出現(xiàn)異常從而抑制成骨細(xì)胞功能并促進(jìn)其凋亡，導(dǎo)致骨質(zhì)疏松的產(chǎn)生[8-9]。LncRNA還能通過Notch通路介導(dǎo)轉(zhuǎn)錄和免疫因子的表達(dá)促進(jìn)破骨細(xì)胞的產(chǎn)生和增加其活性[10]。

骨質(zhì)疏松癥與單核細(xì)胞之間的分子作用機(jī)制受眾多基因調(diào)控[11]，眾多基因的影響增加了該領(lǐng)域的研究難度[12]。而目前，加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)[13]被廣泛報(bào)道用于生物基因研究中，而WGCNA與既往的差異基因表達(dá)分析不同的是，WGCNA關(guān)注基因間的關(guān)聯(lián)，并通過聚類的方式將它們歸納到各類型的模塊部分中。建立起臨床性狀與模塊基因的關(guān)聯(lián)，并確定與感興趣性狀相關(guān)的模塊，之后通過計(jì)算獲得模塊特征值用于選取模塊基因。

在研究中，我們通過使用來自GEO數(shù)據(jù)庫的80個(gè)單核細(xì)胞樣品數(shù)據(jù)構(gòu)建骨質(zhì)疏松的共表達(dá)模塊，將共表達(dá)模塊與骨密度性狀結(jié)合，獲得相關(guān)模塊，尋找與骨密度具有密切關(guān)聯(lián)性的LncRNA，通過數(shù)據(jù)庫獲得LncRNA的靶基因，建立靶基因-LncRNA的互作網(wǎng)絡(luò)，對(duì)這些潛在LncRNA的靶基因進(jìn)行基因的富集分析，分析與模塊基因相關(guān)的生物過程和相關(guān)通路。嘗試為骨質(zhì)疏松在LncRNA層面上的研究指明方向

1 材料和方法

1.1 數(shù)據(jù)獲取

表達(dá)譜數(shù)據(jù)來自GEO數(shù)據(jù)庫，數(shù)據(jù)集編號(hào)為GSE56815[14]。分別獲得40例正常骨密度和40例低骨密度的絕經(jīng)前后婦女的基因數(shù)據(jù)，State組用以區(qū)分絕經(jīng)前女性和絕經(jīng)后女性。所有樣本均采用 GPL96 [HG-U133 A] Affymetrix Human Genome U133 A Array平臺(tái)檢測(cè)。

1.2 數(shù)據(jù)處理和重注釋

采用R軟件中的 oligo軟件[15]通過RMA算法對(duì)下載的基因數(shù)據(jù)進(jìn)行基因校正、標(biāo)準(zhǔn)化和匯總。運(yùn)用Ensemble Gene 97數(shù)據(jù)庫[16]進(jìn)行重注釋，得到基因類型和基因符號(hào)的對(duì)應(yīng)關(guān)系，根據(jù)基因類型獲得LncRNA。當(dāng)多個(gè)探針對(duì)應(yīng)同一個(gè)LncRNA時(shí)，取其平均值作為L(zhǎng)ncRNA的最終表達(dá)值。

1.3 構(gòu)建共表達(dá)模塊

該分析用于對(duì)同一類基因進(jìn)行聚類分群。樣本聚類分析是通過R的 WGCNA[17]軟件對(duì)LncRNA數(shù)據(jù)進(jìn)行處理后獲得的。并對(duì)樣本進(jìn)行構(gòu)建聚類樹，尋找潛在的離群樣本并將其刪除。隨后根據(jù)公式構(gòu)建軟閾值，合適的軟閾值需要具有較好的平均連通性，同時(shí)符合撲擬合指數(shù)曲線并呈現(xiàn)線性關(guān)系。隨后對(duì)LncRNA進(jìn)行共表達(dá)模塊的構(gòu)建，模塊中的最小基因數(shù)設(shè)為3，相關(guān)模塊的顏色可視化通過colors函數(shù)進(jìn)行聚類樹呈現(xiàn)。

1.4 識(shí)別模塊和BMD之間的相關(guān)性

首先，我們對(duì)基因表達(dá)與臨床信息的線性回歸關(guān)系中的P值進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換為log10并將其定義為基因顯著性GS(gene significance)，將模塊成員顯著性(MM)定義為模塊中所有基因與模塊的關(guān)聯(lián)，具有高構(gòu)建性值的模型被認(rèn)為與臨床特征具有顯著相關(guān)性。此外，對(duì)每個(gè)基因模塊進(jìn)行主成分分析。計(jì)算模塊特征基因和BMD之間的皮爾遜相關(guān)性系數(shù)，以確定與骨質(zhì)疏松相關(guān)的模塊。采用t檢驗(yàn)統(tǒng)計(jì)皮爾遜相關(guān)性系數(shù)，當(dāng)P值<0.05時(shí)認(rèn)為模塊與BMD具有顯著相關(guān)性。

1.5 網(wǎng)絡(luò)可視化與LncRNA 功能預(yù)測(cè)

選取與臨床表型顯著相關(guān)的基因模塊，使用starbase[18]和DIAND-LncBase數(shù)據(jù)庫預(yù)測(cè)LncRNA的靶向基因 (http://starbase.sysu.edu.cn/index.php和 http://carolina.imis.athena-innovation.gr/diana_tools/web/)。使用Cytoscape網(wǎng)絡(luò)關(guān)系軟件進(jìn)行共表達(dá)關(guān)系網(wǎng)絡(luò)圖的繪制，同時(shí)結(jié)合LncRNA的靶基因構(gòu)建共表達(dá)網(wǎng)絡(luò)。在 R軟件上采用ClusterProfiler[19]數(shù)據(jù)對(duì)模塊基因進(jìn)行GO(基因本體論)和KEGG(京都基因與基因組百科全書)分析。

2 結(jié)果

2.1 數(shù)據(jù)處理及共表達(dá)模塊的構(gòu)建

預(yù)處理基因數(shù)據(jù)及臨床特征，對(duì)重復(fù)的基因和缺失值進(jìn)行剔除，共獲得127個(gè)LncRNA的表達(dá)矩陣。為確保結(jié)果的準(zhǔn)確性檢測(cè)離群樣本，檢測(cè)未見明顯離群樣本，隨后進(jìn)行樣本聚類分析(圖1)。當(dāng)軟閾值等于4時(shí)，共表達(dá)網(wǎng)絡(luò)接近為無尺度網(wǎng)絡(luò)，此時(shí)的數(shù)值是令曲線趨于平滑的最小閾值，這有利于網(wǎng)絡(luò)的平均連接程度維持在穩(wěn)定的狀態(tài)從而包含足夠多的信息(圖2)。選取軟閾值為4并得到基因聚類樹，并根據(jù)上述方法進(jìn)行切割設(shè)置。最終得到9個(gè)LncRNA模塊，模塊中的基因個(gè)數(shù)從4到29不等(見圖3)。此外， gray模塊不屬于本研究中所定義的模塊(圖3)。

圖1 樣品樹狀圖和性狀熱圖Fig.1 Sample tree diagram and heat map of traits

圖2 軟閾值權(quán)重的網(wǎng)絡(luò)拓?fù)渚仃嘑ig.2 Network topology matrix with soft threshold weight

圖3 基因共表達(dá)網(wǎng)絡(luò)分層聚類樹與共表達(dá)模塊Fig.3 Gene co-expression network hierarchical clustering tree and co-expression module

2.2 基因模塊與臨床性狀關(guān)聯(lián)分析

將患者的臨床信息與模塊數(shù)據(jù)合并后進(jìn)行相關(guān)性分析，尋找特征相關(guān)性最高的模塊(圖4)。我們發(fā)現(xiàn)多個(gè)模塊與骨質(zhì)疏松相關(guān)，進(jìn)一步分析發(fā)現(xiàn)， turquoise模塊與骨質(zhì)疏松呈正相關(guān)性(r=0.31，P=0.005)。針對(duì)這個(gè)模塊，進(jìn)行相關(guān)性分析并繪制散點(diǎn)圖分析模塊內(nèi)基因是否符合線性(圖5)。我們根據(jù)皮爾遜相關(guān)性系數(shù)發(fā)現(xiàn)這些基因不僅是模塊中的重要基因還與臨床性狀密切相關(guān)。

注：左側(cè)為模塊基因與性狀的相關(guān)關(guān)系圖，每個(gè)單元格包含相應(yīng)的相關(guān)性和P值，紅色為上調(diào)，藍(lán)色為下調(diào)。右側(cè)為模塊的重要性柱狀圖，條形柱的高度帶模塊與性狀的重要性越高代表越重要。圖4 性狀模塊相關(guān)關(guān)系圖Fig.4 Correlation diagram of trait modules

圖5 基因模塊表達(dá)水平與基因重要性的散點(diǎn)圖Fig.5 Scatter diagram of gene module expression level and gene importance

2.3 獲得LncRNA靶基因并構(gòu)建可視化網(wǎng)絡(luò)

分析模塊內(nèi)基因的相關(guān)基因，模塊內(nèi)共有29個(gè)LncRNA。登錄STARBASE和DIAND-LncBase網(wǎng)站輸入獲得的29個(gè)LncRNA，最終共獲得204個(gè)靶向基因。將靶基因與LncRNA的關(guān)聯(lián)數(shù)據(jù)導(dǎo)入Cytoscape構(gòu)建靶基因-LncRNA可視化網(wǎng)絡(luò)和核心基因互作用網(wǎng)絡(luò)(圖6)。turquoise模塊的樞紐LncRNA為ACVR2B-AS1、LINC01278、SND1-IT1、C22orf24、ZNF213-AS1、WT1-AS、PLAC4、RHPN1-AS1和LINC01140。

圖6 靶基因-LncRNA相互作用網(wǎng)絡(luò)和核心基因互作網(wǎng)絡(luò)圖Fig.6 Target gene-LncRNA interaction network and core gene interaction network diagram

2.4 GO和KEGG富集分析

對(duì)turquoise模塊中LncRNA的靶向的204個(gè)基因進(jìn)行GO功能分析和KEGG富集分析發(fā)現(xiàn)。分子功能主要集中在核糖核酸、RNA和蛋白質(zhì)的結(jié)合；生物功能主要集中在細(xì)胞和蛋白質(zhì)運(yùn)輸、成骨細(xì)胞的分化和神經(jīng)等調(diào)控上；細(xì)胞功能主要集中在細(xì)胞核、細(xì)胞膜和線粒體中。GO功能主要與細(xì)胞物質(zhì)的轉(zhuǎn)運(yùn)、蛋白和RNA的合成相關(guān)，涉及成骨細(xì)胞分化的調(diào)控。KEGG主要富集于GnRH、Notch、神經(jīng)營(yíng)養(yǎng)和鈣離子信號(hào)通路上(圖7)。

注：MF為分子功能，BP為生物過程，CC為細(xì)胞組成，KEGG為富集信號(hào)通路，顏色的變化與P值相關(guān)，柱圖的長(zhǎng)度與關(guān)聯(lián)基因的數(shù)量呈正相關(guān)。圖7 turquoise模塊的GO和KEGG富集分析Fig.7 GO and KEGG enrichment analysis of turquoise module

3 討論

骨質(zhì)疏松癥作為常見骨病，可被分為原發(fā)型和繼發(fā)型兩種類型，最常見的一種是絕經(jīng)后骨質(zhì)疏松癥。由于患者成骨和破骨的平衡被打破導(dǎo)致成骨減少和破骨增多，最終引起骨質(zhì)疏松癥的發(fā)生[20]。臨床上，骨質(zhì)疏松常伴隨著骨折的發(fā)生[21]。因此，為了尋找靶向治療方法，需要深入了解其相關(guān)機(jī)制及基因調(diào)控的運(yùn)行。既往的研究中，Zhang等[22]基于基因組表達(dá)數(shù)據(jù)，通過整合多組學(xué)數(shù)據(jù)找到與患者BMD相關(guān)的Mrna和Mirna。而本研究則是基于已有的測(cè)序數(shù)據(jù)，進(jìn)行WGCNA分析尋找與骨質(zhì)疏松相關(guān)的LncRNA。

相較于既往的差異基因表達(dá)分類法，WGCNA不僅對(duì)基因進(jìn)行聚類分群，還將基因與臨床性狀信息進(jìn)行關(guān)聯(lián)，有利于發(fā)現(xiàn)某些在統(tǒng)計(jì)學(xué)上無差異但又在人體生物過程中具有重要意義的基因。結(jié)合臨床特征，綜合分析不同臨床特征的重要性，有利于尋找不同臨床特征在基因?qū)用娴南嚓P(guān)性，并提高了芯片數(shù)據(jù)的利用率，從而獲得針對(duì)特定形狀的關(guān)鍵聚類基因。獲得關(guān)鍵模塊的基因后，對(duì)聚類基因進(jìn)GO功能注釋和KEGG富集分析幫助我們更好的了解基因的互作用及潛在通路和機(jī)制。

我們得到9個(gè)有效模塊，而結(jié)果顯示turquoise模塊與骨質(zhì)疏松顯著相關(guān)。提取了29個(gè)LncRNA的相關(guān)靶點(diǎn)基因信息，并進(jìn)行網(wǎng)絡(luò)可視化和明確樞紐基因。LncRNA在體內(nèi)不能翻譯為蛋白質(zhì)，但參與多種細(xì)胞過程，由于具有組織特異性，被廣泛用于分子靶點(diǎn)和標(biāo)志物的研究中。許多LncRNA僅位于細(xì)胞核中，LncRNA調(diào)節(jié)著基因的表達(dá)，它可以通過多種不同的機(jī)制對(duì)基因進(jìn)行干預(yù)。雖然本研究中在turquoise模塊中只獲得29個(gè)LncRNA，但獲得其潛在靶基因204個(gè)，可見LncRNA可通過多種機(jī)制調(diào)控基因的表達(dá)。

這些LncRNA主要參與細(xì)胞間信息交流途徑，骨形成依賴于骨髓單核細(xì)胞供應(yīng)代謝物，通過協(xié)調(diào)骨形成和骨吸收的平衡維持正常骨代謝，防止骨量丟失，這些過程都需要細(xì)胞之間的信息交換[23]。例如，細(xì)胞外纖調(diào)蛋白和BMP2通路可以被LncRNA激活進(jìn)而促進(jìn)BMCs的成骨分化，最終抑制成脂分化[24]。此外，我們的研究也發(fā)現(xiàn)GnRH、Notch、神經(jīng)營(yíng)養(yǎng)因子和鈣離子信號(hào)通路，也參與了LncRNA相關(guān)的骨形成途徑。有研究指出腦垂體分泌的促性腺激素(Gn)刺激破骨細(xì)胞的形成和功能調(diào)控，其可抑制骨轉(zhuǎn)換和骨量變化，與骨質(zhì)疏松的發(fā)病密切相關(guān)[25]。而鈣離子通路的調(diào)控，與鈣離子的吸收和骨代謝密切相關(guān)。神經(jīng)營(yíng)養(yǎng)因子通路參與調(diào)節(jié)骨髓干細(xì)胞的成骨分化。Nahlé的研究發(fā)現(xiàn)了神經(jīng)營(yíng)養(yǎng)因子通過結(jié)合小鼠的MSCs，激活STAT1和STAT3的磷酸化，抑制與成骨相關(guān)的主要基因的上調(diào)，最終阻止成骨細(xì)胞的生成和礦化[26]。值得注意的是，已有研究表明Notch通路在骨質(zhì)疏松中發(fā)揮重要作用[27]。成骨細(xì)胞和破骨細(xì)胞的動(dòng)態(tài)平衡受到Notch通路的調(diào)控，成骨細(xì)胞的表達(dá)會(huì)隨著Notch的敲除而受到抑制[28]。此外，Canalis等[29]的研究發(fā)現(xiàn)Notch可以抑制成骨細(xì)胞的分化。有實(shí)驗(yàn)證明了Notch與Runx2存在負(fù)反饋調(diào)節(jié)作用，Runx2作為OPG的啟動(dòng)抑制劑，當(dāng)Notch被沉默后，Runx2表達(dá)增多，會(huì)引起OPG的轉(zhuǎn)錄被抑制，破骨細(xì)胞形成增加[30]。這些通路解釋了LncRNA靶基因與骨質(zhì)疏松的關(guān)聯(lián)，為骨質(zhì)疏松的未來的研究提供了新的研究思路。

據(jù)我們所知，這是首次將WGCNA與骨質(zhì)疏松相關(guān)LncRNA相結(jié)合的研究。但是，本文仍存在著幾點(diǎn)不足。首先，探針重注釋的方法可以幫助我們從既往的Geo芯片中獲得可靠的LncRNA數(shù)據(jù)，但由于芯片測(cè)序平臺(tái)的不同以及芯片的更新迭代不同，探針重注釋法不能覆蓋所有LncRNA。其次，本研究的結(jié)果仍需要通過細(xì)胞和分子實(shí)驗(yàn)來驗(yàn)證我們的結(jié)論，以了解這些檢測(cè)到的lncRNA是否與骨質(zhì)疏松具有因果關(guān)系。再者，本研究的樣本量相對(duì)較小，需要結(jié)合更大樣本的數(shù)據(jù)進(jìn)行外驗(yàn)證并通過相關(guān)的機(jī)器學(xué)習(xí)算法尋找重要性更高的靶標(biāo)LncRNA。再者，這些樣本的臨床數(shù)據(jù)不完整，無法獲得更多的相關(guān)數(shù)據(jù)，我們無法獲得骨質(zhì)疏松的血液指標(biāo)數(shù)據(jù)和更多臨床相關(guān)數(shù)據(jù)，如果有了這部分的數(shù)據(jù)，我們的結(jié)果會(huì)更有說服力，結(jié)果也會(huì)更有意義。

本研究發(fā)現(xiàn)這些LncRNA與骨質(zhì)疏松的潛在發(fā)病具有密切關(guān)系。但尚未檢索到關(guān)于本研究所發(fā)現(xiàn)的核心LncRNA與骨質(zhì)疏松相關(guān)文獻(xiàn)，因此這些潛在的LncRNA靶點(diǎn)既有對(duì)骨質(zhì)疏松發(fā)病的預(yù)測(cè)意義，也為我們未來的骨質(zhì)疏松研究提供了潛在的靶點(diǎn)。綜上所述，我們通過加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析方法尋找了的與骨質(zhì)疏松相關(guān)的LncRNA，獲得了核心LncRNA，并通過尋找潛在的靶基因，進(jìn)行了富集分析。這些結(jié)果為L(zhǎng)ncRNA在骨質(zhì)疏松癥病理生理機(jī)制的研究中提供新的見解。