喉鳞状细胞癌的转录组数据分析及生物标志物筛选

王娜娜，杨川

1延安大学附属医院耳鼻咽喉科，陕西延安716000

2河北医科大学第二医院耳鼻咽喉科，石家庄050000

3平昌县人民医院耳鼻咽喉科，四川巴中636040

喉癌是耳鼻咽喉-头颈外科第二高发肿瘤，其中喉鳞状细胞癌占全部病理分型的95%。喉鳞状细胞癌的病因复杂，近年来其发病率不断增加，患者的晚期生存率仍较低。随着机器人辅助手术在头颈部肿瘤治疗中的开展及食管、气管一期重建技术的不断成熟，中国临床肿瘤学会（Chinese Society of Clinical Oncology，CSCO）头颈部肿瘤诊疗指南也强调手术治疗对提高生存率的重要性，但仍面临术后功能障碍、生活质量下降等严峻挑战。另一方面，头颈部肿瘤手术仍被认为是外科“皇冠”，手术水平直接关系到患者的生存预后。因此，进一步阐明喉鳞状细胞癌的分子机制，寻找潜在的生物学靶点，对喉鳞状细胞癌进行早期诊断和治疗仍十分必要。信使RNA（messenger RNA，mRNA）、微小RNA（microRNA，miRNA）、长链非编码RNA（long noncoding RNA，lncRNA）在喉鳞状细胞癌的发生、发展中发挥重要作用。本研究通过肿瘤基因图谱（The Cancer Genome Atlas，TCGA）转录组数据挖掘，分别阐明lncRNA、miRNA和mRNA在喉鳞状细胞癌中的作用机制及其潜在联系，现报道如下。

1 资料与方法

1.1 数据获取及差异分析

在TCGA中下载喉鳞状细胞癌转录组数据和相关临床信息，其中肿瘤样本111例，正常样本12例。采用“edgR”包进行差异分析，设|logFC|﹥2，错误发现率（false discovery rate，FDR）﹤0.01，分别得到差异lncRNA、差异mRNA和差异miRNA，采用R软件绘制可视化图形并进行统计分析。

1.2 蛋白相互作用网络分析

采用STRING11.0（https://string-db.org/cgi/input.pl）在线工具将喉鳞状细胞癌Top500差异mRNA构建蛋白相互作用（protein-protein interaction，PPI）网络，相互作用值=0.7；Cytoscape插件CentiScape2.2、cytoHubba0.1用于寻找关键基因，并进行可视化展示。

1.3 功能富集和通路富集分析

喉鳞状细胞癌差异mRNA用Clusterprofiler包进行细胞组分（cellular component，CC）、生物学过程（biological process，BP）、分子功能（molecular function，MF）的基因文本（gene ontology，GO）富集及京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）信号通路富集。

1.4 MYHAS的共表达mRNA（cor-mRNA）

采用Pearson方法计算cor-mRNA，并取相关系数cor≥0.4。

1.5 构建竞争性内源RNA（competitive endogenous RNA，ceRNA）网络和转录因子网络

ceRNA网络：采用本地版miRanda和Targetscan分别预测lncRNA结合的miRNA，并与差异miRNA取交集。再通过miRTarBase、Targetscan、miRBD分别预测miRNA结合的mRNA，并与差异mRNA取交集，建立lncRNA-miRNA-mRNA网络联系，并通过Cytoscape3.7.2软件进行可视化展示。

转录因子网络：将钙信号通路相关的mRNA于DAVID6.8在线网站（https://david.ncifcrf.gov/）中进行转录因子富集，建立lncRNA相关mRNA的转录因子联系，并通过Cytoscape3.7.2软件进行可视化展示。

2 结果

2.1 一般资料

肿瘤患者的年龄、肿瘤T分期、病理分级、饮酒、吸烟情况对肿瘤患者总体生存率的影响不明显。远处转移对患者总生存率的影响显著，临床转移患者的5年生存率为0%，明显低于无临床转移患者的50.33%，差异有统计学意义（χ=12.230，P=0.0005）；病理转移患者的5年生存率为0%，明显低于无病理转移患者的44.83%，差异有统计学意义（χ=18.260，P=0.0001）。女性患者的5年生存率为30.00%，低于男性患者的50.11%，差异有统计学意义（χ=4.454，P=0.0348）。发生病理淋巴结转移患者的5年生存率为48.08%，低于未发生病理淋巴结转移患者的70.14%，差异有统计学意义（χ=4.539，P=0.0331），但临床上发生淋巴结转移与未发生淋巴结转移患者的5年生存率比较，差异无统计学意义（χ=1.150，P=0.2835）。

2.2 lncRNA、miRNA、mRNA差异表达情况

lncRNA、miRNA和mRNA各自的主成分分析结果显示，肿瘤组与对照组间表达谱数据相互独立，组内重复性较好，数据可信度高（图1A、1C和1E）。当|logFC|﹥2，FDR﹤0.01时，得到的差异lncRNA共687个，其中上调561个，下调126个（图1B）；得到的差异miRNA共54个，其中上调28个，下调26个（图1D）；得到的差异mRNA共1818个，其中上调950个，下调868个（图1F）。

图1 lncRNA、miRNA、mRNA的主成分分析、火山图

2.3 差异lncRNA和miRNA表达的喉鳞状细胞癌预后的单因素分析

根据RNA-Seq临床数据计算了差异lncRNA表达的喉鳞状细胞癌预后的单因素分析结果，取

P

﹤0.05，共有64个lncRNA满足条件，根据是否能在NCBI中查到相关lncRNA信息，总生存率是否有统计学意义，筛选了11个lncRNA，其中

LINC00278

、

LINC00689

、

MYHAS

属于低风险基因（HR﹤1），

MNX1-AS1

、

LINC02575

、

HOXB-AS4

、

LSAMP-AS1

、

LINC02086

、

IGFL2-AS1

、

LINC02253

、

CASC20

属于高风险基因（HR﹥1）。同理，计算出了差异miRNA表达的喉鳞状细胞癌预后的单因素分析结果，

P

﹤0.05的miRNA共7个，其中miRNA-383、miRNA-196a-2、miRNA-196a-1、miRNA-100、miRNA-4652是高风险基因（HR﹥1），miRNA-99a、miRNA-301a是低风险基因（HR﹤1）。（表1）

表1 差异lncRNA和miRNA表达的喉鳞状细胞癌预后的单因素分析

2.4 差异mRNA功能分析和关键基因筛选

将1818个差异mRNA用Clusterprofiler进行GO分析和KEGG信号通路分析，取

P

值或

P.adjust

值﹤0.05。分子功能（MF）以通道活性和激活跨膜转运活性为主，细胞组分（CC）以细胞外基质（extracelluar matrix，ECM）包含的胶原纤维和转运复合体为主，生物学过程（BP）以细胞外基质形成和调节离子跨膜转运为主（图2A～2C）。KEGG信号通路中以细胞因子-细胞因子受体间相互作用和钙信号通路为主要作用通路，其中也包括环磷酸腺苷（cyclic adenosine monophosphate，cAMP）信号通路、环磷酸鸟苷（cyclic guanosine monophosphate，cGMP）-cGMP依赖性蛋白激酶（cGMP-dependent protein kinase，PKG）信号通路和ECM受体相互作用等肿瘤发生、发展的经典信号通路（图2D）。将差异倍数前500的mRNA进行String在线分析，构建PPI网络，并通过Centiscape插件进行关键基因筛选，根据Degree进行筛选，绿色﹤10，青色≥10，紫色≥15，红色≥20，其中

ACTN2

、

ACTN3

、

MYL1

、

MYL2

、

MYL3

、

MYH6

、

TTN

网络连接最多，属于关键基因（图2E）。

图2 差异mRNA的GO、KEGG信号通路富集和PPI网络

2.5 MYHAS的cor-mRNA、关键基因筛选及功能富集

MYHAS的共表达mRNA（cor≥0.4）共253个，采用插件cytoHubba分别计算Closeness、Degree、EPC、MCC、MNC和Radiality 6种算法下前10的基因并取交集，得到关键基因：

TNN

、

MYL1

、

MYL3

、

ACTN2

、

MYH6

。GO富集显示cor-mRNA主要参与调节金属离子转运、跨膜转运复合体和肌动蛋白结合等（图3A～3C）。信号通路分析显示，cormRNA主要富集在钙信号通路、cGMP-PKG信号通路和催产素信号通路等（图3D）。其中钙离子信号通路中富集的基因分别是PRKACA、CAMK2A、SLC8A3、CAMK2B、RYR1、PLN、CASQ2、PHKG1、HRC、MYLK3、TNNC1、CACNA1S、MYLK2、CASQ1、ATP2A1、TNNC2。

图3 MYHAS 的cor-mRNA的GO和KEGG信号通路富集

2.6 ceRNA网络和转录因子网络

将单因素分析与生存分析中P﹤0.05的lncRNA取交集，得到21个lncRNA，并分别进行mi-Randa和Targetscan靶基因预测并取交集，得到2363个与lncRNA可以结合的miRNA，然后将预测的miRNA与差异miRNA的成熟体取交集，得到66个miRNA。对66个miRNA分别进行miRTarBase、Targetscan和MiRDB靶基因预测后与差异mRNA取交集，并构建ceRNA网络（图4A）。ceRNA网络中有11个lncRNA（青色菱形），分别为LINC02576、LINC02086、AC020659.1、LINC00528、LINC00689、HOXB-AS4、MNX1-AS1、LINC00278、AC010624.1、AC016773.1、MYHAS；5个 miRNA（绿色三角形），分别为has-miRNA-206、has-miRNA-573、has-miRNA-3662、has-miRNA-133b、hasmiRNA-449a；12个mRNA（其中7个红色6边型为上调mRNA，5个蓝色圆形为下调mRNA），分别为STC2、PAX3、NETO2、EIF5A2、ADPRHL1、SYNM、ACTA1、GPR156、CLIC5、BMP3、HNF4A、FOXL2。将cor-mRNA钙信号通路中的16个基因进行转录因子富集，并绘制转录因子调控网络，网络中绿色椭圆形代表mRNA，蓝色菱形代表转录因子，这些转录因 子 分 别 是 MYOD、RSRFC4、TAL1BETAITF2、YY1、P300、PAX3、PAX5、ZID和OLF1（图4B）。

图4 ceRNA网络和转录因子网络

2.7 DNA甲基化分析

为了进一步阐明关键基因在喉鳞状细胞癌中的调控作用，进一步分析了TCGA和GEO中的喉鳞状细胞癌DNA甲基化数据。设Δβ绝对值≥0.4，P﹤0.05，取二者交集，共得到75个甲基化位点；其中cg15700197和cg17892556分别为甲基化差异最大的位点，cg15700197在正常组织中高甲基化，在喉鳞状细胞癌组织中低甲基化，cg17892556在正常组织中低甲基化，在喉鳞状细胞癌组织中高甲基化。这75个差异甲基化位点主要位于1stExon、TSS1500和TSS200，其次为5'非翻译区（untranslated region，UTR）和基因体区，3'非翻译区和基因间隔区（intergenic region，IGR）分布较少。上述位点对应高甲基化基因56个，低甲基化基因15个，ZNF625是Δβ值最大的高甲基化基因，OR10J1是Δβ值最大的低甲基化基因。

3 讨论

喉鳞状细胞癌是头颈部常见的肿瘤之一，病因较多，影响因素复杂。本研究针对喉鳞状细胞癌患者的一般情况进行分析，结果显示，患者年龄、T分期、病理分级、饮酒、吸烟情况对总生存率的影响不明显，而远处转移、淋巴结转移会降低患者的总生存率。关于肿瘤分期（T分期）和病理分级对肿瘤预后无影响，可能与目前分期系统的局限性有关；因此针对喉鳞状细胞癌的分期，根据肿瘤萌芽活性（tumor budding activity，BA）和细胞巢大小（cell nest size，CNS）进行分期对预后似乎更具指导意义。本研究结果中淋巴结浸润对喉鳞状细胞癌预后的影响与Ho等研究结果一致，尤其是病理组织上的淋巴结转移对指导临床治疗意义重大，当然这也可能与从临床到病理诊断的过程增加了“开始治疗时间（time to treatment initiation，TTI）”，从而降低了预后有关。

近年来有关lncRNA在喉鳞状细胞癌中作用的研究较多，尤其是miRNA“海绵”靶基因是喉癌发生、发展的重要机制之一。作为ceRNA，lncRNA通过靶向结合miRNA间接调控靶基因，是较好的生物标志物和潜在的治疗靶点。生存分析显示，miRNA-210、miRNA-455和miRNA-4326的表达量与喉鳞状细胞癌预后相关（数据未展示），其中有关miRNA-210在肺癌、胶质瘤、口腔癌、骨肉瘤等恶性肿瘤中的作用已广泛报道，并且miRNA-210在组织、血清和外泌体中均可存在。miRNA-455也可通过自噬、p21活化蛋白激酶2（p21-activated kinase 2，PAK2）轴和 Janus激酶 1（Janus kinase 1，JAK1）轴等作用方式抑制肿瘤生长，并且miRNA-455上述功能的发挥有赖于不同lncRNA的作用，这些lncRNA作为ceRNA竞争性结合miRNA-455；当然，也有较多miRNA-455单独发挥作用的报道。同时，Xu等亦报道miRNA-4326通过抑制APC2基因表达促进肺癌细胞增殖。但遗憾的是，目前尚未见上述3个miRNA在喉鳞状细胞癌中被报道，有待于进一步验证。

在具有抑癌作用的lncRNA中，LINC00689、LINC00278已分别在胶质瘤和食管鳞状细胞癌中报道，尤其是LINC00278编码的多肽YY1BM，可抑制YY1与雄激素受体结合从而降低真核生物延伸因子2激酶（eukaryotic elongation factor 2 kinase，eEF2K）的表达，促进肿瘤细胞凋亡。喉结是第二性征标志，雄激素受体高表达，由此可推断LINC00278在喉鳞状细胞癌中是一个较好的潜在靶点。在具有促癌作用的lncRNA中，IGFL2-AS1、LSAMP-AS1和MNX1-AS1已分别在胃癌、前列腺癌和食管鳞状细胞癌中报道，为下一步喉鳞状细胞癌的研究提供了较好的参考意义。

基于目前报道的有关lncRNA-miRNA-mRNA的相互作用关系，本研究利用差异lncRNA、差异miRNA和差异mRNA构建了ceRNA网络，通过该网络发现，网络中miRNA既无生存意义，又无单因素意义；网络中mRNA亦不属于Centiscape中的关键基因；同时，与miRTarBase疾病搜索中喉鳞状细胞癌的miRNA无交集。因此认为lncRNA并不能完全主导miRNA和mRNA发挥作用，lncRNA与miRNA在喉鳞状细胞癌发病过程中各自存在独立的作用机制。为了进一步说明lncRNA的功能，本研究选择了MYHAS作为代表，通过分析MYHAS共表达基因cor-mRNA，发现cor-mRNA得到的5个关键基因全部与差异mRNA得到的关键基因吻合；并且cor-mRNA的通路分析中主要信号通路钙信号通路、cGMP-PKG信号通路和cAMP信号通路亦是差异mRNA的主要信号通路。有关cAMP、cGMP-PKG和钙信号通路在不同肿瘤中的调控作用已有很多报道，有研究报道了cGMP-PKG在头颈部鳞状细胞癌中的作用，认为cGMP-PKG在头颈部鳞状细胞癌中是较好的治疗靶点，并且全球第一个鸟苷酸环化酶激动剂利那洛肽已于2019年初在中国上市。此外，Park和Juhnn亦报道cAMP信号激活可通过增加去乙酰化酶8的表达促进顺铂对肺癌细胞的凋亡作用。而关于钙信号通路在肿瘤中的作用，Raynal等研究表明，钙信号通路中一系列激酶相互作用后，引起钙/钙调素依赖蛋白激酶Ⅱ（calcium/calmodulin-dependent protein kinaseⅡ，CAMKⅡ）激活，使甲基化结合域蛋白甲基化CPG结合蛋白2（methyl-CPG binding protein 2，MeCP2）出核，导致基因组发生去甲基化，使已经甲基化的抑癌基因重新激活，以发挥抑癌作用。MYHAS属于抑癌基因，其通过相关基因cormRNA在喉鳞状细胞癌组织中参与钙信号通路的失活，从而引起下游靶基因高甲基化。基于上述结果，本研究又分析了数据库中喉鳞状细胞癌的DNA甲基化数据，共筛选出了71个基因；本实验室已经对ZNF667基因进行了试验验证，其在喉鳞状细胞癌组织和细胞中均呈高甲基化、低表达，发挥抑制喉鳞状细胞癌细胞增殖、侵袭的作用。由此认为MYHAS可以较好地代表喉鳞状细胞癌中的差异lncRNA，在喉鳞状细胞癌中是一个较好的生物标志物和治疗靶点，具有潜在研究价值。

同时，本研究针对cor-mRNA中参与钙信号通路的16个基因构建了转录因子调控网络，拟进一步说明非编码RNA在喉鳞状细胞癌信号通路中的上下游调控关系。这些转录因子与MYHAS相互作用，直接在转录水平调控靶基因表达。网络中PAX5是B淋巴细胞特异性转录因子，在B细胞淋巴瘤的发生、发展中具有重要作用，但也见于在其他恶性肿瘤中的报道，并且Brazao等报道了PAX5对不同阶段B细胞中lncRNA的调控作用。而转录因子PAX3的功能较为复杂，这是因为PAX3具有7个可变剪切体，不同剪切体在肿瘤中表现的功能不同。P300被认为是转录共激活因子，是一种组蛋白乙酰转移酶，通常被招募到转录增强子中，并通过乙酰化染色质来调节基因表达。Liang等研究报道，P300可在食管鳞状细胞癌中调控LINC00460的表达，这亦为喉鳞状细胞癌的研究提供了可行性。转录因子YY1是研究较多的肿瘤相关转录因子之一，其可调控lncRNA的表达，引起直肠癌、恶性胶质瘤、乳腺癌和肝癌等肿瘤的发生、发展，其中有lncRNANPCCAT1通过上调YY1促进鼻咽癌发生发展的报道，这对研究YY1在喉鳞状细胞癌中的作用具有重要参考意义。MYOD属于肌源性转录因子，在横纹肌肉瘤中与同类转录因子MYF5共同发挥重要作用；此外Dodd等亦报道了MYOD通过直接结合miRNA-182的启动子增加miRNA-182的表达，从而发挥促肿瘤作用。因此，在喉鳞状细胞癌中转录因子MYOD与lncRNA的作用还是值得期待。

综上所述，公共数据分析显示喉鳞状细胞癌中存在众多差异lncRNA和miRNA，其中lncRNA MYHAS通过抑制钙信号通路引起下游靶基因甲基化可能是喉鳞状细胞癌发生、发展的主要机制之一，lncRNA MYHAS是一个较好的生物标志物和治疗靶点，具有潜在研究价值。