Identification of Potential Biomarkers of Prognosis-Related Long Non-Coding RNA (lncRNA) in Pediatric Rhabdoid Tumor of the Kidney Based on ceRNA Networks

Background

Rhabdoid tumors of the kidney (RTK) are rare but invasive renal tumors that occur in infancy and early childhood and account for approximately 1.3% of all renal tumors in children [1]. The majority of RTKs have SMARCB1 (SWI/SNF Related, Matrix Associated, Actin-Dependent Regulator of Chromatin, Subfamily B, Member 1, SMARCB1) double allele-inactivating mutations [2]. Patients with RTK generally have low survival rates, and recurrence, disease progression, and death usually occur early, making it the most lethal malignant solid tumor among children. Studies have reported that patients with RTK aged <6 months have a 2-year event-free survival (EFS) rate of 15%, while those >27 months have a 2-year EFS of 48% [3]. The main modalities of treatment currently include surgery, radiotherapy, and chemotherapy. However, it has been found that high-dose chemotherapy does not improve the prognosis of patients [4]. At present, there are few reports regarding the prognostic molecular markers of RTK. Therefore, we sought to identify some novel markers to screen RTK earlier and predict patient prognosis.

Long non-coding RNA (lncRNA) refers to a class of RNA molecules located in the nucleus and cytoplasm of cells widely found in eukaryotes, which form a complex regulatory network with target genes and are involved in the pathophysiological processes of the organism [5]. Recently, an increasing number of studies have discovered that lncRNA can be used as a diagnostic and prognostic marker for a variety of tumors, such as gastric cancer [6], liver cancer [7], and pancreatic cancer [8]. The competitive endogenous RNA (ceRNA) hypothesis [9] suggests that lncRNAs not only directly regulate the expression of target genes but also likely affect the number of miRNAs by adsorbing the core seed sequences of miRNAs to further affect the mRNA levels of target genes [10]. The ceRNA network is extensively involved in the regulation of vital activities such as inflammation, apoptosis, and differentiation in the body. However, little is known regarding the role of lncRNA and related ceRNA regulatory networks in RTK.

We downloaded the RNA expression profiles of PTK from the TARGET database for analysis and screened for differentially expressed RNAs. Through use of various bioinformatics software, the ceRNA network was constructed and many genes of potential research interest were identified.

Material and Methods

ACQUISITION OF TARGET PEDIATRIC RTK DATA:

The RNA sequencing data and corresponding clinical information of pediatric RTK were downloaded from the TARGET database (https://ocg.cancer.gov/programs/target; last updated: March 31, 2020). The miRNA expression microarray includes 6 normal samples and 66 tumor samples. The mRNA microarray and lncRNA microarray includes 6 normal samples and 65 tumor samples. DElncRNAs, DEmRNAs, and DEmiRNAs were selected based on the “EdgeR” package of R software (3.6.3). For the cutoff criteria, a false discovery rate (FDR) <0.05 and an absolute log 2-fold change >2 were the cutoff criteria for differentially expressed genes. Using the R software for “ggplot2” and the “pheatmap” package, we built volcano maps and heat maps of DEmRNAs, DElncRNAs, and DEmiRNAs.

CERNA NETWORK CONSTRUCTION:

In all, 107 DElncRNAs were used to predict targeted miRNAs by miRcode [11] (http://starbase.sysu.edu.cn/). The targeted miRNAs were further compared with DEmiRNAs; the intersected miRNAs were retained to develop the DElncRNAs-DEmiRNAs ceRNA network. Then, the DEmiRNAs in the DElncRNAs-DEmiRNAs ceRNA network were distributed via MiRTarBase [12] (http://mirtarbase.mbc.nctu.edu.tw), MiRDB [13] (http://www.mirdb.org/miRDB/), and TargetScan [14] (http://www.tar-getscan.org/) to predict miRNA-targeted mRNAs. The targeted mRNAs were further compared with DEmRNAs, and the intersected DEmRNAs were retained to develop the DEmiRNAs-DEmRNAs ceRNA network. Finally, Cytoscape (v3.8.0) [15] was used for constructing the ceRNA network.

FUNCTIONAL ENRICHMENT ANALYSIS:

GO enrichment and KEGG pathway enrichment analysis of the ceRNA network-related DEmRNA were achieved using R software with the “clusterProfiler” package, and the “enrichplot” package was used to draw related functional enrichment and pathway enrichment.

SURVIVAL ANALYSIS:

Survival curves and the Kaplan-Meier estimator for the ceRNA network-related DElncRNA were computed and plotted with the ‘survival’ package, and P<0.05 served as the statistical threshold to screen for survival-associated DElncRNAs in the ceRNA network.

UNIVARIATE AND MULTIVARIATE ANALYSES TO CONFIRM THE PROGNOSTIC-ASSOCIATED LNCRNA FACTORS IN RTK:

Based on the clinical information of RTK patients, univariate analysis was performed for survival-associated DElncRNAs in the ceRNA network using the log rank test. For survival-associated DElncRNAs with a P value of <0.05 in univariate analysis, least absolute shrinkage and selection operator (LASSO) analysis was undertaken using the “glmnet” package to select the minimum error value to eliminate over-fitting data, Then, multivariate analyses were performed using the Cox regression model to identify the prognostic-associated DElncRNA factors in RTK. By using the “survivalROC” package, we constructed time-dependent receiver-operating characteristic (ROC) curves to analyze the predictive accuracy of the risk score model. Based on the “survminer” R package, patients were divided into high- and low-risk groups according to the risk score (median cutoff value), and survival was analyzed using the Kaplan-Meier method. Risk heat maps, risk scores, and survival status maps were plotted to verify the predictive ability of the model.

Results

DIFFERENTIALLY EXPRESSED RNAS IN PEDIATRIC RTK:

A total of 1109 lncRNAs were screened out, of which 421 were upregulated and 688 were downregulated. Of the 215 miRNAs, 170 were upregulated and 45 were downregulated. Of the 3436 differential mRNAs, 1607 were upregulated and 1829 were downregulated. Figure 1 shows the volcano maps and heat maps.

CONSTRUCTION OF CERNA NETWORKS AND GO AND KEGG ANALYSES:

We searched the 1109 DElncRNAs in the miRcode database. Then, 107 DElncRNAs and 21 DEmiRNAs were used to construct the 513 pairs in the DElncRNAs-DEmiRNAs network. Next, we predicted 21 DEmiRNAs in all 3 target-predicting databases; then, 21 DEmiRNAs and 74 DEmRNAs constructed 93 pairs of the DEmiRNA-DEmRNAs network. The lncRNA-miRNA-mRNA ceRNA network was visualized using Cytoscape software (Figure 2A). The GO and KEGG enrichment analyses were further used to display the biological function of the 74 ceRNA network-related DEmRNAs. The GO enrichment analysis showed that DEmRNAs were mostly related to the regulation of DNA-binding transcriptional activator activity, RNA polymerase II specificity, and phospholipase C activity (Figure 2B, Table 1). Further, KEGG enrichment analysis showed that the mRNAs were highly enriched in MAPK signaling pathway (Table 1).

SURVIVAL ANALYSIS AND ESTABLISHMENT OF THE CERNA NETWORK-RELATED DELNCRNAS PROGNOSTIC MODEL:

Survival analysis results showed that ADAMTS9-AS2, DSCR9, EMX2OS, HCG22, HNF1A-AS1, HOTAIR, SNHG6, TPTEP1, and ZNF503-AS2 were correlated with prognosis (P<0.05) (Figure 3). These significant survival-associated DElncRNAs were further submitted to univariate and multivariate Cox regression analyses to screen prognosis-related DElncRNAs for RTK. Subsequently, 6 lncRNAs were obtained after univariate Cox regression analysis by using P<0.05 (Table 2). All 6 lncRNAs underwent LASSO regression analysis (Figure 4A, 4B), and a prognostic model for prognostic-associated lncRNAs was constructed. The 4 prognostic-associated lncRNAs contained HNF1A-AS1, TPTEP1, and SNHG6, and ZNF503-AS2 (Figure 4C, Table 2). The model risk score was calculated as follows: risk score=(0.460038×exp_ HNF1A-AS1)+(0.249016×exp_ ZNF503-AS2)+(−1.557758×exp_) SNHG6)+(0.443621×exp_ TPTEP1), where HNF1A-AS1, TPTEP1, and ZNF503-AS2 were associated with high risk of death and SNHG6 is associated with low risk of death. The 4 prognostic-associated lncRNAs of the distribution of survival risk score (Figure 5A), survival time (Figure 5B), and expression heatmap (Figure 5C) revealed that the high-risk group was more likely to have a worse prognosis. According to the median risk score (median risk score=0.863512), pediatric RTK patients were divided into a low-group and high-risk group. Survival analysis revealed that the low-risk group was had a better prognosis (P<0.0001, Figure 5D). The ROC curves for the 4 lncRNAs in the signature had area under the curve values >0.5 (Figure 5E, 5F). These findings suggested that the 4 identified lncRNAs may serve as meaningful prognostic markers.

Discussion

RTK is a highly malignant and aggressive tumor, and, unfortunately, the overall prognosis for RTK is still only about 20–25% [3]. Therefore, accurate prediction of prognosis is essential for the personalized treatment of these patients. We sought to find novel biomarkers to diagnose RTK earlier and predict the prognosis of patients. We screened our 4 prognosis-associated lncRNAs based on a ceRNA network combined with clinical information of patients to identify prognostic biomarkers of RTK.

At present, the ceRNA hypothesis suggests that lncRNA can competitively bind to miRNA and act as an “miRNA sponge” to indirectly regulate the downstream target genes of the miRNA [9]. A growing body of evidence suggests that lncRNAs can serve as molecular biomarkers in multiple cancer types to predict patient prognosis. However, thus far, there have been very few reports on the ceRNA network related to pediatric RTK. Ren et al. obtained 109 samples of rhabdoid tumor through the GEO database, including 31 rhabdoid tumor and 21 normal kidneys, and through use of the co-expression method, they finally identified RPL5 and RPL10 as potential diagnostic markers for AT/RT [16]. Their research provides an idea for the diagnosis of rhabdoid tumor. However, due to the lack of clinical survival information of patients, the markers they screened cannot be judged to predict prognosis. In our research, we downloaded RNA-seq data from the TARGET database, which is a database with complete clinical information and sequencing information, and we constructed a ceRNA network containing 107 DElncRNAs, 21 DEmiRNAs, and 78 DEmRNAs. To screen out prognosis-related lncRNAs, we conducted survival analysis on 107 DElncRNAs in the ceRNA network and found that 9 were related to survival. These significant survival-associated DElncRNAs were further submitted to univariate and multivariate Cox regression analyses to screen prognosis-related lncRNAs for RTK. Finally, 4 prognosis-related lncRNAs (HNF1A-AS1, TPTEP1, SNHG6, and ZNF503-AS2) were used to establish a prognostic model. The 3- and 5-year survival probability areas under the ROC curves were 0.93 and 0.889, respectively, which fully explained the accuracy of the lncRNA prognostic models.

To study the function of the ceRNA network, we performed KEGG and GO enrichment analysis on DEmRNAs in ceRNA, which are mainly enriched in the MAPK signaling pathway. In cancer, the constitutively active MAPK signaling pathway often leads to promotion of abnormal cell growth and tumorigenesis [17] such as liver cancer [18] and gastric cancer [19]. In colon cancer, HNF1A-AS1 upregulates OTX1 to activate the ERK/MAPK pathway [20]. Among the 4 prognostic-associated lncRNAs found in our study, HNF1A-AS1 was enriched in the MAPK signaling pathway. Therefore, we predicted that HNF1A-AS1 can promote the development of RTK by activating the ERK/MAPK pathway. The high expression of TPTEP1 can predict the prolonged overall survival of patients with lung cancer [21]. TPTEP1 overexpression suppressed cell proliferation and metastasis but did not affect cell apoptosis in hepatocellular carcinoma [22]. However, in our study, the survival rate of patients with TPTEP1 overexpression was low, which likely means that in RTK, TPTEP1 overexpression can promote tumor development, which needs to be further confirmed in future studies. Overexpression of SNHG6 is related to tumor progression and poor prognosis. In rectal cancer, SNHG6 promotes the migration, invasion, and epithelial-mesenchymal transition of colorectal cancer cells through the miR-26a/EZH2 axis [23]. SNHG6 is overexpressed in childhood Wilms’ tumors, and SNHG6 knockout can inhibit proliferation, migration, and invasion [24]. Our study found that patients with high SNHG6 expression have a lower survival rate. It is possible that in patients with RTK, high SNHG6 expression can inhibit tumor growth. ZNF503 is involved in the development of tumors. In hepatocellular carcinoma, ZNF503 promotes migration, invasion, and epithelial-mesenchymal transition (EMT) processes by regulating GATA3 expression [25]. However, very few studies have focused on the relationships among HNF1A-AS1, TPTEP1, SNHG6, and ZNF503-AS2 in RTK. Therefore, we need further experimental research to clarify the molecular mechanism and biological functions of these prognostic-related lncRNAs in RTK.

Our study has some limitations. First, our clinical sample size was not large enough. Second, we were able to predict the relationship between lncRNA and miRNA and between miRNA and mRNA only through online databases; hence, further experimental research should be conducted to address this issue. Last, the 4 prognostic-related lncRNAs selected in this study did not perform any biological functions. However, the strength of our research lies in being one of the first studies to predict the molecular mechanism of RTK from the perspective of the ceRNA network, and at the same time screen out 4 lncRNA prognostic biomarkers based on the ceRNA network combined with clinical prognostic information. This likely provides new ideas for RTK research.

Conclusions

In conclusion, we identified 4 lncRNAs related to RTK patient prognosis that can be used as candidate biomarkers and therapeutic targets. Further large-scale studies are required to validate and strengthen our findings.