于曉慶

(上海應(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ù)。

長(zhǎng)非編碼RNA(long noncoding RNA,LncRNA)是長(zhǎng)度超過200個(gè)核苷酸的非編碼基因。過去很長(zhǎng)一段時(shí)間，這些非編碼基因曾被認(rèn)為是基因組的“暗物質(zhì)”，并不具有生物學(xué)功能。然而，隨著實(shí)驗(yàn)技術(shù)和計(jì)算機(jī)技術(shù)的發(fā)展，越來越多的研究表明，LncRNA通過不同的機(jī)制參與了細(xì)胞的幾乎整個(gè)生命過程，并在細(xì)胞分化、代謝、細(xì)胞增殖和凋亡、轉(zhuǎn)錄、染色質(zhì)水平的表觀遺傳學(xué)狀態(tài)調(diào)控等許多基本和關(guān)鍵的生物過程中都起到了非常重要的作用[4-5]。同時(shí)，很多研究發(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]。

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

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

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

2 方法

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

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

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

單樣本網(wǎng)絡(luò)是基于單個(gè)樣本對(duì)于給定的一組參考樣本的擾動(dòng)分析而構(gòu)建的，它可以精確地刻畫個(gè)體樣本的疾病狀態(tài)[10]。單樣本網(wǎng)絡(luò)的構(gòu)建需要兩組數(shù)據(jù)，一組是參照樣本數(shù)據(jù)，另一組是測(cè)試樣本數(shù)據(jù)。通常，這兩組數(shù)據(jù)有著共同的特征屬性。在本研究中，59個(gè)肺腺癌的癌旁樣本作為參照樣本組，535個(gè)肺腺癌的腫瘤樣本作為測(cè)試樣本組。首先，利用皮爾森相關(guān)系數(shù)(pearson correlation coefficient，PCC)計(jì)算每對(duì)lncRNA基因表達(dá)譜之間的相關(guān)性從而得到一個(gè)相關(guān)性網(wǎng)絡(luò)，稱為參照網(wǎng)絡(luò)(reference network)，記為Nr。然后，再將一個(gè)測(cè)試樣本s添加到該參考樣本組中，利用上述相關(guān)性的計(jì)算方法，可以得到一個(gè)新的相關(guān)性網(wǎng)絡(luò)，稱為擾動(dòng)網(wǎng)絡(luò)(perturbed network)，記為Np。最后，通過比較參照網(wǎng)絡(luò)和擾動(dòng)網(wǎng)絡(luò)之間的差異獲得一個(gè)差異網(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個(gè)腫瘤樣本依次作為測(cè)試樣本添加到參照樣本組中，可以得到535個(gè)單樣本網(wǎng)絡(luò)。容易看出，如果測(cè)試樣本在基因表達(dá)模式上與參照樣本組的表達(dá)模式相似，即使將測(cè)試樣本添加到參照樣本組中，任何2個(gè)分子之間的PCC變化將非常小。相反，如果測(cè)試樣本和參照樣本組在表達(dá)模式上存在明顯差異，則將測(cè)試樣本添加到參照樣本組后將引起明顯的變化?？梢钥闯?，單樣本網(wǎng)絡(luò)是由樣本s決定的，因此它可以用來作為描述樣本s的一個(gè)特定特征。我們對(duì)所有的腫瘤樣本構(gòu)建其對(duì)應(yīng)的單樣本網(wǎng)絡(luò)并進(jìn)行統(tǒng)計(jì)分析，篩選出可能與肺腺癌相關(guān)的lncRNA分子標(biāo)志物。

2.3 生存分析

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

2.4 功能分析

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

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

3 結(jié)果與討論

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

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

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

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

(1)

式中，SDi可通過以下方式計(jì)算：

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

(2)

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

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

3.3 生存分析結(jié)果

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

圖3 基于6個(gè)lncRNA的風(fēng)險(xiǎn)得分對(duì)肺腺癌患者進(jìn)行預(yù)測(cè)的結(jié)果(a)風(fēng)險(xiǎn)得分分布情況和肺腺癌患者的生存狀態(tài)；(b)整體生存情況預(yù)測(cè)的風(fēng)險(xiǎn)得分ROC曲線分析；(c)高風(fēng)險(xiǎn)組(n=256)和低風(fēng)險(xiǎn)組(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個(gè)lncRNATab.1 The 6 lncRNAs significantly associated with prognosis of LUAD patients

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

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

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