基于加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析篩選結(jié)腸癌預(yù)后風(fēng)險(xiǎn)長(zhǎng)鏈非編碼RNA

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

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

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

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

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

1 資料與方法

1.1 數(shù)據(jù)來(lái)源

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

1.2 方法

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

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

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

1.3 統(tǒng)計(jì)學(xué)方法

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

2 結(jié) 果

2.1 差異表達(dá)基因篩選

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

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

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

A：動(dòng)態(tài)分層聚類(lèi)圖；B：特征基因模塊和結(jié)腸癌不同臨床信息相關(guān)性的熱圖

2.3 篩選與結(jié)腸癌患者總體生存率相關(guān)的LncRNA

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

2.4 候選LncRNA 在結(jié)腸癌中的表達(dá)及功能預(yù)測(cè)

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

圖5 GO功能富集分析和KEGG信號(hào)通路富集分析

A：LINC01021在結(jié)腸癌組織及正常組織中的表達(dá)水平；B：LINC01021在原位結(jié)腸癌組織及轉(zhuǎn)移結(jié)腸癌組織中的表達(dá)水平

圖2 結(jié)腸癌患者差異表達(dá)的LINC01021生存曲線(xiàn)分析

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

3 討 論

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

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

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

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