乳腺癌致病microRNA調控網絡的識別與生物信息學分析

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

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

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

本研究整合TCGA、ENCODE、Fantom和GTEx等多組學數據，構建轉錄因子-miRNA-基因調控網絡，篩選出乳腺癌致病miRNA，并分析了這些miRNA的功能富集和生存情況，為探究miRNA調控網絡在乳腺癌中的作用提供參考依據。

1 材料與方法

1.1 材料

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

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

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

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

1.3 乳腺癌致病miRNA調控網絡構建與分析

依據Fantom5識別的miRNA的TSS，定義其上、下游5 kb范圍內有轉錄因子結合確定為存在對應的轉錄因子-miRNA調控關系。為了降低miRNA靶基因預測假陽性，將miRNA預測靶基因數據庫 miRanda、miRDB、TargetScan取交集，取Tarbase與mirTarBase實驗證實的miRNA靶基因數據并集作為miRNA-基因調控關系。轉錄因子-基因調控關系來自BioGrid數據庫中實驗證實的蛋白質與蛋白質相互作用。使用Cytoscape Network Analyzer[18]計算miRNA調控網絡的度中心性與聚類系數。

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

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

2 結果

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

圖1 正常組與癌癥組的變異系數、表達量與主成分分析

由于正常乳腺與乳腺癌中顯著差異表達的miRNA在腫瘤的發生發展中起關鍵作用[21]，篩選具有顯著差異表達的miRNA尤為重要。本研究共得到132個乳腺癌與正常乳腺顯著差異表達miRNA。對這132個miRNA進行變異系數分析（圖1A），以識別乳腺癌中樣本變化差異大的miRNA，這些miRNA往往意味著在乳腺癌中高度失調。最終，篩選出19個具有顯著差異表達且樣本表達變化差異大的miRNA（圖1B），包括乳腺癌上調 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）、下調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個顯著差異表達miRNA在正常樣本與腫瘤樣本中是否具有顯著特征，采用PCA進行評估。結果表明，正常樣本與癌癥樣本各自有明顯的聚集，證明以正常與癌癥為特征量能有效區分這部分miRNA（圖1C）。

表1 19個差異表達miRNA

2.2 乳腺癌致病miRNA調控網絡構建與分析

miRNA可被轉錄因子調控而作用于脆性位點或腫瘤發生相關區域，并且在乳腺癌中致病性表達，對下游基因產生調控作用并影響腫瘤發生、發展與轉移等生物過程。為此，我們將篩選出的19個乳腺癌致病miRNA進行調控網絡構建，并用度中心性與聚類系數評估miRNA在網絡中的貢獻。其中，度中心性代表元件重要程度，聚類系數量化元件聚集程度。

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

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

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

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

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

3 討論

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

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

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