肺腺癌相關(guān)長鏈非編碼RNA生物標(biāo)志物識別

于曉慶

(上海應(yīng)用技術(shù)大學(xué) 理學(xué)院，上海 201418)

肺癌是全球因癌癥導(dǎo)致人類死亡的最大因素。到目前為止，肺癌的5年生存率僅為18%[1]。在肺癌的亞型中，非小細(xì)胞肺癌約占85%。而肺腺癌(lung adenocarcinoma,LUAD)是非小細(xì)胞肺癌中最常見的惡性病理類型。近些年，隨著靶向藥物的出現(xiàn)，肺腺癌的5年生存率不斷在提高，但由于缺乏有效的診斷和預(yù)后生物標(biāo)志物，死亡率仍然很高[2-3]。 因此，為了降低肺腺癌患者的死亡率，迫切需要尋找有效的分子標(biāo)志物，為患者的診斷和臨床治療提供一定的理論依據(jù)。

長非編碼RNA(long noncoding RNA,LncRNA)是長度超過200個核苷酸的非編碼基因。過去很長一段時間，這些非編碼基因曾被認(rèn)為是基因組的“暗物質(zhì)”，并不具有生物學(xué)功能。然而，隨著實驗技術(shù)和計算機(jī)技術(shù)的發(fā)展，越來越多的研究表明，LncRNA通過不同的機(jī)制參與了細(xì)胞的幾乎整個生命過程，并在細(xì)胞分化、代謝、細(xì)胞增殖和凋亡、轉(zhuǎn)錄、染色質(zhì)水平的表觀遺傳學(xué)狀態(tài)調(diào)控等許多基本和關(guān)鍵的生物過程中都起到了非常重要的作用[4-5]。同時，很多研究發(fā)現(xiàn)LncRNA與包括肺腺癌在內(nèi)的多種癌癥的發(fā)展和生存都有關(guān)聯(lián)性，強(qiáng)調(diào)了LncRNA在肺腺癌發(fā)生和發(fā)展中的分子機(jī)制和生物學(xué)特性，表明某些lncRNA可以作為肺腺癌患者預(yù)后的分子標(biāo)志物[6-7]。

個體化醫(yī)療的一個關(guān)鍵問題就是闡明疾病個體的分子機(jī)制。復(fù)雜疾病的發(fā)生，通常并不是由于單個分子的功能障礙導(dǎo)致，而是由于隨著時間和條件的動態(tài)變化，相關(guān)系統(tǒng)或者網(wǎng)絡(luò)的功能紊亂而引起[8-9]。一個樣本個體在特定狀態(tài)下的分子網(wǎng)絡(luò)，稱其為單樣本網(wǎng)絡(luò)。單樣本網(wǎng)絡(luò)是能夠準(zhǔn)確描述疾病個體在特定條件下疾病狀態(tài)的關(guān)鍵個性化特征。對樣本個體在特定條件下的單樣本網(wǎng)絡(luò)進(jìn)行評估和分析，能夠從系統(tǒng)和網(wǎng)絡(luò)水平揭示復(fù)雜疾病發(fā)生發(fā)展的分子機(jī)制。本文基于TCGA數(shù)據(jù)庫中肺腺癌的RNA-seq數(shù)據(jù)和臨床數(shù)據(jù)，利用單樣本網(wǎng)絡(luò)方法構(gòu)建肺腺癌患者的單樣本網(wǎng)絡(luò)模型，分析和識別出6個與肺腺癌相關(guān)的lncRNA分子標(biāo)志物。功能富集分析表明，6個lncRNA可能參與了肺腺癌患者細(xì)胞的代謝過程。本文的研究可以為肺腺癌的診斷、臨床治療以及預(yù)后提供可能的分子標(biāo)志物。

1 數(shù)據(jù)準(zhǔn)備

從TCGA (The Cancer Genome Atlas) 數(shù)據(jù)庫下載了肺腺癌的 RNA-seq 數(shù)據(jù)和臨床試驗數(shù)據(jù)，包括535個肺腺癌的組織樣本(tumor samples)以及59個癌旁的組織樣本(normal samples)。

2 方法

2.1 差異表達(dá)lncRNA的識別

我們使用R軟件包“edgeR”對肺腺癌表達(dá)數(shù)據(jù)中所有的lncRNA進(jìn)行基因的差異表達(dá)分析。該軟件包為多組實驗提供了精確的統(tǒng)計方法，它的一個特殊功能是經(jīng)驗貝葉斯方法，該方法即使在生物復(fù)制水平最低的情況下也可以估算基因特異性的生物變異。在我們的研究中，F(xiàn)DR<0.05且 |FC|>2為篩選差異表達(dá)的lncRNAs的臨界條件。

2.2 單樣本網(wǎng)絡(luò)的構(gòu)建

單樣本網(wǎng)絡(luò)是基于單個樣本對于給定的一組參考樣本的擾動分析而構(gòu)建的，它可以精確地刻畫個體樣本的疾病狀態(tài)[10]。單樣本網(wǎng)絡(luò)的構(gòu)建需要兩組數(shù)據(jù)，一組是參照樣本數(shù)據(jù)，另一組是測試樣本數(shù)據(jù)。通常，這兩組數(shù)據(jù)有著共同的特征屬性。在本研究中，59個肺腺癌的癌旁樣本作為參照樣本組，535個肺腺癌的腫瘤樣本作為測試樣本組。首先，利用皮爾森相關(guān)系數(shù)(pearson correlation coefficient，PCC)計算每對lncRNA基因表達(dá)譜之間的相關(guān)性從而得到一個相關(guān)性網(wǎng)絡(luò)，稱為參照網(wǎng)絡(luò)(reference network)，記為Nr。然后，再將一個測試樣本s添加到該參考樣本組中，利用上述相關(guān)性的計算方法，可以得到一個新的相關(guān)性網(wǎng)絡(luò)，稱為擾動網(wǎng)絡(luò)(perturbed network)，記為Np。最后，通過比較參照網(wǎng)絡(luò)和擾動網(wǎng)絡(luò)之間的差異獲得一個差異網(wǎng)絡(luò)，被稱為樣本s的單樣本網(wǎng)絡(luò)(single-sample network),記為Nssn(Nssn=|Nr-Np|)(見圖1)。

圖1 單樣本網(wǎng)絡(luò)構(gòu)建流程圖Fig.1 Flowchart for constructing a single-sample network

按照上述方法，輪流將535個腫瘤樣本依次作為測試樣本添加到參照樣本組中，可以得到535個單樣本網(wǎng)絡(luò)。容易看出，如果測試樣本在基因表達(dá)模式上與參照樣本組的表達(dá)模式相似，即使將測試樣本添加到參照樣本組中，任何2個分子之間的PCC變化將非常小。相反，如果測試樣本和參照樣本組在表達(dá)模式上存在明顯差異，則將測試樣本添加到參照樣本組后將引起明顯的變化?？梢钥闯觯瑔螛颖揪W(wǎng)絡(luò)是由樣本s決定的，因此它可以用來作為描述樣本s的一個特定特征。我們對所有的腫瘤樣本構(gòu)建其對應(yīng)的單樣本網(wǎng)絡(luò)并進(jìn)行統(tǒng)計分析，篩選出可能與肺腺癌相關(guān)的lncRNA分子標(biāo)志物。

2.3 生存分析

COX比例風(fēng)險回歸模型，是由英國統(tǒng)計學(xué)家D.R.Cox 提出的一種半?yún)?shù)回歸模型。該模型以生存結(jié)局和生存時間為因變量，能夠分析各種影響因素對病人生存期的影響。COX回歸分析是醫(yī)學(xué)領(lǐng)域?qū)τ诩膊☆A(yù)后分析最常用的統(tǒng)計方法之一。我們從TCGA數(shù)據(jù)庫下載了肺腺癌病人的臨床試驗數(shù)據(jù)，同時篩選掉沒有臨床數(shù)據(jù)的樣本，一共得到512個肺腺癌的樣本。利用單變量Cox回歸分析來評估生存率和每個lncRNA表達(dá)水平之間的關(guān)聯(lián)性。

2.4 功能分析

到目前為止，對與肺腺癌相關(guān)的lncRNA的功能機(jī)制的報道還很少。為了進(jìn)一步驗證與肺腺癌相關(guān)的lncRNA的分子功能，利用相關(guān)性分析計算了并識別出lncRNA的共表達(dá)基因，通過對其共表達(dá)基因功能的探索，進(jìn)而研究lncRNA的分子功能。這里，把計算結(jié)果中PCC>0.45，且P<0.01的編碼基因，認(rèn)為是識別出的lncRNA的共表達(dá)基因。

Metascape(http://metascape.org/gp/index.html)是一個免費(fèi)的基因注釋和分析資源，可以幫助研究者了解一個或多個基因的相關(guān)功能[11]。我們通過Metascape對共表達(dá)基因進(jìn)行Go和KEGG功能富集分析。這里，P的臨界值設(shè)定為P<0.01。

3 結(jié)果與討論

3.1 差異表達(dá)分析結(jié)果

從下載的RNA-Seq數(shù)據(jù)中，我們提取了 14 823 個lncRNA的表達(dá)數(shù)據(jù)。 其中，腫瘤樣本未表達(dá)數(shù)量超過10% 的lncRNA數(shù)據(jù)被剔除，共得到了 4 311 個lncRNA。通過edgeR差異表達(dá)分析得到了 1 178 個差異表達(dá)的lncRNA，其中包括851個表達(dá)下調(diào)的lncRNA和327個表達(dá)上調(diào)的lncRNA(見圖2)。圖2中3個部分分別表示表達(dá)下調(diào)、無差異表達(dá)、表達(dá)上調(diào)的lncRNA。

圖2 差異表達(dá)lncRNAs的火山圖Fig.2 The volcano plot of the differentially expressed lncRNAs

3.2 單樣本網(wǎng)絡(luò)分析結(jié)果

為了方便計算，我們把通過2.2得到的單樣本網(wǎng)絡(luò)轉(zhuǎn)化為一個鄰接矩陣 ΔD。容易知道，ΔDi,j表示單樣本網(wǎng)絡(luò)中第i個lncRNA和第j個lncRNA在參照網(wǎng)絡(luò)和擾動網(wǎng)絡(luò)中對應(yīng)邊的PCC的差值，記為ΔPCC?；诩僭O(shè)：如果某個lncRNA可能是肺腺癌相關(guān)的分子標(biāo)志物，則在差異網(wǎng)絡(luò)中，由該lncRNA連接的所有邊對應(yīng)的ΔPCC的總和將比其他lncRNA對應(yīng)的總和更高。這里，用向量SD來表示 1 178 個lncRNA在每個單樣本網(wǎng)絡(luò)中連接的所有邊所對應(yīng)的ΔPCC的總和，記作

SD=(SD1,SD2,…,SDi,…,SD1 178)T

(1)

式中，SDi可通過以下方式計算：

(i=1,2,…,1 178,j=1,2,…,1 178)

(2)

為了找到與肺腺癌相關(guān)的lncRNA分子標(biāo)志物，我們對每個單樣本網(wǎng)絡(luò)中所有l(wèi)ncRNA的SDi值進(jìn)行排序。SDi的值越大，其成為肺腺癌的分子標(biāo)志物的可能性也就越大。對于所有535個單樣本網(wǎng)絡(luò)，我們用一個矩陣M1 178×535來表示

M=

i=1,2,…,1 178,j=1,2,…,535

(3)

式中，SD1,j>SD2,j>…>SDi,j>SD1 178,j，i表示第i個lncRNA，j表示第j個腫瘤樣本。

我們?nèi)∩鲜鼍仃嘙的前10行(即在每個腫瘤樣本排序中的前10個lncRNA)進(jìn)行統(tǒng)計分析。在535個腫瘤樣本中，出現(xiàn)頻率最高的30個lncRNA作為我們通過單樣本網(wǎng)絡(luò)方法識別出的可能與肺腺癌相關(guān)的候選lncRNA分子標(biāo)志物。

3.3 生存分析結(jié)果

從肺腺癌患者的臨床數(shù)據(jù)中選擇整體生存期(overall survival，OS)進(jìn)行單因素cox風(fēng)險分析，結(jié)果顯示在30個lncRNA中有6個與患者的OS顯著相關(guān)(P<0.05)。6個lncRNA的詳細(xì)信息見表1?；谶@6個lncRNA建立了一個風(fēng)險得分評估模型，式為：OS-risk score=0.034×LINC00460+0.057×RP11-284F21.9+0.426×RP11-539E17.5+0.227×CTD-2377D24.6+0.053×RP1-27K12.4+0.007×RP11-783K16.5。式中的6個系數(shù)分別為6個lncRNA在生存分析回歸模型中對應(yīng)的風(fēng)險系數(shù)，基因名代表該lncRNA的表達(dá)數(shù)據(jù)。利用上述公式對512個肺腺癌患者進(jìn)行風(fēng)險得分計算，并根據(jù)所有患者風(fēng)險得分的中位數(shù)將患者分為高低風(fēng)險兩組。這512個肺腺癌患者在預(yù)后模型中的分布以及生存時間的分布見圖3(a)。ROC曲線被用來評估該風(fēng)險得分評估模型的預(yù)測性能，結(jié)果顯示該模型在生存時間為1年、3年和5年時間的AUC值分別為0.68、0.7和0.71,這說明在5年的生存期內(nèi)，隨著生存時間越來越長，該風(fēng)險評估模型對于患者生存預(yù)后的準(zhǔn)確性也越來越好，見圖3(b)。同時也采用KM生存曲線對這些患者的生存情況進(jìn)行了分析，結(jié)果顯示低風(fēng)險組患者的死亡率顯著低于高風(fēng)險組(HR=2.72,95%CI=2.13～3.48,P<0.000 1)，見圖3(c)。

圖3 基于6個lncRNA的風(fēng)險得分對肺腺癌患者進(jìn)行預(yù)測的結(jié)果(a)風(fēng)險得分分布情況和肺腺癌患者的生存狀態(tài)；(b)整體生存情況預(yù)測的風(fēng)險得分ROC曲線分析；(c)高風(fēng)險組(n=256)和低風(fēng)險組(n=256)患者整體生存情況的KM分析Fig.3 The result of predicting the overall survival of patients with LUAD based on the risk score of the 6 lncRNAs (a)The risk score distribution and patients’ survival status,(b)The receiver operating characteristic (ROC) analysis of the risk score for predicting the overall survival,(c)Kaplan-Meier analysis of patients’ overall survival in the high risk (n=256) and low-risk (n=256)

表1 與肺腺癌具有顯著預(yù)后關(guān)系的6個lncRNATab.1 The 6 lncRNAs significantly associated with prognosis of LUAD patients

3.4 功能富集分析結(jié)果

在共表達(dá)基因分析中，與6個lncRNA至少1個相關(guān)的蛋白質(zhì)編碼基因就被認(rèn)為是顯著相關(guān)的共表達(dá)基因。通過皮爾森相關(guān)分析計算，一共得到了52個共表達(dá)基因,這些基因主要與LINC00460、RP11-284F21.9、RP11-539E17.5、CTD-2377D24.6、RP11-783K16.5這5個lncRNA相關(guān)，見表2。采用Metascape對52個共表達(dá)基因進(jìn)行Go和KEGG功能富集分析，見圖4。圖4中每個條形代表富集到的一簇功能，可以看到共表達(dá)基因主要與輔因子代謝過程、類固醇代謝過程以及細(xì)胞酮等代謝過程有關(guān)，而這些代謝過程是疾病產(chǎn)生和發(fā)展過程中的重要影響因素。

表2 5個lncRNA的共表達(dá)基因Tab.2 The co-express genes of the 5 LncRNAs

圖4 共表達(dá)基因功能富集分析(a)共表達(dá)基因功能富集熱圖；(b)共表達(dá)基因富集GO的詳細(xì)網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 Functional enrichment analysis of the co-express genes (a) Heatmap of the functional enrichment analysis of the co-expressed genes,(b)Detailed net structure of enriched GO term of the co-expressed genes

4 結(jié) 語

基于TCGA數(shù)據(jù)庫中肺腺癌的表達(dá)數(shù)據(jù)，利用基因差異分析，單樣本網(wǎng)絡(luò)方法以及生存分析識別出與肺腺癌相關(guān)的6個lncRNA生物標(biāo)志物。利用單樣本網(wǎng)絡(luò)識別lncRNA 生物標(biāo)志物，是對樣本個體特定條件下動態(tài)分子網(wǎng)絡(luò)進(jìn)行分析研究，對于揭示復(fù)雜疾病在網(wǎng)絡(luò)水平的LncRNA分子機(jī)制有著積極的指導(dǎo)意義，也為從系統(tǒng)水平識別人類復(fù)雜疾病相關(guān)LncRNA生物標(biāo)志物提供了一個新視角。