乳腺癌增强子-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调控网络与癌症的关联提供了理论和方法参考。