Check for updates
Alternative Splicing Events and Splicing Factors Are Prognostic in Adrenocortical Carcinoma
Jian Lv1,21, Yuan He1t, Lili Li1* and Zhihua Wang1,2*
1 Central Laboratory, Renmin Hospital of Wuhan University, Wuhan, China, 2 Department of Cardiology, Renmin Hospital of Wuhan University, Wuhan, China
Alternative splicing is involved in the pathogenesis of human diseases, including cancer. Here, we investigated the potential application of alternative splicing events (ASEs) and splicing factors (SFs) in the prognosis of adrenocortical carcinoma (ACC). Transcriptome data from 79 ACC cases were downloaded from The Cancer Genome Atlas database, and percent spliced-in values of seven splicing types were downloaded from The Cancer Genome Atlas SpliceSeq database. By the univariate Cox regression analysis, 1,839 survival-related ASEs were identified. Prognostic indices based on seven types of survival-related ASEs were calculated by multivariate Cox regression analysis. Survival curves and receiver operating characteristic curves were used to assess the diagnostic value of the prognostic model. Independent prognosis analysis identified several ASEs (e.g., THNSL2| 54469| ME) that could be used as biomarkers to predict the prognosis of patients with ACC accurately. By analyzing the co-expression correlation between SFs and ASEs, 188 highly correlated interactions were established. From the protein interaction network, we finally screened six hub SFs, including YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4, whose expression levels were significantly related to the overall survival and prognosis of ACC. Our findings provide a reliable model for predicting the prognosis of ACC patients based on aberrant alternative splicing patterns.
Keywords: alternative splicing, Adrenocortical carcinoma, splicing factor, prognosis, hub gene
INTRODUCTION
Alternative splicing of pre-messenger RNA (mRNA) produces transcript isoforms for 95% of human genes, increases protein diversity, and provides functional diversity at various regulation level (Pan et al., 2008). There are seven types of splicing patterns, including alternate acceptor site (AA), alternate donor site (AD), alternate promoter (AP), alternate terminator (AT), exon skip (ES), mutually exclusive exons (ME), and retained intron (RI), as listed in The Cancer Genome Atlas (TCGA) SpliceSeq database (Ryan et al., 2016). Splicing factors (SFs) are involved in the removal of introns to create mature mRNAs, a process catalyzed by a large complex termed spliceosome (Yan et al., 2019). Alterations in SF expression lead to missplicing of key cancer-associated genes (Anczukow and Krainer, 2016; Lee and Abdel-Wahab, 2016). Aberrant alternative splicing events (ASEs) have been frequently observed in cancers and is recognized as an important signature for tumorigenesis and related pathologies, such as initiation and development of cancer (Oltean and Bates, 2014; Chen and Weiss, 2015; Lee and Abdel-Wahab, 2016),
frontiers in Genetics
OPEN ACCESS
Edited by:
Claes Wahlestedt,
Leonard M. Miller School of Medicine, University of Miami, United States
Reviewed by:
Argyris Papantonis,
University Medical Center Göttingen, Germany
Juw Won Park,
University of Louisville, United States
*Correspondence: Lili Li
Zhihua Wang
t These authors have contributed equally to this work
Specialty section:
This article was submitted to RNA,
a section of the journal Frontiers in Genetics
Received: 19 December 2019
Accepted: 23 July 2020
Published: 03 September 2020
Citation:
Lv J, He Y, Li L and Wang Z (2020) Alternative Splicing Events and Splicing Factors Are Prognostic in Adrenocortical Carcinoma. Front. Genet. 11:918.
doi: 10.3389/fgene.2020.00918
cancer metabolism (Kozlovski et al., 2017), cancer immunotherapy (Frankiw et al., 2019), cancer drug resistance (Siegfried and Karni, 2018), and so on.
Adrenocortical carcinoma (ACC) is a rare aggressive tumor with poor prognosis and less than 40% survival rate in 5 years (Jouinot and Bertherat, 2018; Crona and Beuschlein, 2019). Recent studies highlighted that specific molecular signatures could predict the survival and prognosis of ACC patients, which came from genomic approaches, including transcriptome, exome or whole-genome sequencing, chromosome alterations, methylome, and miRnome (Barreau et al., 2013; Patel et al., 2013; Assie et al., 2014; Szabo et al., 2014; Jouinot and Bertherat, 2018). However, only limited studies have focused on the potential roles of alternative splicing patterns in the pathogenesis of ACC. The present study aims to investigate the relationship between seven types of ASEs and SFs with the prognosis of ACC. Our findings
provide a new path to identify potential targets for the diagnosis and treatment of ACC.
RESULTS
Overview of Alternative Splicing Events in Adrenocortical Carcinoma
Transcriptome data from the TCGA database provide identity information for up to 56,754 transcripts, which represents a key resource for exploring ASEs in tumors concurrently deposited in the database. Individual ASE is assigned a unique annotation in the TCGA SpliceSeq database; for instance, in the term CIRBP| 46443| ES, CIRBP is the official gene symbol, 46443 is the unique ID number of a specific ASE, and ES represents the type of the alternative splicing pattern. Full datasets of splicing events of
A
C
1592
1500
Gene Intersections
1403
Gene Intersections
800
771
822
1000
600
623
400-
646
325
500
457
325 19
200-
250
2098789776456
136
S.
132
07
19 191909 00 7 87 16 75 14 71 69 69 68 66 66 65 63 60 56 50 47 A5 N/A 39 38 35 33 33 33 31 25 23 23 21 .
18
D4338342220
0
0
1010g 1
T
6
B
5
5
L
A
4
3
5
3
2
1
1
1
1
1
ME
.. .
RI
ME
.
·
. .
..
.
AA
..
AD
RI
0
.
AA
AD
·
AP
AP
AT
AT
ES
ES
B
10
SR-ASEs
D
CIRBP-46443-ES
n.s.
BLOC1S1-22229-AP
TRAFD1-24580-AD
w
CIRBP-46445-ES
P Value
-log10(P Value)
UNG-24278-AP
UNG-24277-AP
-
1.5e-07
10
EIF6-59080-RI
METTL 15-14782-AT
1.0e-07
CMC2-37702-AP
-
5.0e-08
CMC2-37703-AP
+
HM13-58892-ES
KLHL7-78950-AT
-log10(P Value)
KLHL7-78952-AT
7.0
2
TECPR2-29416-AT
7.5
TECPR2-29415-AT
DNAJC12-11898-AT
8.0
DNAJC12-11899-AT
8.5
0
DUT-30484-AP
-5
0
5
MPND-46796-ES
DUT-30485-AP
z-score
-6
-3
0
3
6
z-score
79 ACC cases have been deposited in TCGA SpliceSeq, which contains 34,419 ASEs corresponding to 8,994 parent genes. The splicing event data are presented as the percent-spliced-in (PSI), a score to evaluate the level of specific ASEs. The data were firstly filtered by deleting the unreliable ASEs with mean PSI of less than 0.05 and a standard deviation of PSI of less than 0.01. This generated a total of 22,521 ASEs from 8,040 parent genes. As summarized in the UpSet plot (Figure 1A), ES and AT were the top splicing types with high frequency in most genes, whereas ME was the least frequent splicing type.
Identification of Survival-Related Alternative Splicing Events
To explore the relationship between the alternative splicing pattern and the prognosis of ACC, we performed a univariate Cox regression analysis by comparing ASEs with the overall survival of ACC patients. A total of 1,839 ASEs were significantly associated with overall survival in ACC patients (P < 0.01), termed survival-related ASEs (SR-ASEs; Figure 1B). The upSet plot showed that ES, AT, and AP were the most frequent SR-ASEs, whereas only a small number of genes displayed a combination of multiple splicing forms (Figure 1C). Supplementary Figure 1 listed the top 20 most significant SR-ASEs in each splicing pattern according to their Z score and P-value. Taking all types together, splicing events within the parent genes CIRBP, BLOC1S1, TRAFD1, UNG, EIF6, METTL15, CMC2, HM13, KLHL7, TECPR2, DNAJC12, DUT, and MPND were highly related to the overall survival of ACC patients (Figure 1D). Survival curves with a cutoff at the median PSI value of the top four ASEs showed a remarkable difference (Figures 2A-D).
Independent Prognostic Predictors of Survival-Related Alternative Splicing Events in Adrenocortical Carcinoma
To identify independent prognostic indices and to explore the relationship between aberrant types of ASEs and ACC survival outcomes, we performed the multivariate Cox regression analysis for each splicing type to build a prognostic model. Lasso regression analysis was done to select the most significant SR-ASEs by avoiding overfitting. The median value of risk scores was then used to stratify the 79 ACC patients into low- and high-risk groups. Kaplan-Meier method was used to analyze the efficacy of the prognostic indices to predict the overall survival. The plotted survival curves and receiver operating characteristic (ROC) curves were shown in Figure 3. Significant differences in survival curves were observed in individual splicing type as well as all together (ALL), indicating that each alternative splicing type could be recognized as an independent prognostic indicator (Figure 3). The greatest difference in overall survival curves was observed in ES type (P = 6.684e-11), which is the most frequent splicing type among ASEs and SR-ASEs (Figures 1A,C). The area under the curve (AUC) of each ROC curve was more than 0.7, indicating the predictive efficiency of the eight models. We found that the AUC for AA type is 0.978, which is higher than all others (Figure 3), indicating that the prognostic indices based on AA type demonstrated the greatest efficacy in stratifying patients.
A
CIRBP.46443.ES
High
Low
100
AUC=0.627
Survival probability
1.00
80
0.75
Sensitivity (%)
60
0.50
0.25
P=2.894e-07
40
0.00
0
1
2
a
4
6
7
20
8
9
10 11 12
Time(years)
CIRBP|46443|ES
High LOW
36
39 1
15
7
A
a
O
0
2
9
0
2
3
2
9
4
5
6
7
8
9
10
11
12
100
80
60
40
20
0
Time(years)
Specificity (%)
B
BLOC1S1.22229.AP
High
Low
100
1.00
AUC=0.575
0.75
8
60
0.50 0.25 P=3.669e-05 0.00 5 0 1 2 3 4 6 36 10 21 S 3 4 5 8 7 8 1 Time(years) Survival probability . 4
Sensitivity (%)
Ov
7
8
9
10 11 12
2
Time(years)
High
2
6
0
BLOC1S1|22229|AP
Y
0
1
2
9
10
11
12
100
80
60
40
20
0
Specificity (%)
C
TRAFD1.24580.AD
High
Low
100
Survival probability
1.00-
AUC=0.635
0.75
30
Sensitivity (%)
0,50
6
0.25
P=2.587e-06
40
0.00
20
0
1
2
3
4
5
6
7
8
9
10 11 12
Time(years)
High LOW
16
3.
TRAFD1|24580|AD
A
A
G
3
8
5
3
6
6
0
100
80
60
40
20
0
0
1
2
3
4
5
6
7
8
9
10
11
12
Time(years)
Specificity (%)
D
CIRBP.46445.ES
High
Low
100
Survival probability
1.00-
AUC=0.635
0.75
80
Sensitivity (%)
0.50
5
0.25
P=5.392e-07
4
0.00
0
2
3
4
5
6
7
8
9
Time(years)
10 11
12
2
1
High LOW
38 39
800
4
18
โอ
3
3
3
1
1
0
CIRBP|46445|ES
0
1
2
3
4
Time(years)
5
6
7
B
9
10
11
12
100
80
60
40
20
0
Specificity (%)
E
THNSL2.54469.ME
High
Low
1.00
100
Survival probability
AUC=0.673
0.75
00
Sensitivity (%)
0,50
60
0.25
P=2.347e-04
40
0.00
0
1
2
3
4
5
6
7
8
9
10 11 12
2
Time(years)
High LOW
37
Não
15
10
1
2
10
3
6
0
THNSL2.54469.ME
0
1
2
3
4
5
5
1
8
9
10
11
12
100
80
60
40
20
0
Time(years)
Specificity (%)
FIGURE 2 | Kaplan-Meier survival curves and ROC curves for five alternative splicing type: CIRBP| 46443| ES (A), BLOC1S1| 22229| AP (B), TRAFD1| 24580| AD (C), CIRBP| 46445| ES (D), and THNSL2| 54469| ME (E). Seventy-nine ACC patients were divided into high- and low-risk groups based on the median of PSI. Red line indicates the high PSI score group, and the blue line indicates the low PSI score group.
We also performed multivariate Cox regression analysis to evaluate the effect of other clinical parameters, including sex, tumor stage, T and N stages in the tumor-node-metastasis (TNM) classification in Table 1. The significant hazard ratios (HRs) for T stage and tumor stage were 3.349 (95% confidence interval [CI]: 2.075-5.403) and 2.886 (CI: 1.819- 4.579), respectively (Table 1). HRs of AA, AD, AP, AT, ES, RI, and
Survival probability >
Risk - High risk - Low risk
Ture positive rate
1.0
Survival probability w
+
+
1.00
0.8
AUC=0.978
Risk
High risk
Low risk
1.0
1.00
Ture positive rate
0.8
AUC=0.994
+
+
0.75
0.75
+
0.50
0.50
+
0.25
0.4
0.25
0.4
+
+
0.00
p=1.016e-09
+
0.2
0.00
p=1.062e-08
0.2
0
1
2
3
4
5
6
7
8
9
101112
0
1
2
3
4
5
6
7
8
9
101112
Time(years)
AA
0.0
0.0 0.2 0.4 0.6 0.8 1.0
Time(years)
AD
0.0
0.0
False postive rate
0.2 0.4 0.6 0.8 1.0
False postive rate
C
Risk - High risk = Low risk
Survival probability O
1.0
Risk = High risk = Low risk
Survival probability
1.00
Ture positive rate
0.8
AUC=0.972
1.00
Ture positive rate
1.0
0.75
0.8
AUC=0.764
0.75
0.50
0.6
0.50
0.6
+
+
0.25
0.4
0.25
0.4
+
p=3.4e-09
+
p=2.862e-
07
+
0.00
0 1 2 3 4 5 6 7 8 9 101112
0.2
0.00
0 1 2 3 4 5 6 7 8 9 101112
0.2
Time(years)
AP
0.0
Time(years)
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
False postive rate
AT
False postive rate
E
F
Survival probability
Risk
-
High risk - L
Low risk
Ture positive rate
1.0
AUC=0.869
Survival probability
Risk
-
High risk
+
Low risk
1.0
1.00
0.8
1.00
Ture positive rate
0.8
AUC=0.767
+
0.75
+
0.75
+
0.50
0.6
0.50
0.6
+
0.25
0.4
0.25
0.4
+
0.00
p=6.684e-
+
11
0.2
0.00
p=1.599e-04
0 1 2 3 4 5 6 7 8 9 101112
0.0
0 1 2 3 4 5 6 7 8 9 101112
0.2
Time(years)
0.0
ES
Time(years)
0.0
0.2
0.4
0.6
0.8
1.0
False postive rate
ME
0.0
0.2
0.4
1.0
False postive rate
0.6
0.8
G
H
Risk
High risk = Low risk
= High risk
+
Low risk
Survival probability
1.00
Ture positive rate
1.0
Survival probability
Risk
1.0
0.8
AUC=0.928
1.00
Ture positive rate
0.75
0.8
AUC=0.892
0.75
+
+
+
0.50
4 0.6
0.6
+
0.50
+
0.25
0.4
0.25
0.4
+
+
+
+
+
+
0.00
p=1.547e-08
0.2
0.00
p=2.39e-08
0.2
0
1 2 3 4 5 6 7 8 9 101
101112
0
2
3
4
5
6
8 9
101112
Time(years)
RI
0.0
Time(years)
0.0
0.0
0.2
0.4
0
.6
0.8
1.0
ALL
0.0
0.2
0.4
0.6
0.8
1.0
False postive rate
False postive rate
FIGURE 3 | Kaplan-Meier survival curves and ROC curves for each alternative splicing type: AA (A), AD (B), AP (C), AT (D), ES (E), ME (F), RI (G), and ALL (H) model. Seventy-nine ACC patients were divided into high- and low-risk groups based on the median of prognostic indices. Red line indicates the high-risk group, and blue line indicates the low-risk group. Prediction time of ROC curve is 1 year. AA, alternate acceptor site; AD, alternate donor site; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.
ALL models are also shown in Table 1. Among the seven types, the HR for ME type was 1.270 (CI: 1.093-1.476), the highest HR value (Table 1).
Hazard ratios in each of AA, AD, AP, AT, ES, RI, and ALL models with clinical parameters are shown in Figure 4. The HRs for all the AS type were under a significant level (P-value < 0.05), and ME type has the highest HR value of 1.681 (CI: 1.344-2.103) (Figure 4F). Moreover, THNSL2| 54469| ME ranks that most significant event in the ME type (Supplementary Figure 1F). THNSL2| 54469| ME was the top
significant SR-ASE in ME model to predict the prognostic status of ACC cases. Therefore, the most significant SR-ASEs in ME pattern is THNSL2| 54469| ME (Figure 2E). To test the accuracy of ME model of THNSL2| 54469| ME, the survival state of ACC patients could be significantly classified into high- and low-PSI groups according to the median value of PSI scores of THNSL2| 54469| ME (Figure 2E). These results indicate that THNSL2| 54469| ME could be the best independent prognostic indicator to predict the prognosis of ACC cases.
| Term | Hazard ratio | HR.95L# | HR.95H# | P-value |
|---|---|---|---|---|
| Sex | 0.963 | 0.425 | 2.181 | 0.928 |
| Tumor stage | 2.886 | 1.819 | 4.579 | 6.75E-06 |
| T | 3.349 | 2.075 | 5.403 | 7.36E-07 |
| N | 2.088 | 0.781 | 5.578 | 0.142 |
| AA | 1.003 | 1.001 | 1.004 | 0.001 |
| AD | 1.007 | 1.003 | 1.011 | 6.17E-04 |
| AP | 1.018 | 1.010 | 1.026 | 3.29E-06 |
| AT | 1.095 | 1.060 | 1.130 | 2.29E-08 |
| ES | 1.009 | 1.005 | 1.012 | 8.51E-07 |
| ME | 1.270 | 1.093 | 1.476 | 0.002 |
| RI | 1.009 | 1.005 | 1.012 | 3.18E-06 |
| ALL | 1.014 | 1.007 | 1.021 | 4.89E-05 |
HR.95L and HR.95H: 95% confidence interval.
Also, the AUCs for 1, 3, and 5 years of ME model are predicted to be 0.767 (Figure 3F), 0.792 (Figure 5A), and 0.829 (Figure 5B), respectively. The number of deaths increased in the high-risk group, with short survival time (Figures 5C,D). In summary, those results indicated that the constructed ME model had great efficacy in predicting the prognosis of ACC patients.
Enrichment Analysis of Survival-Related Alternative Splicing Events
To explore the biological processes and signal pathways related to alternative splicing in the progression of ACC, enrichment analysis of the SR-ASEs parent genes was performed by gene ontology and pathway analysis in Metascape. The most significantly enriched terms were regulation of mitotic cell cycle, cofactor metabolic process, covalent chromatin modification, cell cycle G2/M phase transition, and mitochondrion organization (Figure 6), pathways that are all involved in tumorigenesis.
Correlation Network of Survival-Related Alternative Splicing Events and Splicing Factors Expression
Alteration of ASEs is largely attributable to changes in SF expression. We then extracted the expression data of 404 SFs, summarized by a previous study (Seiler et al., 2018), from the transcriptome data of the 79 ACC patients. Principal component analysis (PCA) plots could discriminate the distribution pattern of SF expression levels according to survival state, tumor stage, and TNM classification (Supplementary Figure 2), suggesting an impact of altered SF expression on the ACC outcome.
As our current knowledge is incapable of dissecting the sequence specificity for each SF, we could not establish a direct network for SF-regulated ASEs. Thus, we analyzed the relationship between SFs and SR-ASEs based on their co-expression patterns using the Spearman method, which has been widely used in alternative splicing studies (He et al., 2018; Xiong et al., 2018; Zhang et al., 2020). A total of 188 highly correlated interactions between SFs and SR-ASEs were detected with a correlation coefficient larger than 0.65 (Figure 7A and
Table 2). Hub SFs with no less than five SR-ASE connections were HSPA1B, YBX1, SRPK1, SART1, PRCC, ILF2, SNRPG, SNRPE, SF3B4, BUD13, INTS4, and CLK2 (Table 3). Among these SFs, interestingly, SNRPE was exclusively correlated with AT-type ASEs, suggesting a specific role in regulating terminal exon selection (Figure 7A). Moreover, HSPA1B was exclusively correlated with detrimental ASEs showing negative correlations with overall survival, indicating a potential causal role of HSPA1B in the progression of ACC (Figure 7A).
The expression data showed that SF YBX1 was positively correlated with the PSI values of C16orf13| 32916| ES, COX4I1| 37906| RI, FHAD1| 747| AT, MYL6| 22381| AA, PI4K2A| 12728| AP, ASH2L| 83369| AP, PTAR1| 86546| AT, and IRF3| 51012| ES and was negatively correlated with the PSI values of FHAD1| 749| AT, AIG1| 77970| AT, PI4K2A| 12729| AP, STARD3NL| 79286| ES, and ASH2L| 83368| AP; SF SNRFE was positively correlated with the PSI values of VWA8| 25742| AT, MSI2| 42617| AT, and PELI3| 17034| AT and was negatively correlated with the PSI values of VWA8| 25741| AT, MSI2| 42616| AT, and PELI3| 17033| AT (Figure 7A). We then constructed the Kaplan-Meier survival curves for each ASEs related with YBX1 and SNRFE and found that high levels of PSI values for all the negatively correlated ASEs have a better prognostic except FHAD1| 747| AT, whereas the low levels of PSI values for all the positively correlated ASEs have a better prognostic (Figure 8). Interestingly, six pairs of ASEs from three parent genes were observed for these two SFs in the network. For example, ASH2L| 83368| AP and ASH2L| 83369| AP are two alternative splicing sites for ASH2L exon 1 selection. PI4K2A| 12728| AP and PI4K2A| 12729| AP are two alternative splicing sites for ASH2L exon 1 selection (Figures 8A,B). Our results indicate that promoter selections of ASH2L and PI4K2A are important for tumorigenesis of ACC. Similar results were also observed for the SF SNRFE (Figures 8C,D). These results suggest that terminator selections of MSI2 and PELI3 and VWA8 are important for the progression of ACC development (Figures 8C,D).
Considering that SFs would influence each other when regulating the work mode of the spliceosome, we constructed the protein-protein internecion network to illustrate the interactions among SFs using the Search Tool for the Retrieval of Interacting Genes database (Figure 7B). A total of 66 nodes and 485 edges were revealed, with the interaction score of 0.900 (Figure 7B). Module analysis was done to identify hub genes. Two modules were further identified by the app MCODE in Cytoscape. Module 1 has 30 genes with a score of 27, and module 2 has four genes with a score of 8 (Table 3). By combining the ASE correlation results and protein interaction networks, we finally identified six hub SFs, including YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4 (Table 4). Overall, these results indicate that these six hub genes may play an important role in the progression of ACC by regulating the pattern of SR-ASEs.
Analysis of Hub Splicing Factors
To confirm whether these six genes, YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4, are high-risk factors, the overall survival and expression level of these genes were further investigated. The results showed that the survival rate of patients
A
P Value
Hazard Ratio
B
P Value
Hazard Ratio
Gender
0.977
1.013(0.432-2.374)
Gender
0.683
0.826(0.330-2.068)
Stage
0.984
1.010(0.373-2.733)
stage
0.767
1.167(0.420-3.241)
T
0.016
3.046(1.229-7.550)
T
0.007
3.495(1.397-8.744)
N
0.192
2.318(0.656-8.193)
N
0.128
2.712(0.751-9.793)
AA
0.003
1.002(1.001-1.003)
AD
<0.001
1.008(1.004-1.012)
0.50
1.0
2.0
4.0
8.0
0.50
1.0
2.0
4.0
8.0
C
P Value
Hazard Ratio
D
P Value
Hazard Ratio
Gender
0.706
1.181(0.498-2.802)
Gender
0.966
1.020(0.421-2.468)
Sage
0.895
1.078(0.354-3.283)
Stage
0.529
0.706(0.238-2.088)
T
0.057
2.784(0.970-7.993)
T
0.006
4.076(1.507-11.026)
N
0.918
1.082(0.244-4.788)
N
0.076
3.151(0.887-11.194)
AP
0.001
1.014(1.005-1.022)
AT
<0.001
1.079(1.044-1.116)
0.25
0.50
1.0
2.0
4.0
8.0
0.25
0.50
1.0
2.0
4.0
8.0
E
P Value
Hazard Ratio
F
P Value
Hazard Ratio
Gender
0.796
0.892(0.375-2.121)
Gender
0.757
1.146(0.484-2.716)
Stage
0.731
1.192(0.437-3.249)
Stage
0.764
1.167(0.427-3.190)
T
0.041
2.704(1.041-7.022)
T
0.007
3.709(1.421-9.681)
N
0.604
1.373(0.414-4.556)
N
0.198
2.160(0.668-6.982)
ES
<0.001
1.006(1.003-1.009)
ME
<0.001
1.681(1.344-2.103)
0.50
1.0
2.0
4.0
0.50
1.0
2.0
4.0
8.0
G
P Value
Hazard Ratio
H
P Value
Hazard Ratio
Gender
0.453
0.699(0.275-1.780)
Gender
0.947
0.970(0.403-2.339)
Stage
0.467
1.426(0.548-3.714)
Stage
0.806
0.877(0.306-2.508)
T
0.045
2.486(1.022-6.048)
T
0.005
3.881(1.495-10.073)
N
0.396
1.717(0.493-5.974)
N
0.141
2.592(0.730-9.198)
RI
<0.001
1.008(1.004-1.011)
ALL
0.006
1.010(1.003-1.018)
0.25
0.50
1.0
2.0
4.0
0.50
1.0
2.0
4.0
8.0
FIGURE 4 | Multivariate Cox regression analysis for AA (A), AD (B), AP (C), AT (D), ES (E), ME (F), RI (G), and ALL (H) model. Hazard ratio is shown as hazard ratio (95% confidence interval). AA, alternate acceptor site; AD, alternate donor site; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.
A
B
1.0
AUC = 0.792
1.0
AUC = 0.829
True positive rate
0.8
True positive rate
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0 0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0
False positive rate
0.0 0.2 0.4 0.6 0.
0.8
1.0
False positive rate
C
..
8
Risk score
6
4
2
0
20
Patients (increasing risk socre)
40
60
80
Survival time (years)
10 12
8
6
4
2
0
0
20
40
60
80
Patients (increasing risk socre)
FIGURE 5 | Overview of ME model. ROC curves of ME type with prediction time of 3 years (A) and 5 years (B). (C) Risk factors of ME model of 79 ACC cases. Patients were divided into low- and high-risk groups based on the median of prognostic indices. Red indicates high-risk group, and green indicates low-risk group. (D) Survival time and survival status of ME model of 79 ACC cases. Red indicates the deaths of the cases. Green indicates cases who are alive.
with high expression of hub SFs was significantly lower than that in the low expression group (P < 0.001) (Figure 9), and the results were consistent with those from the Gene Expression Profiling Interactive Analysis (GEPIA) database (Supplementary Figure 3).
The correlation between the expression levels of six SFs and clinical information, including clinical stages, tumor stage, and lymph nodes stage, was further analyzed. We found that the expression levels of SFs increased along with the clinical stages. The expression level of six hub SFs in tumor patients has an increasing trend with the tumor progress stage. Among six SFs, the level of SF3B4 in ACC cases with the tumor stage III and the level of PRCC and SF3B4 in ACC cases with stage IV were statistically different from those in stage I (Figure 9), which is consistent with the results from the GEPIA database (Supplementary Figure 3). The levels of the SFs in ACC cases with lymph node stage also increased, and the expression of YBX1 was statistically changed. The results indicate that the expression
levels of these SFs are closely related to the survival time and prognosis of patients with ACC (Figure 9).
Considering the lack of normal control tissue for ACC in the TCGA database, two microarray datasets from Gene Expression Omnibus (GEO) database were used to compare the expression of SFs between tumor tissues and normal controls (Figure 10). The mRNA levels of YBX1 and SNRPE were increased in tumor samples in the two microarrays, although the expression levels only showed statistically significant differences in GSE19750 (Figure 10). Combining the mRNA level of YBX1 and SNRPE in different tissues and the effects of YBX1 and SNRPE on the survival curve of ACC cases (Figure 7), we could conclude that YBX1 and SNRPE could exert positive regulation in the progression of ACC.
We next performed a Cox regression analysis to evaluate the prognostic value of these six hub SFs and other clinical parameters. Results showed that the HRs of these six genes ranged from 1.003 to 1.156; although the HR is lower than the T and N
0
2
4
6
8
10
12
14
-log10 (P)
GO:0007346: regulation of mitotic cell cycle
GO:0051186: cofactor metabolic process
GO:0016569: covalent chromatin modification
GO:0044839: cell cycle G2/M phase transition
GO:0007005: mitochondrion organization
GO:0140014: mitotic nuclear division
GO:0007017: microtubule-based process
GO:0006605: protein targeting
R-HSA-1428517: The citric acid (TCA) cycle and respiratory electron transport
GO:0006402: mRNA catabolic process
R-HSA-8953854: Metabolism of RNA
GO:0051640: organelle localization
R-HSA-3700989: Transcriptional Regulation by TP53
GO:0033044: regulation of chromosome organization
GO:0006914: autophagy
R-HSA-9609507: Protein localization
GO:0045787: positive regulation of cell cycle
R-HSA-556833: Metabolism of lipids
GO:0006412: translation
GO:0006631: fatty acid metabolic process
A
B
Val
-
ISY1
₥
un
GEMIN2
INTS4
TLR
4
SNRPESNRPG
SRRT
A
PRPF39
PRPF38A
SNRPA1
SNRNP40
DNAJC8
DDX39A
.-
SART1
SNRPD1
SF3B6
LY RIS
1
LUC7L
SNRPC
LSM7
RBM8A
PRCC
THOC3
CLK
4
ul
PRPF40A
SF3A3
-
CCAR1
SRSF2
MAGOH
35
TCERG1
SNRNP27
ZC3H11A
SF3B1
SUGP1
YBX1
ALYREF
5
NUO
SF3B4
RBM5
SRSF3
THOC5
un
BUD13
DHX9
2
HNRNPM
HNRNPH1
SRPK1
CURT
TA
TIA1
ILF3
ILF2
PICA
TIAL1
CHDIRE
×
Mª
JUP
A
6
*-
FIGURE 7 | Correlation network of SR-ASEs and SFs. (A) Co-expression networks between parent genes of SR-ASEs and SFs. Yellow square indicates SFs, red circle indicates the poor prognostic SR-ASEs, and blue circle indicates the better prognostic SR-ASEs. Red line indicates that SFs were positively correlated with the SR-ASE, and blue line indicates SFs were negatively correlated with the SR-ASE. (B) Protein-protein interaction network of SFs.
stages, it statistically significant for all the six SFs (Supplementary Figure 4). The AUC of each ROC curve is higher than 70% (Supplementary Figure 4), indicating that all of these six SFs could be used to classify the stages of ACC.
DISCUSSION
Adrenocortical carcinoma is a rare malignancy tumor with a poor prognosis. Currently, surgery is the only available curative treatment option (Crona and Beuschlein, 2019). Recent studies (Barreau et al., 2013; Patel et al., 2013;
Assie et al., 2014; Szabo et al., 2014; Jouinot and Bertherat, 2018) highlighted that genomic approaches derived from TCGA database and GEO datasets could provide specific molecular signatures for the diagnosis and prognosis of ACC. At present, few studies focused on the different isoforms of alternative splicing in ACC (Bie et al., 2019). The present study is to explore the aberrant of ASEs and hub SFs to develop novel diagnostic and prognostic markers for ACC.
Seven types of ASEs were investigated in this study. ES type was the top splicing type with high-frequency ASEs and SR-ASEs
| Splicing factor | ASEs# | Correlation | P-value | Regulation |
|---|---|---|---|---|
| HSPA1B | KCNJ5-19431-AP | -0.84181986 | 8.92E-22 | Negative |
| HSPA1B | UBE3A-29716-AP | -0.82743598 | 1.77E-20 | Negative |
| MAGOH | ZNF131-71936-ES | -0.81209442 | 3.20E-19 | Negative |
| HSPA1B | IKBKB-83691-ES | -0.80904662 | 5.51E-19 | Negative |
| HSPA1B | CYTH1-43892-ES | -0.80356862 | 1.43E-18 | Negative |
| HSPA1B | EYS-76614-AT | -0.79837952 | 3.43E-18 | Negative |
| PRPF38A | ZNF131-71936-ES | -0.78822643 | 1.77E-17 | Negative |
| SF3B4 | ZNF131-71936-ES | -0.788164 | 1.79E-17 | Negative |
| HSPA1B | GNG7-46614-AP | -0.77278859 | 1.83E-16 | Negative |
| PRPF38A | ZNF131-71932-ES | -0.77084258 | 2.42E-16 | Negative |
| BUD13 | TNFRSF12A-33348-ES | -0.77031719 | 2.61E-16 | Negative |
| PRPF38A | F6-23309-ES | -0.75603278 | 1.88E-15 | Negative |
| PRPF38A | HDDC3-32524-AA | 0.751788406 | 3.30E-15 | Positive |
| SNRPA1 | ZNF131-71936-ES | -0.74748153 | 5.77E-15 | Negative |
| HSPA1B | DAGLB-78732-ES | -0.74496838 | 7.96E-15 | Negative |
| SRPK1 | NDUFA12-23740-ES | -0.74060614 | 1.38E-14 | Negative |
| SNRPG | KIF20B-12497-AT | -0.73362818 | 3.23E-14 | Negative |
| SNRPG | KIF20B-12498-AT | 0.733628184 | 3.23E-14 | Positive |
| PRCC | ARHGEF28-72492-AT | -0.73325406 | 3.38E-14 | Negative |
| PRCC | ARHGEF28-72493-AT | 0.733254055 | 3.38E-14 | Positive |
# In terms of KCNJ5-19431-AP, KCNJ5 is an official gene symbol, 19431 is the number of alternative splicing events, and AP is the type of AS pattern.
| Module | Score | Gene symbol |
|---|---|---|
| Module 1 | 27 | SART1, SNRPA1, SNRPG, SF3B1, SF3B6, SNRPD1, RBM8A, MAGOH, SNRPE, SRSF3, LSM7, SF3B4, SF3A3, SRSF2, HNRNPM, SNRPC, SNRNP40, HNRNPH1, SRRT, ALYREF, PRPF38A, PRPF40A, DHX9, YBX1, SNRNP27, RBM5, DNAJC8, CCAR1, PRCC, SUGP1 |
| Module 2 | 8 | THOC5, THOC3, DDX39A, ZC3H11A |
(Figures 1A,C), indicating the ES is the dominant alternative splicing type in ACC.
In this study, we identified 1,839 SR-ASEs by univariate Cox regression analysis. Interestingly, different ASEs in the same gene could exert opposite functions in the overall survival of ACC, indicating that the parent genes of these ASEs may play an important role in ACC development (Figures 1D, 8). Because ME type has the highest HR value (Figure 4F) THNSL2| 54469| ME ranks the most significant event in the ME type (Supplementary Figure 1F). Therefore, THNSL2| 54469| ME could be used as an independent prognostic indicator to predict the prognosis of ACC cases (Figure 2E). THNSL2 encodes a threonine synthase-like protein and has multiple transcript variants. The function of THNSL2 and the alternative splicing of THNSL2| 54469| ME could be further investigated. Enrichment analysis revealed several important pathways, including regulation of the mitotic cell cycle and cell cycle G2/M phase transition,
which could impact the occurrence and development of ACC (Figure 6).
Genes involved in the same biological process or signaling pathway are usually co-regulated in the disease context. Co-expression network analysis has been widely used to dissect the functional gene panels in large datasets, including alternative splicing studies (He et al., 2018; Xiong et al., 2018; Zhang et al., 2020). One hundred eighty-eight highly correlated interactions between SFs and SR-ASEs were identified (Figure 7). SR-ASEs that were positively correlated with SFs have a poor prognostic value, whereas SR-ASEs that were positively correlated with SFs have a poor prognostic value (Figures 7, 8). It provided a new insight for the molecular mechanism of alternative splicing in ACC.
The number of pairs of ASEs from one parent gene was observed for SFs in the network (Figure 7), and six pairs of ASEs related with YBX1 and SNRFE were illustrated in detail (Figure 8). We observed that the pairs of ASEs conferred the opposite function for the ACC progress shown in the overall survival curve (Figure 8), indicating that specific exons were important, such as the alternative promoters’ selection of ASH2L, and alternative terminator selection of MSI2. Moreover, among these six pairs of ASEs, only alternative promoters of ASH2L have been reported in embryonal carcinoma (Alagaratnam et al., 2013).
A previous study on ASEs in endometrial cancer also identified YBX1 as the hub SF. One pair of ASEs that correlated with YBX1 was identified in that study: DNAH9-AT-39292 and DNAH9-AT-39293 (Wang et al., 2019). We could conclude that, firstly, it seems YBX1 more specifically regulates the first exon, as the splicing type is AT type in both studies;
A
YBX1_negatively correlate
ASH2L.83369.AP
High
Low
FHAD1.747.AT
High
Low
PI4K2A.12728.AP
High
Low
Survival probability
1.00-
Survival probability
1.00
Survival probability
1.00 -
0.75
0.75
0.75
0.50
0.50
0.50
0.25
P=4.147e-03
0.25
P=2.285e-01
0.25
P=1.883e-03
0.00
0.00
0
1
2
3
4
5
6
7
8
9
10 11
12
0
1
2
3
4
5
6
7
8
9
Time(years)
10 11 12
0.00
0
1
2
3
4
5
6
7
8
9
Time(years)
10 1
Time(years)
12
C16orf13.32916.ES
High
Low
COX411.37906.RI
High
Low
PTAR1.86546.AT
High
Low
1.00-
Survival probability
Survival probability
1.00
Survival probability
1.00-
0.75
0.75
0.75
0.50
0.50
0.50
0.25
P=1.066e-01
0.25
P=1.046e-01
0.25
P=6.527e-02
0.00
0.00
0.00
0
1
2
3
4
Time(years)
5
6
7
8
9
10 11 12
0
1
2
3
4
5
6
7
8
9
Time(years)
10 11 12
0
1
2
3
4
5
6
7
8
9
Time(years)
10 1
12
B
YBX1_positively correlate
ASH2L.83368.AP
High
Low
FHAD1.749.AT
High
Low
PI4K2A.12729.AP
High
Low
Survival probability
1.00-
1.00
Survival probability
Survival probability
1.00
0.75
0.75
0.75
0.50
0.50
0.50
0.25
P=1.268e-02
0.25
P=1.437e-01
0.25
P=6.481e-04
0.00
0.00
0.00
0
1
2
3
4
5
6
7
8
9
10 11 12
0
1
2
3
4
5
6
7
8
9
10
0 11
12
0
1
2
3
4
5
6
7
8
9
Time(years)
Time(years)
Time(years)
10 11 12
STARD3NL.79286.ES
-
High
Low
AIG1.77970.AT
High
Low
1.00
1.00
Survival probability
Survival probability
0.75
0.75
0.50
0.50
0.25
P=7.559e-02
0.25
P=1.388e-05
0.00
0.00
0
1
2
3
4
6
7
Time(years)
5
8
10 11 12
0
1
2
3
4
5
6
7
8
9
C
Time(years)
10 1
12
SNRFE_negatively correlate
MSI2.42617.AT
High
Low
PELI3.17034.AT
High
Low
VWA8.25742.AT
High
Low
Survival probability
1.00
Survival probability
1.00
Survival probability
1.00-
0.75
0.75
0.75
.50
0.50
0.50
0.25
P=8.748e-04
0.25
P=2.455e-03
0.25
P=5.11e-02
0.00
0.00
0.00
0
1
2
3
4
5
6
7
Time(years)
8
9
10 11 12
0
1
2
3
4
5
6
7
8
9
Time(years)
10 11 12
0
1
2
3
4
Time(years)
5
6
7
8
9
10 11 12
D
SNRFE_positively correlate
MSI2.42616.AT
High
Low
PELI3.17033.AT
High
Low
VWA8.25741.AT
High
Low
Survival probability
1.00
Survival probability
1.00
Survival probability
1.00
0.75
0.75
0.75
0.50
0.50
0.50
0.25
P=3.054e-04
0.25
P=1.765e-03
0.25
P=4.078e-02
0.00
0.00
0.00
0
1
2
3
4
5
6
7
8
Time(years)
9
10 11 12
0
1
2
3
4
5
6
7
8
9
10 11 12
0
1
2
3
4
5
6
7
8
9
10 11 12
Time(years)
Time(years)
| Splicing factor | Count |
|---|---|
| HSPA1B | 14 |
| YBX1 | 13 |
| SRPK1 | 7 |
| SART1 | 7 |
| PRCC | 7 |
| ILF2 | 7 |
| SNRPG | 6 |
| SNRPE | 6 |
| SF3B4 | 6 |
| BUD13 | 6 |
| INTS4 | 5 |
| CLK2 | 5 |
secondly, YBX1, as the hub SF, might regulate specific genes in different cancer types.
Six hub SFs, including YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4, were identified in this study. We found that the expression level of hub SFs was negatively correlated with the survival time and survival state of ACC patients (Figure 9). There was still some evidence that the expression level of hub SFs in tumor tissues was higher than that in normal tissues (Figure 10). The altered expression level of hub SFs has been reported in multiple types of cancer, such as YBX1 (Rossner et al., 2016; Chen et al., 2019), SART1 (Ishida et al., 2000; Sasatomi et al., 2000; Yutani et al., 2001), SNRPE (Tamura et al., 2007; Jia et al., 2011), and SF3B4 (Liu et al., 2018). Other studies have shown that SNRPE (Quidville et al., 2013) and SF3B4 (Shen and Nam, 2018) could develop a new therapeutic agent in cancer. Further study found that the inactivation of SF3B4 inhibited liver tumorigenesis in vitro and in vivo (Shen et al., 2018). In terms of molecular mechanism, SF3B4 triggers SF3b complex to splice tumor suppressor KLF4 transcript to non-functional skipped exon transcripts, downregulates the transcriptional activity of p27Kip1, and upregulates the transcriptional activation of Slug genes (Shen et al., 2018). However, the functions of hub SFs in ACC development need to be further studied.
Our present study performed a bioinformatic analysis of SR-ASEs and hub SFs and provided insight into the function of aberrant ASEs in ACC development and progression. SR-ASEs and hub SFs identified in this study could be potential targets for the diagnosis of ACC patients.
MATERIALS AND METHODS
Data Collection
The transcriptome data of 79 ACC cases and the corresponding clinical information, including survival time, survival status, sex, tumor stage, T stage, and lymph node metastasis, were downloaded from TCGA1. Seven types of ASEs of 79 ACC cases
were downloaded from the TCGA SpliceSeq database2 (Ryan et al., 2016). ASEs were quantified using the PSI, rating from 0 to 1. Splicing event data of 79 ACC cases included 34,419 ASEs, corresponding to 8,994 genes. To better track the AS events, the name of the ASEs contains three parts: the official gene symbol, a unique splicing event ID number, and splicing type. The number of ASEs for different types of ASEs are shown in the UpSet plot (Figure 1A).
Screening of Survival-Related Alternative Splicing Events
The missing value of PSI in seven AS event types was supplemented by impute.knn function using nearest neighbor averaging method with k = 10, rowmax = 0.5, colmax = 0.8, and other default setting in R. Then, the PSI data of 79 ACC cases were filtered by deleting the ASEs with mean PSI of less than 0.05 and standard deviation of PSI of less than 0.01. The filtered PSI data, including 22,521 ASEs from 8,040 genes in 79 ACC cases, were used for visualization and subsequent analysis. The survival time and survival status of 79 ACC cases were integrated with the filtered PSI data. Then, SR-ASEs were selected using univariate Cox regression analysis with a threshold set to a P-value < 0.01 and were visualized by volcano plot and bubble plot.
Construction of Prognosis Model of Survival-Related Alternative Splicing Events
Multivariate Cox regression analysis was performed to calculate the prognostic indices of each for each type of splicing pattern. To prevent over-fitting of the model, seven different types of SR-ASEs and ALL SR-ASEs were analyzed by lasso regression analysis, making the Lambda values at smaller level. The SR-ASEs with high correlation have been removed to guarantee the accuracy of the model. Risk factors were calculated using the following formula BSR-ASE1 x PSISR-ASE1 + BSR- ASE2 × PSISR-ASE2 + … + BSR-ASEn × PSISR-ASEn, where ß corresponded to the regression coefficient. The samples were stratified to high- and low-risk groups based on the median of risk scores. The survival curve and ROC curve of the 79 ACC cases were used to evaluate the accuracy of the model. Univariate and multivariate Cox regression analyses were used to determine whether the prognostic indices could be used as an independent prognostic factor for predicting the prognosis of ACC cases.
Enrichment Analysis of Survival-Related Alternative Splicing Events
To illustrate the biological functions and pathway associated with SR-ASEs of ACC cases, gene enrichment analysis for parent genes of SR-ASEs was performed in the online database Metascape3 (Zhou et al., 2019) with the default setting. Metascape is an online analytical tool for pooled gene annotation, enrichment analysis, and protein interaction analysis.
1https://tcga-data.nci.nih.gov/tcga/
2https://bioinformatics.mdanderson.org/TCGASpliceSeq/
A
Survival probability
YBX1
High
Low
12
*
12
n.s.
**
*
*
1.00
n.s.
n.s
10
n.s.
n.s.
10
n.s.
0.75
n.s.
8
YBX1
n.s
n.s.
..
+.
0.50
×
+·
YBX1
8
”
$
m
8
*1
t
0.25
AT
6
-
+
p=1.138e-05
6
6
9
G
0.00
·
0
1
2
3
4
5
6
7
8
9
101112
4
Time(years)
4
..
..
1
Stage
Stage 2
Stage 3
4
Stage 4
T1 T2 T3 T4
NO
N1
B
Survival probability
SART1 + High
Low
n.s.
n.s.
*
*
1.00
6
**
6
n.s.
n.s.
n.s.
n.s.
5
0.75
n.s.
n.s.
SART1
n.s.
n.s.
.
+#
+#
SAR
SART1
0.50
4
=
R
0.25
p=1.39e-04
4
4
+
·
0.00
-
+1
3
0 1 2 3 4 5 6 7 8 9 101112 Time(years)
3
3
Stage 1
Stage 2
Stage 3
A
Stage
T1 T2 T3 T4
NO
N1
C
Survival probability
PRCC
High
Low
8-
n.s. **
8-
n.s.
*
6
1.00
n.s
n.s.
n.s.
*
n.s.
0.75
6
n.s.
6
n.s.
PRCC
n.s
+·
0
n.s.
+·
PRCC
0.50
2
C
1
0.25
4
4
.
A
0.00
p=2.073e-06
2
0 1 2 3 4 5 6 7 8 9 101112
2
2
Time(years)
…
Stage
1
Stage 2
Stage 3 Stage
..
A
T1 T2 T3 T4
NO
N1
D
Survival probability
SNRPG - High
Low
8
n.s.
8-
n.s
1.00
**
*
n.s.
6
n.s.
n.s.
n.s.
0.75
SNRPG
n.s.
n.s.
6
n.s.
SNRPG
n.s.
SNRPG
0.50
+;
5
J
Y
4
0.25
0.00
p=7.44e-06
4
4
”
3
0 1 2 3 4 5 6 7 8 9 101112 Time(years)
2
..
..
Stage 1
Stage 2
Stage 3
Stage 4
T1
T2 T3 T4
NO
N1
E
Survival probability
SNRPE
High
Low
n.s.
n.s.
7
1.00
8
**
n.s.
n.s.
8
n.s.
n.s
n.s.
n.s.
6
0.75
SNRPE
n.s.
n.s.
n.s.
6
SNRPE
n.s.
SNRPE
5
0.50
.
7
#
4
0.25
C
0.00
p=5.548e-04
4
4
.
3
0 1 2 3 4 5 6 7 8 9 101112
Time(years)
2
2
2
..
Stage 1
Stage 2
Stage 3
Stage 4
T1 T2 T3 T4
NO
N1
F
Survival probability
SF3B4
High
Low
n.s.
n.s.
w
*
1.00
*
*
8
n.s.
9
9-
0.75
*
*
SF3B4
n.s
SF3B4
n.s
+·
+·
SF3B4
0.50
6
V
0.25
=
1
5-
12
-
0.00
p=2.176e-05
5
4
0 1 2 3 4 5 6 7 8 9 101112
-
4
Time(years)
.
.
A
Stage 1
Stage 2
Stage 3
Stage 4
T1
T2
T3
T4
NO
N1
FIGURE 9 | Relationship between hub SFs and survival curve, tumor stage, T stage and lymph node state. (A) YBX1. (B) SART1. (C) PRCC. (D) SNRPG. (E) SNRPE. (F) SF3B4. * P < 0.05; ** P < 0.01; *** P < 0.001 and n.s. indicates no significance.
A
GSE19750
n.s.
13
7
12
8
12.0
6
12
11.5
5
YBX1
11
SART1
7
PRCC
6
SNRPG
SNRPE
11
SF3B4
11.0
10
10.5
4
6
5
9
S
10
10.0
3
8
9.5
Tumor
Normal
Tumor
Normal
Tumor
Normal
Tumor
Normal
Tumor
Normal
Tumor
Normal
B
GSE33371
n.s.
n.s.
n.s.
n.s.
4.1
n.s.
n.s.
3.0
3.5
4.4
3.2-
4.3
2.9
YBX1
4.2
A
SART1
4.2
PRCC
SNRPG
SNRPE
4.0
SF3B4
3.4
3.0
2.8
4.1
4.0
2.7
3.3
2.8
4.0
3.9
2.6
3.8
3.9
3.2
Tumor
Normal
Tumor
Normal
Tumor
Normal
Tumor
Normal
Tumor
Normal
Tumor
Normal
FIGURE 10 | Expression level of hub SFs of ACC or normal samples in the microarray: (A) GSE19750 and (B) GSE33371. * P < 0.05; *** P < 0.001 and n.s. indicates no significance.
Construction of Splicing Correlation Network
Four hundred four SFs summarized from the previous study were used in this study for SF analysis (Seiler et al., 2018). They collected and prioritized the final list of 404 SF genes by compiling and filtering spliceosome and splicing related genes from three sources, which were all experimentally validated in the previous study (Barbosa-Morais et al., 2006; Hegele et al., 2012; Cvitkovic and Jurica, 2013). Source 1 (Hegele et al., 2012) was from a comprehensive yeast two-hybrid study using spliceosome components as bait. Source 2 (Barbosa-Morais et al., 2006) included 254 SFs and splicing related proteins. Source 3 was from SpliceosomeDB (Cvitkovic and Jurica, 2013). Of 404 SFs, 390 were extracted from the transcriptome data of 79 ACC cases. Co-expression network analysis calculated by spearman method was used in the construction of alternative splicing regulation networks to screen the important SFs. Cytoscape visualized a highly correlated interaction network with the threshold of correlation coefficient of 0.65.
Also, protein-protein interaction network of SFs was established to find the important SFs in the progression of ACC, using the online database Search Tool for the Retrieval of Interacting Genes4 (version 11.0). Then, the hub SFs were selected by the MCODE app of Cytoscape. Six prognostic-related hub SFs for patients with ACC were
identified by combining the co-expression networks analysis and protein-protein internecion networks.
Analysis of Hub Splicing Factors
Correlations of the mRNA expression level of hub SFs with clinical data (survival curve, tumor stage, T stage, and lymph node state) of 79 ACC cases were analyzed. The online database GEPIA (Tang et al., 2017)5, which contains the TCGA database with the GTEx database, was further used to evaluate the SFs on clinical features. Also, three gene expression microarray data were downloaded from GEO6 for illustrating the expression levels of hub SFs between normal tissue and tumor tissues. GSE19750 (Demeure et al., 2013) includes 44 ACC samples and 4 control samples. GSE33371 (Heaton et al., 2012) includes 33 ACC samples and 10 control samples.
DATA AVAILABILITY STATEMENT
The transcriptome data of 79 ACC cases and the corresponding clinical information can be downloaded from The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga/). Seven types of alternative splicing events (ASEs) of 79 ACC cases downloaded from TCGA SpliceSeq database (https://bioinf ormatics.mdanderson.org/TCGASpliceSeq/). Gene expression
5http://gepia.cancer-pku.cn/index.html
6http://www.ncbi.nlm.nih.gov/geo
microarray data GSE19750 and GSE33371 were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm. nih.gov/geo).
AUTHOR CONTRIBUTIONS
ZW and LL conceived and designed the project. JL and YH downloaded and analyzed the data. JL and LL drafted the manuscript. All authors contributed to text revision and discussion.
FUNDING
This study was supported by the National Natural Science Foundation of China (no. 81722007) and the
REFERENCES
Alagaratnam, S., Harrison, N., Bakken, A. C., Hoff, A. M., Jones, M., Sveen, A., et al. (2013). Transforming pluripotency: an exon-level study of malignancy-specific transcripts in human embryonal carcinoma and embryonic stem cells. Stem Cells Dev. 22, 1136-1146. doi: 10.1089/scd.2012. 0369
Anczukow, O., and Krainer, A. R. (2016). Splicing-factor alterations in cancers. RNA 22, 1285-1301. doi: 10.1261/rna.057919.116
Assie, G., Letouze, E., Fassnacht, M., Jouinot, A., Luscap, W., Barreau, O., et al. (2014). Integrated genomic characterization of adrenocortical carcinoma. Nat. Genet. 46, 607-612. doi: 10.1038/ng.2953
Barbosa-Morais, N. L., Carmo-Fonseca, M., and Aparicio, S. (2006). Systematic genome-wide annotation of spliceosomal proteins reveals differential gene family expansion. Genome Res. 16, 66-77. doi: 10.1101/gr.3936206
Barreau, O., Assie, G., Wilmot-Roussel, H., Ragazzon, B., Baudry, C., Perlemoine, K., et al. (2013). Identification of a CpG island methylator phenotype in adrenocortical carcinomas. J. Clin. Endocrinol. Metab. 98, E174-E184. doi: 10. 1210/jc.2012-2993
Bie, J., Liu, K., Song, G., Hu, X., Xiong, R., Zhang, X., et al. (2019). ENST00000489707.5 is a preferred alternative splicing variant of PTK7 in adrenocortical cancer and shows potential prognostic value. Med. Sci. Monit. 25, 8326-8334. doi: 10.12659/MSM.919818
Chen, J., and Weiss, W. A. (2015). Alternative splicing in cancer: implications for biology and therapy. Oncogene 34, 1-14. doi: 10.1038/onc.2013.570
Chen, Q. F., Li, W., Wu, P., Shen, L., and Huang, Z. L. (2019). Alternative splicing events are prognostic in hepatocellular carcinoma. Aging 11, 4720-4735. doi: 10.18632/aging.102085
Crona, J., and Beuschlein, F. (2019). Adrenocortical carcinoma - towards genomics guided clinical care. Nat. Rev. Endocrinol. 15, 548-560. doi: 10.1038/s41574- 019-0221-227
Cvitkovic, I., and Jurica, M. S. (2013). Spliceosome database: a tool for tracking components of the spliceosome. Nucleic Acids Res. 41, D132-D141. doi: 10. 1093/nar/gks999
Demeure, M. J., Coan, K. E., Grant, C. S., Komorowski, R. A., Stephan, E., Sinari, S., et al. (2013). PTTG1 overexpression in adrenocortical cancer is associated with poor survival and represents a potential therapeutic target. Surgery 154, 1405-1416. doi: 10.1016/j.surg.2013.06.058
Frankiw, L., Baltimore, D., and Li, G. (2019). Alternative mRNA splicing in cancer immunotherapy. Nat. Rev. Immunol. 19, 675-687. doi: 10.1038/s41577-019- 0195-197
He, R. Q., Zhou, X. G., Yi, Q. Y., Deng, C. W., Gao, J. M., Chen, G., et al. (2018). Prognostic signature of alternative splicing events in bladder urothelial carcinoma based on Spliceseq data from 317 cases. Cell Physiol. Biochem. 48, 1355-1368. doi: 10.1159/000492094
Heaton, J. H., Wood, M. A., Kim, A. C., Lima, L. O., Barlaskar, F. M., Almeida, M. Q., et al. (2012). Progression to adrenocortical tumorigenesis in mice and
National Health Commission of China (no. 2017ZX1030440 2001-008).
ACKNOWLEDGMENTS
We appreciate the generosity of TCGA, TCGA SpliceSeq, GEO, and GEPIA for sharing the huge amount of data. We also thank members of our laboratory for supporting this work.
SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2020.00918/full#supplementary-material
humans through insulin-like growth factor 2 and beta-catenin. Am. J. Pathol. 181, 1017-1033. doi: 10.1016/j.ajpath.2012.05.026
Hegele, A., Kamburov, A., Grossmann, A., Sourlis, C., Wowro, S., Weimann, M., et al. (2012). Dynamic protein-protein interaction wiring of the human spliceosome. Mol. Cell 45, 567-580. doi: 10.1016/j.molcel.2011. 12.034
Ishida, H., Komiya, S., Inoue, Y., Yutani, S., Inoue, A., and Itoh, K. (2000). Expression of the SART1 tumor-rejection antigen in human osteosarcomas. Int. J. Oncol. 17, 29-32. doi: 10.3892/ijo.17.1.29
Jia, D., Wei, L., Guo, W., Zha, R., Bao, M., Chen, Z., et al. (2011). Genome- wide copy number analyses identified novel cancer genes in hepatocellular carcinoma. Hepatology 54, 1227-1236. doi: 10.1002/hep.24495
Jouinot, A., and Bertherat, J. (2018). MANAGEMENT OF ENDOCRINE DISEASE: adrenocortical carcinoma: differentiating the good from the poor prognosis tumors. Eur. J. Endocrinol. 178, R215-R230. doi: 10.1530/EJE-18- 0027
Kozlovski, I., Siegfried, Z., Amar-Schwartz, A., and Karni, R. (2017). The role of RNA alternative splicing in regulating cancer metabolism. Hum. Genet. 136, 1113-1127. doi: 10.1007/s00439-017-1803-x
Lee, S. C., and Abdel-Wahab, O. (2016). Therapeutic targeting of splicing in cancer. Nat. Med. 22, 976-986. doi: 10.1038/nm.4165
Liu, Z., Li, W., Pang, Y., Zhou, Z., Liu, S., Cheng, K., et al. (2018). SF3B4 is regulated by microRNA-133b and promotes cell proliferation and metastasis in hepatocellular carcinoma. eBio Med. 38, 57-68. doi: 10.1016/j.ebiom.2018. 10.067
Oltean, S., and Bates, D. O. (2014). Hallmarks of alternative splicing in cancer. Oncogene 33, 5311-5318. doi: 10.1038/onc.2013.533
Pan, Q., Shai, O., Lee, L. J., Frey, B. J., and Blencowe, B. J. (2008). Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing. Nat. Genet. 40, 1413-1415. doi: 10.1038/ng.259
Patel, D., Boufragech, M., Jain, M., Zhang, L., He, M., Gesuwan, K., et al. (2013). MiR-34a and miR-483-5p are candidate serum biomarkers for adrenocortical tumors. Surgery 154, 1224-1228. doi: 10.1016/j.surg.2013.06.022
Quidville, V., Alsafadi, S., Goubar, A., Commo, F., Scott, V., Pioche-Durieu, C., et al. (2013). Targeting the deregulated spliceosome core machinery in cancer cells triggers mTOR blockade and autophagy. Cancer Res. 73, 2247-2258. doi: 10.1158/0008-5472.CAN-12-2501
Rossner, F., Gieseler, C., Morkel, M., Royer, H. D., Rivera, M., Blaker, H., et al. (2016). Uncoupling of EGFR-RAS signaling and nuclear localization of YBX1 in colorectal cancer. Oncogenesis 5:e187. doi: 10.1038/oncsis. 2015.51
Ryan, M., Wong, W. C., Brown, R., Akbani, R., Su, X., Broom, B., et al. (2016). TCGASpliceSeq a compendium of alternative mRNA splicing in cancer. Nucleic Acids Res. 44, D1018-D1022. doi: 10.1093/nar/gkv1288
Sasatomi, T., Yamana, H., Shichijo, S., Tanaka, S., Okamura, T., Ogata, Y., et al. (2000). Expression of the SART1 tumor-rejection antigens in colorectal cancers. Dis. Colon Rectum. 43, 1754-1758. doi: 10.1007/bf02236863
Seiler, M., Peng, S., Agrawal, A. A., Palacino, J., Teng, T., Zhu, P., et al. (2018). Somatic mutational landscape of splicing factor genes and their functional consequences across 33 cancer Types. Cell Rep. 23, 282-296. doi: 10.1016/j. celrep.2018.01.088
Shen, Q., Eun, J. W., Lee, K., Kim, H. S., Yang, H. D., Kim, S. Y., et al. (2018). Barrier to autointegration factor 1, procollagen-lysine, 2-oxoglutarate 5-dioxygenase 3, and splicing factor 3b subunit 4 as early-stage cancer decision markers and drivers of hepatocellular carcinoma. Hepatology 67, 1360-1377. doi: 10.1002/ hep.29606
Shen, Q., and Nam, S. W. (2018). SF3B4 as an early-stage diagnostic marker and driver of hepatocellular carcinoma. BMB Rep. 51, 57-58. doi: 10.5483/bmbrep. 2018.51.2.021
Siegfried, Z., and Karni, R. (2018). The role of alternative splicing in cancer drug resistance. Curr. Opin. Genet. Dev. 48, 16-21. doi: 10.1016/j.gde.2017.10.001
Szabo, D. R., Luconi, M., Szabo, P. M., Toth, M., Szucs, N., Horanyi, J., et al. (2014). Analysis of circulating microRNAs in adrenocortical tumors. Lab. Invest. 94, 331-339. doi: 10.1038/labinvest.2013.148
Tamura, K., Furihata, M., Tsunoda, T., Ashida, S., Takata, R., Obara, W., et al. (2007). Molecular features of hormone-refractory prostate cancer cells by genome-wide gene expression profiles. Cancer Res. 67, 5117-5125. doi: 10.1158/ 0008-5472.CAN-06-4040
Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017). GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 45, W98-W102. doi: 10.1093/nar/gkx247
Wang, C., Zheng, M., Wang, S., Nie, X., Guo, Q., Gao, L., et al. (2019). Whole genome analysis and prognostic model construction based on alternative splicing events in endometrial cancer. Biomed. Res. Int. 2019:2686875. doi: 10.1155/2019/2686875
Xiong, Y., Deng, Y., Wang, K., Zhou, H., Zheng, X., Si, L., et al. (2018). Profiles of alternative splicing in colorectal cancer and their clinical significance: a study based on large-scale sequencing data. eBio Med. 36, 183-195. doi: 10.1016/j. ebiom.2018.09.021
Yan, C., Wan, R., and Shi, Y. (2019). Molecular mechanisms of pre-mRNA splicing through structural biology of the spliceosome. Cold Spring Harb. Perspect. Biol. 11:a032409. doi: 10.1101/cshperspect.a032409
Yutani, S., Shichijo, S., Inoue, Y., Kawagoe, N., Okuda, K., Kurohiji, T., et al. (2001). Expression of the SART1 tumor-rejection antigen in hepatocellular carcinomas. Oncol. Rep. 8, 369-372. doi: 10.3892/or.8.2.369
Zhang, J., Deng, Y., Zuo, Y., Wang, J., and Zhao, Y. (2020). Analysis of colorectal cancer-associated alternative splicing based on transcriptome. DNA Cell Biol. 39, 16-24. doi: 10.1089/dna.2019.5111
Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10:1523. doi: 10.1038/s41467-019-09234- 9236
Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Copyright @ 2020 Lv, He, Li and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.