乳腺癌致病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曲线