CN110010198A - 一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法及预后模型应用 - Google Patents

Abstract

本发明涉及基因技术及生物医学领域，具体涉及肝细胞癌（HCC）中可变剪切（AS）事件的全转录组系统分析的方法以及基于全转录组分析方法的相关预后模型构建与应用。本发明通过在HCC中发生了剪接变异信号的全局辨识和生物信息学分析的全面描述，以及基于系统的全转录组分析方法的预后模型及剪切网络的构建。本发明提供了HCC中mRNA剪接特征的全局变化的综合图像，并进一步开发了包含存活相关AS事件的预后模型和剪接网络。本发明是对HCC AS事件特征进行全面综合计算研究的第一次尝试，并进一步扩展到预后和分子靶向意义的新领域，为阐明生物发生和HCC的发展的潜在遗传机制提供了新的见解。

Description

一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法 及预后模型应用

技术领域

本发明涉及基因技术及生物医学领域，具体涉及肝细胞癌（HCC）中可变剪切（AS）事件的全转录组系统分析的方法以及基于全转录组分析方法的相关预后模型构建与应用。

背景技术

肝细胞癌（HCC）是全球癌症死亡的第三大原因，每年报告的死亡人数约为50万。在临床中，HCC由不同的亚型组成，具有不同的结果，由于异质性而具有悲观的预后，传统的预后因素对预测HCC患者的个体化预后没有帮助。因此，进一步全面研究HCC中涉及的生物学机制和分子生物标志物是开展靶向治疗和改善HCC患者预后预测的重要步骤。

最新研究主要考虑加速全基因组癌症研究的高通量技术的快速发展。近年来，HCC的许多诊断，预后和治疗分子生物标志物逐渐被鉴定和验证，包括非编码RNA，CpG甲基化位点和拷贝数变异。这些研究尽管结果很有希望，但主要集中在改变HCC的转录组分析水平。然而，转录后修饰的系统分析总是被忽略，而且很大程度上是未知的。

从分子的角度来看，HCC是异质性肿瘤。在过去十年中，已经做出了重大努力来揭示参与HCC发展的基因组谱的分子变化。这些研究有助于确定预后遗传特征，包括基因，微小RNA和非编码RNA等特征促进了相关预后标志物甚至治疗靶标的鉴定。然而，AS作为增加基因组编码能力和扩展蛋白质多样性的主要转录后生物学行为，在癌症生物学中表现出更多的潜在意义。

先前对AS特征的探索已经表明，异常调节的剪接变体可以用作HCC的诊断，预测和预后生物标志物。然而，缺乏关于在HCC中建立具有异常AS的预后特征的详细信息。因此，用于判断HCC预后的全转录组分析方法没有系统描述。

发明内容

鉴于现有技术存在的问题，本发明的目的在于提供了一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法及预后模型应用。本发明通过在HCC中发生了剪接变异信号的全局辨识和生物信息学分析的全面描述，以及基于系统的全转录组分析方法的预后模型及剪切网络的构建。本发明确定了290名HCC患者的1314个基因中的1805个AS事件，这些事件与HCC患者的总体存活率显着相关。最终的预后模型包括26个AS事件被构建。风险评分（全部）的接受者 - 操作者特征（ROC）曲线下面积达到0.835，用于预测5年内HCC患者的生存状态。分类为中位数风险评分（全部）值的高风险和低风险评分组显示出显着不同的存活时间（高与低= 827对3125天，P <2e-16）。此外，剪接调控网络显示出HCC中加香因子作用的明显趋势，这可能是与HCC患者预后相关的潜在生物学机制。此外，Gene Ontology（GO），Kyoto Encyclopedia of Genes and Genomes（KEGG）富集分析以及GSVA揭示了HCC中剪接机制中涉及的某些关键术语。本发明提供了HCC中mRNA剪接特征的全局变化的综合图像，并进一步开发了包含存活相关AS事件的预后模型和剪接网络。

一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法的应用，通过建立HCC中mRNA剪接特征的全局变化的综合图像，并进一步建立包含存活相关AS事件的预后模型和剪接网络。

为了实现上述目的，本发明采用以下技术方案。

一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法，具体步骤如下。

1）数据管理过程：HCC队列的第三级mRNA测序数据和相应的临床信息获自TCGA数据门户http://tcga-data.nci.nih.gov/，其包含374个HCC组织和50个相邻的非肿瘤组织。为了生成每个HCC患者的AS概况，使用一种明确量化每个外显子和剪接点的包含水平的Java应用程序即SpliceSeq来评估HCC患者的mRNA剪接模式。共纳入290名符合以下标准的患者：A.确诊HCC的组织学诊断；B.患者在初始病理诊断后至少30天仍存活；C.患者具有相应的mRNA剪接数据。从0到1的百分比拼接（PSI）评分通常用于评估特定基因与特定剪接模式的转录物比率。为尽可能地生成一组AS事件，实施一系列严格的过滤器，PSI得分≥75的样本百分比，PSI得分平均值≥0.05，并且使用k-最近邻居估算了缺失的表达值（ kNN）算法与impute包。为了准确描述AS事件，为每个AS事件分配了一个唯一的注释，该注释由SpliceSeq中的剪接类型，ID号和匹配的基因符号组成。例如，“RI_C9orf9_ID_87994”表示保留的内含子（RI）表示剪接类型，C9orf9是父基因符号，ID_87994表示SpliceSeq数据库中的特定序列号。

2）HCC患者的生存分析和预后特征的构建：根据中位数切点，将HCC患者AS事件的PSI评分分为低指数和高指数，然后进行适合分类PSI指数和患者总生存（OS）的单变量Cox回归分析，以确定预后AS事件。 UpSetR是一个用于生成镦粗图（ upset plot）的R包，可以可视化七种类型的生存相关AS事件之间的交叉集（如图1d和图1c所示）。前20名（如果可用）最重要的生存相关选择每种AS类型的AS事件作为搜索独立预后因子的候选者，并分别通过多变量Cox回归进一步构建HCC患者的AS特征。每个预后特征选择通过前向逐步选择过程和Akaike信息标准（AIC）[23]进行。然后，将每个AS特征中识别的AS事件组合起来，通过二级前向逐步方法构建最终预测模型，以找到最小的AS事件集。其中，使用以下方法计算OS预测的每个AS特征的风险评分：，此外根据中位风险评分将HCC患者分为低风险和高风险亚组，并进行Kaplan-Meier生存分析以区分这两个亚组。此外，每个预后模型的临床疗效和预测准确性是可比拟的，拟合时间依赖性ROC曲线与survivalROC包，并且分别计算1年，3年和5年的曲线下面积（AUC）。

3）生存相关AS事件的功能丰富和相互作用分析：通过单变量Cox回归分析确定的存活相关AS事件的亲本基因通过设定错误发现率（FDR）<0.05进行基因本体论（GO）和京都基因百科全书和基因组（KEGG）途径分析。上述分析使用clusterProfiler包（版本3.7）进行。此外通过使用Reactome FI Cytoscape PlugIn（版本3.6.1）分析和可视化基因相互作用网络，目的是搜索重要的中枢基因。

4）HCC队列中潜在的SF-AS监管网络：通过筛选文献和SpliceAid2数据库鉴定总共67个SF，并且使用DESeq2包将TCGA HCC样品中的相应表达水平用“方差稳定转化功能”整合和标准化。然后进行单变量Cox危害回归以筛选独立的预后SF。还进行了存活相关SF的表达与预后AS事件的PSI评分之间的Spearman相关性分析。 P值由Benjamini＆Hochberg（BH）相关调整，调整后的P值<0.001被认为是显着的。此外由Cytoscape生成SF-AS监管网络。

5）SF-AS相关对的基因集变异分析（GSVA）：GSVA是一种差异功能基因集富集分析，可通过计算样本方式基因集富集分数来检测异质样本的细微途径活性变化。使用来自bioconductor（www.bioconductor.org）的GSVA包进行GSVA以进一步从HCC组织中搜索显着富集的GO和经典途径（KEGG，Reactome和BioCarta途径数据库），（从Molecular SignaturesDatabase（MSigDB）下载）（www.broadinstitute.org）。然后使用limma包进行肿瘤与邻近正常样本之间的差异基因集分析，并且对于GO术语，| logFC |> 0.4和FDR <0.05，以及|logFC |> 0.2且FDR <0.05用于途径组分别具有统计学意义。将GSVA结果的火山图和热图可视化并显示，如图5所示。

与现有技术比，本发明的有益效果如下。

本发明对290个HCC队列中的AS特征进行了分析，然后与构建的最佳预后模型进行了整合。最终结合所有可用AS类型的最终风险评分（全部）得到包括26个亲本基因的26个AS事件，并增强了临床应用中预测能力的预测能力，用于预测1年，3年和5年AUC上升的HCC患者的存活率分别为0.852,0.753和0.788。本发明是对HCC AS事件特征进行全面综合计算研究的第一次尝试，并进一步扩展到预后和分子靶向意义的新领域。

鉴于肿瘤微环境中广泛的AS改变可能由有限数量的SF广泛协调，因此关注预后SF和SF-AS剪接相关网络。根据本发明的研究，共有165例与生存相关的AS事件与HCC患者的预后SF显着相关，包括SFPQ、HNRNPM、PCBP2、KHDRBS3和SRSF10，而所有这些SF都被认为是不利的预后因素。此外，本发明可以确认大多数有利的预后AS事件基于SF负调节，而不良预后AS事件在HCC中基于SFs正调节异常。此外，这些剪接的功能性GSVA分析也表明剪接模式及其与HCC上游SF的机制关联可能涉及若干致癌过程，例如细胞增殖和免疫逃逸，这为阐明生物发生和HCC的发展的潜在遗传机制提供了新的见解。

附图说明

图1是HCC中AS事件分析概述，其中A为用于在大规模RNA-Seq数据中分析HCC患者的可变剪接事件的流程图；B为七种AS事件的插图，包括外显子Skip（ES），互斥外显子（ME），保留内含子（RI），替代启动子（AP），替代终止子（AT），替代供体位点（AD） 和AlternateAcceptor站点（AA）；C为根据分类的P值0.05描述来自290个HCC患者的AS事件和涉及的基因；D为在HCC中七种类型的存活相关AS事件之间的父基因相互作用的UpSet图。

图2是在HCC中存活相关替代事件的风险比的Forrest图。 其中A-F分别为前20个基因的危险比与存活相关的ES、AA、RI、AT、AP和AE事件；G为显著存活相关ME事件的危险比。（ 此外AT、AP的x轴坐标比例为对数，如果95％CI的上限或下限超出预定义边界，则将置信区间剪切为箭头。）

图3是HCC中7种类型的可变剪接事件的预后信号的Kaplan-Meier图，其中A-G为分别用HCC中的ES、AT、AD、RI、AA、AP和ME建立的预后预测因子的Kaplan-Meier曲线，黄线表示高风险亚组，蓝线表示低风险亚组，分别显示每个特征（虚线）中的P值和中值存活时间。

图4是使用所有类型的预后相关剪接事件，基于预后相关剪接事件构建和评估风险评分。其中A为使用所有类型的存活相关剪接事件建立的预后预测因子的Kaplan-Meier曲线，蓝线表示低风险子组，黄线表示高风险子组；B为HCC中3年和5年的预后信号的ROC曲线；C为基于生存相关剪接事件的风险评分构建与分析，HCC患者分为低风险和高风险组。

图5是HCC中存活相关AS事件的基因网络和功能富集分析，其中A为由Cytocape生成的HCC中存活相关AS的基因网络；（B-E），B为生物过程中的GO术语的前20个途径；C为分子功能；D为细胞组分；E来自存活相关AS事件的亲本基因的KEGG途径分析。 节点的大小和颜色分别代表富集的基因数和FDR值。

图6为HCC中剪接因子预后信号的Kaplan-Meier生存分析。其中A-E分别为HFP患者的SFPQ、HNRNPM、PCBP2、KHDRBS3和SRSF10的Kaplan-Meier生存曲线分为低（绿色）和高（红色）风险组，根据最佳切点（每张图像上方显示） mRNA表达值， 每个模型显示P值，风险比和95％CI。

图7为HCC中的拼接相关网络。其中a为HCC中的拼接相关网络，5个存活相关的剪接因子（蓝点）呈阳性（红线）或阴性（绿线）与86个有利AS事件（绿点）和102个不良AS事件（红点）相关；b和c为从SFs-ASs网络中选择的相关性的代表性点图。

图8是SF-AS相关对的基因集变异分析（GSVA）。 其中A和C为火山图，火山图显示了在HCC癌症和癌旁组织以及血清中用于SF-AS相关对的GO术语和途径的差异富集结果，用文本标记标记； B和D为Heatmaps分别显示了GO项目和SF-AS相关对路径的差异富集结果。

具体实施方式

下面结合附图和实施例详细描述本发明，以下所述仅是本发明优选实施方式，应当指出，对于本技术领域的普通技术人员，在不脱离本发明方法的前提下，还可以做出若干改进和补充，这些改进和补充也应视为本发明的保护范围。

实施例 一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法及预后模型应用。

1.TCGA HCC队列中AS事件概。

具有AS事件概况的290名HCC患者的整合mRNA测序数据包括在本研究中。中位随访期为15.75个月（范围1-122.5个月）。总体而言，基于SpliceSeq分析的结果，检测到来自8986个亲本基因的总共34163个AS事件。如图1B所示，通过不同的机制产生了七种类型的AS事件。此外，AS事件和每种AS类型的基因被分成存活相关组（p <0.05）或存活无关组（p>0.05），如图1C所示。从1314个亲本基因中完全鉴定了1805个存活相关AS事件，其中包含479个基因中的565个外显子跳跃（ES），356个基因中的541个替代终止子（AT），258个基因中的341个替代启动子（AP），137个替代受体位点（ AAs）在128个基因中，104个基因中的112个保留内含子（RI），95个基因中的100个替代供体位点（AD）和9个基因中的9个互斥外显子（ME）。注意到一个基因可能有两个或多个与OS相关的AS事件，因此生成UpSet图来定量分析七种AS类型之间的交互集（如图1D所示）。因此大多数与生存相关的AS事件属于一个亲本基因，而几个基因可能具有多达四种类型的AS事件，这些事件均显着与预后相关。例如基因CIRBP（红点线）中的RI、AA、AP和ES与HCC中的OS显着相关，并且基因BSCL2（绿点线）中的AD，AP和ES与HCC中的OS显着相关。

2. HCC患者的生存分析和预后特征的构建。

为确定独立的预后因素，进行单变量Cox存活分析以评估临床OS与每种AS事件之间的关系。在每个AS事件中，基于PSI中位数分数将HCC患者分成两个亚组。 7种AS中最重要的20个存活相关基因（除了仅有9个基因的ME）在图2中以森林图表的形式呈现。如上图所示，无论AS类型如何，与生存相关的AS事件似乎几乎相同，均分布在有利的预后和反向亚组中。值得注意的是，用11个ES事件，6个AT事件，10个AD事件，10个AA事件，7个AP事件，6个RI事件或7个ME事件构建的AS信号在区分HCC患者的不良或良好临床结果方面都表现出很强的预测能力（ p <0.0001），而AA、AP、AT和ES在7种类型中最为强大，低风险亚组的中位生存期超过3000天，如图3B所示。然后进一步构建HCC的最终预后模型，仅包括26个AS事件以使AIC最小化（高与低= 827对3125天，P <2e-16，如图4A所示）。此外应用1年、3年和5年的ROC曲线来比较不同AS模型之间的效率。如图4B所示，AUC证实具有所有类型AS事件的最终预后模型与其他事件相比表现出更高的预后效率，1年，3年和5年的AUC分别为0.852,0.783和0.788。另外，还显示了患者的存活状态和风险评分的分布，每个特征中包括的AS事件的剪接模式，如图4C所示。

3. 生存相关AS事件的功能丰富和相互作用分析。

将HCC中1805个存活相关AS事件的总共1314个亲本基因用通道输送到Cytoscape软件并用Reactome FI PlugIn进一步分析以产生基因相互作用网络，结果揭示了重要的途径，包括TP53、EGFR、HSPA8和UPF3B（如图5A所示）等中枢基因。 为进一步阐明预后AS事件对其相应蛋白质的潜在影响，进行GO和KEGG富集分析。在生物过程（BP）方面总共鉴定了176个GO术语，包括细胞粘附分子结合和钙粘蛋白结合途径（如图5B所示）。 20个分子功能（MF）术语（如图5d所示）和71个细胞成分（CC）术语（如图5D所示）也显示极大富集，表明嘌呤相关代谢过程，蛋白质靶向以及粘着斑和细胞表面连接途径有显著变化。此外，KEGG富集分析显示共有48条富集的途径，其中大部分与肝癌三部曲，肝癌细胞的侵袭性和远处转移，癌症途径和核苷酸代谢有关。更重要的是，众所周知的一些KEGG途径也被证实可以参与肿瘤的发生和HCC的进展，例如p53信号通路（FDR <0.028），AMPK信号通路（FDR <0.017），HIF-1信号通路（ FDR <0.0059）和EGFR酪氨酸激酶抑制剂耐药（FDR <0.017）。巧合的是，还鉴定了癌症相关途径，例如结肠直肠癌，神经胶质瘤，胰腺癌，子宫内膜癌，非小细胞肺癌等（如图5E所示）。总之发生预后AS事件的亲本基因可能在HCC的致癌，进展和转移中起重要作用。换句话说，这些特异性失调的AS事件可能协调亲本基因的转录后修饰并进一步改变HCC患者的蛋白质特征。

4. HCC队列中潜在的SF-AS监管网络。

获得TCGA HCC队列的3级RNA测序数据，5个SF的表达谱，包括SFPQ（P <0.001，HR =2.02，95％CI：1.25-3.27），HNRNPM（P = 0.015，HR = 1.7， 95％CI：1.06-2.72），PCBP2（P =0.009，HR = 1.77，95％CI：1.09-2.85），KHDRBS3（P = 0.012，HR = 1.93,95％CI：1.22-3.03）和SRSF10（P = 0.019，HR = 1.74，95％CI：1.12-2.7）与OS相关，并且在每一个中发现了对应于相反预后亚组的最佳分离的最佳切割点以进一步验证，如图6所示。此外进行这5个OS相关SF的表达水平（p <0.05）与预后AS事件的PSI评分（p <0.05）之间的spearman相关性分析，以估计潜在的调节关系，然后这些显着的相关对（p <0.001）在图7a中呈现为SF-AS调节网络，其显示165个存活相关的AS事件，包含77个AS事件（HR> 1）表明预后不良，88个不良事件（HR <1）与这5个SFs（P <0.001）显着相关。 值得注意的是，大多数不良生存预后AS事件（红点）与SFs（蓝点）呈正相关（红线），而最有利的生存预后AS事件（绿点）与SFs呈负相关（绿线）。此外，SF和AS事件之间的代表性相关性显示在符合线性回归的散点图中（如图7b和7c所示）。例如SFPQ的表达与ZNF91的AT呈正相关，而与SMAD6的AP呈负相关。

SF-AS相关对的基因集变异分析（GSVA）。

进一步探索SF-AS调节网络的潜在功能，无论是正相关还是负相关。为搜索显着富集的GO和经典途径的基因组，在HCC组织中的KEGG、Reactome和BioCarta途径数据库，进行关于肿瘤和邻接正常样品的GSVA和差异富集分析。此外，结果分别显示为火山图和聚类热图（如图8所示），其显示肿瘤组织表现出增加的细胞增殖活性和与细胞周期相关的蛋白质表达，同时降低免疫应答、细胞粘附活性和正常的生理代谢发生。

Claims (2)

1.一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法，其特征在于，具体步骤如下：

1）数据管理过：HCC队列的第三级mRNA测序数据和相应的临床信息获自TCGA数据门户http://tcga-data.nci.nih.gov/，其包含374个HCC组织和50个相邻的非肿瘤组织；

为了生成每个HCC患者的AS概况，使用一种明确量化每个外显子和剪接点的包含水平的Java应用程序即SpliceSeq来评估HCC患者的mRNA剪接模式，共纳入290名符合以下标准的患者：1）确诊HCC的组织学诊断；2）患者在初始病理诊断后至少30天仍存活；3）患者具有相应的mRNA剪接数据；从0到1的百分比拼接评分通常用于评估特定基因与特定剪接模式的转录物比率，为尽可能地生成一组AS事件，实施一系列严格的过滤器，PSI得分≥75的样本百分比，PSI得分平均值≥0.05，并且使用k-最近邻居估算了缺失的表达值算法与impute包；为准确描述AS事件，为每个AS事件分配了一个唯一的注释，该注释由SpliceSeq中的剪接类型，ID号和匹配的基因符号组成；

2）HCC患者的生存分析和预后特征的构建：根据中位数切点，将HCC患者AS事件的PSI评分分为低指数和高指数，然后进行适合分类PSI指数和患者总生存的单变量Cox回归分析，以确定预后AS事件；UpSetR是一个用于生成镦粗图的R包，可以可视化七种类型的生存相关AS事件之间的交叉集，前20名最重要的生存相关选择每种AS类型的AS事件作为搜索独立预后因子的候选者，并分别通过多变量Cox回归进一步构建HCC患者的AS特征，每个预后特征选择通过前向逐步选择过程和Akaike信息标准（AIC）进行，

然后将每个AS特征中识别的AS事件组合起来，通过二级前向逐步方法构建最终预测模型，以找到最小的AS事件集，其中使用以下方法计算OS预测的每个AS特征的风险评分：，此外根据中位风险评分将HCC患者分为低风险和高风险亚组，并进行Kaplan-Meier生存分析以区分这两个亚组，此外每个预后模型的临床疗效和预测准确性是可比拟的，拟合时间依赖性ROC曲线与survivalROC包，并且分别计算1年，3年和5年的曲线下面积；

3）生存相关AS事件的功能丰富和相互作用分析：通过单变量Cox回归分析确定的存活相关AS事件的亲本基因通过设定错误发现率（FDR）<0.05进行基因本体论（GO）和京都基因百科全书和基因组（KEGG）途径分析，上述分析使用clusterProfiler包-版本3.7进行，此外通过使用Reactome FI Cytoscape PlugIn-版本3.6.1分析和可视化基因相互作用网络，目的是搜索重要的中枢基因；

4）HCC队列中潜在的SF-AS监管网络：通过筛选文献和SpliceAid2数据库鉴定总共67个SF，并且使用DESeq2包将TCGA HCC样品中的相应表达水平用“方差稳定转化功能”整合和标准化；然后进行单变量Cox危害回归以筛选独立的预后SF；还进行存活相关SF的表达与预后AS事件的PSI评分之间的Spearman相关性分析； P值由Benjamini＆Hochberg（BH）相关调整，调整后的P值<0.001被认为是显着的，由Cytoscape生成SF-AS监管网络；

5）SF-AS相关对的基因集变异分析：GSVA是一种差异功能基因集富集分析，可通过计算样本方式基因集富集分数来检测异质样本的细微途径活性变化，使用来自bioconductor-www.bioconductor.org的GSVA包进行GSVA以进一步从HCC组织中搜索显着富集的GO和经典途径，其中经典途径包括KEGG、Reactome和BioCarta途径数据库，从www.broadinstitute.org-Molecular Signatures Database下载，然后使用limma包进行肿瘤与邻近正常样本之间的差异基因集分析，并且对于GO术语，| logFC |> 0.4和FDR <0.05，以及| logFC |>0.2且FDR <0.05用于途径组分别具有统计学意义，将GSVA结果的火山图和热图可视化并显示。

2.一种基于全转录组的肝细胞癌可变剪切事件系统分析的方法的应用，其特征在于，通过建立HCC中mRNA剪接特征的全局变化的综合图像，并进一步建立包含存活相关AS事件的预后模型和剪接网络。