一种新型黑色素瘤的蛋白质预后模型

陈雅静，张桐玮，周 俊，陈舒曼，蒲仕明*

(1.广西师范大学 生命科学学院, 广西 桂林 541006；2.广西高校干细胞与医药生物技术重点实验室(广西师范大学),广西 桂林 541004;3.广西师范大学 生物医学研究中心，广西 桂林 541004;4. 广西珍稀濒危动物生态学 重点实验室(广西师范大学), 广西 桂林 541006)

恶性黑色素瘤是最具侵袭性的癌症之一，其发病率呈不断增加趋势[1-2]。近年来，黑色素瘤相关研究进展很快，例如Nivolumab作为一种抗PD-1抗体，有助于提高对癌细胞的免疫反应，已被FDA批准用于治疗晚期黑色素瘤[3]。其他免疫疗法，例如monalizumab——一种增强NK细胞和CD8+T细胞活性的人源化抗NKG2A抗体，也在开发中[4]。尽管使用了包括靶向疗法和免疫疗法在内的不同治疗方法，但黑色素瘤患者的预后仍较差。

目前，已报道了多个指导黑色素瘤个体化治疗的预后模型[5]。研究表明，肿瘤厚度与患者存活率之间存在极高的相关性[6]。早期识别、积极治疗溃疡和密切随访对于黑色素瘤的预防和改善预后至关重要[7]。有无前哨淋巴结转移是临床Ⅰ期和Ⅱ期黑色素瘤中的独立预后因素[8]。然而，这些根据患者临床特征，包括肿瘤的厚度、是否伴随溃疡和原发肿瘤的部位等基于组织学特征构建的预后模型并不能很好地预测每个个体的肿瘤预后[9-10]。

近年来，众多研究相继发现了一系列与黑色素瘤预后相关的分子标志，例如特异基因的异常表达(包括蛋白质和mRNA水平)与非编码lncRNA和miRNA等。许多研究还发现调节因子m6A通过调节lncRNA来影响癌症的进展和预后[11]。例如，一些lncRNAs可以区分黑色素瘤的亚型并预测其生存情况。过表达的lncRNA HCP5还可以在体外降低黑色素瘤细胞恶变的概率，这可能与RARRES3上调有关[12]；血清微miRNA作为黑色素瘤复发的生物标志物[13]。然而，仅使用一种生物标志物来预测患者预后的敏感性和特异性并不能达到理想效果，使用多种蛋白质生物标志物的预后模型将在黑色素瘤患者的诊断和预后中显示出巨大的潜力。

TCPA(the cancer proteome atlas，https:∥www.tcpaportal.org/)是一个患者肿瘤样本的蛋白质组学数据库，TCGA(the cancer tenome atlas，https:∥portal.gdc.cancer.gov/)则是癌症患者的临床特征数据库。近来，许多研究人员对这2个数据库进行联合分析，建立了不同肿瘤的蛋白质预后模型，发现可以很好地预测肿瘤患者的预后。Fang等[14]利用Kaplan-Meier生存分析、单变量和多变量Cox分析，建立了肺鳞状细胞癌(LUSCC)患者的预后预测模型。根据风险模型中每个蛋白质的系数计算出每个患者的风险值，能有效预测LUSCC患者的预后。 Tian等[15]通过利用TCPA、GEO(gene expression omnibus)和TCGA获得了成人肾上腺皮质癌(adrenocortical carcinoma，ACC)患者的RPPA数据、基因表达谱和临床信息，并结合上海癌症中心免疫组合数据库数据，建立综合预后相关蛋白质(IPRPs)模型。该预后模型在多队列中显示出比其他生物标志物(Ki-67、β-catenin等)更强的预测价值，提示了利用TCPA和TCGA等数据库建立的疾病预后模型，对肿瘤相关蛋白质风险与肿瘤治疗预后具有较好的指导意义。本文通过TCPA和TCGA这2个数据库建立一个基于11种蛋白质表达数据的蛋白质预后模型，并通过B16黑色素瘤细胞系开展相关试验，对部分蛋白质的风险相关性进行验证。该模型的建立拟为黑色素瘤患者的预后判断提供参考。

1 材料和方法

1.1 数据采集

黑色素瘤蛋白质表达的数据集TCGA-SKCM-L4于2020年3月23日从TCPA下载，该数据集由354个样本组成。与之相关的临床特征数据从TCGA下载，并选择TCGA和TCGA-SKCM选项。该数据集由344个样本组成，包含各种临床特征的信息，包括生存状态(0=存活，1=死亡)、生存时间(d)、年龄、性别、分级和TNM分期。

1.2 生物信息学分析

本文使用R语言生存包中Kaplan-Meier(KM)和Cox比例风险分析统计方法，确定蛋白质表达与黑色素瘤预后的关系。简单地说，根据蛋白质中位表达水平将样本分为低表达组和高表达组。比较低表达组和高表达组的生存状态和生存时间，用KM统计方法计算P值，P<0.05被视为有统计学意义。采用Cox法计算P值，P<0.05，差异有统计学意义。使用Cox方法计算危险比(HR)和P值，并用于生成火山图，使用ggplot2和ggrepel软件包将数据绘制在图表上。HR值>1的蛋白质被认为与高风险相关，HR值<1的蛋白质被认为与低风险相关。

采用多因素Cox分析，以生存状态和生存时间为因变量，以蛋白质表达为自变量，生成基于蛋白质的预后模型。采用双向步骤对模型作进一步优化。利用11种蛋白质及其系数建立预测模型。随后，使用这种预后模型计算所有患者的风险评分和中位风险评分，并将患者分组：患者风险评分大于中位风险评分的为高危组，其余患者为低危组[16-17]。

使用survival和survminer软件包进行基于预后模型中的11种蛋白质的蛋白质生存分析。计算中位蛋白质表达量，将蛋白质水平高于中位的患者归为高表达组，其余患者分为低表达组；比较高、低表达组的生存状态和存活时间，使用每个患者的风险评分进行风险生存分析，并比较高危和低危患者的生存状态和生存时间。

风险分析是通过对患者的风险评分进行排序而实现[18]。按照排序，将每个患者的蛋白质表达情况和生存状态显示在热图中。生存包用于单因素Cox和多因素Cox模型的单因素和多因素独立预后分析，这些因素对生存状态和时间的影响用HR值和P值显示，并使用森林图进行可视化。

采用生存ROC软件包绘制ROC曲线。基于KM分析，使用危险度评分和其他因素来预测生存状态和时间，估计时间为1 a。随后，计算假阳性率和真阳性率，绘制ROC曲线，计算AUC值。蛋白质共表达分析采用相关性检验，以预测与模型中蛋白质共表达的其他蛋白质。相关系数和P值滤波器分别为0.4和0.001。对满足过滤条件的蛋白质进行保存并绘制相关图。使用ggplot2和ggalluvial绘制桑基图，以显示蛋白质与预后模型中蛋白质之间的相关性。

1.3 细胞培养和慢病毒感染

使用添加10%胎牛血清的RPMI-1640进行B16黑色素瘤细胞培养。GATA6和X1433ZETA过表达慢病毒购自南京Abmgood公司，并使用MOI为10的比例进行慢病毒感染。在荧光显微镜下监测GFP的表达，确定基因转移效率。

1.4 细胞迁移

细胞迁移实验如文献[19]所述。将感染GATA6、X1433ZETA或对照慢病毒的细胞接种到含有聚乙烯吡咯烷酮涂层的聚碳酸酯过滤器(8 μm孔径)的24孔小室中，并评估其迁移能力。小室下方的孔中加入添加0.1%BSA的RPMI-1640培养基，将细胞(无血清)接种到小室。培养24 h后，固定小室下方的细胞，用结晶紫染色。

1.5 凋亡测定

细胞凋亡检测方法如文献[20]所述。先重悬细胞，随后用5 μL Annexin-V和5 μL PI孵育15 min(避光)。流式细胞术分析各组凋亡早期细胞(Annexin V+/PI-)和晚期细胞(Annexin V+/PI+)的比例。

1.6 肿瘤生长分析

肿瘤生长分析如文献[21]所述。将1.5×106个GATA6或X1433ZETA过表达的B16黑色素瘤癌细胞注射到6～8周龄的雌性C57BL/6J小鼠皮下。一周后，每隔一天用游标卡尺测量肿瘤体积，计算公式为:体积=0.5×宽2×长。

2 结果与分析

2.1 鉴定黑色素瘤预后相关蛋白质

使用Kaplan-Meier (KM)和Cox(比例风险模型)生存包分析与黑色素瘤预后相关的蛋白质，样本分组为高表达或低表达(相对于中位表达)。比较低表达组和高表达组的生存时间和状态。计算P值和HR值，HR值>1表明与高危相关，HR值<1表明与低风险相关，用P值和HR值绘制火山图。如图1所示，与黑色素瘤预后相关的蛋白质有52个，其中蛋白质高表达和低表达的患者预后有显著差异(P<0.05为显著差异，有统计学意义)。

图1 黑色素瘤预后相关蛋白质的火山图Fig. 1 Volcano map of melanoma prognosis-related proteins

2.2 构建蛋白质预测模型

采用生存包中多因素Cox分析构建蛋白质预后模型，鉴定出11个蛋白质及其系数(表1)，其中GATA6、P27、S6、CD20的系数较小，X1433ZETA、P21、YAP的系数较大，表明GATA6与低风险相关，X1433ZETA与高风险相关。

表1 黑色素瘤预后相关蛋白质及其系数Tab. 1 Melanoma prognosis-related proteins and their coefficients

患者风险评分计算方法:风险评分=∑[表达(蛋白)×系数]。根据每个患者的风险评分以及中位风险评分，将患者分为高风险组或低风险组。通过比较蛋白质表达水平高的患者和表达水平低的患者的预后，该模型的有效性得到了验证。图2显示：LCK、P27、S6、SRC_pY416、CD20、CD49B和GATA6与低风险相关，因为这些蛋白质水平高的患者预后较好；CKIT、P21、YAP和X1433ZETA与高危相关，这些蛋白质表达水平高的患者预后较差。

图2 黑色素瘤预后与11种预后模型蛋白质表达的关系Fig. 2 Prognosis of melanoma in relation to protein expression in 11 prognostic models

2.3 利用蛋白质预后模型确定预后

首先，为检验模型的有效性，本文进行了风险分析，结果如图3所示：风险评分高的患者预后较差。低风险评分患者的5年生存率大于75%，而高风险评分患者的5年生存率低于50%(P=5.49×10-6)。这显示了该模型的可靠性。

图3 黑色素瘤预后与风险评分的关系Fig. 3 Prognosis of melanoma in relation to risk score

接下来，使用这个模型分析每个病人的预后。如图4所示，根据风险评分对患者进行排序，风险评分的增加与死亡的增加相关。此外，这11种蛋白在低危患者和高危患者中的表达存在显著差异。高危相关蛋白CKIT、P21、X1433ZETA在高危评分患者中高表达，而低危相关蛋白LCK、P27、S6、SRC_pY416、CD20、GATA6主要在低危评分患者中表达。

然后，我们使用单因素Cox和多因素Cox模型进行单因素和多因素独立预后分析，以确定预后相关因素。在单因素独立预后模型中，肿瘤分期、TNM分型、风险评分与预后相关(P<0.05)。而在多因素独立预后分析中，风险评分和TNM分型与黑色素瘤预后相关；风险评分的HR值为1.759(1.329～2.328)，大于其他因素。采用ROC曲线分析进行预测效率评价，风险评分的AUC得分最大，为0.761分(图4c)，说明风险评分相对于其他特征因素是更好的预测因子。

图4 单因素和多因素独立预后分析及与黑色素瘤预后相关因素的ROC曲线Fig. 4 Univariate and multifactorial independent prognostic analyses and ROC curves of factors associated with melanoma prognosis

因为这11种蛋白质不是分泌蛋白质，它们的水平只能在肿瘤组织中测量，因此，我们作进一步分析，以确定是否有分泌蛋白质与这11个蛋白质共表达。为了避免共表达蛋白质过多，相关系数和P值过滤分别设置为0.4和0.001。上述过滤条件鉴定出了与11个模型蛋白质共表达的65个蛋白质，与11个模型蛋白质高度相关的分泌蛋白质有4个，详见图5。

图5 通过相关分析鉴定了与预后模型共表达的4个分泌蛋白质Fig. 5 Four secreted proteins co-expressed with the prognostic model were identified by correlation analysis

桑基图(图6)显示了蛋白质共表达模式。其中GATA6与最多的蛋白质共表达，而X1433ZETA仅与PEA15共表达。

左侧为预后模型里的蛋白质，右侧是与模型蛋白质共表达的蛋白质。有最多共表达蛋白质的是GATA6，X1433ZETA仅与PEA15共表达图6 共表达蛋白质与预后模型蛋白质的相关性Fig. 6 Correlation of co-expressed proteins with prognostic model proteins

2.4 确定X1433ZETA蛋白质和GATA6蛋白质作用

为了验证GATA6和X1433ZETA蛋白质的作用，采用慢病毒感染方法在B16黑色素瘤细胞中过表达GATA6和X1433ZETA，然后检测细胞凋亡和细胞迁移。结果显示，过表达GATA6可显著诱导细胞凋亡，而X1433ZETA过表达对细胞凋亡无影响；过表达X1433ZETA可促进细胞迁移，但过表达GATA6可抑制细胞迁移(图7)。

图7 GATA6和X1433ZETA过表达对B16黑色素瘤细胞凋亡和迁移的影响Fig. 7 Effect of GATA6 and X1433ZETA overexpression on apoptosis and migration of B16 melanoma cells

我们使用X1433ZETA和GATA6过表达细胞建立了黑色素瘤皮下异种移植小鼠模型，发现X1433ZETA过表达可以促进肿瘤生长，GATA6过表达则抑制肿瘤生长。同时，与对照组相比较，X1433ZETA过表达缩短小鼠存活时间，而GATA6过表达则会延长小鼠的存活时间。结果详见图8。

图8 GATA6和X1433ZETA过表达对B16黑色素瘤小鼠肿瘤生长曲线及存活率的影响Fig. 8 Effect of GATA6 and X1433ZETA overexpression on tumor growth curve and survival rate in B16 melanoma mice

3 讨论

本文建立了一个基于11种蛋白质表达水平的预后模型，并用来计算每个患者的风险评分。分析表明，与低风险评分患者相比，高风险评分患者预后明显较差。多因素独立预后分析显示，风险评分与黑色素瘤预后相关。ROC曲线分析显示，风险评分AUC值为0.761，表明风险评分优于年龄、肿瘤分期、TNM分型等其他预后指标。通过相关性分析，本文鉴定了多个与模型蛋白质共表达的分泌蛋白质。体外和体内验证实验表明，与高风险相关的X1433ZETA过表达促进肿瘤生长，而与低风险相关的GATA6过表达抑制肿瘤生长。据估计，全球黑色素瘤发病率在过去20年中翻了一番，其发病率被认为高于任何其他类型的癌症(3.1%/年)[22-23]。鉴于其高死亡率，需要一种有效的黑色素瘤预后模型来改善预后。

本文预后模型包括11种蛋白质，其中LCK、P27、S6、SRC_pY416、CD20和GATA6在预后较好的患者中高表达，提示这些蛋白质与低风险相关。这些结果与之前的研究结果一致，例如，Deepak等[24]发现GATA6是一个星形细胞瘤肿瘤抑制基因。对于高风险蛋白质CKIT、P21、X1433ZETA，相关研究也提示这些蛋白质与较差的预后相关。例如CKIT在某些癌症的发生中起至关重要的作用，抑制CKIT激酶活性可以成为癌症治疗的靶点[25]。这些结果表明，本文模型与文献[24-25]的研究结果一致，模型是可靠的。

为了验证该模型，进行了单因素和多因素独立预后分析，结果显示风险评分较高的患者预后明显较低。此外，使用多因素独立预后分析发现风险评分与黑色素瘤预后相关。有研究报道了黑色素瘤的预后模型，例如，一些预后模型是根据患者临床特征建立的，包括肿瘤厚度、有无溃疡和原发肿瘤的位置[15]。然而这些模型没有考虑蛋白质表达，所以这些模式无法代表患者的分子生物学特征。同时，与疾病分级或TNM分期这些用于预后分析的模型相比，我们的风险评分AUC值高于TNM分期或分级(图4)，因而，本文所建立模型的预测能力优于现有的基于TNM分期或疾病分级的预后模型。

蛋白质预后模型是近些年开发的一种新的预后模型，其特点是可以用每个患者的某些蛋白质的表达量来预测患者的预后，而非基于患者的临床特征。因为基于数据库，我们可以获取每个患者上万种蛋白质的表达量，从中可以挑选若干跟预后显著相关的蛋白质来建立预后模型，与传统的基于TNM分期等模型相比，更好地代表了每个患者的肿瘤细胞分子生物学特征，在更深入的分子层次建立患者的预后模型，有更好的预测能力。本文研究中，比较了本文模型与TNM分期等在预后预测中的价值，发现蛋白质预后模型明显优于传统的基于临床特征的模型。正是基于上述优点，蛋白质预后模型已广泛用于肿瘤患者的预后预测[26]。因此，本文模型比过去的模型更有优势。相关性分析发现，本文模型中有多个蛋白质与这11个蛋白质共表达。应用严格的相关系数和P值过滤(0.4和0.001)鉴定出许多分泌蛋白质与11个模型蛋白质共表达。当无法获得肿瘤组织时，这些分泌蛋白质可能特别重要，因为它们可以作为11种模型蛋白质的替代品来预测患者的预后。

综上所述，我们得到了一种基于蛋白质的黑色素瘤预后模型，该模型比基于患者临床特征的传统预后模型表现更好，这有助于人们找到一个更优的个体化治疗方案。