基于加权网络的卵巢癌可变剪接预后特征分析

顾云婧,潘以红,朱东月,朱 平

(江南大学理学院,中国江苏无锡214122)

卵巢癌是女性常见的癌症之一,居于女性癌症死因的第四位,2012年全球共诊断出新增卵巢癌患者约238 700例[1]。研究显示,由于早期症状不明显,约70%的卵巢癌患者被诊断时已是晚期,晚期患者的5年生存率低于30%[2]。目前,卵巢癌的治疗手段主要是手术切除和化疗,但预后情况并不理想[3],大多数晚期患者会在18个月内复发[4]。因此,迫切需要更深刻地了解卵巢癌的预后分子机制,寻找新型治疗靶标,以开发更有效的治疗方法来改善卵巢癌的预后效果。

已有研究显示,基因表达失调对肿瘤的发生和发展有重要影响,通过分析基因表达模式可以发现具有预后价值的癌症生物标志物[5]。目前,这些研究主要集中在对编码基因表达水平的分析,而对mRNA可变剪接(alternative splicing,AS)的系统研究还较为缺乏。可变剪接是单个前体RNA分子产生多种成熟mRNA的机制,mRNA的可变剪接是高等真核生物扩大蛋白质组多样性的重要基因调节机制之一[6～7]。因而,可变剪接模式的变化可能导致蛋白质功能的变化并影响基于这些功能的重要活动。近年来,大量的基因组学和功能性研究表明可变剪接与癌症的发生发展密不可分,可变剪接事件参与了一系列的致瘤过程,包括侵袭、增殖、血管生成、免疫和转移[8]。据报道,含有外显子v8-10的CD44变体与卵巢癌的转移和预后有关[9]；剪接变体AIMP2-DX2可能是治疗卵巢癌化疗耐药性的有效辅助靶点[10]。癌症研究的新趋势表明,可变剪接在癌症治疗方面具有临床潜力[11～14],可以作为癌症诊断和治疗的生物标志物,并作为新型治疗药物的潜在靶标。

本文根据卵巢癌患者的全基因组可变剪接数据及临床数据,筛选出与总生存率(overall survival,OS)相关的AS事件,并融合其对应基因的蛋白质相互作用(protein-protein interaction,PPI)数据构建卵巢癌预后加权网络；对预后加权网络进行拓扑分析,筛选有意义的AS事件作为卵巢癌预后特征；利用预后特征建立预后预测模型。结果显示模型可用于预测卵巢癌患者的临床结果,显示了可变剪接在卵巢癌中的预后价值。

1 材料与方法

1.1 数据来源

卵巢癌患者的可变剪接数据下载自TCGASpliceSeq数据库(http://projects.insilico.us.com/TCGASpliceSeq/)[15]。组数据库(The Cancer Genome Atlas,TCGA)中mRNA可变剪接模式的公开数据库资源。本文从中选取了393个具有临床数据的卵巢癌患者样本,生存时间均大于30 d,且样本的可变剪接表达估计值(percent spliced in,PSI)百分比大于75%。同时,从TCGA数据库中额外下载了16个卵巢癌患者样本的RNA-seq数据,并利用SUPPA软件[16]计算样本中可变剪接事件的PSI值。

1.2 预后相关AS事件的筛选

对于临床数据,由于每个病例随访时间不同,发生终点事件的可能性也不相等,因此本文选用Cox单因素回归分析筛选预后相关的AS事件(prognostic-associated alternative splicing events,PASEs)。

Cox单因素回归模型[17]形式:

其中,h(t|X)是具有协变量X的样本在t时的风险函数,t为生存时间；h0(t)是协变量为0时的风险函数；β为模型中待估计的回归系数。

回归系数β＞0时,协变量的取值越大,风险函数h(t|X)的值越大,表示病人死亡的风险越大；回归系数β=0时,表示协变量对风险函数h(t|X)没有影响；回归系数β＜0时,协变量的取值越大,风险函数h(t|X)的值越小,表示病人死亡的风险越小。

对于每个可变剪接事件,将患者根据该事件PSI值的中位数进行分类,分为高、低表达两组进行单因素回归分析,计算危险比(hazard radio,HR)及95%置信区间(confidence interval,CI),以P＜0.05为标准筛选出与卵巢癌预后相关的AS事件。

1.3 卵巢癌预后加权网络的构建

考虑PASEs之间的调控关系,计算预后相关AS事件之间的相关性,并融合其对应基因的PPI数据建立加权网络。对于每一个PASE,计算它与其余所有PASEs之间的皮尔逊相关性系数(Pearson correlation coefficient,PCC)。选取PCC值排名前5%,且相应P值小于0.01的AS事件对构建网络。相关性矩阵MAS由PASEs中符合条件的PCC构成,矩阵元素M(i,j)代表第i和j个事件之间PCC的绝对值。按如下方式对相关性矩阵进行标准化,标准化后的矩阵为

其中EAS(i,i)为矩阵MAS第i行元素的和。

除了相关性关系外,本文在网络构建中融入蛋白质相互作用数据。将预后相关的AS事件对应基因(genes of prognostic-associated alternative splicing events,PASEGs)提交至STRING在线数据库(https://string-db.org)[18],以得分大于0.4为标准,选取PASEGs之间的PPI数据。随后,根据PASEGs之间PPI的置信度构造矩阵MPPI,对其进行标准化:

其中EPPI(i,i)为矩阵MPPI第i行元素的和。

将AS事件对的相关性关系和其对应基因的PPI相结合,构成网络权重:

1.4 预后加权网络的拓扑分析及PASEGs的功能富集分析

对预后加权网络进行网络拓扑特征分析,以挖掘拓扑重要性强的网络节点。网络特征分析包含网络节点的度(degree)、特征路径长度(characteristic path length,CPL)、聚类系数(clustering coefficient,CC)和小世界指数(small-world index,SW)。节点的度为该节点与其他节点直接连接的边数。网络的CPL为网络中所有节点对之间的最短路径长度的平均值。节点的CC为与该节点直接连接的节点中实际存在的边与所有可能存在的边的比例,网络的CC为网络中所有单个节点的CC平均值。通过与相同规模的ER(Erdos-Rényi)随机网络比较可定义小世界网络[19],即该网络需同时满足以下条件:

其中,CCrandom和CPLrandom分别表示1 000个与目标网络相同规模的ER随机网络的CC、CPL平均值。

进一步地,网络小世界指数[20]被定义为:

当网络的小世界指数бSW大于1时,该网络为小世界网络。这意味着,该网络的节点高度聚集,网络的平均最短路径长度几乎与随机网络一样低,该生物网络具有高效的信息传递能力。

通过DAVID在线工具[21]对PASEGs进行GO(gene ontology)功能富集分析。GO分析包括分子功能(molecular function,MF)、生物过程(biological process,BP)和细胞组分(cellular component,CC)3个部分。同时,运用KOBAS数据库(http://kobas.cbi.pku.edu.cn/)[22]对PASEGs进行KEGG(kyoto encyclopedia of genes and genomes)通路富集分析,P＜0.05的通路被认为是显著富集的通路。

1.5 预后特征的筛选及预测模型的建立

卵巢癌预后加权网络的拓扑性质显示了AS事件在卵巢癌中的预后价值。为了预测卵巢癌患者的临床结果,本文根据网络节点的度筛选预后特征,并构建预后预测模型。首先选择预后加权网络中节点度排名前20的AS事件作为预后特征,随后通过Cox比例风险模型(proportional hazards model)构建卵巢癌预后预测模型,预后指数(prognosis index,PI)的计算公式为:

其中,x和β表示预后特征及其回归系数,n表示预后特征的个数。

根据预后指数PI的中值将患者分为高、低风险两组,运用Kaplan-Meier图表估计生存函数,分析两组之间的预后差异。此外,ROC(receiver operating characteristic)曲线[23]和 AUC(area under curve)值被用来评估预后预测模型的效率,AUC值越大,预测模型的正确率越高。

2 结果及分析

2.1 预后相关的可变剪接事件

可变剪接事件一般分为7种[24],分别为外显子跳跃(exon skip,ES)、互斥可变外显子(mutually exclusive exon,ME)、内含子保留(retained intron,RI)、可变 3′剪接位点(alternative acceptor site,AA)、可变 5′剪接位点(alternative donor site,AD)、可变初始外显子(alternative promoter,AP)和可变终末外显子(alternative terminator,AT)。

利用Cox回归对393个样本的所有剪接事件分别进行单因素生存分析,共筛选出290个与卵巢癌OS显著相关的AS事件。其中,AA事件16个,AD事件26个,AP事件69个,AT事件33个,ES事件116个,ME事件1个,RI事件29个。P值排名前10的AS事件见表1。

2.2 预后加权网络的拓扑分析

从STRING数据库中获取了118个PASEGs的158对互作关系,与相应AS事件对的PCC结合后构建了卵巢癌预后加权网络。该网络共包含281个节点,2 284条边。图1是利用Cytoscape软件[25]绘制的权重大于0.1的部分网络图。网络拓扑分析表明,预后加权网络的聚类系数为0.255,特征路径长度为2.498。1 000个相同节点数和边数的ER随机网络的平均聚类系数为0.058,平均特征路径长度为2.309。与ER随机网络相比,预后加权网络满足公式(5)和公式(6),因此预后加权网络为小世界网络,由公式(7)可得其小世界指数为4.064。

表1 排名前10的预后相关AS事件Table 1 Top 10 prognostic-associated alternative splicing events

图1 卵巢癌预后加权网络的代表性结果节点的大小代表该AS事件度的大小,边的粗细代表两AS事件间的权重大小。Fig.1 Representative results of prognostic weighted network for ovarian cancerThe size of the node represents the degree of the AS event,and the thickness of the edge represents the weight between the two AS events.

2.3 PASEGs的功能富集分析

通过DAVID在线工具对PASEGs进行GO富集分析,结果表明,PASEGs主要富集于细胞核、核质、RNA聚合酶II启动子转录调控、DNA转录调节、细胞凋亡、蛋白质结合、核酸结合等29个生物学过程(图2A)。显然,这些生物学过程与细胞中遗传物质的转录、蛋白质的合成等具有密切的相关性,可对细胞的活动起到一定的调控作用,从而影响卵巢癌的发生发展。

同时,运用KOBAS数据库对PASEGs进行KEGG通路富集分析,结果显示,PASEGs主要富集在肠道免疫网络、p53信号通路、凋亡途径、核糖体相关途径等24条通路上。这与GO富集分析结果一致,表明细胞中遗传物质的转录、蛋白质的合成及细胞凋亡等功能的异常在卵巢癌发生发展的过程中具有重要作用,可能是影响卵巢癌预后的重要途径。排名前10的显著富集通路见图2B。

2.4 可变剪接预后特征及预后模型

根据卵巢癌预后加权网络中节点的度对网络节点进行排序,选择排名前20的AS事件作为可变剪接的预后特征(表2),并运用Cox回归模型建立预后预测模型。预后模型的预后指数由20个预后特征的表达值和Cox回归分析得出的回归系数线性组合构建。根据公式(8)计算每个患者的预后指数,并将患者根据预后指数的中位数分为高风险组和低风险组。图3A的Kaplan-Meier生存曲线显示,高风险组的患者生存率明显低于低风险组(logrankP＜0.000 1),表明预后特征能够显著区分不同预后的患者。同时,为了评估模型的灵敏度和特异性,我们进行了时间依赖性ROC曲线分析,结果如图3B所示,预后模型的AUC值为0.767。用此模型对TCGA中的16个样本进行测试,测试集的AUC值为0.785,显示了模型良好的预测性能。

3 讨论

随着近几年对可变剪接的深入研究,已有文献证明异常可变剪接对于卵巢癌的发生发展具有重要影响[9～10]。本文确定了一批与卵巢癌预后相关的AS事件,并通过融合PPI数据构造了卵巢癌预后加权网络。此外,通过可变剪接预后特征建立了预后模型,模型显示,这些预后特征可以较好地预测卵巢癌患者的预后情况。因此,这些预后特征可能是卵巢癌患者潜在的预后因素,可作为癌症预后的生物标志物和治疗靶点。

本研究的拓扑分析表明卵巢癌预后加权网络为小世界网络,小世界指数为4.064。这意味着,一方面,网络中的节点高度聚集:当一个节点连接到另外两个节点时,后两个节点也倾向于彼此直接连接；另一方面,网络中的平均最短路径长度几乎与随机网络一样低,说明预后加权网络具有高效的信息传递能力。因此,本文利用网络节点的拓扑性质筛选了可变剪接预后特征,并以此构建了卵巢癌的预后预测模型,模型显示了良好的预测性能。

表2 用于构建预后模型的AS预后特征Table 2 Prognostic signatures of AS events involved in prognostic model

图2 PASEGs的富集分析(A)PASEGs的GO富集分析；(B)PASEGs的KEGG富集分析中排名前10的通路。Fig.2 The enrichment analysis of PASEGs(A)GO enrichment analysis of PASEGs；(B)KEGG enrichment analysis of PASEGs(top 10).

图3 可变剪接预后特征的Kaplan-Meier生存曲线及ROC曲线(A)预后特征的Kaplan-Meier生存曲线；(B)预后模型的ROC曲线。红色曲线表示训练集,紫色曲线表示测试集。Fig.3 Kaplan-Meier plot and ROC curve of alternative splicing prognostic signatures(A)Kaplan-Meier plot of prognostic signatures；(B)ROC curve of prognostic model.Red curve and purple curve indicate training set and test set,respectively.

在可变剪接预后特征所对应的基因中,AP2B1[26]、ATF5[27]、BCL2L11[28]、ING3[29]、SOD2[30]、USP36[31]等已被证实在卵巢癌细胞增殖和生长中具有重要功能,这表明拓扑学上重要的节点往往在疾病中发挥关键作用。例如:与化学敏感性细胞或组织相比,AP2B1在卵巢癌化学抗性细胞或组织中始终过表达[26]；ATF5在大多数上皮性卵巢癌样本中高表达,且与晚期临床分期和肿瘤分化差异显著相关[27]；BCL2L11介导的AKT磷酸化水平可影响卵巢癌细胞中顺铂的耐药性[28]。虽然现阶段其余基因在卵巢癌中的功能尚未得到证实,但已被证明与癌症相关。例如:DUSP28在胰腺癌中高度表达并发挥关键作用,靶向DUSP28可以抑制恶性胰腺癌的发展[32]；IL17RC缺失的结直肠肿瘤细胞表现为侵袭、生长和转移,提示IL17RC可作为结直肠癌的预后标志物[33]。因此,这些基因可能是卵巢癌的潜在疾病基因,值得进一步研究。

在PASEGs的GO功能富集分析中,我们发现了几类一致的GO注释。例如:预后相关基因主要参与细胞核、核质等细胞组分,与RNA聚合酶II启动子转录调控、DNA转录调节、蛋白质结合、核酸结合等生物学过程一致。显然,这些生物学过程与细胞中遗传物质的转录、蛋白质的合成等具有密切的相关性,可对细胞活动起到一定的调控作用,进而影响卵巢癌的发生发展。PASEGs还参与了细胞凋亡过程,这与KEGG通路富集分析的结果相互印证。KEGG通路富集分析表明,PASEGs主要显著富集在肠道免疫网络、p53信号通路、凋亡途径、核糖体相关途径等,提示这些通路可能是与卵巢癌AS预后相关的重要途径。已有研究显示,通过激活MDM2(或HDM2)和MDM4的可变剪接,可诱导p53通路发生改变[34]。此外,卵巢癌细胞的凋亡程序受到p53家族转录因子的调控[35],而p53的促凋亡功能是癌症治疗的有效靶标[36]。因此,在卵巢癌的发生发展机制方面深入研究这些相互作用的AS预后相关基因将是非常有意义的。

总的来讲,本文在全基因组水平上对卵巢癌的可变剪接事件进行了综合挖掘,揭示了AS事件在卵巢癌中的预后价值。特别是在预后预测中,对mRNA可变剪接的分析可能揭示了新的癌症驱动因素,并以此提出了理想的卵巢癌预后预测特征。