Nanopore测序揭示m6A修饰铁死亡相关基因在中晚期头颈部鳞状细胞癌中的表达及临床意义

庄黎明,丁雅婷,蔡成福,叶辉,蔡耿明

(1.福建医科大学附属泉州市第一医院 耳鼻咽喉头颈外科,福建 泉州 362000;2.泉州市医学高等专科学校附属人民医院 耳鼻咽喉科,福建 泉州 362000;3.厦门医学院附属海沧医院 耳鼻咽喉头颈外科,福建 厦门 361026;4.福建医科大学临床医学部,福建 福州 350122)

头颈癌(head neck cancer,HNC)是世界上第8种最常见的癌症。根据2020年全球癌症统计,HNC的发病人数约为870 000例,占所有恶性肿瘤的4.5%。HNC的发病机制是一个多步骤的过程,涉及分子变化的渐进积累。进一步阐明HNC发展中涉及的分子事件,可能有助于识别潜在和有效的生物标志物,并为靶向治疗提供新的方向[1]。头颈部鳞状细胞癌(head neck squamous cell carcinoma,HNSCC)是其中发生率最高的一类。标准疗法如保留器官的微创手术方法、放疗和多模式治愈性疗法等的改进已经可以更好地保留患者的头颈部器官的功能,并提高患者的5年生存率[2-4],但是复发性或转移性患者的疗效及预后仍差强人意。基于T细胞的免疫治疗,已被证明能提高复发性或转移性HNSCC患者的总生存率[5],但只有很低比例的患者对此疗法有响应,还存在免疫耐受和免疫逃逸等问题,因此,转移性、复发性和晚期HNSCC患者的预后仍然不理想。对HNSCC进展和转移形成的机制已经有深入的研究,但仍需要开发新的治疗策略[2]。

RNA修饰,例如N6-甲基腺苷(N6-methyladenosine,m6A),在表观遗传基因调控和细胞功能中发挥着不可忽视的作用,与许多人类疾病密切相关。m6A是在真核细胞中记录的最常见的转录后RNA修饰,这些修饰大多发生在5’-和3’-非翻译区(UTRs)以及内部长外显子的信使核糖核酸终止密码子附近,在许多生理和病理过程中起着重要作用,包括产生免疫反应、微小RNA的编辑和各种癌症的进展[6]。目前已有研究表明,m6A修饰失调与胶质母细胞瘤、乳腺癌症、胃癌、结直肠癌及甲状腺乳头状癌[7]等肿瘤的发生发展密切相关。有研究表明[8],喉鳞状细胞癌(laryngeal squamous cell carcinoma,LSCC)患者的预后,可以通过m6A相关的LncRNA表达量构建的预后模型预测。然而,m6A修饰在HNSCC进展中的作用还需要进一步研究[9]。

1 材料与方法

1.1 分析流程

分析过程流程见图1。

图1 数据收集分析过程 注:HNSCC(头颈部鳞状细胞癌);GO(基因本体);KEGG(京都基因与基因组);PPI(蛋白相互作用);GSVA(基因集变异分析)。

1.2 数据处理

1.2.1 RNA提取、文库制备和测序 RNA提取、文库制备和高通量测序的数据分析均由中国有限公司赛格健康科技有限公司进行。按照Chomczynski等[10]的方法,使用TRIzol(InvitrogenTM,Cat.No.15596018),从临床收集的3对癌及癌旁组织样品中提取总RNA。在DNaseI(Thermo ScientificTM, Cat.No.EN0521)提取RNA后进行DNA消化。RNA质量通过用NanodropTMOneC分光光度计(Thermo ScientificTM, Cat.No.ND-ONEC-W)检测A260/A280来测定。RNA经证实完整性及合格后通过QubitTM3荧光计(InvitrogenTM, Cat.No.Q33216)与QubitTM RNA HS测定试剂盒(InvitrogenTM, Cat.No.Q32855)进行定量。对应于250～500 bp的文库产物进行富集、定量,并最终在NovaSeq(Illumina®)上测序。

1.2.2 RNA序列数据分析 原始测序数据首先通过SOAPnuke(版本1.6.0)进行过滤,去除低质量及被污染的读数。Clean Reads首先根据唯一分子标识符(unique molecular identifier,UMI)序列进行聚类,相同UMI序列分组到同一聚类。生成所有子聚类后,进行多序列比对,得到每个子聚类的一个一致序列。

1.2.3 新基因和转录物的预处理、比对和分析 使用Guppy软件(版本5.0.16)[11]对原始读取进行基础调用,评估读取质量,Guppy还对适配器序列进行了修剪。Nanofilt(2.7.1版)过滤低质量(Q值<7)和短长度读数(<50bp)[12]。使用Fclmr2(版本0.1.2)和短长度RNA-Seq数据校正剩余的干净读数[13]。然后,使用Minimap2(2.17-r941版)将每个文库的干净读数与其基因组进行比对[14]。使用Samtools(1.10版)计算干净读数与参考基因的比对率[15]。为了减少结果的冗余,合并仅在5’端外显子差异的比对,并使用StringTie软件(版本:2.1.4;参数:-保守性-L-R)获得非冗余转录序列。Gffcompare软件(版本:0.12.1;参数:-R-C-K-M)用于将转录物与基因组的已知转录物进行比较,并发现新的转录物和新的基因。为了获得新转录物的全面功能信息,根据序列相似性和基序相似性对7个数据库的转录物进行了注释,包括NR、Pfam、UniProt、京都基因与基因组(kyoto encyclopedia of gengs and genomes KEGG)、基因本体(gene ontology,GO)、KOG/COG和pathway。

1.2.4 基因转录物的结构分析 将转录本与已知的基因组转录本进行比较,以校正边界。

1.2.5 测序质量控制 对牛津纳米孔测序技术(oxford nanopore technologies, ONT)[16]原始测序数据进行质量控制。

1.2.6 数据过滤 根据二代及三代数据过滤标准,对数据进行过滤,从而得到有效数据。

1.2.7 功能注释和通路富集分析 对上述对照组鉴定的差异表达基因(differentially expressed gene,DEG)进行功能注释,GO富集分析使用了GO术语的注释和可视化和元景观。GO术语的不同表达基因列表之间的重叠通过富集分析圆图进行。然后将DEG引入FunRich(功能富集分析工具)(http://www.funrich.org/)用于KEGG通路分析。GENEMANIA 用于构建DEG的基因-基因相互作用网络,以评估这些基因的功能。

1.2.8 蛋白-蛋白相互作用(protein protein interaction,PPI)网络的构建 交集的基因通过string V11.0进行基因之间的相互作用预测,得到它们之间的关系,筛选标准为combined_score>=0.4。并通过Cytoscape软件(V3.5.1)进行网络图的构建可视化。 共表达的hub基因通过对PPI中的基因进行共表达分析,从而发现它们之间的表达关系,进而发现它们在此实验中的调控关系。

1.2.9 基因集变异分析(gene set variation analysis,GSVA) 使用软件R(版本3.6)的GSVA包分析m6A。

2 结果

2.1 转录组差异基因与差异m6A的差异基因的交集

首先通过癌与癌旁组织对比得到转录组差异基因筛选出来,标准为P<0.05且log2>1或者log2<-1。将m6A的差异基因筛选出,标准为P<0.05且log2>1或者log2<-1。通过两组差异基因筛选出转录组上调和m6A上调的基因以及转录组下调和m6A下调的基因。我们得到96个上调基因(转录组上调m6A上调)和159个下调基因(转录组下调m6A下调)。

2.2 交集基因的富集分析

首先将交集基因按照上调和下调基因来做富集分析,从而得到它们各自的方向,具体见表1。

表1 交集基因的富集分析

对这部分的基因进行功能以及pathway的富集分析,对显著性差异基因功能分别做了柱状图(图2),显著性筛选的标准为P<0.05。

图2 部分基因功能及pathway富集分析中显著性差异的基因 注:纵坐标为差异基因功能名称, 横坐标为P值的负对数(-LgP);图中仅显示按照-LgP排序最大的各前10项。collagen-containing extracellular matrix(含胶原蛋白的细胞外基质);secretory granule lumen(分泌颗粒腔);cytoplasmic vesicle lumen(细胞质囊腔);vesicle lumen(囊泡腔);blood microparticle(血液微粒子);tertiary granule lumen(三级颗粒腔);microvillus membrane(微滤膜);endoplasmic reticulum lumen(内质网腔);apical part of cell(细胞顶端);apical plasma membrane(顶端质膜);extracellular matrix structural constituent(细胞外基质结构成分);glycosaminoglycan binding(糖胺聚糖结合);heparin binding(肝素结合);extracellular matrix structural constituent conferring compression resistance(细胞外基质中具有抗压性的结构成分);oxidoreductase activity(氧化还原酶活性);collagen binding(胶原蛋白结合);endopeptidase inhibitor activity(内肽酶抑制剂活性);peptidase inhibitor activity(肽酶抑制剂活性);endopeptidase regulator activity(内肽酶调节器活性);sulfur compound binding(硫化合物结合);complement activation, alternative pathway(补体激活,替代途径);O-glycan processing(O型糖加工);iron ion homeostasis(铁离子稳态);body fluid secretion(体液分泌);secretion by tissue(组织分泌);protein O-linked glycosylation(蛋白质 O 链接糖基化);transition metal ion homeostasis(过渡金属离子稳态);iron ion transport(铁离子运输);regulation of ossification(骨化调节);lung development(肺发育)。下同。

2.3 差异基因信号通路分析

根据差异基因通路中每一通路的显著性水平(P值)和所包含的差异基因数据(表2)构建显著性差异基因通路的分布图。对部分显著性差异基因通路做了气泡图(图3)。从通路的结果中可以看出:在癌组织中PPAR通路下调,免疫相关通路下调,以及代谢通路的上调,从而影响了肿瘤的发展和增殖。

表2 差异基因数据分析

图3 显著性差异基因通路气泡图 注:纵坐标为差异基因通路名称,横坐标为富集程度(enrich factor);mucin type O-glycan biosynthesis(黏蛋白型 O-糖的生物合成);complement and coagulation cascades(补体和凝血级联);renin secretion(肾素分泌);PPAR signaling pathway (PPAR信号通路);glycosphingolipid biosynthesis-lacto and neolacto series(糖磷脂的生物合成-内酯和新内酯系列);staphylococcus aureus infection(金黄色葡萄球菌感染);ferroptosis(铁死亡);fatty acid degradation(脂肪酸降解);other types of O-glycan biosynthesis(其他类型的 O-糖生物合成);phosphonate and phosphinate metabolism(膦酸盐和膦酸盐代谢);intestinal immune network for IgA production(产生 IgA 的肠道免疫网络);malaria(疟疾)。下同。

2.4 PPI的建立

通过string V11.0对上述交集的基因进行基因之间的相互作用预测,得到它们之间的关系,并通过Cytoscape软件进行网络图的构建可视化构建,点用来表示基因,线用来表示基因之间的关系,其中颜色表示上下调的关系和程度,红色表示上调基因,而绿色表示下调基因,颜色越深表示变化越显著。

通过构建差异基因的作用网络,得到基因之间的作用关系。基因的度排序结果见表3。通过癌与癌旁组织进行差异基因分析发现,DCN、GLU、ARG2、C3均呈现下调趋势,而CDK1、APOE、POSTN、RPL8、YWHAZ、MMP13均呈现上调的趋势(表3)。

表3 基因的度与作用趋势

2.5 m6A_GSVA

通过癌和癌旁组织的6个样本的转录组测序技术(RNA sequencing,RNAseq),其每千个碱基的转录每百万映射的Transcript(Transcript per Kilobase per Million mapped reads,TPM)值进行GSVA分析(GSVA:用来体现每个样本的在某个通路的变化情况)。运用t检验方法筛选,通过基因的信号值预测每个样本的富集得分值,计算两组之间的富集差异性,得到有差异富集的通路。从图4可见PI3K/AKT信号通路为在肿瘤中上调,PPAR通路在肿瘤中为抑制通路。

图4 6个样本在不同分型的样本富集得分值 注:图中颜色是由蓝向红来变化,表示富集分数数值的升高,表示通路被激活。

2.6 共表达的hub基因

通过对PPI中的基因进行共表达分析,从而发现它们之间的表达关系,进而发现它们在此实验中的调控关系。通过计算PPI中的基因之间的共表达关系,具体见图5。从图5中PPI共表达关系,我们发现AGR2、CDK1在m6A调控下,与周围的基因联系很密切,具体结果见表4。

图5 PPI中基因的共表达关系 注:深蓝色表示基因的度越大,而浅色表示基因的度越小。

表4 共表达基因的度

3 讨论

由于肿瘤的异质性、对包括免疫治疗在内的多种治疗效果不佳,以及耐药性等问题的存在。HNSCC患者,特别是复发性和中晚期患者的治疗效果和预后并不令人满意,适当的个体化治疗方案需要进一步去探索和建立,从而有效提高和改善HNSCC患者的治疗效果和预后。既往的研究表明,m6A调节因子的异常表达可能与多种癌症的发生发展有关,包括乳腺癌症、膀胱癌症、胶质瘤和结直肠癌等。与正常组织相比,某些m6A相关基因在HNSCC组织中的表达上调,与免疫细胞浸润以及免疫检查点抑制剂(immune checkpoint inhibitors,ICIs)治疗的反应和结果有关[6]。从PPI共表达关系图中,我们发现在m6A调控下,共表达基因与周围的基因联系很密切,它们处于网络的中心位置。

细胞周期蛋白依赖性激酶1(cell cycle protein-dependent kinase 1, CDK1)是一种丝氨酸/苏氨酸激酶,是细胞周期的主调节因子[17],它对细胞周期调节至关重要。CDK1在控制DNA复制、修复及mRNA转录等功能作用中均能起到一定的作用[18]。一般情况下,细胞周期控制的失调主要引起肿瘤生长[19],CDK1的上调与肿瘤的增殖和进展相关。一项研究表明[20],铁死亡相关基因与CDK1密切相关。并且有研究已验证了CDK1 可使LA 相关蛋白 1(LARP1)磷酸化[21]。几项研究均显示了LARP1在几种恶性肿瘤中显著上调[22-25]。与之相关的是,LARP1可作为CDK1的底物在部分癌症中上调,从而提高异常翻译存活率,导致癌细胞增殖增加。载脂蛋白E(apolipoprotein E,ApoE)是一种参与脂质转运和脂蛋白代谢的多功能蛋白,介导组织和细胞中的脂质分布/再分配。它还可以调节炎症和免疫功能,保持细胞骨架的稳定性,改善神经组织功能。近年来,ApoE频繁出现在肿瘤研究中,并逐渐成为肿瘤的生物标志物。研究发现ApoE在许多肿瘤细胞,包括神经胶质瘤、卵巢肿瘤、口腔鳞状细胞癌症、胰腺导管细胞癌等大多数实体瘤组织中高表达。ApoE可以调节巨噬细胞的极化变化,参与肿瘤免疫微环境的构建,调节肿瘤炎症和免疫反应,在肿瘤侵袭和转移中发挥作用。许多研究已经证实ApoE促进肿瘤的发展、增殖和转移,但最近的研究也发现ApoE具有保护作用[26]。Periostin(POSTN)属于基质细胞蛋白,分泌到细胞外环境中的非结构ECM蛋白,在大多数成年组织中以低水平表达。这些蛋白质与细胞表面受体相互作用,介导细胞和细胞外通讯。Periostin在各种实体上皮肿瘤中过表达,其与细胞表面受体整合素的相互作用调节细胞内信号通路,对癌症有直接影响。它在转移中上调,并可影响转移病灶的大小和数量[27]。编码核糖体蛋白(ribosomal protein L8,RPL8)是60S亚基的组成部分,该蛋白属于核糖体蛋白的L2P家族,位于细胞质中。RPL8是一种铁死亡相关基因,它的表达随胰腺肿瘤分级而显著改变,与胰腺癌症预后相关,与胶质瘤癌症进展有关[28-29]。血小板参与肿瘤血管生成和癌症进展,是癌症生物标记物的潜在来源。研究发现,RPL8在癌症患者的血小板蛋白质组中上调[30]。

簇蛋白是一种进化高度保守的糖蛋白,簇蛋白通过调节多种信号通路介导各种癌症(如前列腺癌、乳腺癌、肺癌、肝癌、结肠癌、膀胱癌和癌症)中的癌症进展,在调节几个基本生理过程中发挥了突出作用,包括程序性细胞死亡、转移、侵袭、增殖和细胞生长[31]。簇蛋白介导的线粒体自噬可以抑制口腔癌症的干性和自我更新能力[32]。它在神经保护和癌症以及化疗耐药性方面具有重要意义,被认为既是一种生物标志物,也是一种非常可探索的新治疗靶点,适用于多种疾病[33]。但簇蛋白在恶性肿瘤中的作用一直存在争议,它可以促进体外的肿瘤发生和对化疗药物的耐药性。然而, 簇蛋白在小鼠体内也可能会限制神经母细胞瘤和癌症的发展,研究发现,它在抑制细胞死亡、调节促生存信号传导和增强肿瘤细胞化疗耐药性方面至关重要。簇蛋白高表达与膀胱癌患者生存率低明显相关,被定义为膀胱癌中的一种癌基因。它通过抑制细胞凋亡来增强癌症对顺铂的耐药性[34]。蛋白聚糖由铁死亡的细胞释放,然后作为警报信号触发先天和适应性免疫反应。铁死亡期间蛋白聚糖的早期释放是一个活跃的过程,涉及分泌性大自噬/自噬和溶酶体胞吐作用。一旦释放,细胞外蛋白聚糖与巨噬细胞上的受体高级糖基化终产物特异性受体(AGER)结合,以NFKB/NF-κB依赖的方式触发促炎细胞因子的产生[35]。YWHAZ(也称为14-3-3ζ)是许多信号转导途径的中心枢纽蛋白,在肿瘤进展中起着重要作用。越来越多的证据表明,YWHAZ在多种类型的癌症中经常上调,并在细胞生长、细胞周期、凋亡、迁移和侵袭等广泛的细胞活动中起致癌基因的作用。此外,YWHAZ可能是几种癌症的诊断、预后和化疗耐药性的潜在生物标志物[36]。临床上, YWHAZ高表达提示胃癌、膀胱尿路上皮癌、乳腺癌、前列腺癌、结肠癌、非小细胞肺癌等多种癌症患者预后差有关[37-42]。经典的铁死亡诱导剂erastin促进了软骨细胞中基质金属蛋白酶13的表达,同时抑制了Ⅱ型胶原的表达[43-44]。AGR2基因位于染色体7p21.3内,是一种蛋白质二硫键异构酶(protein disulfide isomerase,PDI),PDI家族成员作为氧化还原酶而发挥作用,负责催化二硫键的形成或重排,这是发生于细胞内质网中的一种关键蛋白质修饰。许多研究已证实,AGR2高表达于含有黏液分泌功能细胞的组织。肿瘤微环境(tumour microenvironment,TME)中成分复杂,包含肿瘤细胞、免疫细胞、肿瘤相关炎症细胞和成纤维细胞等细胞,此外还有组织间质、血管、细胞因子和趋化因子等。颈癌研究中发现,AGR2的表达于与ALDH1、Sox2和Oct4之间存在相关性,表明AGR2可能与头颈癌有关[45]。此外,AGR2在各种上皮肿瘤中表达上调,通过激活包括p53激活途径在内的多种途径促进肿瘤的侵袭性生长、血管生成、侵袭、转移和耐药性等恶行生物学行为[46]。

有趣的是,与已有研究相反,测序结果发现CLU、AGR2在中晚期HNSCC癌组织中表达,我们认为:① 实验所收集的标本为中晚期HNSCC标本,与上述研究的标本有所不同。②它们可以逃离这个位置,在细胞表面或细胞外基质中发挥作用。③转化生长因子-β(transforming growth factor-β,TGF-β)的缺失可以抑制AGR2。因为除了通过关键信号系统诱导AGR2外,还可以通过TGF-β的介导调节AGR2,即TGF-β/AGR2轴。TGF-β的缺失可以抑制AGR2从而可以抑制p53肿瘤抑制活性[20]。④测序显示 FOXA1低表达,而功能研究表明,FOXA1在前列腺癌中存在两种相悖的功能:FOXA1通过雄激素受体(androgen receptor,AR)途径诱导前列腺癌细胞繁殖,同时FOXA1也可以通过不依赖AR途径抑制癌细胞迁移和上皮-间皮转化(epithelialto-mesenchymal transition,EMT)。AR及其辅因子FOXA1与一个特定的内含子雄激素响应元件结合,通过MBOAT2调节前列腺癌的铁死亡[47]。FOXA1是Yes相关蛋白(yes-associated protein,YAP)依赖性的转录因子,在转录因子CP2(transcription factor CP2,TFCP2)作为YAP调控转录因子的情况下,可以于YAP一起调节FOXA1,诱导铁蛋白轻链(ferritin light chain,FLT)转录,从而抑制铁死亡[35,48]。

我们测序的数据结果表明,m6A修饰与FOXA1、RPL8、DCN等铁死亡相关基因紧密相关。铁死亡是铁依赖性的,以脂质过氧化为主要特征的一种区别于凋亡细胞死亡的一种独特的程序性死亡类型。已在包括癌症等的自身免疫疾病中被探索研究,并被作为新的治疗策略[49-50]。因此我们推测,m6A修饰铁死亡相关基因在中晚期HNSCC中起重要作用,这需要进一步的实验深入研究证实。本研究探讨了m6A介导的中晚期HNSCC表观转录组学,对中晚期HNSCC的免疫逃逸进行初步的探讨,为HNSCC的个体化精准治疗提出新的观点。