可变剪切联合免疫分析构建膀胱癌预后模型的研究

王睿昊 成丁财 潘钰 王逢民

膀胱癌是泌尿生殖系统中最常见的恶性肿瘤之一，根据全国501个登记处上报国家癌症中心的1998—2015年间中国膀胱癌发病率统计结果，其发病率在我国恶性肿瘤中位居第十三[1]。膀胱癌通常分为非肌层浸润性膀胱癌(non muscle-invasive bladder cancer, NMIBC)和肌层浸润性膀胱癌(muscle-invasive bladder cancer, MIBC)两类。其中，MIBC临床表现为侵袭性病程，易发生局部复发和远处转移。由于其高复发率和预后差，现已成为威胁人类健康的主要疾病之一[2]。

可变剪切(alternative splicing, AS)是指在mRNA前体到成熟mRNA的过程中，不同的剪切方式使得同一个基因可以产生多个不同成熟mRNA，最终产生不同的蛋白质的过程。AS事件包括可变受体位点(alternate acceptor site, AA)、可变供体位点(alternate donor site, AD)、可变启动子(alternate promoter, AP)、可变终止子(alternate terminator, AT)、外显子跳跃(exon skipping, ES)、外显子互斥(mutually exclusive, ME)和保留内含子(intron retention, RI)7种亚型[3]。AS事件通过影响mRNA前体成熟过程导致了不同的mRNA剪接异构体的产生和蛋白质变体的翻译，这是转录组和蛋白质组多样性的关键调控机制之一[4]。目前有研究指出AS相关的遗传变异会导致膀胱癌风险增加[5]，但其对于临床应用的相关研究却不多。本研究旨在构建基于AS的膀胱癌预后模型，并探讨其与免疫治疗的关系。

对象与方法

一、数据的获取与处理

在癌症基因组图谱(The Cancer Genome Atlas, TCGA)数据库(https://portal.gdc.cancer.gov/)中提取405例膀胱癌患者基因测序数据及相关临床信息。并在TCGA Splice-Seq(https://bioinformatics.mdanderson.org/TCGASpliceSeq/)中获取AS拼接百分比(percent-spliced-in, PSI)，其是可以量化AS的指标，可以实现某个外显子是否纳入转录本定量，从而用于单个样本或组间多个样本的比较。

二、预后相关AS事件的筛选

在R软件(版本4.0.5，https://www.r-project.org/)中通过UpSetR、survival、ggplot2包进行单变量Cox回归分析以筛选预后相关AS事件，以PSI值>0.75 作为筛选标准，并排除PSI 值标准偏差<0.01 的AS 事件。绘制AS UpSet图、火山图、气泡图。AS 事件多数是通过剪切类型、Splice-Seq中的 ID 号和对应的基因名称来注释的。例如，在“TRMU|62711|AA”中，TRMU表示对应的基因名称，62711表示剪接变异的ID，AA表示剪接类型。

三、预后模型的构建及验证

根据RS将患者分为低风险组和高风险组，通过survival包进行生存分析，并通过多元Cox回归分析验证RS为膀胱癌预后独立危险因素。借助survminer、timeROC、ggpubr包绘制受试者工作特征(receiver operating characteristic curve, ROC)曲线，验证模型的预后价值及准确性，并分析RS与临床变量的相关性。为进一步使模型可视化，基于rms包建立包含RS和其他临床特征的预后列线图，并绘制校准曲线。

四、RS与肿瘤免疫微环境相关性分析

为了进一步检验RS是否可以作为免疫指标，基于estimate、limma、ggpubr包通过CIBERSORT、ssGSEA算法分析RS与肿瘤微环境的相关性。并基于ggExtra、corrplot包分析RS与免疫检查点的相关性。

五、RS相关靶基因的筛选与验证

基于构成RS的AS相关基因，通过limma包进行差异分析，并对筛选出的靶基因进行生存分析。通过人蛋白质图谱在线数据库(https://www.proteinatlas.org/)进一步验证。

结 果

一、纳入对象的基本情况

病例总数为405例，其中男299例(73.83%)，女106例(26.17%)。患者基本特征见表1。

表1 膀胱癌患者基本临床信息[例(%)]

二、预后相关AS事件

经过单变量Cox回归分析共筛选出2 828例与生存显著相关的AS事件(P<0.05)，其中ES 831例(29.38%)，AP 699例(24.72%)，AT 540例(19.10%)，RI 367例(12.98%)，AA 228例(24.72%)，AD 156例(5.52%)，ME 7例(2.48%)。由UpSet图(图1)可知，生存相关的AS事件中ES为主要亚型，且AA与ES的交集基因最多。

注：横坐标代表AS事件相关基因数目，纵坐标代表交集基因数目

三、预后模型的构建与验证

利用Lasso回归(图2A、2B)及单因素、多因素Cox回归分析构建了一个包含8个AS事件的预后模型(表2)，其中ANK3|11845|AP、AK9|77203|AT、PCSK5|86634|AT、PTGER3|3415|AT为危险性因素。多因素Cox回归分析也验证了RS可以作为预测膀胱癌患者预后的独立危险因素(P<0.05)(图2C、2D)。ROC曲线显示，相较于年龄、肿瘤分期等，RS具有更好的预后预测能力，且RS 1、2和3年生存率的ROC曲线下面积(area under curve, AUC)均>0.7，表明其对生存预测具有较高的敏感性和特异性(图2E、2F)。根据RS值将患者分为高、低风险组，高风险组患者的总生存期(overall survival, OS)显著缩短(P<0.05)。

A：用于Lasso回归中调整参数选择的10次交叉验证；B：所有AS事件的Lasso系数分布；C：Cox单因素回归森林图；D：Cox多因素回归森林图；E：RS用于1、2、3年生存率预测的ROC曲线；F：RS及不同临床因素的ROC曲线

表2 模型构成剪切事件

为进一步探讨其临床应用价值，我们完善了RS与临床变量的相关性分析，由图3A、3B可知，随着肿瘤分期增加，RS显著上升(P<0.05)，表明RS与肿瘤进展呈正相关，但是对于高分期肿瘤其意义并不大。我们建立了包括RS、年龄、分期在内的列线图来可视化预后模型(图3C)，排除性别和分级是因为其AUC<0.6，且在多因素Cox分析中无意义。校准曲线近似对角线，表明列线图具有较好的 1、2 和 3 年生存率预测能力(图3D～3F)。

A：RS与肿瘤分期的相关性分析；B：RS与肿瘤T分期的相关性分析；C：整合了RS、年龄、分期的列线图；D～F：校准曲线(红线与对角线越接近则预测与实际越相符合)

四、RS与肿瘤免疫微环境及免疫检查点相关基因相关性分析

RS与肿瘤微环境相关性分析结果显示，高风险组伴随着更高的免疫细胞评分及基质细胞评分(P<0.05)(图4A～4D)，提示肿瘤的发展可能与免疫细胞浸润相关。进一步的免疫细胞浸润分析显示RS与以下免疫细胞显著相关：未刺激活化的巨噬细胞(Macrophages M0)(r=0.33，P=1.9e-06)、替代活化型巨噬细胞(Macrophages M2)(r=0.24，P=0.000 69)、浆细胞(Plasma cells)(r=-0.29，P=4.7e-05)、CD8+T细胞(T cells CD8)(r=-0.22，P=0.002 2)、调节性T细胞(T cells regulatory, Tregs)(r=-0.32，P=5e-06)(图4E、5A～5E)。

免疫检查点相关基因相关性分析提示RS与表面抗原分化簇274(cluster of differentiation 274, CD274)(r=0.26，P=1.2e-07)、甲型肝炎病毒细胞受体2(hepatitis A virus cellular receptor 2, HAVCR2)(r=0.31，P=3.5e-10)、细胞程序性死亡蛋白1配体2(recombinant programmed cell death protein 1 ligand 2, PDCD1LG2)(r=0.42，P=2.2e-16)呈正相关(图4F～4G、5F～5H)。

A：RS与免疫细胞评分相关性分析小提琴图；B：RS与基质细胞评分相关性分析小提琴图；C：RS与综合评分相关性分析小提琴图，综合评分为免疫细胞评分与基质细胞评分之和；D：RS与肿瘤纯度相关性分析小提琴图；E：免疫细胞浸润分析(横坐标代表免疫细胞名称，纵坐标代表免疫细胞含量)；F：免疫检查点相关基因相关性分析(红色代表正相关，颜色越深即相关性越高)；G：免疫检查点相关基因差异分析(纵坐标代表基因表达量，红色为高风险组)(***P<0.001；**P<0.01；*P<0.05)

A～E：RS与免疫细胞相关性分析；F～H：RS与免疫检查点相关基因相关性分析(黄色和蓝色的峰代表细胞的集中分布)

五、靶基因筛选与验证

将构成预后模型的8个基因进行差异分析，有意义的为TRMU和PTGER3(P<0.05)。差异分析提示TRMU在肿瘤组织中高表达(图6A)，PTGER3在肿瘤组织中低表达(图6B)，生存分析也验证了PTGER3对于膀胱癌而言为保护性因素(P<0.05)(图6C、6D)。在人蛋白质图谱在线数据库(http://www.proteinatlas.org)中检索正常组织和膀胱肿瘤组织中TRMU及PTGER3免疫组织化学染色图像，进一步验证了上述结论。

A：TRMU差异分析；B：PTGER3差异分析；C：TRMU生存分析；D：PTGER3生存分析

讨 论

在过去的30年里，膀胱癌的治疗进展甚微。大约25%的膀胱癌是治疗难度较大的MIBC[6]。其中85%的患者在未接受治疗的情况下在2年内死亡[6]。目前临床上需要能准确预测膀胱癌患者预后的生物标志物，以利于临床决策和个性化预后监测。近年来免疫疗法在肌层浸润性和转移性膀胱癌的治疗中显示出一定的临床应用潜力[7]。

AS不仅对造血、肌肉收缩等机体的正常生理过程有着重要影响，而且在多种肿瘤的发生、发展和肿瘤微环境的形成中也起着至关重要的作用[8-9]。AS事件的预后价值已在多种肿瘤中得到证实，包括肾癌、肝癌、头颈部肿瘤等[10-12]。对AS事件的综合分析有望为膀胱癌的预后分析提供潜在生物标志物。

在目前的研究中，我们基于2 828例与预后显著相关的AS事件构建了一个包含8个AS事件的预后模型，并通过生存分析、ROC曲线及Cox多因素回归证实了该模型的预后准确性。其中，本模型1、2和3年生存率的AUC均>0.7，表明其对生存预测具有较高的敏感性和特异性。此外，相较于年龄、肿瘤分期等临床特征，本模型有着更好的预后预测能力。我们进一步通过绘制包含RS和年龄、分期等临床特征的列线图来可视化该模型，从而提高其实用性。

为了揭示AS事件在免疫浸润中的作用，我们进行了一系列的免疫相关分析。结果表明，高风险组免疫细胞浸润程度高，RS与抑制性免疫细胞(巨噬细胞 M0/M2)之间存在显著正相关，与抗肿瘤免疫细胞(CD8+T细胞)呈负相关。肿瘤相关巨噬细胞(tumor-associated macrophage, TAM)是有助于免疫逃逸的抑制性免疫细胞之一，可分化为M1(抗肿瘤)和M2(促肿瘤)表型[13]。由单核细胞形成的M0巨噬细胞尚未极化为肿瘤内的M1或M2亚型，虽然M1巨噬细胞参与了炎性活动并发挥了抗肿瘤活性，但是TAM倾向于获得M2致癌表型。M2巨噬细胞常通过基质重塑、新血管生成和抑制局部免疫支持肿瘤生长和转移[14]。有研究指出，M0巨噬细胞浸润增加与OS降低和肿瘤分期增加有关[15]。因此，M0巨噬细胞可能在膀胱癌的发展中发挥重要作用，可能具有治疗靶向性。此外，TAM可以表达细胞因子和酶，这些细胞因子和酶可以抑制T细胞的募集和激活，从而增强对免疫治疗的抵抗性[16]。本研究还显示，RS与CD274、HAVCR2、PDCD1LG2等免疫检查点相关基因呈正相关。CD274也被称为细胞程序性死亡-配体1(programmed cell death-1, PD-1)，最近有研究发现，TAM可以导致PD-1的上调以实现免疫抵抗[17]。因而TAM在免疫逃逸及免疫抵抗中所发挥的重要作用可能是RS成为膀胱癌预后标志物的内在因素。

TRMU是一种对线粒体DNA翻译至关重要的核基因，通过编码tRNA 5-甲基氨基甲基-2-硫尿苷酸甲基转移酶，使线粒体tRNA硫醇化[18]。PTGER3是一种G蛋白偶联受体，其对肿瘤相关血管生成和肿瘤生长至关重要[19]。本研究显示TRMU及PTGER3在正常膀胱组织和膀胱肿瘤组织中差异表达，且预后分析提示TRMU高表达患者的预后较差，PTGER3高表达患者的预后较好，但与之相关的研究甚少。TRMU及PTGER3可能成为新的研究靶点。

综上所述，本研究基于AS构建了一个具有可行性和准确性的膀胱癌预后模型，在探寻模型背后免疫机制的同时挖掘了潜在的靶基因。但是受限于客观条件，本研究的发现仍需进一步的实验和临床验证。