劉俊彤，王培培，葉明全，黎青青

(1. 皖南醫(yī)學(xué)院醫(yī)學(xué)信息學(xué)院，安徽 蕪湖 241002；2. 皖南醫(yī)學(xué)院健康大數(shù)據(jù)挖掘與應(yīng)用研究中心，安徽 蕪湖 241002)

乳腺癌是全球范圍內(nèi)女性常見的癌癥之一，約占全球新發(fā)癌癥病例的11.7%，也是造成女性患癌死亡的主要原因之一[1]。近年來，乳腺癌的發(fā)病率呈逐年上升的趨勢[2]。而早期的乳腺癌并沒有典型的臨床表現(xiàn)，一旦確診往往已是中晚期[3]，因此，通過尋找高特異性的乳腺癌生物標(biāo)志物，對改善患者的預(yù)后及生存率具有十分重要的意義。

隨著RNA高通測序技術(shù)的發(fā)展，轉(zhuǎn)錄組數(shù)據(jù)呈爆發(fā)式增長，信使RNA(messageRNA，mRNA)一直是主要研究重點(diǎn)，是蛋白質(zhì)的合成模板[4]。非編碼蛋白質(zhì)的RNA主要包括微小RNA(microRNA，miRNA)和長鏈非編碼RNA(long non-coding RNA，lncRNA)，二者可以通過與蛋白質(zhì)、DNA、RNA相互作用，參與基因表達(dá)的調(diào)控[5]，已有研究表明microRNA和lncRNA可以通過多種方式影響癌基因的表達(dá)[6-8]，但是，目前腫瘤預(yù)后相關(guān)生物標(biāo)志物的挖掘主要基于mRNA、miRNA、lncRNA中的一種分子標(biāo)簽，而忽略了mRNA-miRNA-lncRNA之間的相互調(diào)控作用以及它們對癌癥的共同作用。

本研究利用對癌癥公共數(shù)據(jù)庫(TCGA)的挖掘，分別篩選乳腺癌預(yù)后相關(guān)的mRNA、miRNA和lncRNA分子標(biāo)簽，綜合考慮每種標(biāo)簽的可靠性和標(biāo)簽之間的關(guān)聯(lián)，構(gòu)建基于mRNA-miRNA-lncRNA的競爭性內(nèi)源RNA(competing endogenous RNA，ceRNA)網(wǎng)絡(luò)，同時探究調(diào)控網(wǎng)絡(luò)中的mRNA與免疫系統(tǒng)的相關(guān)性，提出的核心ceRNA調(diào)控網(wǎng)絡(luò)(3個mRNA、1個lncRNA、2個miRNA)為乳腺癌的預(yù)后研究提供了新的研究方向，同時本研究的模型也可用于其他腫瘤預(yù)后相關(guān)標(biāo)志物的篩選參考。

1 材料與方法

1.1 TCGA數(shù)據(jù)下載和預(yù)處理 從TCGA數(shù)據(jù)庫(https://portal.gdc.cancer.gov/)獲得乳腺癌的轉(zhuǎn)錄組測序序列(RNA-seq)、miRNA測序序列(miRNA-seq)以及臨床信息，包括1 078例BRCA樣本和104個正常對照。乳腺癌的項(xiàng)目ID為TCGA-BRCA。據(jù)檢索得到的數(shù)據(jù)按以下條件篩選：①有完整且可供分析的生存數(shù)據(jù)(包括生存狀態(tài)、總生存時間)；②每一個樣本的mRNA數(shù)據(jù)、lncRNA數(shù)據(jù)以及miRNA數(shù)據(jù)完整。最后對符合標(biāo)準(zhǔn)患者的測序數(shù)據(jù)和生存信息下載和整合，以供進(jìn)一步分析。

1.2 差異基因篩選及可視化 利用R包“Deseq2”對基因數(shù)據(jù)的差異表達(dá)進(jìn)行分析，篩選出 |Log2 Fold change| >1和P值< 0.05的基因。使用R包“pheatmap”分別繪制差異表達(dá)mRNA、lncRNA、miRNA的聚類熱圖，使用R包“EnhancedVolcano”分別繪制差異表達(dá)mRNA、lncRNA、miRNA的火山圖。

1.3 ceRNA網(wǎng)絡(luò)構(gòu)建及可視化 Cytoscape是一個基于Java的可視化平臺，可以整合分子交互網(wǎng)絡(luò)與基因表達(dá)數(shù)據(jù)和其它分子的狀態(tài)信息。利用miRcode數(shù)據(jù)庫對特征miRNA調(diào)控的lncRNA進(jìn)行分析預(yù)測，利用miRTarBase、starbase和miRDB等數(shù)據(jù)庫對miRNA調(diào)控的靶基因進(jìn)行分析預(yù)測，篩選符合ceRNA機(jī)制的基因，篩選標(biāo)準(zhǔn)如下：利用Pearson相關(guān)系數(shù)計(jì)算mRNA與lncRNA的相關(guān)性，篩選出表達(dá)正相關(guān)的pairs(P值<0.01)；根據(jù)ceRNA網(wǎng)絡(luò)理論，lncRNA和mRNA競爭性結(jié)合miRNA，二者之間必須共享大量的miRNA，因此，通過超幾何檢驗(yàn)來判斷基因數(shù)據(jù)間是否存在顯著水平的多miRNA共享，設(shè)定超幾何檢驗(yàn)的P值<0.01滿足篩選條件。將得到的網(wǎng)絡(luò)關(guān)系導(dǎo)入Cytoscape中繪制ceRNA網(wǎng)絡(luò)調(diào)控關(guān)系圖。

1.4 ceRNA網(wǎng)絡(luò)核心RNA預(yù)后分析 構(gòu)建ceRNA網(wǎng)絡(luò)后，為了進(jìn)一步明確網(wǎng)絡(luò)中的差異mRNA、lncRNA和miRNA與乳腺癌預(yù)后的關(guān)系，我們采用R包“Survival”進(jìn)行單因素COX回歸分析，作出Kaplan-Meier生存曲線，并設(shè)定P值<0.05具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 乳腺癌預(yù)后相關(guān)的mRNA分子標(biāo)簽 從TCGA下載乳腺癌患者的RNA測序序列(RNA-seq)以及臨床信息，挑選其中RNA-seq數(shù)據(jù)中“group”列為protein_coding的數(shù)據(jù)作為編碼蛋白的mRNA表達(dá)數(shù)據(jù)，使用“Deseq2”進(jìn)行差異表達(dá)分析，篩選log2FC絕對值>1和Pvalue<0.01認(rèn)為是差異表達(dá)基因。與癌旁正常組織比較，共篩選出2 767個編碼蛋白差異表達(dá)mRNA，其中，上調(diào)1 125個，下調(diào)1 642個。使用R包“pheatmap”和“EnhancedVolcano”對已篩選差異表達(dá)基因的熱圖和火山圖分別進(jìn)行繪制，結(jié)果見圖1，圖1A為熱圖，紅色表示高表達(dá)，藍(lán)色表示低表達(dá)，橫軸代表樣本，樣本聚類顯示在橫軸上方，縱軸代表mRNA。圖1B為火山圖，Log2值>1的基因以藍(lán)色標(biāo)記；Log2值>2的基因以紅色標(biāo)記。為進(jìn)一步探究差異基因可能涉及的生物學(xué)功能和通路，我們進(jìn)行了通路富集分析。富集分析結(jié)果如圖2所示，乳腺癌組織中上調(diào)mRNA富集通路的前5名是細(xì)胞周期、細(xì)胞分裂、含CENPA核小體在著絲粒上的沉積、細(xì)胞周期過程調(diào)控、癌癥中的視網(wǎng)膜母細(xì)胞瘤基因；乳腺癌組織中下調(diào)mRNA富集通路的前5名是血管生成、肌肉結(jié)構(gòu)發(fā)育、細(xì)胞成分運(yùn)動的正向調(diào)節(jié)、循環(huán)系統(tǒng)過程、阿米巴型細(xì)胞遷移。利用單因素COX比例風(fēng)險回歸模型分析與生存期顯著相關(guān)的mRNA，Log-rank檢驗(yàn)結(jié)果的P值<0.05的共306個基因，其中，與生存期正相關(guān)103個，負(fù)相關(guān)203個。表1顯示與生存期相關(guān)性最高的前20個mRNA。

圖1 差異表達(dá)的mRNA可視化

圖2 乳腺癌差異表達(dá)上調(diào)和下調(diào)mRNA通路富集結(jié)果

表1 與乳腺癌生存期相關(guān)性最高的前20個mRNA

2.2 乳腺癌預(yù)后相關(guān)的LncRNA分子標(biāo)簽 與癌旁正常組織比較，共篩選出167個差異表達(dá)的lncRNA，其中上調(diào)61個，下調(diào)106個。差異表達(dá)lncRNA的熱圖和火山圖見圖3。利用單因素COX比例風(fēng)險回歸模型分析與生存期顯著相關(guān)的lncRNA，Log-rank檢驗(yàn)結(jié)果的P值<0.05的共18個基因，其中，與生存期正相關(guān)2個，負(fù)相關(guān)16個。表2顯示與生存期顯著相關(guān)性的前18個lncRNA。

圖3 差異表達(dá)的lncRNA可視化

表2 與乳腺癌生存期具有顯著相關(guān)性的18個lncRNA

2.3 乳腺癌預(yù)后相關(guān)的miRNA分子標(biāo)簽 與癌旁正常組織比較，共篩選出158個差異表達(dá)的miRNA，其中上調(diào)71個，下調(diào)87個。差異表達(dá)miRNA的熱圖和火山圖見圖4。利用單因素COX比例風(fēng)險回歸模型分析與生存期顯著相關(guān)的mRNA，Log-rank檢驗(yàn)結(jié)果的P值<0.05的共14個基因，其中，與生存期正相關(guān)4個，負(fù)相關(guān)10個。表3顯示與生存期顯著相關(guān)性的前14個miRNA。

圖4 差異表達(dá)的miRNA可視化

表3 與乳腺癌生存期具有顯著相關(guān)性的14個miRNA

2.4 乳腺癌ceRNA調(diào)控網(wǎng)絡(luò)的構(gòu)建 為了更好地理解與生存期密切相關(guān)的差異RNAs在乳腺癌中的作用，本研究通過ceRNA機(jī)制進(jìn)一步確認(rèn)這些RNAs的相互關(guān)系和生物學(xué)功能。借助miRcode數(shù)據(jù)庫對lncRNA和miRNA之間的關(guān)系進(jìn)行分析預(yù)測，miRNA和mRNA之間的作用關(guān)系通過starBase數(shù)據(jù)庫完成預(yù)測，整理結(jié)果后對網(wǎng)絡(luò)進(jìn)行可視化。將構(gòu)建的mRNA-miRNA-lncRNA關(guān)系結(jié)果導(dǎo)入Cytoscape軟件，繪制調(diào)控網(wǎng)絡(luò)圖，如圖5所示。其中圖5A為利用差異且生存相關(guān)mRNA為中心構(gòu)建的ceRNA，圖5B為利用差異且生存相關(guān)的mRNA和差異miRNA、lncRNA為中心構(gòu)建的ceRNA，圖5C為利用差異且生存相關(guān)的mRNA、miRNA和lncRNA為中心構(gòu)建的ceRNA。隨著篩選條件的提高，核心子網(wǎng)中的基因越來越少，最終為3個mRNA、1個lncRNA和2個miRNA。

圖5 乳腺癌的ceRNA網(wǎng)絡(luò)

2.5 乳腺癌ceRNA網(wǎng)絡(luò)核心RNA功能分析 我們使用R包“survival”對ceRNA核心網(wǎng)絡(luò)3個mRNA(KLF11、EDA、STXBP1)、1個lncRNA(XIST)和2個miRNA(hsa-miR-130a-3p、hsa-miR-195-5p)分別進(jìn)行生存分析，生存曲線結(jié)果見圖6。結(jié)果表明，hsa-miR-130a-3p、hsa-miR-195-5p兩個miRNA高表達(dá)，總生存期長；mRNA和lncRNA高表達(dá)，總生存期短。乳腺癌根據(jù)病理亞型分為Luminal A型、Luminal B型、Her2陽性和三陰性乳腺癌，因同時具有明確分子分型、基因測序數(shù)據(jù)、臨床數(shù)據(jù)的樣本量較少，部分分型甚至出現(xiàn)無樣本情況，我們將預(yù)后最差的三陰性乳腺癌作為一組，其他三種分型合并為一組，對核心ceRNA子網(wǎng)絡(luò)中的6個RNA進(jìn)一步做生存分析，生存曲線結(jié)果見圖7。結(jié)果能表明在不同乳腺癌分子分型中，已篩選核心子網(wǎng)RNA對生存周期的影響具有一致性。因?yàn)閙iRNA和lncRNA在體內(nèi)主要是調(diào)控mRNA作用，接下來我們對核心ceRNA網(wǎng)絡(luò)中的mRNA進(jìn)行了通路分析，如表4。結(jié)果表明，KLF11和EDA與免疫系統(tǒng)功能相關(guān)，而STXBP1主要與突觸和神經(jīng)遞質(zhì)功能相關(guān)。利用GEPIA工具，我們對3個mRNA(KLF11、EDA、STXBP1)組成的分子標(biāo)簽與免疫系統(tǒng)功能相關(guān)性進(jìn)行探究見圖8。結(jié)果表明KLF11、EDA、STXBP1組成的分子標(biāo)簽與效應(yīng)T細(xì)胞和駐留記憶T細(xì)胞存在顯著的正相關(guān)性。效應(yīng)T細(xì)胞分泌的穿孔素可以溶解腫瘤細(xì)胞膜或者被感染的細(xì)胞膜，進(jìn)而導(dǎo)致靶細(xì)胞形成孔洞并最終解體。而記憶T細(xì)胞保留了一些對抗原的免疫記憶，當(dāng)抗原再次入侵時，可迅速增值分化產(chǎn)生大量的效應(yīng)T細(xì)胞和記憶T細(xì)胞，再次破壞靶細(xì)胞并釋放出其中的抗原。

圖6 乳腺癌ceRNA核心子網(wǎng)RNA的生存分析

圖7 三陰性乳腺癌/非三陰性乳腺癌ceRNA核心子網(wǎng)RNA的生存分析

表4 乳腺癌ceRNA核心子網(wǎng)mRNA的通路分析

圖8 乳腺癌ceRNA核心子網(wǎng)mRNA標(biāo)簽與免疫細(xì)胞的相關(guān)性分析

3 討論

隨著乳腺癌發(fā)病率的逐年增加，尋找具有高特異性的乳腺癌生物標(biāo)志物，對改善乳腺癌患者的預(yù)后及生存率具有重要的意義[9-10]。近年來，基于公共癌癥數(shù)據(jù)庫，挖掘乳腺癌預(yù)后相關(guān)分子標(biāo)簽的研究逐漸成為國內(nèi)外研究熱點(diǎn)[11-13]。

在該研究中，通過對TCGA的挖掘，使用R包“Deseq2”進(jìn)行差異表達(dá)分析，篩選出差異表達(dá)基因(2767個mRNA、167個lncRNA、158個miRNA)，利用單因素COX比例風(fēng)險回歸模型對差異表達(dá)基因進(jìn)行生存分析，篩選出差異表達(dá)且生存相關(guān)基因(306個mRNA、18個lncRNA、14個miRNA)，借助miRcode與starBase數(shù)據(jù)庫，設(shè)定超幾何檢驗(yàn)的P值<0.01滿足篩選條件，將得到的網(wǎng)絡(luò)關(guān)系導(dǎo)入Cytoscape中繪制ceRNA網(wǎng)絡(luò)調(diào)控關(guān)系圖，最終篩選出3個mRNA(KLF11、EDA、STXBP1)、1個lncRNA(XIST)和2個miRNA(hsa-miR-130a-3p、hsa-miR-195-5p)參與構(gòu)建ceRNA網(wǎng)絡(luò)，使用R包“survival”對ceRNA核心網(wǎng)的6個RNA進(jìn)行生存分析，結(jié)果表明mRNA和lncRNA高表達(dá)、miRNA低表達(dá)，會促使乳腺癌進(jìn)一步增值和侵襲，相互調(diào)控關(guān)系見圖9，對核心ceRNA網(wǎng)絡(luò)中mRNA的進(jìn)一步通路分析和免疫分析也驗(yàn)證了篩選出的生物標(biāo)志物是有意義的。

圖9 乳腺癌ceRNA核心子網(wǎng)RNA的在乳腺癌增值、侵襲的作用

本研究獲得了一個有意義的乳腺癌ceRNA核心子網(wǎng)，其中心分子相互調(diào)節(jié)并在乳腺癌中發(fā)揮重要作用，有望成為潛在的治療靶點(diǎn)，同時本研究的策略為其它腫瘤預(yù)后相關(guān)標(biāo)志物的篩選提供了新的思路。進(jìn)一步擬通過生物學(xué)實(shí)驗(yàn)探索和驗(yàn)證這些RNA在乳腺癌各個分子分型發(fā)展的機(jī)制。