蒙俊樺，郎梅，康冉，張茵，唐飛，郭志云

西南交通大學(xué) 生命科學(xué)與工程學(xué)院，四川 成都 610031

2018年的癌癥統(tǒng)計(jì)分析報(bào)告顯示，乳腺癌發(fā)病率為24.2%，病亡率為15%，已經(jīng)成為女性頭號殺手[1]。microRNA（miRNA）是一段保守的非編碼RNA片段，長度約為22 nt，通過結(jié)合靶基因的3'非翻譯區(qū)（UTR）而降解或抑制基因表達(dá)[2]。先前大量研究表明，miRNA失調(diào)顯著參與乳腺癌的發(fā)生與發(fā)展[3]。后續(xù)研究發(fā)現(xiàn)，轉(zhuǎn)錄因子可以與miRNA的啟動(dòng)子區(qū)結(jié)合，從而參與miRNA的調(diào)控。因此，解析乳腺癌中轉(zhuǎn)錄因子、miRNA與其靶基因間的調(diào)控網(wǎng)絡(luò)關(guān)系，對于從系統(tǒng)生物學(xué)角度分析該疾病的發(fā)病機(jī)理，以及篩選乳腺癌相關(guān)致病miRNA有重要意義。

本研究整合TCGA、ENCODE、Fantom和GTEx等多組學(xué)數(shù)據(jù)，構(gòu)建轉(zhuǎn)錄因子-miRNA-基因調(diào)控網(wǎng)絡(luò)，篩選出乳腺癌致病miRNA，并分析了這些miRNA的功能富集和生存情況，為探究miRNA調(diào)控網(wǎng)絡(luò)在乳腺癌中的作用提供參考依據(jù)。

1 材料與方法

1.1 材料

miRNA表達(dá)量數(shù)據(jù)與臨床數(shù)據(jù)下載自TCGA（The Cancer Genome Atlas）數(shù)據(jù)庫[4]；miRNA的轉(zhuǎn)錄起始位點(diǎn)（transcription start site，TSS）數(shù)據(jù)下載自 Fantom（Functional Annotation of the Mouse）數(shù)據(jù)庫[5]；轉(zhuǎn)錄因子的ChIP-seq數(shù)據(jù)下載自UCSC（University of California Santa Cruz）數(shù)據(jù)庫[6]Uniform track和Txn track；miRNA靶基因數(shù)據(jù)下載自實(shí)驗(yàn)證實(shí)的miRNA靶基因數(shù)據(jù)庫Tarbase[7]、mir-Tarbase[8]和預(yù)測數(shù)據(jù)庫miRanda[9]、miRDB[10]、TargetScan[11]；包括轉(zhuǎn)錄因子在內(nèi)的基因表達(dá)量下載自 GTEx（Genotype-Tissue Expression）[12]、HPA（Human Protein Altas）[13]和 Encyclopedia of DNA Elements（ENCODE）[14]數(shù)據(jù)庫。

1.2 乳腺癌致病miRNA的識別與分析

從TCGA數(shù)據(jù)庫得到104個(gè)正常樣本和1103個(gè)腫瘤樣本的miRNA表達(dá)量數(shù)據(jù)[15]，通過計(jì)算CPM（count per million）進(jìn)行數(shù)據(jù)歸一化。將60%以上的樣本中均有表達(dá)且CPM＞0.5的miRNA作為乳腺癌中有效表達(dá)miRNA[16]，有效表達(dá)miRNA正常與腫瘤兩者表達(dá)量相差2倍以上認(rèn)為是存在顯著差異表達(dá)的miRNA，并根據(jù)下式，在腫瘤樣本中計(jì)算得到變異系數(shù)（CV）：

其中，M為乳腺癌miRNA在各樣本間表達(dá)的標(biāo)準(zhǔn)差，N為乳腺癌miRNA在各樣本間表達(dá)的平均值。取變異系數(shù)排名前15%的miRNA作為顯著差異表達(dá)的乳腺癌致病miRNA，并應(yīng)用SPSS軟件對這部分miRNA正常樣本與腫瘤樣本的表達(dá)量做主成分分析（principalcomponentsanalysis，PCA）[17]。

1.3 乳腺癌致病miRNA調(diào)控網(wǎng)絡(luò)構(gòu)建與分析

依據(jù)Fantom5識別的miRNA的TSS，定義其上、下游5 kb范圍內(nèi)有轉(zhuǎn)錄因子結(jié)合確定為存在對應(yīng)的轉(zhuǎn)錄因子-miRNA調(diào)控關(guān)系。為了降低miRNA靶基因預(yù)測假陽性，將miRNA預(yù)測靶基因數(shù)據(jù)庫 miRanda、miRDB、TargetScan取交集，取Tarbase與mirTarBase實(shí)驗(yàn)證實(shí)的miRNA靶基因數(shù)據(jù)并集作為miRNA-基因調(diào)控關(guān)系。轉(zhuǎn)錄因子-基因調(diào)控關(guān)系來自BioGrid數(shù)據(jù)庫中實(shí)驗(yàn)證實(shí)的蛋白質(zhì)與蛋白質(zhì)相互作用。使用Cytoscape Network Analyzer[18]計(jì)算miRNA調(diào)控網(wǎng)絡(luò)的度中心性與聚類系數(shù)。

1.4 乳腺癌致病miRNA功能富集與生存分析

采用DAVID對上述篩選出的miRNA靶基因進(jìn)行 GO（Gene Ontology）[19]和 KEGG（Kyoto Encyclopedia of Genes and Genomes）分析[20]，將 GO 與KEGG篩選條件為P＜0.05作為有效結(jié)果。獲取TCGA乳腺癌miRNA臨床數(shù)據(jù)，提取隨訪起止周期及存活數(shù)據(jù)，匹配miRNA表達(dá)量樣本，取乳腺癌樣本中miRNA表達(dá)中位數(shù)區(qū)分高低表達(dá)，構(gòu)建Cox回歸模型做生存曲線（Kaplan-Meier曲線），并確定P＜0.05為顯著性閾值。

2 結(jié)果

2.1 乳腺癌致病miRNA的識別與分析

圖1 正常組與癌癥組的變異系數(shù)、表達(dá)量與主成分分析

由于正常乳腺與乳腺癌中顯著差異表達(dá)的miRNA在腫瘤的發(fā)生發(fā)展中起關(guān)鍵作用[21]，篩選具有顯著差異表達(dá)的miRNA尤為重要。本研究共得到132個(gè)乳腺癌與正常乳腺顯著差異表達(dá)miRNA。對這132個(gè)miRNA進(jìn)行變異系數(shù)分析（圖1A），以識別乳腺癌中樣本變化差異大的miRNA，這些miRNA往往意味著在乳腺癌中高度失調(diào)。最終，篩選出19個(gè)具有顯著差異表達(dá)且樣本表達(dá)變化差異大的miRNA（圖1B），包括乳腺癌上調(diào) miRNA（hsa-mir-508、hsa-mir-509-3、hsamir-509-1、hsa-mir-184、hsa-mir-1248、hsa-mir-577、hsa-mir-153-2、hsa-mir-20b、hsa-mir-153-1、hsa-mir-9-1、hsa-mir-9-2、hsa-mir-9-3）、下調(diào)miRNA（hsa-mir-1-1、hsa-mir-1-2、hsa-mir-133a-1、hsa-mir-133a-2、hsa-mir-451a、hsa-mir-204、hsa-mir-144（表1）。為探究這19個(gè)顯著差異表達(dá)miRNA在正常樣本與腫瘤樣本中是否具有顯著特征，采用PCA進(jìn)行評估。結(jié)果表明，正常樣本與癌癥樣本各自有明顯的聚集，證明以正常與癌癥為特征量能有效區(qū)分這部分miRNA（圖1C）。

表1 19個(gè)差異表達(dá)miRNA

2.2 乳腺癌致病miRNA調(diào)控網(wǎng)絡(luò)構(gòu)建與分析

miRNA可被轉(zhuǎn)錄因子調(diào)控而作用于脆性位點(diǎn)或腫瘤發(fā)生相關(guān)區(qū)域，并且在乳腺癌中致病性表達(dá)，對下游基因產(chǎn)生調(diào)控作用并影響腫瘤發(fā)生、發(fā)展與轉(zhuǎn)移等生物過程。為此，我們將篩選出的19個(gè)乳腺癌致病miRNA進(jìn)行調(diào)控網(wǎng)絡(luò)構(gòu)建，并用度中心性與聚類系數(shù)評估m(xù)iRNA在網(wǎng)絡(luò)中的貢獻(xiàn)。其中，度中心性代表元件重要程度，聚類系數(shù)量化元件聚集程度。

為了識別調(diào)控核心miRNA，我們篩選了在網(wǎng)絡(luò)中與轉(zhuǎn)錄因子、基因均有調(diào)控的miRNA，最終得到由5個(gè)miRNA、8個(gè)轉(zhuǎn)錄因子、130個(gè)基因構(gòu)成的262個(gè)乳腺癌核心調(diào)控網(wǎng)絡(luò)（圖2A）。度中心性是網(wǎng)絡(luò)中元件與元件相關(guān)聯(lián)程度的體現(xiàn)，其值越大表明元件處于網(wǎng)絡(luò)越中心位置，元件處于中心位置則大概率為hub元件，hub元件所形成的網(wǎng)絡(luò)稱為hub網(wǎng)絡(luò)。圖2B是5個(gè)miRNA的度中心性與元件個(gè)數(shù)的關(guān)系，本網(wǎng)絡(luò)度中心性自大到小依次為 hsa-mir1-1、hsa-mir-1-2、hsa-mir-184、hsa-mir-1248、hsa-mir-144。依據(jù)冪律定義，度中心性與元件數(shù)的自然對數(shù)分布滿足線性關(guān)系。圖2B中紅色線滿足線性關(guān)系，表明本網(wǎng)絡(luò)服從冪律分布，是無標(biāo)度網(wǎng)絡(luò)，即網(wǎng)絡(luò)中元件與元件調(diào)控關(guān)系分布不均勻，少數(shù)元件可調(diào)控多個(gè)元件。

聚類系數(shù)可反映相關(guān)元件間調(diào)控關(guān)系的聚集程度，圖2C是相鄰元件數(shù)的自然對數(shù)與聚類系數(shù)平均值的關(guān)系，聚類系數(shù)平均值越小，表明兩元件間產(chǎn)生調(diào)控關(guān)系的平均概率越大。同時(shí)，聚類系數(shù)平均值也可以反映哪些元件可能會形成模塊。元件數(shù)的自然對數(shù)與平均聚類系數(shù)服從冪律分布，hsa-mir-1-1、hsa-mir-1-2是網(wǎng)絡(luò)中的hub元件并形成hub網(wǎng)絡(luò)，hsa-mir-1248、hsa-mir-144形成獨(dú)立于hub網(wǎng)絡(luò)的閉合網(wǎng)絡(luò)。結(jié)合圖2B、C，確定hsa-mir-1248、hsa-mir-144構(gòu)成網(wǎng)絡(luò)中最簡單的調(diào)控模式（僅由miRNA、轉(zhuǎn)錄因子與基因組成）。雖然這2個(gè)miRNA獨(dú)立于hub網(wǎng)絡(luò)，但其自身也是顯著差異表達(dá)的miRNA，暗示可能對研究乳腺癌的發(fā)生與發(fā)展具有重要意義。

2.3 乳腺癌致病miRNA功能富集與生存分析

為進(jìn)一步了解上述5個(gè)miRNA靶點(diǎn)在乳腺癌中參與的生物進(jìn)程與信號通路，采用DAVID進(jìn)行GO與KEGG分析。GO分析表明這些miRNA靶基因顯著參與細(xì)胞周期、細(xì)胞分化、細(xì)胞生長、轉(zhuǎn)移等轉(zhuǎn)錄后調(diào)控的腫瘤相關(guān)生物進(jìn)程（圖3A）；KEGG分析發(fā)現(xiàn)這些miRNA的靶基因顯著參與癌癥轉(zhuǎn)錄失調(diào)信號通路，如FoxO信號通路、p53信號通路、基因監(jiān)測通路等（圖3B）。

為了篩選出臨床上顯著影響生存的miRNA，對顯著差異表達(dá)的miRNA建立Cox回歸模型并分析，發(fā)現(xiàn)hsa-mir-144、hsa-mir-133a-2與乳腺癌患者生存顯著相關(guān)（P＜0.05）。由 hsa-mir-144與hsa-mir-133a-2生存曲線可知表達(dá)與生存率成反比（圖3C、D），暗示乳腺癌患者h(yuǎn)sa-mir-144與hsa-mir-133a-2高表達(dá)意味著生存率顯著下降。

3 討論

本研究整合miRNA表達(dá)量、變異系數(shù)與PCA篩選出乳腺癌中19個(gè)顯著差異表達(dá)的miRNA；構(gòu)建乳腺癌中轉(zhuǎn)錄因子-miRNA-基因調(diào)控網(wǎng)絡(luò)，獲得miRNA乳腺癌調(diào)控網(wǎng)絡(luò)262個(gè)，涉及miRNA 5個(gè)、轉(zhuǎn)錄因子8個(gè)、基因130個(gè)；結(jié)合功能富集分析可知這5個(gè)miRNA靶基因與細(xì)胞周期、細(xì)胞分化、細(xì)胞生長、轉(zhuǎn)移等轉(zhuǎn)錄后調(diào)控生物進(jìn)程相關(guān)，并且與癌癥轉(zhuǎn)錄失調(diào)信號通路、FoxO信號通路和p53信號通路等癌癥信號通路高度相關(guān)；由Cox回歸模型得知hsa-mir-144與hsa-mir-133a-2顯著與乳腺癌生存相關(guān)。

乳腺癌中各元件相互作用關(guān)系是探究乳腺癌發(fā)病機(jī)制與治療的關(guān)鍵，元件本身特性對網(wǎng)絡(luò)的影響有重要意義，網(wǎng)絡(luò)中的hub元件是影響最大的因素之一。本研究為探討乳腺癌發(fā)病機(jī)制提供了理論依據(jù)與數(shù)據(jù)基礎(chǔ)，且隨著miRNA調(diào)控網(wǎng)絡(luò)的不斷完善，本研究所提供的方法可擴(kuò)展到其他疾病研究中，可為了解復(fù)雜疾病的發(fā)生提供參考。

圖3 5個(gè)miRNA的GO（A）、KEGG（B）分析以及hsa-mir-144（C）、hsa-mir-133a-2（D）的Kaplan-Meier曲線