3.1 BAP1 acts as a tumor suppressor in KIRC
A flow diagram of ceRNA construction and analysis is shown in Figure 1. We conducted a pan-cancer analysis to evaluate the RNA expression level of BAP1 using data from UCSC XENA (https://xenabrowser.net/datapages/), and protein expression level through Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). We found that the RNA expression level of BAP1 was significantly downregulated in KIRC (Fig. 2A), and interestingly BAP1 protein level was the lowest in renal cancer (Fig S1). Immunohistochemistry (IHC) staining obtained from HPA also validated the downregulation of BAP1 in tumor tissue (Fig. 2B).
To investigate the relationship between BAP1 and the prognosis of KIRC, we used Kaplan-Meier survival curves to compared the survival of high- and low-BAP1 expression group. Our result indicated that low expression of BAP1 was associated with poor OS in KIRC (Fig. 2C).
Moreover, we utilized cBioPortal (http://www.cbioportal.org/) to explore the potential mechanisms underlying the abnormally lower expression of BAP1 in KIRC [31]. Figure 2D showed the genetic alteration rate of BAP1 is 13%. Gene deletions (deep deletion and shallow deletion) account for more than half of the copy number alterations in KIRC samples (Fig. 2E). In addition, the mRNA expression level of BAP1 is found to be positively correlated with the copy number value (Fig. 2F).
As listed above, BAP1 is down-regulated in KIRC, and deletions of copy number may be one of the potential mechanisms underlying the lower expression level of BAP1 in KIRC.
3.2 Identification of different expressed lncRNAs, miRNAs, and mRNAs
We aimed to construct a prognostic-related ceRNA network based on the expression of BAP1. Therefore, we performed differentially expressed genes (DEGs) analysis comparing BAP1high and BAP1low group in KIRC samples. Finally, a total of 3425 DElncRNA, 84 DEmiRNA, and 2753 DEmRNA were screened according to the cut-off value (Fig. 3).
3.3 Construction of lncRNA-miRNA-mRNAs networks
After taking intersection of the predicted miRNAs with the obtained DEmiRNAs, four miRNAs were selected. Then, we used miRBD database and Targetscan databases to identify the downstream targeting mRNAs of the four predicted miRNAs. Similarly, by taking intersections with DEmRNAs, 387 predicted mRNAs were identified. Finally, 4 lncRNAs, 4 miRNAs, and 387 mRNAs were selected.
Then, we put these genes into Cytoscape (version 3.6.1) for ceRNA network construction (Fig. 4A), and plug-in “cytoHubba” was utilized for hub genes network construction. Finally, lncRNAs (XIST, HELLPAR, PURPL), miRNAs (miR-10a-5p, miR-508-3p, miR-135a-5p), and mRNAs (IRS1, SERPINE1, KCAN1, TRIM2, RORB, SIX4) were identified (Fig. 4B).
We conducted functional enrichment analysis of DEmRNAs to investigate the potential function of these ceRNA networks through Metascape. The result demonstrated that these DEmRNAs were involved in mesenchyme development, transmembrane receptor protein tyrosine kinase signaling pathway, regulation of cell adhesion (Fig. 4C).
3.4 Construction of prognostic-related ceRNA in KIRC
To construct a prognostic-related ceRNA, we first analyzed the expression of hub genes comparing tumor and normal tissue, and the results showed that differences were significant in all hub genes (Fig. 5). Then, survival analyses were performed to select prognostic-related genes as shown in Figure 6. In total, two DElncRNAs (XIST, HELLPAR), one DEmiRNA (miR-10a-5p), and four DEmRNAs (SERPINE1, TRIM2, RORB, SIX4) were found to be prognostic-related genes.
We utilized lncLocator (www.csbio.sjtu.edu.cn/bioinf/lncLocator/) to determine the subcellular localization of two DElncRNAs, and found that XIST was mainly distributed in cytoplasm (Fig. 7A), which indicated that XIST may act as a ceRNA in improving the expression of SERPINE1 through sponge of miR-10a-5p. Finally, a prognostic-related XIST/miR-10a-5p/SERPINE1 ceRNA network was constructed (Fig. 7B). In addition, the target sites of XIST and SERPINE1 to miR-10a-5p were predicted by Targetscan (Fig. 7C).
3.5 Clinical relevance of the XIST/has-miR-10a-5p/SERPINE1 axis in KIRC patients
We evaluated the expression level of XIST, has-miR-10a-5p and SERPINE1 with different clinical factors (T stage, N stage, M stage, pTNM stage, tumor grade, gender, and age) to determine the clinical relevance of this ceRNA axis. The results showed that the expression of XIST was correlated with tumor grade, and sex (P < 0.05, Fig. S2A). Lower expression level of has-miR-10a-5p was found to associated with higher T, M stages, pTNM stages, tumor grade and gender (P < 0.05, Fig. S2B). In addition to age and M stage, the expression level of SERPINE1 was strongly correlated with T stage, N stage, pTNM stage, tumor grade and gender (P < 0.05, Fig. S2C). These findings are generally consistent with the results of survival analysis.
Moreover, we performed univariate and multivariate Cox regression model analysis and ROC (receiver operating characteristic curve) analysis to determine the OS-related prognostic significance of this ceRNA axis. In multivariate Cox regression analysis, SERPINE1 (HR=1.456, P=0.015) and has-miR-10a-5p (HR=0.681, P=0.014) were found to be independent prognostic factor in KIRC, while XIST was not (Table. S1-3). The AUC (area under curves) of ROC analysis suggested a good prognostic performance of SERPINE1 (AUC=0.789) and has-miR-10a-5p (AUC=0.892). These results were showed in Fig. S3.
In addition, we performed a pan-cancer analysis of SERPINE1 mRNA expression and found that SERPINE1 was highly expressed in kidney cancer (Fig. S4A). Immunohistochemical analysis revealed that SERPINE1 was located in the cytoplasmic/membranous and was only detected in the renal tubules (Fig. S4B).
3.6 DNA methylation analysis of SERPINE1
To elucidate the mechanism of the abnormally high expression of SERPINE1 in KIRC, we performed a series of methylation analyses of SERPINE1. Firstly, the results of co-expression analysis suggested that the expression of SERPINE1 was positively correlated with the expression levels of the DNMT1, DNMT3A and DNMT3B (P < 0.05, Fig. 8A). Secondly, we obtained the same result with DiseaseMeth version 2.0, that the methylation level of SERPINE1 in normal tissues were much higher than in KIRC (P < 0.001, Fig. 8B). Thirdly, we found that the methylation level of SERPINE1 correlated with OS. Moreover, we also identified 12 methylation sites on DNA sequences that were negatively correlated with SERPINE1 expression (Fig. 8C).
3.7 Immune infiltration analysis of SERPINE1 in KIRC
To further investigate the relationship between the expression level of SERPINE1 and the immune microenvironment in KIRC, we performed an immune infiltration analysis with TIMER. The “SCNA” module analysis indicated that the immune infiltration level of CD 4+ T cell was associated with altered copy numbers of SERPINE1 (Fig. 9A). Moreover, the “Gene” module analysis showed that the expression level of SERPINE1 was positively related with immune infiltration level of CD 4+ T cell, CD 8+ T cell, macrophages, dendritic cells, neutrophils (Fig. 9B).