<BRIEFINGS IN BIOINFORMATICS> ISSN: 1467-5463;
2018_IF = 9.101;
2019_IF = 8.99;
中科大类: 生物 2区;
中科小类: 生化研究方法 1区, 数学­­与计算生物学 1区;
JCR分区: Q1

文章利用TCGA SpliceSeq 数据库对肾癌细胞（renal cell carcinoma，RCC）分类，并系统分析RCC各个亚型与临床病程、遗传变异、剪接因子、通路激活和免疫异质性的相关性。研究人员发现，选择性剪接对于RCC亚型的分类很有帮助的。接下来文章会简要介绍这个思路是如合展开的。

背景概述|

肾癌或肾细胞癌（RCC）在美国2020年新增确诊病例73750例，死亡14830例，在美国男性肿瘤中排名第六，女性肿瘤中排名第八。RCC包括KIRC（75%）,KIRP（10%-16%）和KICH（~5%）不同的病理亚型。有文献指出，与近端肾癌相比，远端肾癌的KICH线粒体基因组拷贝数增加，促进了Krebs循环和电子传递链途径的激活并增加了大量的ATP。代谢异常的KICH组表现出更差的总体生存周期。另外KIRP研究发现，存在2种临床亚型，一种为惰性型，一种为侵袭型。两种亚型分别以嗜碱性细胞质的小细胞丰富和嗜酸性细胞质的大细胞丰富，并且其分子水平也存在差异。KIR异质性强，并且出现较多的遗传变异，染色质重塑现象。

 可以看出RCC作为研究对象，其本身多样性和异质性就已经增加了研究难度。这里考虑顺式作用效应（cis-acting effect）和反式作用效应（trans-acting effect）的异常可变剪接对于RCC产生的影响。

研究数据|

TCGA SpliceSeq database (http://bioinformatics.mdanderson.org/TCGASpliceSeq/index.jsp)

从上述数据库获取三种肾癌可变剪接数据，拼接率（PSI）从0到1表示外显子变化。同时获取RNA-seq表达数据（FPKM-UQ），UCSC Xena database (https://xenabrowser.net)数据库获取GDC泛癌队列临床数据，cbioPortal database (http://cbioportal.org)数据库获取TCGA RCC队列的病理报告，molecular signatures database (MSigDB version 7.0) (https://www.gsea-msigdb.org/gsea/msigdb)获取hallmark (H)、position (C1)、curated canonical pathways (C2-CP)和gene ontology terms (C5)信息，最后从CancerSplicingQTL database (http://www.cancersplicingqtl-hust.com/#/)这个数据库获取剪接定量性状位点(sQTLs)数据。

方法总述|

经过上面所有的数据准备，基本上可以进行RCC亚型的可变剪接情况分析。作者这里使用基于t-SNE方法，进行聚类分析然后进行亚型分类，将TCGA RCC样本分为了KIRC, KIRP, KICH 和 KIRN四类。接着，对RCC样本进行选择性剪接，通路激活喝mRNA表达的差异分析，这里对PSI进行卡方检验和RCC亚型与KIRN差异选择性剪接事件（DASE）进行组间Wilcoxon秩和检验。然后使用最大熵模型得到log odds得分，TCGA SpliceSeq database获取TCGA基因结构文件，根据外显子剪接情况，确定候选片段的起始位置和结束位置，进而分别得到不同剪接类型中剪接位点区域的边界。参考序列是GRCh37.p13，3‘ss表示3’端剪接位点MaxEnt分数，5‘ss端剪接位点MaxEnt分数。采用非参数无监督方法进行基因集变异分析（GSVA），目的是样本中进行富集分析。然后对KICH、KIRC、KIRP和KIRN的GO项进行比较。考虑到异质性存在，这里使用了剪接相关抗原(SRAs)分析，理解PSI在KIRN样本种存在差异是否会产生不规则剪接相关的新蛋白异构体或已知异构体的变异，即SRA。

结果分析|

RCC亚型的剪接聚类与临床特征

整合TCGA队列KICH、KIRC和KIRP的剪接集，889例RCC肿瘤和129例邻近正常组织(KIRN)，共1018例889患者纳入本研究。包括典型的五种剪接类型-交替受体位点(AA)，交替供体位点(AD)，外显子跳变(ES), 互斥外显子(ME)和保留内含子(RI)，在KICH队列中8281个基因中检测到27303个ASE，KIRC组8630个基因中检测到28274个ASE，KIRP中的7963个基因检测到24757个ASE，7454个基因中的22,839个ASE在3种亚型中都存在。

经过t-SNE分析，如图1所示，同一TCGA亚型样本发生聚类，RCC样本大致分为4类，分别对应KIRC(类群1)、KIRP(类群2)、KICH(簇3)和KIRN(簇4)。但从图1A和1B中，可以看到，仍有一些样本落入错误的聚类集，这里重新根据病理报告，文献等，重新定义肾癌样本的亚型分类如图1B和1C所示，除了RCC亚型和KIRN的分离外，剪接分类还得到了包括代谢紊乱在内的两个小类别KICH亚型KICH- md(6个样本均不同于KICH聚类)和KIRP的CIMP样亚型KIRP- cl(10个中有9个在前期研究中报道为CIMP亚型)。与其他亚型相比，这两种亚型的样本的总生存率明显较低分别为KICH和KIRP，在本研究中Kaplan-Meier plot也得到了证实 (图1E)。此外，KIRC的样品进一步划分为4个亚型，包括KIRC-T1、KIRC-T2A、KIRC-T2B、KIRC-T2C(图1B-C) 。Kaplan-Meier图显示，这些KIRC亚型之间的临床结果有显著差异(图1D)。

Fisher’s exact检验的临床相关性分析显示，几个侵袭性RCC亚型具有明显的特征，病理肿瘤大小较大(T)，远端转移(M)和整体TNM晚期，晚期Fuhrmann分级和肉瘤样改变可在图1F的KIRCT2B中看到。更多的非裔美国人具有KIRC-T2B的特征和KIRC-T2C的特征并且在KIRC的临床效果最差（图1F-L）。这些数据表明选择性剪接在不同的RCC组织病理类型中存在差异，且临床特征与患者生存期高度相关，因此适合分型。

RCC亚型DASEs剪接特点

差异选择性剪接分析表明，所有RCC亚型中，上调的DASEs比下调的DASEs要多(图2A)。在前文提到的22,839个ASE中，RI的比例升高，ES的比例降低。在下调的DASE的RCC亚型中，ES占比增加，AA和AD占比相对减少（图2B）。通过重叠分析，尽管上调的DASE的数量变化很大，但是在上调组和下调组中，不同KIRC亚型的DASE依然重叠(图2C)。因此认为KIRC-T2B是由于剪接效率低下而保留了大量的内含子和外显子区。

为了描述不同KIRC亚型DASEs的拼接特征，剪接位点评分、GC含量和改变长度，这里对于所有的DASEs分别进行计算。如图2D-F所示，下调的PSI诱导外显子跳跃的DASE明显表现为3 ' ss较弱，GC含量较高，外显子长度改变较短，对于AA、AD和RI上调的PSI则表现出相反的现象。脾酪氨酸激酶（SYK）参与信号转导，SYK的长剪接体，在KIRC所有亚型中均显著过表达，并在RT-PCR中对KIRC组织进行验证（图2G,H）。在KIRC亚型中，ICAM3 mRNA表达上调，PSI值下调，二者显著相关。总的来说，KIRC-T1和KIRC-T2B中，肿瘤相关通路剪接的异质性和均一性发生了变化，这可能表明选择性剪接在肿瘤发生和维持癌细胞生存中的作用不同。在以低效剪接为特征的侵袭性KIRC亚型中，非常规剪接可能通过调节癌症生物通路中基因的亚型切换和mRNA表达促进癌症发展。

RCC亚型富集到的基因集和不同的剪接调控模式

从MSigDB数据库对所有RCC样本中的12515个基因集进行GSVA分析，基因集的大小一般设置为10-500个基因，如图3A所示，与KIRN相比，KIRC-T2B、KIRCT2A和KICH-MD亚型样本的GSVA评分相对较低，KIRC-T1亚型样本的GSVA评分较高。差异分析在KIRC-T1中鉴定出1554个上调和365个下调DESs，而KIRC-T2B中有410个上调DESs和 2087个下调DESs(图3B)。

为了广泛研究通路和剪接之间的潜在关系，分别对KIRC、KIRP、KICH和KIRN样本进行了Spearman相关分析。如图3C所示，在KIRC中，相关性大多为负，它们大多是保留内含子和较低的PSI值。因此，进一步分析这些相关基因集和ASEs在不同样本类型之间的差异。研究发现，KIRC和KIRP的基因集中于对刺激反应和cellular anatomical entity， KICH集中于转录后调控基因表达和ATPase活性而KIRN则侧重于调节的生物质量和转运蛋白活性(图3D)。这些结果部分说明了RCC亚型与正常组织的细胞内和细胞外环境的差异以及相应的不同的剪接调控机制。

接着，图3E-G中可以看出，大多数ASEs与剪接通路相关（SP）相关性小于-0.5，与剪接因子（SF）相关性大于0.5。这些ASEs通常在KIRC-T2A、-T2B和-T2C中表达上调，而在KIRC-T1中不表达，通常与相关SFs、SPs或RTS通路的相应变化相关。此外，维恩图显示SP相关的ASEs主要包含在ASEs相关RTS通路中，KIRC-T2B中的DASEs有67%的上调和17%的下调（图3H-I）。因此，KIRC-T2B中SFs通路的激活和表达等亚型特异性变化可能与其严重的异常剪接相关。

基因组变异导致RCC的剪接缺失

这里重点分析RCC亚型中剪接位点的基因组结构的变化情况。使用301个KIRC样本的VHL变异和剪接数据和149个VHL变异数据纳入分析。VHL变异区位置分为剪接区、外显子区、侧翼区和基因区4个相互独立的区域。拼接区域从5‘到3’，分为先5‘ss,先3’ss,再5‘ss 和再3’ss(图4A)。3‘和5’拼接区域包括典型的拼接位点和相邻序列Ent得分推导出的3’ss 和5’ss区。与VHL野生样本相比，VHL变异样本在第一个3‘ss，第二个5’ss和交替外显子显著促进了第二个外显子的跳跃，这与符合MaxEnt评分的下降和较弱的剪接强度一致。然而对于第二个3‘ss变异，剪接位点受损和减弱不会导致VHL显著的外显子跳变（图4B-D）。因此，具有剪接相关VHL突变的KIRC患者可能具有更强的肿瘤侵袭性，需要进行更多的个体化评估和治疗。

另外图4E所示，剪接数量性状位点划分为与前面提到类似的染色体、基因、外显子和剪接区域，通过比较不同区域每百万碱基中与剪接相关的SNPs数量，发现外显子SNP和剪接区SNP更有可能影响KIRC的选择性剪接（图4 F-H）。

RCC中肿瘤微环境和免疫异质性与不规则剪接相关

血管生成被证实是免疫细胞浸润的先决条件之一，而KIRC具有血管生成能力是由组成性假缺氧表型和缺氧诱导因子(HIF-1α)通路的激活所决定的。在本研究中，与其他RCC亚型相比(图5A)，KIRC患者的免疫细胞浸润率较高。有趣的是，KIRC- t2b的血管生成、缺氧和HIF-1α信号明显低于其他KIRC簇。

不良临床结果的KIRC患者中存在剪接缺失的情况，如图5B-C所示。考虑到KIRC具有复杂的肿瘤间和肿瘤内异质性，确定特定的治疗性生物标志物以提高患者的总体生存率和减少经济负担仍然具有挑战性和必要性。这些免疫特征在KIRC剪接亚型之间的差异表明了对ICIs和TKIs可能有不同的治疗反应，因此，不规则剪接可以作为潜在的生物标志物，并指导晚期RCC患者选择个性化治疗方案。

 简要回顾|

一般而言，我们对于TCGA数据的样本情况分析使用上，很难想到从亚型的基因结构出发，进行新的临床分型。这个文章提供了很好的思路，并且关键性做图都很漂亮，下一期将出一篇专门的文章来学习一下他们的基因剪接分析方法，敬请期待吧。