基于生物信息學(xué)的腦膠質(zhì)瘤預(yù)后風(fēng)險(xiǎn)長鏈非編碼RNA篩選

俞盛健，麻懷露，陳王洋，丁曉飛，陳光，梁勇

（1.溫州醫(yī)科大學(xué) 第一臨床醫(yī)學(xué)院，浙江 溫州 325035；2.臺(tái)州學(xué)院 醫(yī)學(xué)院，浙江 臺(tái)州 318000；3.河北北方學(xué)院 研究生院，河北 張家口 075000）

膠質(zhì)瘤是頭部最常見的原發(fā)性惡性腫瘤[1]。目前，主要的腦膠質(zhì)瘤治療方式有手術(shù)、化療、放療、靶向治療等，但是膠質(zhì)瘤患者的預(yù)后仍較差，其中膠質(zhì)母細(xì)胞瘤平均生存時(shí)間不超過18個(gè)月[2]。因此，探討腦膠質(zhì)瘤發(fā)生發(fā)展的潛在分子機(jī)制來發(fā)現(xiàn)有效的診斷和治療靶點(diǎn)非常重要。長鏈非編碼RNA（long non-coding RNA，lncRNA）是一類長度超過200個(gè)堿基的非編碼RNA。雖然它不參與蛋白質(zhì)的編碼，但是近來許多研究發(fā)現(xiàn)lncRNA參與了體內(nèi)多種生命活動(dòng)，如胚胎干細(xì)胞的多樣性、細(xì)胞周期的調(diào)節(jié)，以及癌癥的發(fā)生發(fā)展等[3]。癌組織中l(wèi)ncRNA的差異表達(dá)可以通過調(diào)控DNA的突變，如改變?cè)鰪?qiáng)子和啟動(dòng)子的活性或者染色質(zhì)狀態(tài)廣泛影響轉(zhuǎn)錄來實(shí)現(xiàn)[4]。同樣，lncRNA的差異表達(dá)也可以用來預(yù)測(cè)因組織DNA突變而導(dǎo)致的癌癥發(fā)生。如今越來越多的研究證實(shí)lncRNA可以作為癌癥的診斷標(biāo)志物，甚至作為一種治療選擇[5]。例如在早期手術(shù)切除的乳腺癌中，lncRNA HOTAIR的過表達(dá)對(duì)于腫瘤轉(zhuǎn)移和整體生存率具有很高的預(yù)測(cè)作用[6]。lncRNA PCA3作為前列腺癌患者尿中的一類生物標(biāo)志物，與血清前列腺特異性抗原相比具有更高的特異性和預(yù)測(cè)價(jià)值，現(xiàn)已成為在前列腺癌診斷中一個(gè)非常有效可靠的指標(biāo)[7]。

癌癥基因圖譜（the Cancer Genome Atlas，TCGA）是開放的公共數(shù)據(jù)平臺(tái)，其包括30種癌癥11 000例患者信息，對(duì)于這些巨大數(shù)據(jù)的挖掘分析，有助于研究者全面研究各種癌癥發(fā)生的分子機(jī)制，確定新的治療靶點(diǎn)和與腫瘤相關(guān)的分子標(biāo)志物，從而提高對(duì)多種人類惡性疾病的診斷能力和治療效果[8]?；蛐徒M織表達(dá)（GTEx）數(shù)據(jù)庫研究來自449名生前健康的人類捐獻(xiàn)者的7 000多份尸檢樣本，GTEx對(duì)幾乎所有轉(zhuǎn)錄基因的基因表達(dá)模式進(jìn)行了觀察，從而能夠確定基因組中影響基因表達(dá)的特定區(qū)域。在本研究中，我們通過對(duì)聯(lián)合GTEx和TCGA數(shù)據(jù)庫中的腦膠質(zhì)瘤樣本進(jìn)行基因轉(zhuǎn)錄水平的差異表達(dá)分析，構(gòu)建了可用于預(yù)測(cè)患者預(yù)后的風(fēng)險(xiǎn)模型，并提取了可能與膠質(zhì)瘤惡性進(jìn)展有關(guān)的lncRNA靶點(diǎn)。

1 材料和方法

1.1 數(shù)據(jù)收集 從GDC Data Portal（https://portal.gdc.cancer.gov/）中下載人類腦膠質(zhì)瘤的RNA-Seq數(shù)據(jù)和臨床數(shù)據(jù)，時(shí)間截至2019-01-14。RNA-Seq數(shù)據(jù)包含lncRNA和mRNA，包括700個(gè)腦膠質(zhì)瘤樣本。按照TCGA官方指導(dǎo)說明中的方法，使用Data Transfer Toll進(jìn)行數(shù)據(jù)下載。從UCSC Xena（https://xena.ucsc.edu/）中下載GTEx對(duì)應(yīng)正常腦部組織表達(dá)數(shù)據(jù)，得到1 152個(gè)正常樣本。

1.2 差異基因的篩選 利用edgeR包對(duì)腦膠質(zhì)瘤以及相對(duì)應(yīng)的正常樣本進(jìn)行分析，篩選出差異表達(dá)的lncRNA和mRNA。篩選標(biāo)準(zhǔn)為P＜0.05，F(xiàn)old Change＞2。

1.3 生存分析以及多因素Cox模型的建立 根據(jù)lncRNA表達(dá)量的中位值將RNA表達(dá)量分為高表達(dá)組和低表達(dá)組，結(jié)合臨床信息，通過survival包對(duì)lncRNA做Kaplan-Meier生存分析并畫出生存曲線圖。使用survival包對(duì)lncRNA做單因素Cox分析來挑選生存相關(guān)基因（survival related genes，SRGs），并且根據(jù)P值排序篩選出P值最小的前20個(gè)lncRNA，對(duì)挑選出的20個(gè)lncRNA進(jìn)行多因素Cox模型的建立，其中剔除有共表達(dá)關(guān)系的lncRNA。根據(jù)風(fēng)險(xiǎn)值的中位數(shù)將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組并畫出風(fēng)險(xiǎn)生存曲線。最后通過ROC曲線評(píng)估多因素Cox模型的鑒別能力。

1.4 共表達(dá)基因預(yù)測(cè) 通過使用雙側(cè)Pearson相關(guān)系數(shù)檢驗(yàn)lncRNA與蛋白編碼基因的表達(dá)關(guān)系（|皮爾遜相關(guān)系數(shù)|＞0.4，P＜0.001）。符合條件的蛋白編碼基因被認(rèn)為是該lncRNA的共表達(dá)基因。

1.5 富集分析 根據(jù)前述所得的差異表達(dá)基因，通過DAVID在線軟件（https://david.ncif-crf.gov/）和京都基因與基因組百科全書KEGG通路數(shù)據(jù)庫（https://www.kegg.jp/）進(jìn)行功能注釋和通路富集分析。

1.6 蛋白相互作用網(wǎng)絡(luò)構(gòu)建 根據(jù)預(yù)測(cè)共表達(dá)基因中重復(fù)出現(xiàn)基因進(jìn)行蛋白相互作用網(wǎng)絡(luò)的構(gòu)建。利用在線工具STRING11（https://string-db.org/）構(gòu)建蛋白質(zhì)的相互作用網(wǎng)絡(luò)，統(tǒng)計(jì)互作網(wǎng)絡(luò)各基因鄰接節(jié)點(diǎn)數(shù)目。利用Cytoscape對(duì)獲得的蛋白互作網(wǎng)絡(luò)進(jìn)行可視化，并通過MCODE插件對(duì)蛋白相互作用網(wǎng)絡(luò)中各節(jié)點(diǎn)進(jìn)行關(guān)聯(lián)度分析（K值＞5）。

1.7 患者標(biāo)本和臨床評(píng)估 患者標(biāo)本和臨床評(píng)估來源于具有腦膠質(zhì)瘤組織學(xué)診斷，且在手術(shù)前未曾接受放療和化療的患者。從臺(tái)州學(xué)院附屬臺(tái)州市立醫(yī)院和臺(tái)州市中心醫(yī)院獲取5對(duì)腦膠質(zhì)瘤癌組織及癌旁組織。所有的樣本均在手術(shù)切除后立即冷凍在液氮中。研究通過臺(tái)州學(xué)院醫(yī)學(xué)院倫理委員會(huì)批準(zhǔn)，并獲得患者的知情同意。

1.8 熒光實(shí)時(shí)定量PCR TRIzol試劑提取總RNA，Takara試劑盒完成反轉(zhuǎn)錄和qPCR檢測(cè)，HOXA-AS2引物序列為5'-CCCGTAGGAAGAACCGATGA-3'（正）和5'-TTTAGGCCTTCGCAGACAGC-3'（反），GAPDH引物序列為5'-ACCACAGTCCATGCCATC-3'（正）和5'-TCCACCCTG TTGCTGTA-3'（反）。PCR擴(kuò)增循環(huán)：95 ℃ 10 min；40循環(huán)的95 ℃ 10 s，60 ℃ 60 s。通過ΔΔCT法半定量HOXA-AS2相對(duì)水平。

1.9 統(tǒng)計(jì)學(xué)處理方法 應(yīng)用R3.3.5軟件進(jìn)行相關(guān)圖形制作，差異顯著性的篩選閾值為log2foldchange（log2FC）不小于2，采用Kaplan-Meier法繪制生存曲線并應(yīng)用log-rank方法比較生存率，多變量因素利用Cox回歸進(jìn)行生存分析，雙側(cè)Pearson相關(guān)系數(shù)檢驗(yàn)lncRNA與蛋白編碼基因的表達(dá)關(guān)系。采用雙側(cè)檢驗(yàn)，P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 膠質(zhì)瘤中差異表達(dá)的lncRNA和mRNA 通過比較腦膠質(zhì)瘤樣本和正常腦組織樣本，獲得差異表達(dá)的lncRNA 424個(gè)（211個(gè)lncRNA上調(diào)，213個(gè)lncRNA下調(diào)；見圖1A），mRNA 3 827個(gè)（1 618個(gè)mRNA上調(diào)，2 209個(gè)mRNA下調(diào)；見圖1B），差異有統(tǒng)計(jì)學(xué)意義（P＜0.05）。

圖1 腦膠質(zhì)瘤和正常腦組織中差異表達(dá)RNA的火山圖

2.2 生存分析以及多因素Cox模型的建立 運(yùn)用單因素Cox分析探究差異表達(dá)的lncRNA與膠質(zhì)瘤患者生存期的關(guān)系，挑選出統(tǒng)計(jì)學(xué)最顯著的20個(gè)lncRNA（見表1）。利用這20個(gè)lncRNA，我們進(jìn)行了多因素Cox分析?？紤]到臨床應(yīng)用的經(jīng)濟(jì)性，我們剔除了有共表達(dá)關(guān)系的lncRNA，只保留其中的一個(gè)，最終得到包含9個(gè)lncRNA（見表2）的生存預(yù)后模型，并且進(jìn)行生存曲線分析（見圖2A）。為了測(cè)試模型的準(zhǔn)確性，我們計(jì)算了模型的AUC值并進(jìn)行ROC曲線分析（見圖2B），AUC值達(dá)0.907，說明該模型能夠相對(duì)準(zhǔn)確地預(yù)測(cè)患者的預(yù)后。

表1 經(jīng)單因素Cox模型挑選的lncRNA

2.3 膠質(zhì)瘤中顯著差異表達(dá)的lncRNA富集分析 對(duì)挑選出的9個(gè)lncRNA中的CRNDE、LINCOO994、LBX2-AS1、HOXA-AS2進(jìn)行蛋白編碼基因共表達(dá)關(guān)系的預(yù)測(cè)，以相關(guān)系數(shù)大于0.6為臨界點(diǎn)，并對(duì)所得基因進(jìn)行篩選，挑選出重復(fù)出現(xiàn)的基因進(jìn)行富集分析。GO功能和KEGG通路富集分析結(jié)果顯示（見圖3），對(duì)其富集注釋結(jié)果提示（見表3-4），這些基因功能主要富集在通道蛋白活性調(diào)控，分子黏附等方面；介導(dǎo)的信號(hào)轉(zhuǎn)導(dǎo)通路主要集中于免疫缺陷病毒感染，神經(jīng)信號(hào)轉(zhuǎn)導(dǎo)相關(guān)通路。

表2 經(jīng)多因素Cox模型挑選的lncRNA

2.4 蛋白相互作用網(wǎng)絡(luò)的構(gòu)建和關(guān)鍵基因的篩選 在預(yù)測(cè)的共表達(dá)基因中，構(gòu)建重復(fù)基因的蛋白相互作用網(wǎng)絡(luò)，并與差異表達(dá)基因比對(duì)，發(fā)現(xiàn)33個(gè)上調(diào)基因，35個(gè)下調(diào)基因（見圖4A）。對(duì)蛋白相互作用網(wǎng)絡(luò)中各節(jié)點(diǎn)的鄰接數(shù)進(jìn)行統(tǒng)計(jì)（見圖4B），其中KIF2C與CDC20的鄰接節(jié)點(diǎn)最多，分別是16個(gè)和13個(gè)。使用MCODE插件篩選出網(wǎng)絡(luò)中的關(guān)聯(lián)度最緊密的2個(gè)基因簇（見圖4C），其中一個(gè)基因簇包括KIF1A、KIF3A、KIF3C、SPTBN2、CAPZA1、KDELR2、KIF2C、KIF21B、TMED9、KDELR3、KDELR1。對(duì)其蛋白相互作用網(wǎng)絡(luò)關(guān)鍵基因簇功能進(jìn)行注釋，該基因簇功能主要富集在逆行小泡介導(dǎo)轉(zhuǎn)運(yùn)，微管轉(zhuǎn)運(yùn)，ATP結(jié)合和細(xì)胞胞質(zhì)中。另一個(gè)基因簇包括GABBR1、ADCY5、GNG12、HTR1A、ANXA1、SSTR2、GNG5，功能主要富集在細(xì)胞胞膜上。

圖2 Cox分析篩選出的差異基因以及構(gòu)建的生存模型

圖3 預(yù)測(cè)共表達(dá)mRNA的GO（A）及KEGG（B）通路富集結(jié)果

2.5 HOXA-AS2在腦膠質(zhì)瘤組織及正常腦組織中的表達(dá) 為進(jìn)一步驗(yàn)證生物信息學(xué)分析結(jié)果，我們收集5例臨床腦膠質(zhì)瘤樣本，抽提RNA，進(jìn)行Real time PCR定量檢測(cè)。結(jié)果如圖5所示，HOXA-AS2在腦膠質(zhì)瘤組織中表達(dá)水平明顯高于正常腦組織（見圖5）。

表3 GO富集注釋結(jié)果

表4 KEGG通路富集注釋結(jié)果

3 討論

膠質(zhì)瘤是最常見的腦腫瘤之一，在早期腦腫瘤中占50%～60%[9-10]。在WHO分類中，神經(jīng)膠質(zhì)瘤被歸類為I-IV級(jí)。I級(jí)和II級(jí)星型細(xì)胞瘤和少突神經(jīng)膠質(zhì)瘤是低級(jí)別膠質(zhì)瘤，III級(jí)和IV級(jí)星型細(xì)胞瘤和少突神經(jīng)膠質(zhì)瘤是高級(jí)別膠質(zhì)瘤[11]。目前有證據(jù)表明lncRNA作為一類新的非編碼基因調(diào)節(jié)因子，可能在膠質(zhì)瘤的廣泛的生物學(xué)過程中發(fā)揮重要作用，例如，通過充當(dāng)癌基因或腫瘤抑制基因，它有助于膠質(zhì)瘤的發(fā)生、進(jìn)展和其他惡性表型的發(fā)展[12]。在本研究中，我們使用TCGA膠質(zhì)瘤轉(zhuǎn)錄水平數(shù)據(jù)和臨床數(shù)據(jù)，分析篩選出了424個(gè)差異表達(dá)的lncRNA，再對(duì)這些差異表達(dá)lncRNA進(jìn)行單因素Cox分析，挑選出20個(gè)代表性的lncRNA，其中有研究發(fā)現(xiàn)DGCR9在胃癌組織中高表達(dá)，并且有利于胃癌細(xì)胞系A(chǔ)GS，MGC803的增殖、侵襲和糖代謝[13]。

HOTAIRM1是HOXA基因簇的反義鏈的轉(zhuǎn)錄產(chǎn)物，最近有研究顯示其與急性髓性白血病、結(jié)直腸癌、多形性膠質(zhì)母細(xì)胞瘤有關(guān)[14-16]。然而，HOTAIRM1在胃癌中又通過靶向miR-17-5p抑制胃癌細(xì)胞系的增殖和遷移并誘導(dǎo)其凋亡[17]。我們還通過多因素Cox分析構(gòu)建了與腦膠質(zhì)瘤患者生存相關(guān)的生存模型。通過繪制出ROC曲線和計(jì)算AUC值來確定該模型的準(zhǔn)確性。最終AUC值達(dá)0.907，說明該模型能夠相對(duì)準(zhǔn)確地預(yù)測(cè)患者的預(yù)后，并且，我們對(duì)篩選出的CRNDE、LINCOO994、LBX2-AS1、HOXA-AS2這4個(gè)lncRNA進(jìn)行了mRNA共表達(dá)的預(yù)測(cè)，利用這些mRNA的GO功能和KEGG通路富集結(jié)果來預(yù)測(cè)lncRNA可能參與的功能和通路。我們發(fā)現(xiàn)，這些基因主要在控制通道蛋白活性，分子黏附等功能富集。其中，有10個(gè)基因富集在控制鈣通道活性的功能里。鈣黏蛋白是鈣依賴性細(xì)胞-細(xì)胞黏附分子，在神經(jīng)系統(tǒng)發(fā)育成熟過程中起至關(guān)重要的作用，根據(jù)其結(jié)構(gòu)和功能，可分為12個(gè)家族，包括經(jīng)典鈣黏蛋白、橋粒鈣黏蛋白、T-鈣黏蛋白、7D-鈣黏蛋白、原鈣黏蛋白、CDH15和23、脂肪、Dachsous、Flamingo、Celsr、鈣調(diào)素和Ret[18]。在癌癥中，E-鈣黏蛋白作為侵襲抑制因子，常常被下調(diào)，而作為侵襲起始因子的N-鈣黏蛋白則被上調(diào)[19]。原鈣黏蛋白PCDH10甲基化抑制了腫瘤細(xì)胞的生長、遷移、侵襲和克隆形成[20]?；蜻€在免疫缺陷病毒感染，神經(jīng)信號(hào)轉(zhuǎn)導(dǎo)相關(guān)通路上富集。研究表明，HIV-1胞膜糖蛋白120可誘導(dǎo)膠質(zhì)瘤細(xì)胞分泌基質(zhì)金屬蛋白酶[21]。另外，已經(jīng)有報(bào)道稱CRNDE在非小細(xì)胞肺癌、肝癌、宮頸癌和腦膠質(zhì)瘤中起到促癌作用[22-25]。對(duì)中國膠質(zhì)瘤基因組圖譜數(shù)據(jù)庫的RNA測(cè)序數(shù)據(jù)進(jìn)行分析，發(fā)現(xiàn)LBX2-AS1與膠質(zhì)瘤患者的預(yù)后和惡性進(jìn)展有顯著關(guān)系[26]。HOXA-AS2不僅在急性髓系白血病細(xì)胞對(duì)阿霉素的耐藥性中起重要作用，也參與了肝癌的惡性進(jìn)展過程[27-28]。

圖4 預(yù)測(cè)共表達(dá)基因中重復(fù)基因的蛋白相互作用網(wǎng)絡(luò)及關(guān)鍵基因簇分析

圖5 HOXA-AS2在腦膠質(zhì)瘤中的表達(dá)情況

綜上所述，本研究利用生物信息學(xué)方法分析TCGA數(shù)據(jù)，篩選潛在的膠質(zhì)瘤中表達(dá)差異顯著的基因，通過Cox分析構(gòu)建生存模型，并篩選出了可能為膠質(zhì)瘤分子致病機(jī)制和分子靶向治療提供參考的幾個(gè)lncRNA和mRNA，包括CRNDE、LINCOO994、LBX2-AS1和HOXA-AS2。并在有限例數(shù)臨床患者膠質(zhì)瘤樣本HOXA-AS2表達(dá)水平和病理進(jìn)程相關(guān)性的分析中，進(jìn)一步證實(shí)HOXA-AS2可能為診斷預(yù)警分子，但其潛在與膠質(zhì)瘤分子機(jī)制的關(guān)系仍需臨床及基礎(chǔ)實(shí)驗(yàn)驗(yàn)證。