范敬炎，韓欣桐，丁家安，王樹急，杜晶晶，朱玉峰

(1.錦州醫(yī)科大學；2.錦州醫(yī)科大學附屬第一醫(yī)院，遼寧 錦州 121000)

結腸癌(colon cancer)是一種異質(zhì)性疾病，給患者和醫(yī)療系統(tǒng)造成了巨大危害[1]。因此，在分子水平上更深入地了解其調(diào)控機制，確定可靠的預后生物標志物對于結腸癌患者的生存預測及個體化治療具有重要的作用。

長鏈非編碼RNA(LncRNA)的長度超過200個核苷酸，并且在特定組織以及不同類型腫瘤中差異表達，盡管lncRNA不編碼蛋白質(zhì)，但它們?nèi)匀痪哂性S多功能，與腫瘤發(fā)生發(fā)展密切相關[2]。而且lncRNA在體液中易于檢測，使其成為理想的臨床預后判斷的生物標志物[3]。

隨著各種生物信息學分析工具以及公共數(shù)據(jù)庫的出現(xiàn)，從高通量數(shù)據(jù)中鑒定關鍵基因變得越來越容易。其中，加權基因共表達網(wǎng)絡分析(WGCNA)是尋找高度相關基因模塊的有力工具，已被廣泛用于鑒定候選生物標志物。在這項研究中，我們從公共數(shù)據(jù)庫下載結腸癌患者的RNA表達數(shù)據(jù)和臨床信息，并建立共表達網(wǎng)絡以挖掘可以判斷結腸癌患者預后的lncRNA。使用獨立驗證數(shù)據(jù)集，包括結腸癌樣本組織和臨床生存信息，驗證lncRNA表達水平和結腸癌患者預后的關系。同時深入研究候選lncRNAs在CRC中的表達和功能，以提高我們對結腸癌發(fā)生發(fā)展的分子機制的理解，并為CRC診療提供候選靶標。

1 資料與方法

1.1 數(shù)據(jù)來源

結腸癌的基因表達原始數(shù)據(jù)和臨床信息從TCGA數(shù)據(jù)庫中獲取，該數(shù)據(jù)庫中包含了467份結腸癌樣本，包括268份I～II期和199份III～IV期結腸癌樣本。根據(jù)臨床信息整理相應患者的生存狀態(tài)和生存時間。本研究以TCGA數(shù)據(jù)集為訓練集，構建加權基因共表達網(wǎng)絡，識別結腸癌預后相關基因。使用數(shù)據(jù)集GSE39582作為測試集進行生存分析來驗證我們的結果。數(shù)據(jù)集GSE39582下載自GEO數(shù)據(jù)庫，該數(shù)據(jù)集包括566例結腸癌樣本，其中TNM分期I期33例，II期264例，III期205例，IV期60例，總生存數(shù)據(jù)556例，無復發(fā)生存數(shù)據(jù)519例。使用GSE20916數(shù)據(jù)集驗證候選預后分子在結腸癌中的表達水平，該數(shù)據(jù)集包括101例結腸癌及35例正常結腸組織。使用GSE41568數(shù)據(jù)集觀察候選預后分子在原位結腸癌和轉移性結腸癌組織中的表達水平，該數(shù)據(jù)集包括39例原位結腸癌和94例轉移性結腸癌。數(shù)據(jù)集GSE39582、GSE20916和GSE41568的芯片均在同一平臺(GPL570)上進行檢測。

1.2 方法

1.2.1 篩選差異表達的LncRNA 利用R語言軟件的Limma函數(shù)包對原始表達式數(shù)據(jù)進行預處理。從Gencode數(shù)據(jù)庫(https://www.gencodegenes.org)獲取人類基因組(hg38)和相關注釋文件(31版)。該注釋文件用于識別LncRNA?；蝾愋蜑椤發(fā)incrna”、“antisense”、“processed transcript”、“sense_intronic”、“TEC”、“bidirectional promoter lncRNA”、“sense_overlapping”、“macrolncRNA”或“non coding”的分子定義為LncRNA[4]。當多個探針與一個相同的LncRNA匹配時，我們?nèi)”磉_值的平均值[5]。選擇FDR<0.05，Log2(FC)>2，P<0.05作為篩選差異基因的標準。

1.2.2 加權基因共表達網(wǎng)絡分析 剔除TCGA數(shù)據(jù)集中異常離群樣本，然后逐步進行網(wǎng)絡構建和模塊聚類。為了構建無尺度基因共表達網(wǎng)絡，我們使用了R語言中的WGCNA函數(shù)包[6]。首先對所有的基因對進行皮爾遜相關矩陣分析。然后利用冪函數(shù)amn=|cmn|β(其中amn是m和n基因之間的鄰接，cmn是m和n基因之間的皮爾遜相關性)構造加權鄰接矩陣。參數(shù)β是軟閾值參數(shù)，經(jīng)過計算選擇β=14(無尺度r2=0.91)以確保構建無尺度網(wǎng)絡[7]。然后，將鄰接矩陣轉化為拓撲重疊矩陣，由此得到的拓撲重疊是基于兩個基因間共表達關系的一種有生物學意義的基因相似性度量[8]。采用動態(tài)樹切割法進行模塊識別，用顏色命名模塊，采用Pearson相關分析計算各個模塊與TNM分期、總生存期、無復發(fā)生存期等臨床特征的相關性和P值，相關系數(shù)最高的基因模塊內(nèi)的LncRNA作為候選預后分子標志物納入下一步分析。

1.2.3 篩選預后風險LncRNA及功能預測 使用GSE39582數(shù)據(jù)集對候選LncRNA進行生存分析，繪制Kaplan-meier生存曲線，所有P值小于0.05的LncRNA作為潛在的預后分子標志物。在此基礎上，利用另外的獨立數(shù)據(jù)集GSE20916進行分析，研究正常結腸標本和結腸癌標本中LncRNA表達水平差異。利用GSE41568在原位結腸癌和轉移結腸癌樣本之間驗證LncRNA表達水平差異。為了探索這些LncRNA影響相關臨床特征的潛在機制，將候選LncRNA上傳到LncACTdb2.0數(shù)據(jù)庫中，提取lncRNA-mRNA共表達系數(shù)絕對值大于0.6且P<0.05的網(wǎng)絡，并將其中的mRNA定義為該lncRNA的靶基因。對靶基因進行GO功能和KEGG信號通路富集分析，將P<0.01和FDR<0.01設置為篩選標準[9]。

1.3 統(tǒng)計學方法

使用R語言的 survival 函數(shù)包進行Kaplan-Meier生存分析，生存率比較采用log-rank檢驗法，P< 0.05 表示差異具有統(tǒng)計學意義。

2 結 果

2.1 差異表達基因篩選

在TCGA數(shù)據(jù)庫中下載結腸癌數(shù)據(jù)集，共發(fā)現(xiàn)差異表達基因4227個，其中上調(diào)基因1899個，下調(diào)基因2328個。差異表達基因中，mRNA及LncRNA的數(shù)量分別為1817個(上調(diào)755個，下調(diào)1062個)及2410個(上調(diào)1144個，下調(diào)1266個)。差異基因的表達數(shù)據(jù)譜及其臨床信息將用于進行接下來的加權基因共表達網(wǎng)絡分析。

2.2 加權基因共表達網(wǎng)絡分析

基于加權基因共表達網(wǎng)絡分析將所有差異表達的基因分成12個模塊，見圖1A。為了推測這些基因的臨床意義，將基因模塊與結腸癌患者的臨床信息相結合。腫瘤組織分期及患者生存資料作為選擇功能模塊的重要評估指標。結果顯示粉色模塊(pink module)的特征值與結腸癌患者的總生存期具有較高的相關性(r=0.68，P=0.0001)，見圖1B。為了探索預后分子標志物，本研究選擇粉色模塊中的LncRNA進行下一步的生存分析。

A：動態(tài)分層聚類圖；B：特征基因模塊和結腸癌不同臨床信息相關性的熱圖

2.3 篩選與結腸癌患者總體生存率相關的LncRNA

位于粉色模塊中的9個LncRNA(LINC00473，F(xiàn)AM228B,CTC-277H1.6，LOC642846，CHKB-AS1，LRRC75A-AS1，LINC00299，LOC400684，LINC01021)被視為候選的結腸癌預后分子標志物。下載GEO數(shù)據(jù)庫中的GSE39582數(shù)據(jù)集，利用 R 軟件 survival 軟件包及GSE39582數(shù)據(jù)集中的臨床資料對9個候選LncRNA進行生存分析。根據(jù)每個候選lncRNA的表達值的中位數(shù)將結腸癌患者分成兩組， Kaplan-Meier生存分析結果顯示LINC01021高表達的患者的總體存活率更低(P=0.001)，復發(fā)風險更高(P=0.001)，見圖2。其余LncRNA表達水平與結腸癌患者臨床預后的關系，差異無統(tǒng)計學意義。

2.4 候選LncRNA 在結腸癌中的表達及功能預測

使用另外的獨立數(shù)據(jù)集GSE20916對候選lncRNA在結腸癌中的表達水平進行驗證，LINC01021在結腸癌中表達水平高于正常結腸組織，兩組LINC01021表達水平比較，差異有統(tǒng)計學意義(t=9.549，P<0.01)，見圖3A。此外，在數(shù)據(jù)集GSE41568中，LINC01021在轉移結腸癌中的表達水平高于原位結腸癌，差異有統(tǒng)計學意義(t=4.927，P<0.01)，見圖3B。使用LncACTdb2.0數(shù)據(jù)庫預測LINC01021下游靶基因，共137個基因，見圖4。為了進一步推測LINC01021的功能，對下游靶基因進行功能富集分析和信號通路富集分析，GO分析發(fā)現(xiàn)lncRNA-mRNA共表達網(wǎng)絡中的mRNA主要參與細胞增殖負調(diào)控，調(diào)節(jié)細胞周期素依賴蛋白激酶活性，調(diào)節(jié)細胞周期等，KEGG分析發(fā)現(xiàn)lncRNA-mRNA共表達網(wǎng)絡中的mRNA主要參與PI3K/AKT信號通路激活和CTNNB1磷酸化級聯(lián)反應等，結果見圖5。

圖5 GO功能富集分析和KEGG信號通路富集分析

A：LINC01021在結腸癌組織及正常組織中的表達水平；B：LINC01021在原位結腸癌組織及轉移結腸癌組織中的表達水平

圖2 結腸癌患者差異表達的LINC01021生存曲線分析

紅色代表上調(diào)，綠色代表下調(diào)；圓形代表 mRNA，三角形代表 lncRNA，圖形越大表示共表達系數(shù)絕對值越大

3 討 論

在本研究中，利用TCGA數(shù)據(jù)庫的數(shù)據(jù)集，我們首先通過加權基因共表達網(wǎng)絡分析篩選了與結腸癌預后相關的基因模塊，然后利用GEO數(shù)據(jù)庫中的獨立數(shù)據(jù)集對模塊中的基因進行了生存分析，確定了長鏈非編碼RNA(LINC01021)與結腸癌預后相關，高表達LINC01021的結腸癌患者預后差。通過一系列的生物信息學分析，發(fā)現(xiàn)LINC01021下游的靶基因主要參與結腸癌細胞周期的調(diào)控，以及通過激活PI3K/AKT信號通路和參與CTNNB1磷酸化級聯(lián)反應等過程影響結腸癌的發(fā)生發(fā)展，表明LINC01021可能在結腸癌的發(fā)展中發(fā)揮致癌作用。

LINC01021是長度為1112個堿基的長鏈非編碼RNA，定位于人染色體5p14.1。既往研究報道lncRNA LINC01021是p53的靶基因，能夠與p53直接結合[10]。使用CRISPR /Cas9方法敲除p53結合位點的啟動子序列可以降低結腸癌HCT116細胞系中LINC01021表達，增加結腸癌細胞對多柔比星和5-氟尿嘧啶的化療敏感性[11]。國內(nèi)研究報道LINC01021在食管癌組織和食管癌細胞系中的高表達，LINC01021 基因通過影響上皮間質(zhì)轉化促進食管癌細胞侵襲轉移[12]。我們的結果發(fā)現(xiàn)LINC01021在轉移結腸癌中的表達水平高于原位結腸癌，而結腸癌的侵襲轉移與上皮間質(zhì)轉化密切相關，這是一個值得探索的研究方向。我們的研究結果還發(fā)現(xiàn)LINC01021高表達的結腸癌患者的總體存活率更低，復發(fā)風險更高。這些結果說明LINC01021表達的定量可能具有重要的臨床預后價值。

我們通過生物信息學分析發(fā)現(xiàn)LINC01021下游靶基因發(fā)揮的作用機制主要是影響細胞周期調(diào)控，以及通過激活PI3K/AKT信號通路和參與CTNNB1磷酸化級聯(lián)反應。PI3K/AKT信號通路有廣泛的生物活性，可促進細胞增殖，在惡性腫瘤內(nèi)往往呈過度激活的狀態(tài)，是最主要的抑制細胞凋亡的信號途徑[13]。CTNNB1磷酸化的降低往往是由于Wnt信號通路活化，隨著信號通路活化和CTNNB1蛋白在核內(nèi)積累，將導致結腸癌發(fā)生[14]。因此，本研究結果對未來探索LINC01021調(diào)控結腸癌發(fā)生的分子機制奠定了一定的理論基礎，具有非常重要的指導意義。當然，需要進一步的實驗來驗證LINC01021的臨床和生物學功能。

總之，我們采用WGCNA等生物信息學方法從TCGA數(shù)據(jù)庫中研究了結腸癌患者的RNA序列和臨床數(shù)據(jù)。我們得出結論是LINC01021與結腸癌患者的預后相關，高表達LINC01021的結腸癌患者預后差。LINC01021有可能成為一種新的預后指標，有助于結腸癌患者個性化治療及臨床預后判斷。