基于生物信息学分析的子宫内膜异位症miRNAmRNA调控网络的初步构建*

姜丽，章蒙蒙，王思雪，曹乐，周雨馨，方小玲

（中南大学湘雅二医院 妇产科，湖南 长沙 410011）

子宫内膜异位症是育龄期妇女的多发病、常见病，其临床表现多种多样，且具有一定异质性，常见症状有疼痛、月经异常、不孕等。该疾病组织学分类虽然为良性，但是在临床行为学上却具有类似恶性肿瘤的特点，如种植、侵袭、远处转移、容易复发等[1]。子宫内膜异位症的发病机制至今仍未完全阐明，目前主要的理论学说包括异位种植学说、体腔上皮化生学说、诱导学说、遗传因素、免疫与炎症因素、在位内膜决定论等，但是其中任何一种学说都不能完全解释该疾病的所有方面[2]。

近年来基因芯片技术的发展和各大生物信息数据库的建立，为子宫内膜异位症的研究提供有利工具并开拓新的思路和方向。本研究将GEO 数据库中的mRNA 芯片和公共数据库中挖掘的miRNA 进行联合分析，构建miRNA-mRNA 网络，寻找子宫内膜异位症相关的miRNA 及其主要作用的靶基因，并探索子宫内膜异位症的分子发病机制。

1 材料与方法

1.1 基因表达谱数据来源

本研究从GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）检索并下载子宫内膜异位症相关的表达谱数据集GSE7305[3]，该数据集包含从10 例卵巢子宫内膜异位症患者获取的10 个卵巢异位内膜样本和10 个配对的在位内膜样本。所有样本均采用GPL570 芯片平台（Affymetrix human genome U133 Plus 2.0）进行差异表达基因的检测。

1.2 差异表达基因的分析

首先使用R 语言软件中Affy[4]程序包，读取芯片原始数据，然后采用RMA（robust multi-array average）法[5]对数据进行标准化预处理（包括背景矫正、归一化处理和表达值计算）。随后利用GPL570 平台注释文件对探针进行注释，去除没有匹配到基因名称的探针。当不同探针映射到同一基因时，取不同探针的平均值作为这个基因最终的表达值。最后，采用LIMMA[6]程序包（linear models for microarray analysis）筛选差异表达基因，将P≤0.05 和|log2 FC|≥2 作为筛选差异表达基因的阈值。筛选出差异表达基因后，使用R 软件中的Gplots、Pheatmap 等程序包绘制相应的热图和火山图。

1.3 差异表达基因的功能和通路富集分析

DAVID 在线网站（https://david.ncifcrf.gov/home/jsp）整合大规模的生物数据，是为大量基因提供功能注释的有效工具。采用DAVID 在线网站对筛选到的差异表达基因进行基因本体论（gene ontology, GO）注释（即基因功能富集分析）和KEGG（kyoto encyclopedia of genes and genomes）信号通路富集分析。采用Fisher检验和FDR 方法校验，将P＜0.05 作为富集分析的阈值，并按照P值大小升序排序。

1.4 miRNA 数据来源

人类microRNA（miRNA）疾病数据库HMDD v3.0（http://www.cuilab.cn/hmdd）是经过实验验证的人类miRNA 和疾病关联的数据库。以“endometriosis”为关键词在HMDD 中查找已有文献报道的和该疾病相关的miRNA，选取在组织中表达有差异且文献研究报道频次≥2 的相对高频miRNA 共13 个（异位内膜组织中高表达miRNA 6 个，低表达miRNA 7 个）作为本研究对象。

1.5 miRNA 靶基因预测

miRWalk 2.0（http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/）是目前可在线预测miRNA 靶基因的最大综合性数据库之一，提供来自人类、小鼠和大鼠的miRNA 靶基因信息。该数据库内容可分为预测模块和验证模块，其中预测模块提供了综合各大数据库预测miRNA靶基因的功能。本研究使用miRWalk 2.0数据库的预测模块获取miRNA 的预测靶基因，采用miRWalk、MicroT4、miRanda、RNA22、Target Scan 多种预测方法[7]完成预测。为降低假阳性率，同时被≥4 种预测工具标记的靶基因最终被选为miRNA 的预测靶基因。

1.6 构建miRNA-mRNA 网络

将miRNA 的预测靶基因和从GSE7305 分析获得的差异表达基因取得交集，通过综合分析得到miRNA 和mRNA 相互作用关系对。使用可视化软件Cytoscape v3.6.1 对miRNA-mRNA 调控网络进行绘图，并使用CytoHubba插件进行前20个枢纽基因及miRNA的筛选，构建相应的子网络。

2 结果

2.1 筛选出差异表达基因

通过R 语言软件对数据集GSE7305 进行芯片数据处理及差异分析，共得到655 个差异表达基因，其中包括343 个上调基因和312 个下调基因。依据上述结果在R 语言软件中绘制双向聚类热图，横坐标为组织样本，纵坐标为差异表达基因。根据聚类热图（见图1A）及火山图（见图1B）所示，不难发现来自异位内膜组织的样本和来自在位内膜组织的样本形成了相应的两大聚类群。

2.2 Go 和KEGG 富集分析结果

为明确差异表达基因的功能，本文通过DAVID 在线网站对筛选的差异表达基因进行GO 注释和KEGG信号通路富集分析，其中GO 注释主要分为生物进程、细胞组成和分子功能3 个部分。在生物进程中，基因主要富集于cell adhesion, extracellular matrix organization,mitotic nuclear division, cell division。在细胞组成，基因主要富集在extracellular space, extracellular matrix,proteinaceous extracellular matrix, extracellular exosome。在分子功能中，基因主要富集于heparin binding, integrin binding, sequence-specific DNA binding, serine-type carboxypeptidase activity。KEGG 信号通路富集分析提示，该差异表达基因主要富集于complement and coagulation cascades, p53 signaling pathway, cell adhesion molecules 等通路中。见图2。

图1 GSE7305 中异位内膜和在位内膜的mRNA 表达差异

图2 差异表达基因GO 和KEGG 富集分析

2.3 子宫内膜异位症相关miRNA

通过HMDD v3.0 查询经实验证实的与子宫内膜异位症相关的miRNA，初步发现在异位内膜组织中差异表达的miRNA 共有77 个，而有≥2 篇文献报道的miRNA 共有13 个。其中在异位内膜组织表达上调的miRNA 有6 个： hsa-miR-29c-3p、hsa-miR-145-5p、hsa-miR-100-5p、hsa-miR-125a-5p、hsa-miR-143-3p、hsa-miR-202-3p；而在异位内膜组织中表达下调的miRNA 有7 个： hsa-miR-141-3p、hsa-miR-200a-3p、hsa-miR-200b-3p、hsa-miR-424-5p、hsa-miR-20a-5p、hsa-miR-34c-5p、hsa-miR-449b-3p。运用miRWalk 2.0 预测模块，获得上述13 个miRNA 的预测靶基因。miRNA 查找及预测靶基因筛选结果见表1。

2.4 子宫内膜异位症miRNA-mRNA 调控网络

由于目前已知miRNA 对靶基因的调控作用主要是负向调控关系，因此将7 个下调的miRNA 预测靶基因同GSE7305 上调的差异表达基因（343 个）取得交集，共获得332 个靶基因；将6 个上调的miRNA预测靶基因同GEE7305 下调的差异表达基因（312 个）取得交集，共获得114 个靶基因。利用Cytoscape v3.6.1分别绘制下调及上调miRNA 及相应靶基因mRNA 的调控网络图（见图3A、B）。其中hsa-miR-20a-5p 调控的靶基因最多，共73 个，且下调的miRNA 调控的靶基因数均在30 个以上，推测在异位内膜组织中以miRNA 表达下调、相应靶基因表达上调为主。

表1 HMDD v3.0 及miRWalk 2.0 数据库中与子宫内膜异位症相关的miRNA

图3 异位内膜组织中差异表达的miRNAs-mRNAs 调控网络

使用CytoHubba 插件对下调miRNA-上调mRNA和上调miRNA-下调mRNA 2 个网络进行前20 个核心基因及核心miRNAs 的筛选，并构建相应的子网络（见图4）。下调miRNA-上调mRNA 网络的子网络中的核心miRNA 为hsa-miR-20a-5p、hsa-miR-424-5p、hsamiR-34c-5p、hsa-miR-141-3p、hsa-miR-200b-3p、hsamiR-200a-3p、hsa-miR-449b-3p，其调控的核心基因为PLXNA4、SGCD、ANKRD52、NFASC、DENND5B、CELF2、MYOCD、PRELP、TNS1、GATA6、FOXP1、PTGDR、NTRK2。上调miRNA-下调mRNA 网络的子网络中的核心miRNA 为hsa-miR-125a-5p、hsa-miR-29c-3p、hsa-miR-145-5p、hsa-miR-143-3p、hsa-miR-202-3p，其调控的核心基因为ADAM12、PGR、MMP11、FREM2、SHISA6、FXYD3、ITGA6、FRAS1、PAX8、HMGB3、FGFR3、ESR1、ERBB3、DIO2、CPM。通过观察不难发现，图4A 示miRNA 与靶基因的联系更为密切，而图4B 示miRNA 与靶基因之间及靶基因与靶基因之间联系较少，网络较为稀疏。

图4 异位内膜组织中差异表达的miRNAs-mRNAs 调控网络的子网络

3 讨论

本研究利用GEO 数据集GSE7305 分析子宫内膜异位症异位内膜组织和在位内膜组织的mRNA 差异表达情况，共筛选出655 个差异表达基因（343 个上调基因和312 个下调基因）。通过GO 和KEGG 信号通路富集分析发现，这些差异表达基因主要和细胞黏附、细胞增殖分裂、细胞外基质重构等生物学过程和功能有关，同时富集于补体信号通路、P53 信号通路、细胞黏附相关信号通路中。

既往研究表明，子宫内膜异位症是一种自身免疫性疾病，而补体激活导致的体液免疫异常是发病机制中的重要环节[8]。在子宫内膜异位症患者的异位病灶、外周血和腹腔液中出现过各种非器官特异性的自身抗体，如抗多核苷酸类抗体、抗组蛋白抗体、抗磷脂抗体等；还出现特异性自身抗体，如抗子宫内膜抗体和抗卵巢抗体。尤其是抗子宫内膜抗体，可以和异位或在位内膜产生抗原-抗体结合沉积于子宫和异位病灶中，通过激活补体，使患者血清及腹腔液中C3 补体相关蛋白水平升高，并通过激活一系列的免疫反应，导致患者产生更为广泛的细胞免疫和体液免疫异常[9]。本研究也发现，差异表达基因中补体相关蛋白C7、C3、C1R、C1S、CFB、CFH、CFD、PROS1、SERPINE1 的升高，尤其C7、C3 和CFH 表达差异倍数均在log2 FC ＞4 倍以上。

P53作为抑癌基因，在正常人体细胞活动中发挥着促进细胞凋亡、促进DNA 修复、避免受损DNA 堆积以维持基因组稳定的功能[10]。P53基因突变或失活，将通过影响下游增殖、凋亡等信号通路从而导致肿瘤。2000年，VOGELSTEIN 和LEVINE 等学者提出P53基因网络的概念，认为许多基因形成网络协同效应而调节细胞的生命活动，不应该孤立地观察单个基因的生物学功能[11]。本研究中富集在P53信号通路中的基因有CCNB1、CCNE2、CCNB2、CDK1、RRM2、CYCS、SERPINE1、PMAIP1、PERP、THBS1、GTSE1。其中，CCNB1、CCNB2、CCNB2、CDK1 也是影响细胞周期的重要分子。有研究发现，CCNB1 可通过与CDK1 结合形成成熟生长因子MPF，该Cyclin B-CDK1 复合物主要在G2 末期发挥作用，可引导细胞进入有丝分裂期。TANG 等[12]研究表明，子宫内膜异位症患者异位内膜CCNB1 表达量高于在位内膜，说明CCNB1 升高促进异位内膜增殖是重要的致病机制之一。但是其研究表明CDK1 表达无差异。同时YAMAMOTO 等[13]研究发现在与子宫内膜异位症相关的卵巢肿瘤中，细胞周期相关蛋白Cyclin A 和Cyclin E 的表达均升高。

细胞黏附分子是一类介导细胞与细胞、细胞与细胞外基质间黏附作用的膜表面糖蛋白，参与细胞的生长分化、组织修复、肿瘤的侵袭与转移等过程[14]。已有研究表明多种细胞黏附分子与子宫内膜异位症发生、发展密切相关，其通过调节异位内膜细胞与腹膜的黏附，参与异位内膜的种植[15]。本研究中细胞黏附分子NFASC、CLDN11、CDH3、VCAM1、NCAM1、CD22、SELE、NEGR1 表达均上调。其中血管细胞黏附分子VCAM1 在子宫内膜异位症中的研究最多，在子宫内膜异位症患者的在位内膜[16]、异位内膜[17]、腹膜组织[16-17]及血清中[17-18]均有VCAM1 的表达升高，提示VCAM1 参与子宫内膜异位症形成的多个病理过程，并有希望作为子宫内膜异位症的诊断标志物。

miRNA 通过作用于靶基因的mRNA 3'-非编码区（3'-untranslated region, 3'-UTR）使得靶基因mRNA 表达沉默，在转录水平上参与基因的表达调控。已有研究表明，子宫内膜异位症患者的异位内膜、在位内膜和血清中均存在miRNA 的差异表达[19]，且这些差异表达的miRNA 通过调控缺氧、增殖、凋亡、血管新生、细胞外基质重塑等病理过程介导子宫内膜异位症的发生、发展[20]。本研究采用差异表达基因分析、miRNA靶基因预测等生物信息学方法，对筛选的异位内膜组织中表达下调和上调的miRNA 构建相应的miRNAmRNA 调控网络，并从中筛选出2 个最重要的子网络。

在下调miRNA、上调mRNA 网络的子网络中靶基因GATA6在子宫内膜异位症相关研究中被报道最多，属于GATA 转录因子家族的一员，在调节细胞生长分化及其功能的维持等方面有重要作用。IZAWA等[21]研究发现GATA6 在异位内膜组织呈高表达状态，同时去甲基化可增强GATA6 的表达水平；DYSON等[22]在异位病灶基质细胞中也发现了GATA6 的去甲基化状态，并和激素抵抗有关。本研究中hsamiR-20a-5p、hsa-miR-141-3p、hsa-miR-200b-3p、hsa-miR-449b-3p 均可以靶向作用于GATA6，其中hsa-miR-141-3p、hsa-miR-200b-3p 均属于miR-200家族,可在肿瘤侵袭转移过程中起到抑制上皮间质转化作用。DONG 等[23]通过荧光素酶基因报告验证GATA6是miR-141 的靶基因，并且在人脐静脉内皮细胞中验证miR-141 可通过靶向作用于GATA6等多个基因从而抑制血管生成。FAN 等[24]在人主动脉内皮细胞中发现miR-200a-3p 的上调，会抑制TNF-α 引起的GATA6 的表达。本研究均说明miR-200 家族和GATA6 之间存在靶向调控关系，但是该调控关系是否也在子宫内膜异位症发病机制中起重要作用值得进一步研究。

在上调miRNA、下调mRNA 网络的子网络中，值得注意的是PGR、ESR1 均受到hsa-miR-125-5p的调控而处于下调状态。子宫内膜异位症是一种激素依赖性疾病，临床常使用孕激素来拮抗雌激素的作用，抑制异位内膜生长或使之萎缩，达到缓解症状的目的，但是会出现孕激素抵抗而影响其疗效。现有研究推测这和子宫内膜异位症患者体内孕激素受体表达失调，以及孕激素受体亚型分布异常有关，本研究中异位内膜组织中孕激素受体下调与现有研究的结果一致，并且其下调有可能受到miRNA 的调控。另外，雌激素是强有丝分裂原，既可直接刺激异位病灶的细胞增殖，又可以促进新生的血管内皮细胞增生，同时还可以刺激多种炎症相关因子参与子宫内膜异位症的发病过程[25]。雌激素受体有2 种亚型： ESR1（ERα）和ESR2（ERβ），已有研究表明在异位病灶及异位内膜基质细胞中，ESR1 表达下调，而ESR2 表达呈优势状态[26]，同时TRUKHACHEVA 等[27]还证明在异位内膜基质细胞中，高表达的ESR2 可以抑制ESR1 的表达。本研究中ESR1 表达也出现下调，与前人研究结果相符。2016年ZHOU 等[28]曾发现，miR-196a 通过MEK/ERK 信号通路改变孕激素受体亚型在子宫内膜异位症中的表达，因此本研究中的miR-125-5p 是否也可以通过调节雌、孕激素受体而导致子宫内膜异位症的发生、发展值得进一步研究。

综上所述，本研究通过构建miRNA-mRNA 调控网络，研究与子宫内膜异位症中的关键miRNA 和差异表达基因的关系，为子宫内膜异位症发病机制及诊断标志物的研究提供了思路。