摘要
基于公共数据库探究食管鳞癌（Esophageal squamous cell carcinoma，ESCC）中差异表达的类泛素化调节因子及潜在的功能。从TCGA数据库中获取1个食管鳞癌数据集并进行类泛素化调节因子的差异分析。用单因素Cox回归和套索回归（LASSO）筛选与预后相关的关键基因。利用Cox回归构建预后模型评估关键基因的独立预后价值。通过收集食管鳞癌患者癌组织与癌旁组织行RT-qPCR分析，初步验证筛出的差异基因。利用Gene MANIA 数据库对关键基因绘制蛋白质互作图（PPI），对网络中的相关基因进行基因功能和通路富集分析。基于TIMER算法进行关键基因和免疫细胞浸润的相关性研究。结果得到2个与食管鳞状细胞癌预后相关的类泛素化调节因子（DESI2和SMC6）。RT-qPCR结果提示 DESI2和SMC6在ESCC中高表达（P＜0.05）。PPI筛出的核心基因功能和富集分析表明，其可能涉及非同源末端结合和同源性重组加速DNA修复。最后通过TIMER算法发现DESI2和SMC6与肿瘤免疫浸润相关。结果提示DESI2和SMC6在ESCC中的表达具有明显的异质性，可能参与肿瘤进展，为ESCC发病机制、临床诊断和治疗靶点的研究提供理论依据。
Abstract
To investigate differentially expressed class ubiquitin-like modifier (Ub-like) regulatory factors and their potential functions in esophageal squamous cell carcinoma (ESCC) based on public databases, we conducted a study. We obtained one ESCC dataset from the TCGA database for Ub-like regulatory factor differential analysis. We used univariate Cox regression and LASSO regression to screen key genes related to prognosis. We used Cox regression to construct a prognostic model to evaluate the independent prognostic value of the key genes. We used RT-qPCR to initially validate the differential genes in ESCC patients by collecting cancer tissues and adjacent tissues. We used Gene MANIA database to draw a protein-protein interaction (PPI) map for the key genes and performed gene function and pathway enrichment analysis on the network's related genes. We studied the correlation between key genes and immune cell infiltration using the TIMER algorithm. Our study identified two Ub-like regulatory factors, DESI2 and SMC6, that were related to the prognosis of ESCC. The RT-qPCR results showed that DESI2 and SMC6 were highly expressed in ESCC(P＜0.05). The core genes identified by PPI and their functional and enrichment analysis suggested that they may be involved in non-homologous end joining and homologous recombination-accelerated DNA repair. Finally, we found that DESI2 and SMC6 were related to immune cell infiltration using the TIMER algorithm. The results suggest that the expression of DESI2 and SMC6 in ESCC has significant heterogeneity and may be involved in tumor progression, providing theoretical basis for the study of the pathogenesis, clinical diagnosis, and therapeutic targets of ESCC.
2020年新增食管癌病人604 100，死亡病人544 076，在胃肠道肿瘤中排名第五[1]。食管腺癌在西方国家发病率较高，而ESCC是东南亚患者的主要类型，在世界癌症相关死亡的常见原因中排名第六，我国是ESCC的高发区[2]。ESCC早期通常无症状，在诊断时许多患者已经是疾病晚期，治疗选择有限[2]。目前，一些参与ESCC发展的特异性标志物已经被确定，可能有助于早期诊断和预测预后，比如ESCC组织中梭状杆菌含量高。食管癌组织中miR-25、miR-424、miR-151呈高表达、miR-100、miR-99a、miR-29c、miR-140呈低表达，苯丙氨酸可作为ESCC的诊断标志物，热休克蛋白70，并被列为潜在血清标志物[3]。类泛素化（Small ubiquitin-like modifier，SUMOylation）是一种广泛存在的蛋白质翻译后修饰，该修饰存在于几乎所有的真核生物中[4]。SUMOylation可以调控DNA损伤修复、免疫反应、癌变、细胞周期进程和凋亡等生物学过程[5]。SUMOylation失调可导致肿瘤的发生和进展，这在前列腺癌、乳腺癌、卵巢癌等肿瘤研究中已有报道，主要参与诱导肿瘤细胞增殖、抗凋亡和转移潜能[6]。目前其在食管癌中的作用机制尚不明确，已有报道指出，HSP27的SUMOylation修饰调控PKM2促进食管鳞状细胞癌的增殖、侵袭能力[7]。
目前，ESCC中SUMOylation的研究相当有限，本研究利用公共数据库筛选ESCC组织中差异表达的类泛素化调节因子，并构建SUMOylation调节因子的预后模型，探索预后SUMOylation调节因子在ESCC中的潜在功能。
1 资料来源和方法
1.1 提取数据
基于公共数据库TCGA（https: //gdc-portal.nci.nih.gov）获得食管癌患者RNA-seq数据和临床资料，从UCSC数据库（https: //gdc-portal.nci.nih.gov）中下载表达数据，选择ESCA的数据集（https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Esophageal%20Cancer%20（ESCA）&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443），根据临床资料选择ESCC患者的数据。本研究纳入 2017 年 1 月至 2019 年 12 月在新疆医科大学第一附属医院胸外科接受手术的30例ESCC患者。纳入标准： ①病理诊断均为鳞状细胞癌；②无心血管、肝肾及造血等系统严重疾病；③无肿瘤病史。剔除标准：①采集标本前有放化疗等抗肿瘤治疗；②有其他恶性肿瘤者；③无法确切进行病理分期[JP2]者；由本院病理科医生确定癌旁组织（食管癌组织≥5 cm）为无癌正常组织。得到本院伦理委员会批准（编号：NKEBN/274/2009，审核通过日期：2009年10月5日）。
1.2 筛选差异SUMOylation调节因子
采用R软件“Limma”包[8]对TCGA数据库中食管鳞癌样本（n=82）和正常组织（n=11）进行校正，将已知的16个SUMOylation调节因子与TCGA数据库得到的DEGs进行比对，获得SUMOylation调节因子的表达，以（|Cor|＞0.7，P＜0.05）的阈值为显著相关。绘制相应DEGs和SUMOylation调节因子韦恩图和热图。
1.3 筛选与食管鳞癌预后相关的SUMOylation调节因子
利用差异SUMOylation调节因子的表达数据（FPKM），将具有生存信息的81例患者，按照7∶3比例随机分为训练集（57例）和验证集（24例）。在训练集中采用单因素Cox回归分析评估这些基因是否为危险因素。通过预后相关的基因（P＜0.05）构建LASSO回归，使用 R 软件“Glmnet”包，实现LASSO 逻辑回归，选择强相关特征，得到基因系数图形和交叉验证的误差图。当Lambda 值取最小时，即为最优模型。然后用多因素Cox分析和优化的风险评分模型，以中位风险评分为界分为高危组和低危组。最后，通过使用R 软件“Survminer”绘制高低风险组整体生存曲线，评估各组患者的生存情况；采用 ROC 计算模型的 AUC 面积以评估模型的有效性。
1.4 基于预后模型的独立预后分析
为了评估预后模型在ESCC患者的临床病理特征（T、N、性别、风险评分）的独立性，采用单因素和多因素Cox回归分析。以总生存期作为因变量计算危险比（HR）和95%可信区间以及P值，临床特征作为自变量。
1.5 建立与评估用于预测ESCC患者生存的Nomogram图
基于多因素Cox分析筛选出的独立预后因子，使用R软件“cph”包，以患者性别、N分期、Stage、RiskScore作为因素，构建一个评估ESCC患者1、3、5年OS概率的Nomogram图。同时采用校准曲线验证Nomogram图的准确性，由于食管鳞癌患者3~5年的存活的样本较少，因此在分析列线图的校准曲线的时候没办法进行估计，因此这里仅仅展示了1、2年的校准曲线。
1.6 利用GEPIA2数据库验证预后基因的表达和基因的生存分析
基于TCGA数据库分析预后基因在ESCC肿瘤组织和正常组织中的表达，绘制箱式图，并在外部数据集GSE53622数据集进行验证，选取30例ESCC患者临床标本对预后基因进行RT-qPCR实验。用Trizol裂解液裂解食管鳞癌和匹配的食管正常组织，利用RT Kit（Thermofisher，美国）和PCR技术，通过thermofisher 7500 PCR定量系统对SMC6、DESI2、GAPDH基因进行RT-qPCR。严格按照制造商说明书操作，随机将500 ng总RNA逆转录cDNA，总体积为20 ng，循环设定为94℃、15 s，94℃、5 s，58℃、34 s，40次循环。SMC6和DESI2的表达水平由GAPDH进行归一化处理，使用2^-△△CT方法计算。最后进行单基因表达情况与ESCC患者预后相关性分析。
1.7 预后基因在I、II、III期ESCC的ROC曲线下面积和不同临床阶段的表达差异
使用R软件“survivalROC”包分析SMC6和DESI2基因在ESCC中I、II、III期中ROC曲线，以评估SMC6和DESI2基因诊断特征的准确性。通过Wilcoxon检验分析SMC6和DESI2基因在不同临床阶段（T分期、N分期、TNM分期）的表达差异。
1.8 预后基因的共基因表达与免疫浸润的关系
利用GeneMANIA 数据库对预后基因及其 20 个互作基因进行了 PPI分析，利用R软件“clusterProfiler”包对共基因进行GO功能和KEGG途径富集分析，探索关键基因参与的功能通路条目[9]。为了了解关键基因对ESCC免疫浸润的影响，基于TIMER在基因模块中分析了2个关键基因的表达与免疫浸润丰度的相关性，在SCNA中分析体细胞CNAs与免疫浸润丰度的相关性。
1.9 统计分析
采用R语言（3.6.0）进行统计分析，用单变量Cox分析计算SUMOylation调节因子的HR和95%CI。建立套索回归模型，进一步筛选与预后相关的SUMOylation调节因子。以风险评分中位数设定阈值，构建高低风险组，评估不同风险组对生存时间的影响。采用Cox回归分析影响ESCC患者预后的临床特征因素。利用多变量Cox回归分析，构建了Nomogram图，以评估ESCC患者1、3 和5年OS的概率。采用t检验或单因素方差分析比较预后基因的表达差异。使用Kaplan-Meier生存分析法判断预后基因与ESCC患者预后的关系，采用Log-rank检验比较不同表达水平的生存差异，P＜0.05差异有统计学意义。
2 结果
2.1 SUMOylation调节因子的筛选
基于TCGA数据库获得ESCC样本82例和正常组织样本11例的基因表达数据。差异表达基因共有5 331个，上调基因1 487个，下调基因3 844个（P ＜0.05和|Log2FC|＞0.58）。将已知的16个SUMOylation调节因子映射到DEGs数据中，获得16个表达差异的SUMOylation因子（|Cor|＞0.7，P＜0.05）（图1）。
图1差异SUMOylation因子的表达
Fig.1Expression of differential SUMOylation genes
注：（a）韦恩图；（b）热图.
2.2 识别具有预后价值的SUMOylation调节因子并构建其预后模型
采用LASSO法筛选出与预后风险相关基因（DESI2、SMC6）构建预后模型（P＜0.05，图（2（a）~2（b））。通过构建预后基因模型进行风险评分，高危组患者生存期明显缩短，患者死亡率也相应升高（P＜0.05，图（3（a）~3（c））。此结果在验证组同样得到证实（P＜0.05，图（3（e）~3（g））。2个基因热图也显示了2个预后基因与风险评分的相关性（图3（d）、3（h））。
2.3 独立预后分析
为了进一步探究预后模型在ESCC患者中的独立预测价值。训练数据集中，单变量Cox回归显示风险评分、N分期、TNM分期等具有预后价值（P＜0.05，图（4（a）），多因素Cox回归分析提示风险评分、N分期、TNM分期等均不能成为ESCC的独立预后因素（P＜0.05，图（4（b））。总之，结果证实DESI2、SMC6的风险评分是ESCC患者预后的危险因素。
图2LASSO回归分析筛选预后基因
Fig.2LASSO regression analysis for screening of prognostic genes
注：（a）基因系数图形;（b）交叉验证误差图.
图3风险模型评分对患者生存的影响
Fig.3Impact of risk model scoring on patient survival
注：（a）生存分析;（b）风险曲线;（c）散点图;（d）热图;（e）生存分析;（f）风险曲线;（g）散点图;（h）热图.
图4Cox回归分析风险评分及临床病理特征的独立预后价值
Fig.4Cox regression analysis was used to analyze the independent prognostic value of risk score and clinicopathological features
注：（a）单因素COX回归分析;（b）多因素COX回归分析.
2.4 Nomogram图展示ESCC患者的预后
基于TCGA队列数据构建预测1年、3年和5年OS概率的Nomogram图（图5（a））。由于食管鳞癌患者3~5年的存活的样本较少，仅仅展示了1年、2年的校准曲线。结果提示1 年校准曲线斜率为0.783，2 年校准曲线斜率为0.448，1年、2年生存概率与实际结果一致性较好（图5（b））。
图5ESCC患者独立预后模型的展示
Fig.5Presentation of an independent prognostic model for ESCC patients
注：（a）独立预后模型列线图;（b）ESCC患者1年、2年的校准曲线.
2.5 预后基因在ESCC组织中的表达特征
利用TCGA数据、外部数据集GSE53622及30例ESCC患者临床组织标本分析了预后基因的表达水平。相较正常组织，预后基因在肿瘤组织中表达上调（P＜0.05，图6（a）~6（d））。预后基因在食管鳞癌不同临床阶段（T、N、M、TNM分期）的表达水平，结果发现SMC6表达水平在T分期之间差异有统计学意义；SMC6表达水平在N分期之间差异有统计学意义；SMC6在ESCC患者Ⅰ期和Ⅱ期，Ⅱ期和Ⅲ期之间均有差异（P＜0.05，图7（a））。DESI2的表达水平在T分期之间差异有统计学意义；DESI2在ESCC患者Ⅰ期和Ⅱ期之间差异有统计学意义（P＜0.05，图7（b））。
2.6 预后基因在ESCC中的诊断预后分析
通过ROC分析评估SMC6和DESI2基因对ESCC患者早期和中晚期的诊断效果，SMC6和DESI2基因的AUC值是0.618、0.609，AUC值均大于0.5（图8（a））。用Survival算法中的最优截断值，评估了SMC6和DESI2基因高、低表达分组情况与实际生存预后信息之间的关联性，关键基因在ESCC患者中的高表达与总生存率无明显相关（P＞0.05图8（b））。
图6DESI2，SMC6在ESCC中的表达
Fig.6Expression of DESI2 and SMC6 in ESCC
注：（a）DESI2和SMC6在GSEA53622数据集中的表达;（b）DESI2和SMC6在TCGA数据集中的表达;（c）SMC6在食管鳞癌肿瘤组织中的mRNA表达水平;（d）DESI2在食管鳞癌肿瘤组织中的mRNA表达水平.
图7SMC6和DESI2在ESCC不同临床阶段的表达
Fig.7Expression of SMC6 and DESI2 in different clinical stages of ESCC
注：（a）SMC6在食管鳞癌不同临床病理阶段（T、N、M、TNM分期）的表达水平;（b）DESI2在食管鳞癌不同临床病理阶段（T、N、M、TNM分期）的表达水平.
2.7 预后基因的共基因表达及富集分析
利用GeneMANIA数据库分析SMC6和DESI2基因的功能，确定20个与SMC6和DESI2基因显著相关的基因，分别是DESI1、SMC5、NSMCE1、NSMCE2、NSMCE3、NSMCE4A、EID3、RAD50、RPS7、RAD18、TUFT1、SLF2、TTC19、PKMYT1、YIF1、BREXO2、TP53BP2、ATRX、FANCC（图9）。对20个共基因进行GO功能和KEGG信号通路富集分析，共筛选得到了144项显著相关的生物学过程、22项细胞成分、38项分子功能和3条KEGG信号通路（图10（a）~10（d））。
图8预后基因在ESCC患者中ROC曲线和K-M生存分析
Fig.8ROC curve and K-M survival analysis of prognostic genes in ESCC patients
注：（a）基因SMC6、DESI2的ROC分析曲线，（b）预后基因的K-M生存分析.
图9利用GeneMANIA数据库分析SMC6和DESI2基因的基因-基因相互作用网络
Fig.9Gene-gene interaction network of SMC6 and DESI2 genes was analyzed using the GeneMANIA database
图1020个共基因的GO功能和KEGG信号通路富集分析
Fig.10GO function and KEGG signaling pathway enrichment analysis of 20 co-genes
注：（a）G0分析生物学过程气泡图;（b）G0分析生物学过程气泡图;（c）G0分析细胞成分气泡图;（d）KEGG分析条形图.
2.8 预后基因与免疫浸润
为了解关键基因对ESCC免疫浸润的影响，基于TIMER数据库分析2个关键基因的mRNA表达与免疫浸润水平的相关性。结果提示DESI2的水平与巨噬细胞浸润正相关，DESI2的水平与CD4+T细胞、中性粒细胞及树突状细胞浸润负相关（图11（a））；SMC6的CNV与CD4+ T细胞、中性粒细胞及树突状细胞密切相关（图11（b））。
3 讨论
食管鳞状细胞癌是世界上最高发的恶性肿瘤之一，我国是食管鳞癌的高发区，50%以上的患者确诊时都已经处于肿瘤晚期，治疗方法的选择有限[10]。寻找有效的早期诊断标志物和治疗靶点仍是我们继续探索的方向。SUMOylation作为翻译后修饰的关键之一，涉及多种细胞动态调节，调控几乎所有真核生物的基因组完整性、转录调控和细胞内信号转导，受到越来越多的关注[11]。SUMOylation途径的修饰是通过SUMO或泛素分子与蛋白质底物的可逆动态结合而发生的一系列生化级联反应来调节蛋白质功能，涉及的关键酶主要有E1（激活酶）、E2（偶联酶）、E3（连接酶）和去SUMOylation蛋白酶[12]。SUMOylation在人类肿瘤进展当中也扮演了相当重要的角色，参与肿瘤的增殖、抗凋亡、转移潜能等生物学功能[13]。SUMOylation途径的E1激活酶也称为SAE1/2、E1（SAE1和SAE2）在肝癌肿瘤组织中上调，与其较短生存期有关[14]。也有研究表明环状RAPGEF5通过稳定SAE1促进肝内胆管癌肿瘤增殖、转移和SUMOylation。目前，UBC9是唯一一种可进行SUMOylation的E2连接酶，研究表明UBC9能稳定PFKFB3，促进有氧糖酵解和胶质母细胞瘤细胞的增殖[15]。底物蛋白的SUMOylation可以被SENP家族可逆地逆转，SENPs已被证明在正常细胞和包括癌症在内的各种疾病的各种病变细胞中发挥重要作用[16]。前列腺癌的相关研究指出，SENP1通过介导HIF1α信号通路调控骨重塑蛋白的表达，进而影响骨转移[17]。肺癌的相关研究指出，SENP2靶向核Dbf2相关蛋白2（NDR2）去泛素化，提高NDR2激酶活性，进而导致下游靶点p21的不稳定，加速G1/S细胞周期转换，表明SENP2是未来治疗肺癌较有前途的治疗靶点[18]。此外，SUMOylation是肿瘤相关蛋白功能障碍的一个关键点，比如，在乳腺癌的研究中指出，ZNF131是SUMOylation修饰的靶点，ZNF131的SUMOylation修饰减弱了乳腺癌细胞系中雌激素诱导的细胞增殖[19]。结肠癌相关研究指出，MAFB在lysine32位点被SUMO1 类泛素化，进而调控结肠癌中细胞周期的关键，提示MAFB及其SUMOylation过程可能成为结肠癌治疗的靶点[20]。但是，目前在食管癌中关于SUMOylation的研究非常有限。
图11基于TIMER数据库分析SMC6和DESI2基因和6种免疫细胞之间的相关性
Fig.11[JP3]Correlation between SMC6 and DESI2 genes and six types of immune cells were analyzed based on TIMER database
注：（a）DESI2 mRNA表达与B细胞、CD8+T细胞、CD4+T细胞、巨噬细胞、中性粒细胞、树突状细胞;（b） SMC6 mRNA表达与B细胞、CD8+T细胞、CD4+T细胞、巨噬细胞、中性粒细胞、树突状细胞;（c）SMC6基因没有CNV结果.
本研究利用公共数据库的ESCC数据和30例临床组织标本分别验证了DESI2、SMC6的mRNA表达水平在ESCC肿瘤组织中上调，且SMC6、DESI2和ESCC肿瘤分期密切相关，尤其是早期和中晚期之间。ROC分析提示DESI2、SMC6的AUC值均＞0.5。K-M生存分析提示SMC6、DESI2在ESCC患者中的高表达与较短生存期正相关。近期研究报道，去类泛素化异肽酶（DESI2）在DNA损伤的早期反应中被激活，而且参与肿瘤的发生和进展[21]。在胰腺癌的研究中指出，DESI2在肿瘤组织中低表达，上调DESI2可抑制细胞增殖、细胞周期和侵袭等恶性表型[22]。在肺腺癌组织中DESI2低表达，上调DESI2可使细胞组织在S期，阻滞细胞从S期向G2/M期循环，促进细胞凋亡[23]。此外，在肺腺癌研究中指出，DESI2可显著增强免疫反应诱导剂的抗肿瘤活性[24]。鉴于此，本研究也通过TIMER数据库探索DESI2与食管癌肿瘤免疫的潜在联系，发现DESI2的水平与CD4+T细胞、巨噬细胞、中性粒细胞及树突状细胞密切相关（P＜0.05）。在SMC复合物家族中，SMC5/6在调节DNA复制和修复中具有独特的作用[25]。目前，关于SMC5/6在肿瘤中独特的分子调节机制的研究十分少。近期在人肉瘤的研究中指出，SMC5、SMC6的mRNA水平高表达，但预后生存分析提示SMC5、SMC6的表达对肉瘤的OS、DFS无显著影响[26]。在肝癌的研究中指出，SMC6在肝癌组织中mRNA和蛋白水平均上调，且其表达水平与B细胞、CD4+T细胞、巨噬细胞、中性粒细胞密切相关[27]。而本研究指出SMC6的CNV与CD4+ T细胞、中性粒细胞及树突状细胞密切相关（P＜0.05）。
此外，通过GeneMANIA为SMC6、DESI2构建了一个网络，其中包括与其密切相关的20个其他基因，进行了GO功能富集分析，发现它们可能与细胞周期生物功能相关。KEGG富集分析，发现他们可能参与非同源末端结合和同源性重组。比如，在神经胶质瘤细胞中，发现SUMOylation调节因子参与mRNA拼接、DNA复制、ATP酶活性和剪接体[28]。所以针对类泛素化调控因子SMC6、DESI2在食管癌中的具体作用及调控机制值得进一步深究。
本研究还有不足之处。比如SMC6和DESI2基因高、低表达分组情况与实际生存预后信息之间的关联性不显著，考虑到本研究中其与食管鳞癌组织中的表达特点，且食管癌在新疆地区的区域特异性，有必要探索其在不同民族食管鳞癌患者组织中的特异性。未来我们将扩大临床组织样本，收集新疆地区不同民族、不同病理组织，探索SMC6、DESI2的表达水平与临床特征及预后的相关性，通过体外实验证明SMC6和DESI2在食管鳞癌中的生物学功能及具体分子作用机制。
