Comprehensive analysis of aberrant alternative splicing related to carcinogenesis and prognosis of papillary thyroid cancer

Abstract

As a key mechanism, alternative splicing (AS) plays a role in the cancer initiation and development. However, in papillary thyroid cancer (PTC), data for the comprehensive AS event profile and its clinical implications are lacking. Herein, a genome-wide AS event profiling using RNA-Seq data and its correlation with matched clinical information was performed using a 389 PTC patient cohort from the project of The Cancer Genome Atlas (TCGA). We identified 1,925 cancer-associated AS events (CASEs) by comparing paired tumors and neighboring healthy tissues. Parent genes with CASEs remarkably enriched in the pathways were linked with carcinogenesis, such as P53, KRAS, IL6-JAK-STAT3, apoptosis, and MYC signaling. The regulatory networks of AS implied an obvious correlation between the expression of splicing factor and CASE. We identified eight CASEs as predictors for overall survival (OS) and disease-free survival (DFS). The established risk score model based on DFS-associated CASEs successfully predicted the prognosis of PTC patients. From the unsupervised clustering analysis results, it is found that different clusters based on AS correlated with prognosis, molecular features, and immune characteristics. Taken together, the comprehensive genome-wide AS landscape analysis in PTC showed new AS events linked with tumorigenesis and prognosis, which provide new insights for clinical monitoring and therapy for PTC.

Keywords: The Cancer Genome Atlas; alternative splicing; carcinogenesis; papillary thyroid cancer; prognosis.

Figure 1

Overview of AS event profiling in the TCGA PTC cohort. (A) Number of each AS event type and their parent genes in PTC patients. Blue bars represent the AS events, while red bars represent their parent genes. (B) Interactive sets among seven types of AS events (n = 34,773) shown in an UpSet plot. (C) Circos plot shows the details of each AS event and their parent genes in the chromosome. The outer circle represents the chromosome ideogram. The intermediate circle represents the genes with filtered AS events. The inner circle shows the genes with differentially expressed AS events between tumor and normal tissues. The ribbons represent the potential interaction between AS events and their parent genes. ES, exon skipping; AP, alternate promoter; AT, alternate terminator; AD, alternate donor; AA, alternate acceptor site; ME, mutually exclusive exons; RI, retained intron. Mulity-AS, gene contains multiple types of AS.

Figure 2

Identification of CASEs in PTC. (A) Volcano plot of CASEs identified in PTC (log₂^FC > 0.1, adjusted P < 0.05). (B) Heatmap of the CASEs between matched tumor and adjacent normal tissues of 50 PTC patients. (C) Interactive sets among seven AS types of CASEs (n = 1,925) shown in an UpSet plot. (D) The representative CASEs, derived from the same parent gene and exhibited the opposite preference between tumor and adjacent normal tissues, were shown. Student’s t-test was used. (E) Venn diagram of CASEs and DEGs. (F) Interactive sets of AP/AT events and DEGs shown in UpSet plot. ES, exon skipping; AP, alternate promoter; AT, alternate terminator; AD, alternate donor; AA, alternate acceptor site; ME, mutually exclusive exons; RI, retained intron.

Figure 3

Signature enrichment by CASEs in PTC. (A) GO analysis of CASEs. (B) KEGG analysis of CASEs. (C, D) GSEA analysis of all AS events.

Figure 4

Interaction analysis of CASEs. (A) PPI network analysis of CASEs generated by Cytoscape. Nodes indicate parent genes with CASEs, while edges represent the potential interactions between the corresponding proteins. The shape, size, and color of nodes denote AS types, the value of log₂^FC, and change patterns, respectively. ES, exon skipping; AP, alternate promoter; AT, alternate terminator; AD, alternate donor; AA, alternate acceptor site; ME, mutually exclusive exons; RI, retained intron. (B) Hub genes ranked by MCC. (C) Module 1 was correlated with tumorigenesis. (D) Module 2 was correlated with extracellular matrix.

Figure 5

Representative plots of regulatory splicing network in PTC. (A) The correlation analysis between the expression levels of 71 SFs and the PSI values of CASEs. The shape, color, and size of node denote AS types, changes in the pattern (upregulated or downregulated), and the value of log₂^FC, respectively. The breadth of each line represents the extent of interaction strength. (B) Representative dot plots indicate the correlations between the expression of SFs and PSI values of CASEs.

Figure 6

Survival-associated CASEs in PTC. (A, B) Forest plots of hazard ratios for eight CASEs simultaneously associated with OS (A) and DFS (B). (C, D) Kaplan-Meier curves of representative genes associated with OS (C) and DFS (D).

Figure 7

Selection of PTC feature genes and construction of a prognosis model. (A) Identification capability of the prognostic model of the low-risk and high-risk groups. (B) Kaplan-Meier curves of prognostic model for DFS of PTC patients. (C) ROC curves of prognostic model for the 3- and 5-year DFS of PTC patients. (D) Nomogram for predicting PTC patient’s 3- and 5-year DFS.

Figure 8

AS-based clusters significantly associated with prognosis, molecular characteristics, and immune features. (A) Consensus clustering analysis identified four clusters. The white (consensus value = 0, samples never clustered together) and blue (consensus value = 1, samples always clustered together) heatmap display sample consensus. (B) Kaplan-Meier curves show the OS for four AS-based clusters. (C) Kaplan-Meier curves show the DFS for four AS-based clusters. (D) Heatmap shows the molecular characteristics associated with clinical, molecular, and immune features among the four clusters.

Figure 9

AS-based clusters significantly associated with immune features. (A) Bar plots of the relationship between AS-based clusters and infiltrate immune cell types. (B–D) Immune score, stromal score, and cytolytic activity between AS-based clusters. Data were analyzed using Kruskal-Wallis tests. (E) Expression level of immune inhibitory molecules among the four clusters.