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

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

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

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

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

1 材料和方法

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

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

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

采用R軟件中的 oligo軟件[15]通過RMA算法對下載的基因數(shù)據(jù)進行基因校正、標準化和匯總。運用Ensemble Gene 97數(shù)據(jù)庫[16]進行重注釋，得到基因類型和基因符號的對應關系，根據(jù)基因類型獲得LncRNA。當多個探針對應同一個LncRNA時，取其平均值作為LncRNA的最終表達值。

1.3 構建共表達模塊

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

1.4 識別模塊和BMD之間的相關性

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

1.5 網絡可視化與LncRNA 功能預測

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

2 結果

2.1 數(shù)據(jù)處理及共表達模塊的構建

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

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

圖2 軟閾值權重的網絡拓撲矩陣Fig.2 Network topology matrix with soft threshold weight

圖3 基因共表達網絡分層聚類樹與共表達模塊Fig.3 Gene co-expression network hierarchical clustering tree and co-expression module

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

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

注：左側為模塊基因與性狀的相關關系圖，每個單元格包含相應的相關性和P值，紅色為上調，藍色為下調。右側為模塊的重要性柱狀圖，條形柱的高度帶模塊與性狀的重要性越高代表越重要。圖4 性狀模塊相關關系圖Fig.4 Correlation diagram of trait modules

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

2.3 獲得LncRNA靶基因并構建可視化網絡

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

圖6 靶基因-LncRNA相互作用網絡和核心基因互作網絡圖Fig.6 Target gene-LncRNA interaction network and core gene interaction network diagram

2.4 GO和KEGG富集分析

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

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

3 討論

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

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

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

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

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

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