癌症已成为危害人类健康的最大因素之一^[1]，口腔癌是最常见的癌症之一，吸烟和饮酒被认为是其最主要的致病因素^[2]，1990—2017年世界口腔癌发病率和病死率增长了一倍^[3]. 目前, 口腔癌主要的治疗手段是通过手术治疗，而且早期的口腔癌不容易诊断，降低了患者的生存率，且预后较差，晚期口腔癌患者的5年生存率仍然低于63%^[4]. 近年来，大量研究表明，基因的异常表达是导致口腔癌发生与发展的重要因素，因此，从分子水平对口腔癌组织进行研究，探索口腔癌的发病机制和可能的治疗靶点具有重要意义.

微小核糖核酸(MicroRNAs, miRNAs)是由内源性发夹结构转录产物衍生而来的一种具有22个左右核苷的单链RNA^[5]，是基因表达的重要调节因子. miRNAs通过靶向mRNAs的翻译抑制或切割来抑制基因表达^[6]. 据估计，30%的人类基因可能受到miRNAs的调控^[7]. 越来越多的miRNAs与人类疾病的研究表明，miRNAs与人类疾病尤其是癌症的发生密切相关^[8].

目前，生物信息学及基因芯片技术已广泛应用于基因变化所致肿瘤的研究中. 借助现有的miRNAs基因芯片数据，挖掘与口腔癌相关的miRNAs可以探索口腔癌的发生机制，并对口腔癌的早期诊断及治疗提供数据及基础. 本研究利用生物信息学方法对GEO数据库中的口腔癌miRNAs表达数据进行分析，筛选差异表达的miRNAs，并通过基因富集分析分析差异表达的miRNAs的功能，为口腔癌的靶向治疗提供一定的理论依据.

1.   材料与方法

1.1   数据来源

口腔癌miRNAs芯片表达数据集(GSE124566和GSE113956)从GEO数据库^[9](https://www.ncbi.nlm.nih.gov/geo)下载，GSE124566含有10个口腔癌组织样本和10个正常组织样本，GSE113956含有25个口腔癌组织样本和15个正常组织样本. 平台分别为GPL18402和GPL18058，芯片研究类型为Non-coding RNA profiling by array，种属为homo sapiens，都是标准的口腔组织，具有代表性与典型性.

1.2   数据处理与筛选

利用GEO2R在线分析工具(http://www.ncbi.nlm.nih.gov/geo/geo2r/)分析口腔癌与正常组织之间的差异表达miRNAs. 针对差异表达miRNAs采用阈值标准来计算基因表达的Fold Change和校正P值, P≤0.05, |log₂ FC|≥1作为有统计学意义为标准，每个数据集的所有差异用miRNAs火山图表示.

为了使结果更加规范、准确，需要使用Venn Diagram在线工具(http://bioinformatics.psb.uget.be/webtools/Venn)进行筛选上述2个数据集共同的差异表达miRNAs.

1.3   差异miRNAs功能富集分析

基因本体论分析^[10](Gene Ontology，GO)是对基因和蛋白功能进行分类和描述的分析，包括细胞学组分(Cellular Components, CC)、生物学过程(Biological Process, BP)、分子生物学功能(Molecular Function, MF). 京都基因和基因组百科全书^[11](KEGG)用于注释参与其中的基因列表和信号通路网络.

FunRich(http://funrich.org/)是一个独立的软件工具，主要用于基因和蛋白质的功能和相互作用网络分析^[12]，可用于绘制韦恩图，其最大的优点是可以进行基因ID转换、基因富集和miRNAs的富集分析并可直接进行可视化等. 本研究利用FunRich3.1.3对GO和KEGG通路进行富集分析来确定差异表达miRNAs的功能.

1.4   差异表达miRNAs靶基因预测

为了探索上述筛选出来的差异表达miRNAs作用的靶标基因，方便后期对靶标基因进行研究. 本研究使用miRWalk3.0(http://mirwalk.umm.uni-heidelberg.de/)对差异表达的miRNAs进行靶标基因的预测.

1.5   蛋白相互作用分析

利用FunRich寻找靶基因的蛋白相互作用关系，筛选连接度前10位的基因作为网络核心基因.

2.   结果

2.1   差异表达miRNAs筛选

本研究从GEO(Gene Expression Omnibus)数据库下载了2个不同的口腔癌miRNAs芯片表达数据集(GSE124566和GSE113956)，GSE124566含有10个口腔癌组织样本和10个正常组织样本，GSE113956含有25个口腔癌组织样本和15个正常组织样本. 用GEO2R在线分析工具，调整后P≤0.05, |log₂ FC|≥1作为截断值的标准，以分析口腔癌和正常组织之间的差异表达miRNAs.

在GSE124566和GSE1139566中分别筛选了总共109、1 079个差异表达miRNAs，分别包括41、673个上调基因和68、406个下调基因，图 1为这2个数据集的火山图，横坐标表示log₂ FC(癌症样本比正常样本)，纵坐标表示-log₁₀ P. 其中黄色的点表示不显著的差异表达miRNAs，红色的点表示上调的差异表达miRNAs，蓝色的点表示下调的差异表达miRNAs.

图  1  口腔癌miRNAs数据集表达水平火山图

Figure  1.  The volcano map of miRNA expression level in the oral cancer data set

下载: 全尺寸图片 幻灯片

2.2   共同差异表达miRNAs

为了得到更有可信度的差异表达miRNAs，分别对2个数据集上调和下调的miRNAs分别绘制韦恩图取交集(图 2)，结果得出16个上调的共表达差异miRNAs(表 1)和14个下调的共表达差异miRNAs(表 2).

图  2  共同差异表达miRNAs Venn图

Figure  2.  The common differential expression miRNAs Venn map

下载: 全尺寸图片 幻灯片

表  1  16个差异表达上调的miRNAs的P及log₂ FC

Table  1.  The P value and log₂ FC of 16 miRNAs with up-regulated differential expression

miRNA	P	log2 FC
hsa-miR-513b	6.330×10-9	2.750 163
hsa-miR-744-5p	2.233×10-8	2.284 359
hsa-miR-375	7.553×10-16	4.135 624
hsa-miR-1471	1.173×10-3	1.197 684
hsa-miR-598	1.253×10-4	1.090 587
hsa-miR-4324	1.623×10-5	1.356 905
hsa-miR-378a-5p	3.103×10-8	3.065 295
hsa-miR-29c-5p	8.863×10-12	3.915 925
hsa-miR-1224-5p	2.123×10-4	1.316 750
hsa-miR-3911	1.783×10-21	2.475 252
hsa-miR-4647	4.793×10-10	2.891 719
hsa-miR-3188	1.823×10-18	4.603 102
hsa-miR-24-1-5p	4.683×10-3	1.151 63
hsa-miR-204-5p	1.333×10-3	1.574 574
hsa-miR-338-3p	1.023×10-6	2.736 121
hsa-miR-513c-5p	2.153×10-17	2.510 716

下载: 导出CSV 

| 显示表格

表  2  14个差异表达下调的miRNAs的P及log₂ FC

Table  2.  The P value and log₂ FC of 14 miRNAs with down regulated differential expression

miRNAs	P	log2 FC
hsa-miR-4778-5p	2.593×10-14	-2.410 98
hsa-miR-31-5p	1.933×10-2	-1.342 41
hsa-miR-299-3p	1.743×10-6	-2.222 96
hsa-miR-4419a	1.453×10-11	-3.270 17
hsa-miR-223-3p	1.583×10-6	-1.155 63
hsa-miR-142-5p	1.293×10-4	-2.564 01
hsa-miR-23a-5p	1.753×10-12	-1.997 68
hsa-miR-374c-5p	7.603×10-4	-1.400 69
hsa-miR-454-3p	2.263×10-4	-1.954 80
hsa-miR-625-5p	6.983×10-20	-3.165 28
hsa-miR-765	1.823×10-3	-1.113 14
hsa-miR-21-5p	1.523×10-4	-1.207 61
hsa-miR-142-3p	9.313×10-4	-2.640 01
hsa-let-7i-5p	6.863×10-8	-1.426 81

下载: 导出CSV 

| 显示表格

LU等^[13]将6种口腔癌细胞系与正常角质形成细胞作比较，发现在口腔癌细胞系中miR-31表达下调，这与本研究结果一致.

2.3   共同表达差异基因功能分析

为了进一步明确30个共同差异表达miRNAs的生物学功能，利用软件FunRich对差异表达miRNAs进行了GO与KEGG生物途径富集分析.

2.3.1   GO分析

GO分析结果显示：16个上调的miRNAs的细胞组分主要在细胞核(Nucleus)、细胞质(Cytoplasm)、细胞质膜(Plasma membrane)、溶酶体(Lysosome)和核仁(Nucleolus)(图 3A)；分子功能主要集中在转录因子活性(Transcription factor activity)、转录调节活性(Transcription regulator activity)、RNA结合(RNA binding)、丝氨酸/苏氨酸激酶活性(Protein serine/thereonine kinase activity) 和受体信号蛋白丝氨酸/苏氨酸蛋白激酶活性(Receptor signacling protein serine/threonine kinase activity)(图 3B)；生物学过程主要富集在在转运(Transport)，信号转导(Signal transduction)，碱基、核苷、核苷酸和核酸代谢的调节(Regulation of nucleobase, nucleoside, nucleotide and nuclelc acid metabolism)，基因表达和表观遗传的调控(Regulation of gene expression，epigenetics)，囊泡介导的转运(Vesicle-mediated transport)(图 3C).

图  3  共同差异表达上调miRNAs GO分析

Figure  3.  The GO analysis of up-regulated miRNAs in common differential expression

下载: 全尺寸图片 幻灯片

14个下调的miRNAs的细胞组分主要富集细胞核(Nucleus)、细胞质(Cytoplasm)、细胞质膜(Plasma membrane)和溶酶体(Lysosome)和高尔基体(Golgl aparatus)(图 4A)；分子功能主要集中转录因子活性(Transcription factor activity)、转录调节活性(Transcription regulator activity)、RNA结合(RNA binding)、泛素蛋白特异性蛋白酶活性(Ubiquitin-specific protease activity)和转运蛋白活性(Transporter activity)(图 4B)；生物学过程主要富集在转运(Transport)，信号转导(Signal transduction)，碱基、核苷、核苷酸和核酸代谢的调节(Regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism)，基因表达和表观遗传的调控(Regulation of gene expression, epigenetics)和细胞间通讯(Cell communication)(图 4C).

图  4  共同差异表达下调miRNAs的GO分析

Figure  4.  The GO analysis of down-regulated miRNAs in common differential expression

下载: 全尺寸图片 幻灯片

2.3.2   KEGG分析

KEGG/生物途径分析显示：上调的miRNAs主要参与蛋白多糖多配体多糖介导的信号(Proteoglycan syndecan-mediated signaling events)、ErbB受体信号网络(ErbB receptor signaling network) 和磷脂酰肌醇蛋白聚糖通路(Glypican pathway)、质膜雌激素受体信号转导(Plasma membrane estrogen receptor signaling)和LKB1信号(LKB1 signaling events)(图 5A). 另外，下调的miRNAs主要参与蛋白多糖多配体多糖介导的信号(Proteoglycan syndecan-mediated signaling events)、ErbB受体信号网络(ErbB receptorsignaling network)和磷脂酰肌醇蛋白聚糖相关通路(Glypican pathway)、β整合素细胞表面相互作用(β integrin cell surface interactions)，和整合素家族细胞表面相互作用(Integrin family cell surface interactions)(图 5B).有研究表明miR-23、miR-142-3p以及miR-223参与了口腔癌变的磷脂酰肌醇激酶参与的相关通路^[14]，与本研究结果一致.

图  5  共同差异表达miRNAs的KEGG分析

Figure  5.  The KEGG analysis of common differentially expressed miRNAs

下载: 全尺寸图片 幻灯片

2.4   差异表达miRNAs靶基因预测

经过前面的筛选，获得了30个共同差异表达的miRNAs，打开靶标基因预测工具miRWalk3.0，点击Target Ming里的miRNAs，导入30个差异表达的miRNAs进行靶基因预测，最后得到SMG7、NFIA、RIMS1、POTEA 等13 796个靶基因.

2.5   基因蛋白相互作用分析

利用FunRich对靶基因进行相互作用分析，获得了蛋白相互作用网络，经筛选并在蛋白相互作用网络中寻找连接数前10位的作为核心基因. 经过筛选发现，编码生长因子受体结合蛋白2基因(GRB2)、编码表皮生长因子受体基因(EGFR)、编码酪氨酸3单加氧酶基因(YWHAG)、编码调控细胞分裂增殖的Tp53蛋白基因(TP53)、酪氨酸激酶基因(ABL1) 、肿瘤坏死因子受体相关因子(TRAF6)、编码核蛋白的癌基因(MYC)、转化生长因子TGF-β信号转导相关基因(SMAD2)、B细胞k轻肽基因增强子抑制因子(IKBKE)和编码非受体酪氨酸激酶基因(SRC)是靶基因相关蛋白相互作用网络中的核心基因.

3.   讨论

目前，对口腔癌的治疗来说，大多数还是放疗、化疗和手术切除病变位置，这些方法的治疗效果不理想，预后较差，导致口腔癌患者的生存率不高，对于口腔癌在分子层面的发病机制和预后机制的研究也相对较少. 所以，研究发现有关口腔癌在分子层面的发病机制有利于对口腔癌有更清楚的认识，以及更有利于研发癌症靶向治疗的药物.

根据上文可以发现：(1)在分子功能方面，差异表达的miRNAs最显著集中于转录调节活性和转录因子活性. 可见，差异表达的miRNAs对某种口腔癌基因的表达有显著或轻微的抑制作用，是普遍转录因子还是特异性转录因子发挥作用尚不明确. (2)在细胞组分方面，差异表达下调的miRNAs与差异表达上调的miRNAs不同的是下调的miRNAs富集于高尔基体，其作用主要是将蛋白质加工和运输后送到细胞特定的部位或者细胞外，差异表达下调的miRNAs这一独特的作用也在分子功能富集上得到了验证，差异表达下调的miRNAs区别于上调的miRNAs不同的是下调的miRNAs主要富集在蛋白转运. 根据这一特点可知，下调的miRNAs对某些蛋白的转运起作用，后期如果能找到生物标志物，那么下调miRNAs很有可能成为新的治疗靶点. (3)在生物学过程方面，上调和下调miRNAs主要差别在于转运和细胞间通讯过程，可推测它们分别参与物质传递和信息交流相关过程. (4)在差异表达miRNAs参与的生物通路方面，不同之处在于上调的miRNAs主要参与质膜雌激素受体信号转导和LKB1信号通路，下调的miRNAs主要参与整合素家族细胞表面相互作用相关通路，这些都可能对研究口腔癌患者发生恶性肿瘤的机制和预后机制提供一些线索.

本研究也存在一定的局限性，差异表达的miRNAs以及涉及到的生物信号通路，尚未进行分子水平的验证，这将是下一步的研究目标.

4.   结论

本研究对口腔癌与正常组织间的差异表达miRNAs进行分析，并对差异miRNAs进行GO功能注释、KEGG信号通路分析，使我们能够识别出一些miRNA，通过分析预测并通过PPI互作分析筛选出了联系最紧密的靶基因. 这些数据会为口腔癌检测提供一些有用的信息和方向.