成釉细胞瘤患者lncRNA-mRNA共表达网络分析及作用研究

游 婧，傅 瑜，单 杰，张 杰，林晓虎，江裕程，卞成玥，罗颐辰

成釉细胞瘤(ameloblastoma)是多发于颌骨内的较为常见的牙源性上皮性肿瘤，主要临床表现为颌骨进行性、无痛性膨大，由于其具有局部侵袭性，能够导致患者颌面部畸形[1-2]。目前，对于成釉细胞瘤的主要治疗方式为外科手术，但术后复发率高并存在恶变风险，影响患者的生存质量[3-4]。因此，深入研究成釉细胞瘤的发生发展机制，探索有效的治疗手段具有重要意义。

lncRNA是一类长度大于200个核苷酸的非编码RNA，不具有或很少具有蛋白编码功能，近年来研究显示lncRNA在细胞周期、细胞分化、表观遗传等生物学过程中发挥重要调控作用[5-7]。lncRNA及其调控基因的异常表达可能与多种疾病的发生发展有关[8]。Davanian等[9]的研究首次报道成釉细胞瘤与正常组织中lncRNA的表达差异以及lncRNA的表达与成釉细胞瘤瘤体大小的关系，但关于lncRNA及相关mRNA在成釉细胞瘤发生发展中的作用尚未深入分析。本研究通过对公共数据库基因表达数据库(gene expression omnibus, GEO)中的两组成釉细胞瘤患者组织的芯片数据集GES38494及GES132472进行合并分析，筛选差异表达的基因及lncRNA并构建lncRNA-mRNA共表达网络，对异常表达的lncRNA及mRNA进行生物信息学分析，寻找与成釉细胞瘤发生发展相关的靶基因，为成釉细胞瘤的分子机制研究及临床治疗提供理论基础。

1 资料与方法

1.1 研究对象

本研究所选成釉细胞瘤患者组织芯片表达谱数据集GES38494及GES132472下载自美国国立生物信息中心(National Center for Biotechnology Information, NCBI)公共数据库基因表达数据库(Gene Expression Omnibus, GEO)，芯片分析均取材于手术中的相关组织，芯片数据GES38494由Heikinheimo等研究者上传[10]，采用Affymetrix Human Genome U133 Plus 2.0 Array芯片平台，芯片数据包括15个成釉细胞瘤组织样本及4个正常组织样本；GES132472由Kondo等研究者上传[11]，采用Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381芯片平台，芯片数据包括8个成釉细胞瘤组织样本及8个正常组织样本。

1.2 分析方法

1.2.1 数据预处理 使用R-studio软件对芯片数据进行分析。应用R语言“affyPLM”及“affy”程序包进行背景矫正，RMA算法对数据进行标准化预处理，合并两数据集并应用“sva” 程序包的ComBat算法对合并的数据集进行处理以矫正因实验仪器、人员、试剂等因素造成的批次效应。

1.2.2 筛选差异表达基因及lncRNA 应用R语言“limma”程序包对经过预处理的数据集进行分析，以改变倍数在2倍以上(|Fold Change| ≥2)且矫正后P值(adj.P. Value)< 0.05为筛选条件，得到差异表达基因。应用R语言“ggplot2”程序包绘制火山图，选取排名前100的差异表达基因应用“pheatmap”程序包绘制聚类分析图。在GENCODE数据库(https://www.gencodegenes.org)下载lncRNA注释包，应用R语言“AnnotationHub”程序包进行ID转换，获得差异表达基因中的lncRNA，即差异表达lncRNA。

1.2.3 lncRNA-mRNA共表达网络构建 对差异表达lncRNA与mRNA的表达量进行皮尔森相关分析，选取皮尔森相关系数(Pearson correlation coefficient，PCC)的绝对值(|PCC|)>0.85且P<0.05的lncRNA-mRNA对，采用Cytosacape软件绘制lncRNA-mRNA共表达网络图。

1.2.4 GO及KEGG通路富集分析 应用R语言“clusterProfiler”程序包对lncRNA-mRNA共表达网络中的mRNA进行GO及KEGG通路富集分析，以P<0.05为标准进行筛选。

1.2.5 GESA信号通路富集分析 以应用R语言“clusterProfiler”程序包对差异表达基因进行GSEA分析，以标准化P<0.05、FDR<0.25为标准筛选阳性基因集并使用“ridgeplot”及“gseaplot2”程序包进行可视化分析。

1.2.6 PPI网络分析 为评估lncRNA-mRNA共表达网络中的mRNA之间的交互关系，采用STRING(https://string-db.org/cgi)数据库以confidence≥0.4 为阈值条件构建蛋白质相互作用(Protein-protein interaction, PPI)网络，利用Cytosacape软件进行可视化分析，应用cytoHubba插件获得核心基因。

2 结 果

2.1 差异表达基因的筛选结果

通过对合并数据信息进行处理，共有17 631个基因被纳入后续分析。以|Fold Change|≥2且adj.P. Value<0.05为标准，共筛选出801个差异表达的基因，其中416个在成釉细胞瘤组织中表达上调，385个表达下调(图1A)。取排名前100的差异基因绘制热图(图1B)。

图1 差异表达基因的火山图(A)及聚类分析图(B)

2.2 lncRNA-mRNA共表达网络分析

通过对差异表达基因进行lncRNA重注释发现其中包含6个有意义的lncRNA，分别为SOX2-AS1、HCG22、LOC100506098、MAGI2-AS3、MEG3、MIR100HG，筛选与差异表达lncRNA有关联的mRNA，发现283个mRNA与6个lncRNA具有高度共表达关系。构建lncRNA-mRNA共表达网络进行分析，方框节点代表lncRNA，圆形节点代表mRNA，6个lncRNA与283个mRNA共组成了305个相关链接，其中160个呈正相关，145个呈负相关。MAGI2-AS3居于共表达网络的核心位置，与236个mRNA具有共表达关系，且关联度高的mRNA多呈正相关。23个mRNA分别与2个lncRNA相关联(图2)。

方框代表lncRNA，圆形代表mRNA，直线两端连接的mRNA与lncRNA存在共表达关系。实线表示正相关，虚线表示负相关。红色代表在成釉细胞瘤中高表达，绿色代表低表达。方框的大小代表lncRNA在成釉细胞瘤中的表达程度，圆形的大小代表共表达关联度

2.3 GO及KEGG通路富集分析

GO富集分析结果显示lncRNA-mRNA共表达网络中的mRNA参与的生物过程主要富集于上皮发育(GO:0008544)、皮肤发育(GO:0043588)、表皮细胞分化(GO:0009913)、角化细胞分化(GO:0030216)、角化作用(GO:0031424)等，细胞组成主要定位于含胶原的细胞外基质(GO:0062023)、角化细胞套膜(GO:0001533)、膜的锚定成分(GO:0031225)、顶端结合复合体(GO:0043296)、细胞桥粒(GO:0030057)等，显著富集的分子功能包括肽酶调节因子活性(GO:0061134)、肽链内切酶抑制因子活性(GO:0004866)、肽酶抑制因子活性(GO:0030414)、肽链内切酶调节因子活性(GO:0061135)等(图3A)。

KEGG通路富集分析发现lncRNA-mRNA共表达网络中的mRNA主要参与干细胞多能性调节通路(hsa04550)、胰腺分泌(hsa04972)、视黄醇代谢(hsa00830)、急性骨髓性白血病(hsa05221)、糖酵解和糖异生(hsa00010)等通路(图3B)。

图3 Go分析(A)和KEGG通路富集分析(B)

2.4 GSEA富集分析

采用GSEA对差异表达基因进行KEGG通路富集分析，发现共有13个通路与成釉细胞瘤相关，选取富集分数排名前十位的通路分析(图4A)，其中代谢通路(Metabolic pathways)与成釉细胞瘤成正相关(P=0.001)，其他通路则为负相关。按照富集分数进行排序，轴突导向(Axon guidance)、细胞外基质-受体相互作用(ECM-receptor interaction)以及黏着斑(Focal adhesion)相关通路(P=0.003)为排名前三位的通路(图4B)。

A：KEGG通路富集；B：上图：正相关通路；下图：以富集分数排名前三位的负相关通路

2.5 蛋白质-蛋白质相互作用(protein-protein interaction, PPI)分析

基于String数据库中的信息进行PPI分析，包含280个节点，280条连接线，平均每个节点的度为2。利用Cytoscape进行可视化分析(图4)，对关键基因进行筛选排序，筛选出节点连接度(degree)排名前10位的核心基因，其中连接最广泛的是外膜蛋白(involucrin，IVL)(表1)。

图5 PPI分析图

表1 连接度Top 10核心基因

3 讨 论

lncRNA是在真核生物中发现的一类长度大于200个核苷酸的RNA，多不具备编码蛋白质功能[12]，广泛参与染色质的修饰、转录、表观遗传调控以及免疫监控等生物学过程，其转录与功能的失调可能是许多疾病发生发展的重要因素[13]。lncRNA-mRNA共表达网络的建立是研究lncRNA功能和调控机制的重要方法[14]。近年来越来越多的研究者将对成釉细胞瘤的研究聚焦于lncRNA对其生物学特征的调控[15-16]，但是关于lncRNA及相关mRNA在成釉细胞瘤发生发展中的作用尚未见深入分析。在本研究中，我们通过生物信息学方法分析研究了lncRNA-mRNA共表达网络的功能，为后续实验研究提供基础。

研究表明lncRNA的调控机制主要包括以下四个方面：通过碱基互补配对的方式调控mRNA；吸附miRNA间接抑制miRNA对靶基因的调控作用；与转录因子结合形成复合体调控基因转录活性；改变DNA/RNA甲基化状态、染色体结构和修饰状态控制相关基因表达[12-13，17]。本研究通过对差异表达基因进行注释筛选出6个差异表达的lncRNA，其中SOX2-AS1在成釉细胞瘤中表达上调最明显。根据皮尔斯相关分析(|PCC|> 0.85)在差异基因表达谱中筛选出283个与差异表达lncRNA高度相关的mRNA，构建lncRNA-mRNA共表达网络。在此网络中我们发现一个lncRNA可能调控多个mRNA，而一个mRNA也可能与多个lncRNA相关，lncRNA及与其关联度高的mRNA在成釉细胞瘤中的表达趋势大多一致，推测lncRNA可能通过调控与其关联的mRNA参与成釉细胞瘤的发生发展。然而，共表达网络中lncRNA与mRNA的具体作用关系及相关机制仍有待进一步研究。

目前研究认为，成釉细胞瘤是一种上皮来源的牙源性肿瘤[18-19]。本研究中GO分析发现lncRNA-mRNA共表达网络中的差异mRNA主要影响上皮发育、皮肤发育、表皮细胞分化等GO生物学功能，上皮发育等生物学功能与成釉细胞瘤的发生发展密切相关。KEGG通路富集分析显示共表达网络中的差异mRNA在干细胞多能性调节通路、胰腺分泌、视黄醇代谢等通路富集程度较高，而这些通路在成釉细胞瘤的发生机制中研究较少，为发病机制研究提供了新的方向。基因集富集分析(gene set enrichment analysis，GSEA)是一种基于基因集的富集分析方法，其相对于GO/KEGG富集分析的优势在于能够避免遗漏部分差异表达不显著但却具有重要生物学意义的基因[20-21]。本研究中，我们通过GSEA分析共表达网络中的差异mRNA的KEGG通路富集情况，发现了13条与KEGG通路富集分析结果不同的通路信息，其中12条通路与成釉细胞瘤呈负相关，只有代谢通路(hsa01100)的基因集在成釉细胞瘤中表达上调。GSEA结果为成釉细胞瘤中异常表达基因富集通路分析进行补充以及完善。

另外，本研究对lncRNA-mRNA共表达网络中差异表达的基因进行蛋白质互作网络分析，并筛选出互作网络中的核心基因，发现核心基因外膜蛋白(involucrin，IVL)连接度最高。IVL是角质形成细胞终分化的早期指标，在表皮角化包膜形成过程中发挥支架作用。有研究显示，IVL与角化过度的手湿疹[22]、角化病[23]、接触性皮炎[24]有关。该结果进一步为成釉细胞瘤的治疗提供了研究方向。

综上所述，本研究利用生物信息学分析工具，筛选出成釉细胞瘤患者异常表达的基因及lncRNA，构建lncRNA-mRNA共表达网络，对网络中的mRNA参与的生物过程及通路进行分析，并根据网络中的差异表达基因谱分析相关的蛋白相互作用关系，列出连接度最高的10个核心基因，为成釉细胞瘤发生发展的机制研究及临床治疗靶基因的选择提供理论基础。