杨菁菁, 郭怀娟, 茅静贤, 王佳欣, 王 颖, 严雪冰, 潘晓萍

(扬州大学附属医院, 1. 肿瘤科, 2. 健康管理中心, 江苏 扬州, 225000)

结直肠癌(CRC)是人类常见的消化道恶性肿瘤,在所有肿瘤中发病率排名第3位,死亡率排名第2位[1]。尽管近年来风险筛查、外科技术和靶向治疗方面取得了巨大进步,但仍有相当一部分CRC患者预后不佳。CRC的发生发展是一个多步骤、多基因的复杂生物学过程,而如何选择合适的分子标志物帮助判断CRC预后尤为重要。近年来研究[2]表明, N6-甲基腺苷(m6A)修饰不仅参与CRC发生发展,而且在调控抗癌药物敏感性过程中也发挥着重要作用。例如,甲基转移酶样蛋白14(METTL14)作为催化m6A修饰主要甲基化酶,不仅参与调控CRC细胞的生长和转移,而且其表达水平可能影响免疫及靶向治疗的敏感性[3]。在临床实践中, CRC组织特异性免疫细胞浸润也是一个至关重要的预后因素[4]。例如, CD3+T细胞和细胞毒性CD8+T细胞浸润与CRC患者的术后复发及总体生存密切相关[5]。尽管m6A修饰及免疫细胞浸润均与CRC预后相关,但是两者相关分子能否共同用于预后判断尚不清楚。本研究拟利用生物信息学方法筛选m6A修饰及免疫细胞浸润相关分子,构建预后模型并进行验证,为CRC患者的精准预后评估提供新的策略,现报告如下。

1 资料与方法

1.1 CRC患者相关数据来源

从癌症基因组图谱(TCGA, https://portal.gdc.cancer.gov/)下载612例结直肠癌组织的mRNA表达谱分析结果和临床数据。从基因表达综合数据库(GEO, https://www.ncbi.nlm.nih.gov/geo/)中下载归一化基因表达数据和临床数据,数据集登录号为GSE14333(n=226)、GSE17536(n=177)、GSE39582(n=585)和GSE38832(n=122)。

1.2 研究流程

首先,利用m6A调控因子的mRNA表达谱,基于共识聚类方法对TCGA数据库中612例CRC患者的预后进行分层,在同一队列中,基于单样本基因组富集分析(ssGSEA)检测24种免疫细胞的丰度,并对患者预后进行分层; 其次,分析队列中m6A调节因子与免疫细胞的相关性,选择612例CRC组织和51例正常组织之间的差异表达基因,基于加权基因共表达网络分析(WGCNA)筛选与m6A修饰及免疫细胞浸润相关的基因; 最后,基于Lasso回归分析建立预后模型,一方面利用TCGA数据库分析预后模型在不同亚组中的预测功效,另一方面选取4个GEO数据库队列作为外部验证队列,进一步明确预后模型的临床价值。见图1。

1.3 CRC组织免疫细胞浸润分析

利用ssGSEA计算24种免疫细胞的丰度并进一步聚类,解析样本中不同类别。根据上一步得到的样本类别,利用生存分析解析不同类别在预后方面是否具有显著差异。根据不同免疫细胞丰度,结合样本的生存情况,利用单因素Cox回归分析解析免疫细胞类型对CRC预后的影响。

1.4 WGCNA筛选m6A修饰及免疫细胞浸润相关预后基因

1.5 统计学分析

采用t检验对各组间m6A调控因子表达水平进行比较; 采用Kaplan-Meier模型构建生存曲线,并使用Log-rank检验比较组间生存差异; 采用单变量Cox回归分析评估免疫细胞比例与患者预后的相关性; 采用Lasso回归分析筛选出与预后相关的基因; 根据基因表达水平和Lasso回归系数建立预后模型; 采用受试者工作特征(ROC)曲线对所建预后模型的性能进行评估;P<0.05为差异有统计学意义。

2 结 果

2.1 m6A分子表型及免疫细胞浸润在CRC患者预后中的意义

首先,对612个CRC患者组织样本进行聚类并确定最佳聚类K值为4(图2A)。基于总生存期(OS)将612例CRC患者分为亚组1和亚组2(P=0.015, 图2B), 发现14个m6A调节因子[如RNA结合序列蛋白(RBM15)、RBM15B和KIAA1429]在2个亚组间的表达差异显著(图2C)。同时,基于ssGSEA方法计算上述患者肿瘤组织中24个免疫细胞的丰度,进而依据免疫微环境(TME)将患者划分为TME1组(n=180)和TME2组(n=432)(图3A); 生存分析表明, TME1组患者的OS率低于TME2组,差异有统计学意义(P=0.049, 图3B); 单变量Cox回归分析表明,免疫细胞浸润比例与OS无相关性(图3C); 进一步研究发现10种m6A调节因子(如RBM15、RBM15B和WTAP)在TME1组与TME2组间的表达存在显著差异(图3D); 相关性分析表明,大部分m6A调节因子与部分免疫细胞的浸润存在显著相关性(图3E), 例如CRC组织中Alk B同源蛋白5(ALKBH5)表达水平与所有免疫细胞浸润比例均呈正相关。

2.2 筛选m6A修饰和免疫细胞浸润相关预后基因

从TCGA数据库中筛选出612个CRC样本和51个正常组织样本的差异表达基因,纳入标准为: ① |logFC|>1; ②P<0.01。与CRC组织样本相比,对照样本中分别有1 209个基因和922个基因显著上调和下调(图4A)。基于2 131个基因构建了加权基因共表达网络,最佳β值定为6。对基因进行聚类,每个模块的基因数不少于30个。经过动态树切割后,对紧密相邻的模块进行聚类,最终得到9个模块。基于本研究发现m6A分子表型与OS的相关性高于免疫细胞浸润,因此选择与m6A调节因子相关的模块进行下一步研究。根据纳入标准(MM>0.4; GS>0.2)筛选出25个核心基因(图4B)。同时,相关性研究发现部分核心基因表达水平(如ATOH1)与免疫细胞浸润相关(图4C、4D)。

2.3 基于m6A修饰及免疫细胞浸润相关基因构建预后模型

基于Lasso回归分析从25个核心基因中筛选出4个预后基因: 内凝集蛋白1(ITLN1)、淋巴细胞抗原6复合位点G6D(LY6G6D)、无调性同族体1(ATOH1)和基质金属蛋白酶(MMP28)。基于4个预后基因在GEPIA数据库中的表达差异分析,基因LY6G6D与MMP28在结肠癌与直肠癌组织中的表达较癌旁组织有差异性(图5A、5B)。基于4个基因的Lasso分析构建了预后模型: 风险评分=EXPITLN1×(-0.050 23)+EXPLY6G6D×(-0.022 89)+EXPATOH1×(-0.141 97)+EXPMMP28×(-0.045 74)。使用公式计算训练队列每个样本的风险评分,并根据中位值将样本分为高风险组和低风险组。每个样本风险评分分布及对应患者生存信息见图5C、5D。生存分析表明,高风险组的OS低于低风险组(图5E)。ROC曲线分析表明,预后模型曲线下面积(AUC)为0.634, 优于年龄和TNM分期等其他临床因素(图5F)。

2.4 预后模型在训练队列不同亚组中的预测效能

为了进一步探讨预后模型的临床价值,将TCGA数据库队列根据年龄、性别及肿瘤位置分为不同亚组进行分析。结果表明,在>60岁患者中,高风险组的OS显著低于低风险组的OS(P=0.003); 在≤60岁患者中,高风险组与低风险组OS无显著差异(P=0.100)。在女性患者及男性患者中,高风险组OS率均显著低于低风险组(P女性=0.014,P男性=0.047)。结肠癌高风险组OS率显著低于低风险组,但直肠癌高风险组与低风险组OS率无显著差异(P结肠癌=0.003,P直肠癌=0.210)。见图6。

2.5 预后模型在验证队列中的预测效能

从GEO数据库中下载了4个CRC队列作为验证队列以进一步明确预后模型的临床价值。在GSE17536、GSE39582和GSE38832队列中,高风险组OS显著低于低风险组(PGSE17536=0.013,PGSE39582<0.001,PGSE38832=0.045)。在GSE14333队列中,高风险组的无病生存率显著低于低风险组(P=0.016)。见图7。

3 讨 论

研究[6-7]表明,分子标志物能够精准指导CRC患者的早期诊断、预后评估及治疗药物选择。m6A修饰是CRC发生发展中的重要分子事件之一,其调控因子(包括书写蛋白、阅读蛋白及识别蛋白)在CRC组织中表达失调,且可作为患者诊断及预后的潜在标志物[8]。本课题组团队前期基于m6A修饰及上皮细胞间质转化(EMT)筛选预后分子并构建相关临床模型,该模型可在多个CRC队列中识别预后较差的高风险人群[9]。本研究基于m6A修饰及免疫细胞浸润筛选与CRC预后相关的分子标志物(ITLN1、LY6G6D、ATOH1和MMP28), 并构建相应的预后模型,该模型不仅在不同临床特征的亚组中能识别预后不良的高风险人群,而且在多个外部验证队列中表现出较好的临床结局预测能力。

本研究首先从TCGA数据库中提取612例CRC患者临床资料及匹配的肿瘤组织转录组信息,利用聚类分析将患者分为不同预后差异的2个亚组。进一步分析表明,包括RBM15、ALKBH5和脂肪量与肥胖相关蛋白(FTO)在内的14个调控因子在高风险组中较低风险组显著上调或下调,提示上述调控因子可能参与CRC发生发展。在机制研究中,敲除RBM15可影响下游髓样分化因子88的m6A甲基化修饰水平,进而抑制CRC细胞的增殖及转移能力[10]。低表达ALKBH5与CRC患者不良预后相关,其可能通过m6A修饰途径抑制PHD锌指蛋白20的mRNA稳定性而发挥抑癌效应[11]。FTO通过介导锌指蛋白687的m6A修饰进而激活Wnt信号通路,最终促进CRC的增殖、转移及血管生成[12]。在临床研究中, m6A修饰识别蛋白如YTH结构域家族蛋白2和3高表达水平提示CRC患者不良预后,而胰岛素样生长因子2 mRNA结合蛋白1高表达水平则提示相反结局[13]。值得注意的是, m6A修饰书写蛋白甲基转移酶样蛋白3(METTL3)及METTL14的表达水平在高风险组与低风险组间无显著差异,可能与既往部分研究[14-15]结果相悖,因而需更多的回顾性研究来确证。

本研究同时发现多个免疫细胞浸润比例能够共同显著区分TCGA数据库中CRC患者的OS, 提示其在CRC预后评估中存在潜在应用价值。在既往文献中, CD3+和CD8+T细胞浸润比例与CRC患者的良好预后呈正相关,其机制可能与Toll样受体介导的抗肿瘤免疫系统激活相关[16]。在微卫星稳定型伴肺转移的CRC患者中,基于上述2种T细胞的免疫评分可以预测患者术后5年的OS状态[17]。自然杀伤细胞浸润虽然与CRC患者OS无显著相关性,但是与无疾病间隔时间及无复发生存时间呈正相关[18]。值得注意的是,单因素Cox回归分析表明24种免疫细胞与预后均无显著相关性,这与既往部分研究结果相似。例如,一项生物信息学研究[19]表明,在GEO数据库CRC队列中大多数免疫细胞浸润比例与CRC患者的OS无显著相关性。因此,一方面,免疫细胞浸润在CRC中的预后价值仍需更多的回顾性研究来确证; 另一方面,笔者研究结果提示多种免疫细胞类型组合比单一免疫细胞类型更有利于患者预后的精准评估。最后,相关性研究发现m6A调控因子与多数免疫细胞浸润比例相关,提示这些因子在抗肿瘤免疫中可能扮演着重要角色。例如,最新研究[20]发现ALKBH5可抑制细胞核因子κB, 进而减少趋化因子配体5表达水平并促进CD8+T细胞浸润,最终抑制CRC细胞转移。METTL3可以促进髓源性抑制细胞的积聚,进而抑制CD4+和CD8+T细胞的激活和增殖,最终介导CRC细胞对免疫检查点抑制剂的抵抗[21]。

为了进一步探究m6A修饰和免疫细胞浸润在CRC中的预后价值,笔者团队基于WGCNA筛选出与两者相关的4个预后基因(ITLN1、LY6G6D、ATOH1和MMP28)用于构建模型。在既往文献中, ITLN1高表达提示CRC患者更有利的预后,相关机制研究[22]表明其可能通过抑制表皮生长因子受体、细胞外调节蛋白激酶及蛋白激酶B发挥促癌效应。LY6G6D在CRC组织中选择性表达并且可以作为T细胞免疫治疗的潜在靶点[23]。生物信息学研究[24]表明ATOH1低表达为影响CRC患者总体生存的独立不良预后因素。MMP28作为CD36的下游靶基因可通过介导E钙黏蛋白的裂解促进CRC细胞的转移[25]。在进一步验证工作中,一方面基于上述4个预后基因的临床模型可以在不同临床特征亚组中识别潜在预后不良人群,另一方面该临床模型在4个GEO验证队列中仍然可有效区分患者的临床结局。因此,该临床模型有望成为帮助临床医师精准评估患者预后的潜在辅助工具。

本研究存在如下局限性有待改进: 首先,本研究所用验证队列均来源于GEO数据库,缺少本单位及合作单位独立回顾性队列进一步验证。其次, 4个预后基因如何通过调控m6A修饰及抗肿瘤免疫反应影响CRC发生发展仍不清楚,有待于后续基于生命组学的机制研究来阐明。再次,在<60岁及直肠癌亚组分析中,尽管低分险组OS有优于高风险组OS的趋势,但是统计学上无显著差异。因而,预后模型在上述人群中的临床价值仍需更多的回顾性研究来确证。最后,辅助化疗是CRC患者术后的重要治疗手段之一,本预后模型能否识别辅助化疗不敏感的风险人群也是未来值得研究的方向。

综上所述,本研究基于m6A修饰及免疫细胞浸润构建了包含4个预后基因的风险模型,该模型能够在训练队列及验证队列中有效识别预后不良的潜在风险人群。上述研究成果不仅推动相关基础研究领域成果的临床转化,也为基于分子标志物的CRC精准诊疗提供了新的思路。