乳腺癌增強子-miRNA調控網絡識別與分析

錢明明，譚政堂，王文著，李昶鎣，邱正良，郭志云*

(1.西南交通大學生命科學與工程學院，中國四川成都610031；2.中國人民解放軍軍事醫學科學院實驗動物中心，中國北京100071)

增強子(enhancer)是真核生物基因組中的一類順式作用元件，它通過組織特異性的方式招募轉錄因子(transcription factor，TF)及其輔因子來正向調控基因表達。增強子作為一種非編碼序列，在組織中高度保守，增強子失活會導致包括腫瘤在內的多種疾病的發生[1]。微RNA(microRNA，mi-RNA)是一種長約22個核苷酸的非編碼RNA，在大約60%的人類基因表達中起關鍵作用，有些miRNA可能通過調控癌基因表達而成為致癌miRNA[2]。先前研究表明增強子參與了miRNA的合成與調控，并可通過促進Drosha/DGCR8募集和pri-miRNA加工來促進細胞特異性miRNA的產生[3]。研究證明，基因周圍富集增強子的miRNA往往與癌癥預后不良相關[3]。

乳腺癌是女性中最常見的癌癥，約占所有癌癥的16%[4]。目前，增強子與miRNA在乳腺癌中形成何種調控網絡尚不清楚。為此，我們通過癌癥基因組圖譜(The Cancer Genome Atlas，TCGA)數據庫[5]獲取了乳腺浸潤癌(breast invasive carcinoma，BRCA)相關miRNA-seq數據，識別了一系列在BRCA中由增強子調控的差異表達miRNAs，通過對這些miRNAs的靶基因進行基因本體論(Gene Ontology，GO)和京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)富集分析、蛋白質-蛋白質相互作用(protein-protein interaction，PPI)網絡分析以及生存分析，找到了BRCA中關鍵的增強子-miRNA-靶基因調控關系，這些調控關系中的增強子、miRNA、基因有望作為乳腺癌潛在的生物標志物，為乳腺癌遺傳機制研究提供了理論基礎。

1 材料與方法

1.1 材料

人類BRCA和正常組織的miRNA表達量數據(miRNA-seq)來源于在線數據庫TCGA(https://portal.gdc.cancer.gov/)。增強子位點數據(hg19)來源于Chen等[6]的研究。miRNA位點數據(hg19)來源于FANTOM5[7]。基因位點數據(hg19)來源于GENCODE[8]。

1.2 差異表達miRNA的識別

從112個BRCA樣本和104個正常組織樣本的miRNA-seq數據中，分別篩選在至少10%的癌癥/正常樣本中有表達量的miRNA為候選mi-RNA。以|log2FC|＞1(FC:fold change)和 P＜0.05(t檢驗)為閾值，篩選在癌癥中(相較于正常組織)顯著差異表達的miRNA。

1.3 增強子與miRNA相互作用的識別

參照先前關于增強子與miRNA調控關系的研究[3]，識別增強子與miRNA調控的公式如下:S=(EG-EM)/(EM+EG)，其中EG為增強子與其最近基因的距離，EM為增強子與其最近miRNA的距離。如果S＜0.2，則認為增強子調控其最近的miRNA。

1.4 生存分析

利用TCGA中BRCA患者生存時間的數據，使用Survival包(R 3.6.0)分析與癌癥患者生存相關的miRNA。采用Kaplan-Meier生存曲線[9]分析miRNA表達水平與BRCA患者預后的相關性。利用GEPIA[10](http://gepia.cancer-pku.cn/detail.php)分析與BRCA患者生存相關的基因。以P＜0.05為差異具有統計學意義。

1.5 miRNA靶基因的識別

利用miRWalk2.0[11](http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/)的Validated Target Module，得到增強子調控的差異表達miRNA的靶基因。Validated Target Module用于分析已被前人研究證實的miRNA-靶基因調控關系。

1.6 miRNA靶基因的功能分析

為了理解關鍵miRNA靶基因的生物學機制，利用clusterProfiler包[12](R 3.6.0)對靶基因進行GO和KEGG功能富集分析。以P＜0.05為差異具有統計學意義。

1.7 癌癥驅動靶基因的識別及PPI網絡構建

從DriverDBv3[13]下載BRCA的驅動基因數據集。利用Python 3.7編程獲取miRNA靶基因與BRCA驅動基因的交集，交集里的基因被認為是miRNA驅動靶基因。

使用STRING在線數據庫[14](https://string-db.org/)建立miRNA驅動靶基因的PPI網絡，以置信分數＞0.4作為篩選標準。生成的PPI網絡在Cytoscape軟件[15]中可視化。利用Cytoscape App cyto-Hubba計算網絡中各節點基因的得分。

2 結果與討論

2.1 BRCA中差異表達的miRNA

為了獲得在BRCA中增強子調控的差異表達miRNA，我們首先對BRCA中差異表達的miRNA進行了識別。通過分析BRCA miRNA-seq數據，共識別出322個顯著差異表達的miRNAs，其中包括227個上調miRNAs和95個下調miRNAs(|log2FC|＞1，P＜0.05)。上調倍數最大的 miRNA 是hsa-mir-105-1，其在BRCA中上調約67倍；下調幅度最大的miRNA是hsa-mir-486-1，其在BRCA中的表達量約為正常組織的1/26(圖1)。

圖1 差異表達的miRNAs紅色圓點代表在BRCA中上調的miRNAs，其中hsa-mir-105-1上調倍數最大，約67倍；紫色圓點代表在BRCA中下調的miRNAs，其中hsa-mir-486-1下調幅度最大，其在BRCA中的表達量約為正常組織的1/26；灰色圓點代表表達無顯著差異的miRNAs。Fig.1 Differentially expressed miRNAsRed dots represent the up-regulated miRNAs in invasive breast cancer.The most up-regulated miRNA is hsa-mir-105-1(～67 fold).Purple dots represent the down-regulated miRNAs.The most down-regulated miRNA is hsa-mir-486-1(～1/26 of the expression level in the normal tissue).Gray dots represent the miRNAs without significantly differential expression.

2.2 增強子-差異表達miRNA的調控關系

增強子的功能主要通過調控下游的基因(包括miRNA基因)來實現，因此增強子與腫瘤失調mi-RNA的調控關系可能在腫瘤中發揮重要作用。為此，根據先前關于增強子與miRNA調控關系的研究(見方法)，在1 170個BRCA增強子和322個差異表達miRNAs中識別了220對增強子-差異表達miRNA調控關系，涉及220個增強子和56個差異表達miRNAs，平均1個miRNA受到4個增強子的調控。通過對56個受增強子調控的差異表達miRNAs進行生存分析，發現其中6個miRNAs的表達水平顯著與BRCA患者生存相關(P＜0.05，圖2)。除了hsa-mir-195下調(2.90倍)外，其余5個 miRNAs:hsa-mir-671(2.17 倍)、hsa-mir-3619(5.80 倍)、hsa-mir-4664(11.05 倍)、hsa-mir-5003(2.71 倍)、hsa-mir-5691(3.47 倍)在 BRCA 中均上調(P＜0.05)。其中，miR-195、miR-671、miR-3619和miR-5003已經被證實與乳腺癌高度相關。miR-195通過調節乳腺癌中胰島素受體底物1(insulin receptor substrate 1，IRS1)的水平，抑制腫瘤生長和血管生成[16]。在福爾馬林固定石蠟包埋(formalinfixed paraffin-embedded，FFPE)的乳腺癌組織中，研究人員發現miR-671與其化學耐藥性相關，并認為miR-671可作為化學療法中的有效靶標以及乳腺癌的潛在生物標志物[17]。miR-3619在MCF-7乳腺癌細胞中的高表達可以阻礙磷脂酶D2(phospholipase D2，PLD2)的翻譯，發揮腫瘤抑制作用[18]。已知腫瘤抑制因子miR-34也在一些癌癥中上調，包括腎細胞癌、結腸癌和肝細胞癌，但其中的機制有待進一步研究[19]，miR-3619的情況可能與之類似。研究表明，在轉移性乳腺癌細胞中，如果miR-5003-3p受到抑制，那么癌細胞的遷移和侵襲也可能受到抑制[20]。因此，這6個miRNAs被識別為增強子-miRNA調控網絡中的關鍵miRNAs。

圖2 受增強子調控的6個差異表達miRNAs的生存分析Fig.2 Survival analysis of the six differentially expressed miRNAs regulated by enhancers

2.3 關鍵miRNA的靶基因識別與功能富集分析

miRNA的功能主要是通過結合下游靶基因的mRNA來實現。為此，我們通過miRWalk分析了這6個受增強子調控的關鍵miRNAs(hsa-mir-195、hsa-mir-671、hsa-mir-3619、hsa-mir-4664、hsa-mir-5003、hsa-mir-5691)的靶基因，共識別了1 172個實驗證實的靶基因。GO分析的結果表明，miRNAs的靶基因顯著富集在Wnt途徑(圖3A)。Wnt途徑是腫瘤關鍵信號通路，通常可介導和調節一系列生物進程，包括增殖、遷移、黏附和上皮-間質轉化(epithelial-mesenchymal transition，EMT)等，在許多人類癌癥中均異常激活，對于腫瘤的發生、生長和轉移非常重要[21]。KEGG富集分析的結果表明，miRNAs靶基因參與的通路最顯著富集在PI3K-Akt信號通路(PI3K-Akt signaling pathway)、癌癥中的蛋白多糖(proteoglycans in cancer)、乳腺癌(breast cancer)(圖 3B)。其中，PI3K-Akt信號通路與線粒體介導的乳腺癌細胞凋亡有關[22]。這些結果進一步說明，篩選出的這6個關鍵mi-RNAs與乳腺癌的發生發展高度相關。

圖3 6個miRNAs的靶基因功能富集分析(A)6個miRNAs靶基因的GO分析結果；(B)6個miRNAs靶基因的KEGG分析結果。Fig.3 Functional enrichment analysis of target genes of the six miRNAs(A)The GO analysis result；(B)The KEGG analysis result.

2.4 增強子調控miRNAs驅動靶基因識別及PPI網絡構建

為了進一步分析6個miRNAs的靶基因在BRCA中的潛在臨床作用，我們對這些靶基因進行了癌癥驅動基因分析。首先，通過搜索癌癥驅動基因數據庫DriverDBv3，獲取了BRCA驅動靶基因；為了進一步識別驅動靶基因中的關鍵基因，通過STRING獲取了驅動靶基因所對應蛋白質產物的PPI調控網絡；最后，根據Cytoscape App cytoHubba計算，識別了排名前20%的網絡關鍵節點:TP53(tumor protein P53)、CCND1(cyclin D1)、KRAS(KRAS proto-oncogene，GTPase)、CDKN1A(cyclin dependent kinase inhibitor 1A)、RPS6KB1(ribosomal protein S6 kinase B1)和CCNE2(cyclin E2)(圖4A)。這6個節點的基因已被證明與乳腺癌高度相關。研究表明，所有BRCA1缺陷型乳腺癌都含有TP53基因突變，這表明TP53與乳腺癌發展密切相關[23]。CCND1是細胞周期蛋白D家族的成員，研究顯示CCND1基因的擴增導致乳腺癌預后較差[24]。在乳腺癌細胞中，KRAS/MAPK信號通路的激活可以調節細胞增殖、分化、遷移和侵襲[25]。CDKN1A是RUSC1-AS1的靶基因，是參與RUSC1-AS1介導的乳腺癌進展的腫瘤抑制基因[26]。相關研究報道，RPS6KB1基因的表達降低可能與絕經前肥胖婦女乳腺癌風險降低有關[27]。另外，細胞周期蛋白E2(CCNE2)通過促進Hippo通路下游介質YAP的核定位和活性來調節乳腺癌細胞的運動性和侵襲性[28]。進一步的生存分析顯示，在這6個網絡關鍵節點中TP53、CCNE2的基因表達水平與BRCA患者生存時間顯著相關(P＜0.05，圖4B)。這表明在6個網絡關鍵性節點中，相較于其他節點，TP53和CCNE2可能在乳腺癌下游調控網絡中具有更重要的作用，因此TP53、CCNE2基因被進一步識別為增強子-miRNA調控的關鍵驅動靶基因(圖4A)。以上結果表明，增強子-mi-RNA調控網絡可能通過關鍵miRNAs靶向重要的癌癥驅動基因，從而在癌癥中發揮關鍵作用。

圖4 驅動靶基因蛋白質產物的PPI網絡及TP53和CCNE2的生存分析(A)驅動靶基因蛋白質產物參與的PPI網絡及增強子-miRNA-關鍵驅動靶基因調控關系。橢圓表示PPI網絡中的節點，節點紅色越深，其在網絡中的重要性越大，菱形表示調控相應基因的miRNA，長方形表示調控miRNA的增強子；(B)TP53、CCNE2基因在BRCA中的生存分析。Fig.4 PPI network of driver target genes and survival analysis of TP53 and CCNE2(A)The PPI network of the protein product of driver target genes and the regulatory relationship of enhancer-miRNA-key driver target gene.The ellipses represent the nodes.The redder the node，the greater its importance in the network.The diamonds represent the miRNAs that regulate the genes.The rectangles represent the enhancers that regulate miRNAs；(B)Survival analysis of TP53 and CCNE2 in invasive breast cancer.

3 結論

本研究通過差異表達分析，識別了在BRCA中受增強子調控的322個差異表達miRNAs，其中227個顯著上調，95個顯著下調。為了探究乳腺癌增強子與差異表達miRNA的調控網絡，識別了220對增強子-差異表達miRNA調控關系，包括220個增強子和56個差異表達miRNAs，其中6 個 miRNAs(hsa-mir-195、hsa-mir-671、hsa-mir-3619、hsa-mir-4664、hsa-mir-5003、hsa-mir-5691)的表達水平顯著與BRCA患者生存相關。對6個miRNAs的靶基因進行GO/KEGG功能富集分析，結果顯示其顯著富集在癌癥相關通路上，暗示這6個miRNAs可能通過調控關鍵靶基因參與癌癥相關通路。因此，我們進一步對6個miRNAs靶基因進行了癌癥驅動基因與PPI網絡分析，共識別了6個網絡關鍵基因，生存分析表明其中的TP53、CCNE2基因顯著與BRCA患者生存相關。綜上所述，通過差異表達分析與增強子-miRNA調控關系識別，本文構建了BRCA中增強子-miRNA調控網絡，同時通過進一步分析網絡中6個重要miRNAs的關鍵靶基因，獲得了可能在BRCA中有重要作用的增強子-miRNA-基因調控關系，即增強子chr17:6 925 152～6 925 653可能通過調控hsa-mir-195影響CCNE2基因；hsa-mir-5003可能受到5個增強子調控，從而影響其靶向的TP53基因。這些結果為進一步探究乳腺癌中增強子-miRNA調控網絡與癌癥的關聯提供了理論和方法參考。