Integrated Bioinformatics Analysis for the Identification of Key lncRNAs, mRNAs, and Potential Drugs in Clear Cell Renal Cell Carcinomas

Introduction

Renal cell carcinoma (RCC) is the second leading cause of death related to urinary system malignancy, and about 75% of RCCs are clear cell RCC (ccRCC).^1–3 Histologically, ccRCC cells have a clear cytoplasm, with cell clusters arranged in nests, surrounded by a fine vascular network.^4,5 Currently, patients with ccRCC can use a variety of clinical management methods in the early stage, including surgery and radiotherapy. The 5-year survival rate after diagnosis has improved to some extent. However, the overall survival (OS) is still very poor, especially for patients in the advanced stage, mainly due to the lack of more effective treatments.⁶ Therefore, it is urgent to understand the underlying pathological mechanism of ccRCC more comprehensively and to further develop new therapeutic strategies.

The rapid development of genomics and transcriptomics has highlighted the important role of non-coding RNA in tumors. Long non-coding RNAs (lncRNAs), RNA molecules with more than 200 nucleotides in length, do not encode proteins.⁷ LncRNA plays important roles in cell growth, apoptosis, migration, invasion and autophagy, and other cellular processes.^8–10 Compelling data suggest specific lncRNA promotes the development of human cancer. For example, lncRNA HOTAIR is highly expressed in breast cancer samples, and the increased expression is associated with metastasis and death of patients. LncRNA TUG1 acts as a potential biomarker and therapeutic target in bladder urothelial carcinoma.¹¹ Xie et al¹² indicated lncRNA SPRY4-IT1 was up-regulated in esophageal squamous cell carcinoma, which is significantly related to clinicopathological staging and OS rate.

The database represented by TCGA provides an ideal opportunity for researchers to comprehensively explore the molecular mechanisms of tumorigenesis of various cancers, identify new therapeutic targets and prognostic molecular markers. In order to reveal the potential pathogenesis of ccRCC, RNA-seq data of 530 ccRCC cases in TCGA were downloaded, and a comprehensive analysis was carried out using bioinformatics tools. The machine learning methods were used to construct three machine learning models of ccRCC, survival analysis of related lncRNA was carried out, and the co-expression network of lncRNA-mRNA was constructed to identify potential therapeutic targets or prognostic biomarkers in ccRCC. In addition, potential drugs were screened through DGIdb database, and molecular docking was conducted to explore the interaction between candidate drugs and targets to further develop new therapeutic strategies of ccRCC.

Materials and Methods

Patients and Samples

RNA-seq data and clinical information were collected from the TCGA database. The inclusion criteria are as follows: (1) patients were histopathologically diagnosed with ccRCC; (2) patients had complete demographic data, including age, life status, race, pathological stage, TNM classification, and OS time. A total of 530 ccRCC samples and 72 adjacent non-tumor samples were included. Demographic data of included ccRCC patients is listed in Table 1.

Table 1 Demographic Data of 530 ccRCC Cases

RNA-Seq Data Analysis

The original data of lncRNA and mRNA were processed. More than 20% of the samples with expression value of 0 were filtered out for mRNAs. The lncRNA and mRNA were normalized through Voom and TNM normalization methods. All analyses were performed in R software (version 3.5.2). The Limma package (version 3.36.5) was applied to analyze differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) in ccRCC and adjacent non-tumor samples. Statistical significance was set as log2 fold change (FC) >3, and mRNA/lncRNA with the false discovery rate (FDR) <0.05.

Functional Enrichment Analysis

In addition to the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway analysis, gene ontology (GO) functional enrichment analysis based on biological process, cellular component, and molecular function was also carried out. Functional enrichment analysis on DEmRNAs was performed using the David database. The ggplot2 package in the R software was used to visualize the results.

Immune Cell Infiltration Analysis and Screening of Key Immune Cells

The single sample gene set enrichment analysis (ssGSEA) algorithm was used to quantify the relative abundance of each cell infiltration in the tumor immune microenvironment (TIME) of ccRCC. From Charoentong’s research, the gene set of labeled infiltrating immune cell type in TIME was obtained, which is rich in a variety of human immune cell subtypes.^13,14 Enrichment score calculated by ssGSEA analysis was used to represent the relative abundance of TIME infiltrating cells in each sample. The difference of immune cell infiltration between the normal tissue group and the ccRCC group was analyzed by the Wilcoxon test hypothesis test. To reduce the number of immune cell variables, the key immune cells were screened from the obtained differential immune cells, and the operator LASSO in minimum absolute contraction and selection regression analysis of the R software “glmnet” package was used.

Feature Selection and Modeling

Machine learning algorithm modeling was used to identify lncRNA biomarkers that were most suitable for ccRCC diagnosis. Through the regression analysis of the least absolute contraction and selection operator (LASSO) of the glmnet package (version 2.0.18) in R,¹⁵ key mRNAs were screened from the DElncRNAs. The DElncRNAs were selected by LASSO in descending order according to the value of the average decreasing accuracy through the random forest algorithm in the R package randomForest (version 4.6.14). The random forest algorithm and the top-down method of 10 cross-validations were used, and one DElncRNAs was added each time to obtain the optimal number of features. Finally, classification models of random forest, support vector machine (e1071, version 1.1.7), and decision tree (rpart, version 4.1.15) were established, respectively. The diagnostic capabilities of the three models and each lncRNA biomarker were evaluated by the receiver operating characteristic (ROC) curve. The area under the curve (AUC) was the evaluation indicator of model performance.

RNA Extraction and Quantitative Real-Time PCR (qRT-PCR)

To further verify the expression of selected lncRNAs, another 14 tissue samples, including 7 ccRCC tissues and 7 paired adjacent tissue samples, were included. All ccRCC samples were pathologically confirmed and had no other malignant tumors or autoimmune diseases. These patients did not receive radiotherapy, chemotherapy, immunotherapy, or molecular targeted therapy before the diagnosis. A total of 50 to 100 mg of tumor tissues were lysed by TRIzol reagent. Total RNA was extracted and purified with the ethanol protocol reported previously.¹⁶ After checking the purity and integrity of RNA, a total of 1 μg RNA was reverse transcribed using the FastKing cDNA reverse transcription kit (#KR116, TIANGEN, China) according to the manufacturer’s instruction. The expression was quantitatively reversed by reverse transcription PCR (qRT-PCR) by SYBR Green SuperReal PreMix Plus and the SuperReal enhanced fluorescence quantitative premixed reagent systems (#FP205, TIANGEN, China). GADPH and ACTB served as control. The data obtained were calculated by the formula relative quantification (RQ) = 2^{− ΔΔ CT} method by ABI 7300 fluorescence quantitative PCR instrument (Thermo Fisher Scientific, USA). The used primers are listed in Supplementary Table 1. The reaction conditions were 95°C for 15 min, 95°C for 10s, 55°C for 30s and 72°C for 32s for 40 cycles.

The obtained tissue sample was approved by the ethics committee of Taizhou Municipal Hospital Affiliated with Taizhou University (LW113). The informed consent was obtained from each patient and their families for clinical data use and publication.

Correlation Analysis Between the Optimal Diagnostic lncRNAs and Clinical Features

To explore the correlation between clinical characteristics and the optimal diagnostic lncRNAs, the rank-sum test was used to analyze their expressions in adjacent non-tumor samples and tumors with different stages and grades. Finally, a violin chart was used to show the relationship between the optimal diagnostic lncRNAs and clinical features. In addition, the Kaplan–Meier plotter database was used to conduct a survival analysis of the optimal diagnostic lncRNAs.

Co-Expression Analysis of the Optimal Diagnostic lncRNAs and DEmRNAs

The correlation analysis between the optimal diagnostic lncRNAs and DEmRNAs was carried out to construct the lncRNA-mRNA co-expression network. The Pearson correlation coefficient was calculated. Coefficient with cor > 0.6 and P < 0.01 was considered as significant correlation. The bioinformatics software Cytoscape (version 3.7.2) was used to construct the co-expression network, and the degree of centrality was used to reveal the relative importance of mRNA/lncRNA in network analysis.

Functional Analysis of Co-Expressed DEmRNAs with the Optimal Diagnostic lncRNAs

GO and KEGG pathway enrichment on the co-expressed DEmRNAs with the optimal diagnostic lncRNAs was performed. The David database was used for enrichment analysis, and the results were visualized using the ggplot2 package of R software.

Prediction of ccRCC-Associated Drugs and Molecular Docking

The drug prediction of ccRCC was based on Drug-Gene Interaction Database (DGIdb), which is a web source that helps discover drug–gene interactions. The key genes were imported into DGIdb (https://dgidb.org/) database to screen out potential drugs. In this study, the final drug list included only drugs that were approved by the Food and Drug Administration (FDA).

Molecular docking is a computational approach that “dock” small molecules into the structures of macromolecular targets and “score” their potential complementarity to binding sites. RSCBPDB is a database, which has Countless gene targets. The target was imported into the RSCBPDB (http://www.rcsb.org/pdb/home/home.do) database to download the 3D structures of target protein, and the 3D structures of the drugs were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The protein and drug molecules were pretreated in AutoDockTools-1.5.6, a molecular docking software based on the principle of semi-flexibility,¹⁷ respectively, and molecular docking is calculated. The binding energy was used as a reference for docking results to screen the most active ligand molecules and target genes. The smaller the binding energy, the more stable the binding between ligand and receptor, and the binding energy less than 0 means that ligand and receptor can spontaneously bind. Via PyMOL, a software for pre-docking of small-molecule components and proteins,³ the high-quality 3D structures of small molecules and proteins were created, and the corresponding protein residues and binding bonds were displayed. It is generally believed that the binding capacity is stronger when the binding energy was less than −5.0 kJ/mol.¹⁸

Results

Identification of DElncRNAs and DEmRNAs

Analysis results showed that a total of 58 (27 up-regulated and 31 down-regulated) DElncRNAs (Figure 1A) and 660 (207 up-regulated and 453 down-regulated) DEmRNAs (Figure 1B) were identified in ccRCC. The top 10 DElncRNAs and DEmRNAs are listed in Table 2.

Table 2 Top 10 DElncRNAs and Top 10 DEmRNAs in ccRCC

Figure 1 Top100 DEmRNAs and all DElncRNAs in ccRCC. (A) All DElncRNAs demonstrated by heat map (left) and volcano plot (right); (B) the top 100 DEmRNAs demonstrated by heat map (left) and volcano plot (right).

GO Enrichment and KEGG Pathway Analysis of DEmRNAs

GO enrichment analysis showed that DEmRNAs were enriched in multiple biological processes, including the immune response, metal ion transport, cell adhesion, the adaptive immune response, inflammatory response, and T cell activation (Figure 2A). In addition, GO analysis of cell component demonstrated DEmRNAs were mainly distributed in the extracellular area and the plasma membrane. GO analysis of molecular function indicated that these DEmRNAs were mainly involved in the transport of metal ions and the transport of proteins across membranes. KEGG pathway analysis showed that DEmRNAs were significantly involved in cytokine and receptor interactions, primary immunodeficiency, PPAR signaling pathway, and arachidonic acid metabolism pathway (Figure 2B).

Figure 2 GO enrichment analysis and KEGG pathway analysis of DEmRNAs in ccRCC. (A) GO enrichment analysis; (B) KEGG pathway analysis.

Identification of the Optimal Diagnostic lncRNAs

Totally, 12 DElncRNAs were obtained after LASSO selection. The correlation of these 12 DElncRNAs is shown in Figure 3A. Twelve DElncRNAs were sorted according to the average decreasing precision (Figure 3B). The result of 10-fold cross-validation showed that when the number of lncRNA was 4, the accuracy was the highest (Figure 3C). Therefore, the top 4 DElncRNAs were selected as the optimal diagnostic lncRNAs. Hierarchical clustering of 4 lncRNAs is shown in Figure 3D. The classification model was built using a random forest package, support vector machine, and decision tree of the 4 optimal diagnostic lncRNAs. The AUC of the random forest model was 0.991, with a specificity of 95.8% and a sensitivity of 99.4%. The AUC of the decision tree model was 0.951, with a specificity of 91.7%, and a sensitivity of 98.3%. The AUC of the support vector machine model was 0.983, with a specificity of 97.2% and sensitivity of 99.6% (Figure 4).

Figure 3 Determining process of the optimal diagnostic lncRNAs for ccRCC. (A) Heat map of DElncRNAs selected by LASSO; (B) average decreasing precision sorting of 12 DElncRNAs; (C) accuracy changes with increasing number of lncRNA; (D) heatmaps for the 4 optimal diagnoses of lncRNA.

Figure 4 ROC curve of decision tree (DT) model (A), random forest (RF) model (B) and support vector machine (SVM) model (C).

Correlation Analysis Between the Optimal Diagnosis lncRNAs and Clinical Data

In ccRCC tissues, the expression of 3 lncRNAs, FAM242C, LINC01762, and LINC02038 (Figure 5A) was significantly lower than that of adjacent non-tumor tissues, while the expression of PVT1 was significantly higher than that of adjacent non-tumor tissues, whether in the early or late stage of ccRCC. These results were consistent with the analysis results of our clinical samples to verify (Figure 5B). The “maftools” function in R was used to analyze the mutation data of ccRCC, and the mutation of 4 lncRNA and top 10 key mRNA is displayed in Figure 5C. Results showed that only SIGLEC8 and CA9 had mutations in a small number of samples. Survival analysis (Figure 5D) showed that PVT1 was significantly correlated with the survival rate of ccRCC (P= 0.00019), revealing that PVT1 could be used as a prognostic indicator of ccRCC.

Figure 5 The expression of 4 lncRNAs in TCGA data set and collected samples. (A) expressions of FAM242C, LINC01762, LINC02038 and PVT1 in different stages; (B) expression of FAM242C, LINC01762, LINC02038 and PVT1 in clinical samples (up: fold change; down: log2 fold change); (C) waterfall plot of mutation analysis of 4 lncRNAs and top 10 mRNAs; (D) survival analysis of different PVT1 expressions in TCGA data set. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. Abbreviation: ns, not significant.

Functional Analysis of Co-Expressed DEmRNAs with the Optimal Diagnosis lncRNAs

Totally, 353 DEmRNAs were co-expressed with 4 optimal diagnosis DElncRNAs, resulting in 583 co-expression relationship pairs. The interacting network is shown in Supplementary Figure 1. Among which, PVT1 was the most interacted mRNA (degree = 325), followed by LINC01762 (degree = 220, Supplementary Table 2). These 353 DEmRNAs were annotated by GO and KEGG analysis. GO enrichment analysis showed that co-expressed DEmRNAs were enriched in different biological processes, including the immune response, ion transport, cell proliferation, the adaptive immune response, regulation of lymphocyte activation, and T cell activation (Figure 6A). GO analysis of cell component showed that DEmRNAs were mainly distributed in the extracellular area and the plasma membrane. GO analysis of molecule function indicated that these DEmRNAs were mainly involved in the transport of metal ions, the transport of proteins across membranes, and the binding of immunoglobulins. KEGG enrichment analysis showed that DEmRNAs were significantly involved in cytokine and receptor interactions, primary immunodeficiency, and other metabolic pathways (Figure 6B and Table 3).

Table 3 KEGG Enrichment Analysis of DEmRNAs in the Co-Expression Network

Figure 6 Functional enrichment of co-expressed DEmRNAs with the optimal diagnostic lncRNAs in ccRCC. (A) GO enrichment analysis of co-expressed DEmRNAs; (B) KEGG enrichment analysis of DEmRNAs.

Correlation Between Key lncRNAs, Key mRNAs and Key Immune Cells

Immune infiltration analysis showed that there was no significant difference between CD56bright natural killer cell and type 2 T helper cell in the normal group and the ccRCC group. Eosinophil, immature dendritic cell, and type 17 T helper cell had a high degree of infiltration in the normal group, and the other 19 immune cells had a high degree of infiltration in the ccRCC group (Figure 7A). Six immune cells were further screened by LASSO regression, including activated dendritic cell, eosinophil, immature dendritic cell, myeloid-derived suppressor cells (MDSCs), plasmacytoid dendritic cell and type 17 T helper cell. Further analysis explored the correlation between 4 key lncRNAs, 10 key mRNAs and 6 key immune cells. Among the 4 key lncRNAs, PVT1 was associated with 4 key immune cells. Other 3 lncRNAs were related to one or two immune cells, and the correlation between PVT1 and MDSC was the highest with a value of 0.33. Among the 10 mRNAs, SIGLEC8 was correlated with 6 key immune cells, and the correlation between SIGLEC8 and MDSC was the highest (value = 0.68). The correlation results were demonstrated as Cytoscape network (Figure 7B and C).

Figure 7 Correlation between 4 lncRNAs, 10 mRNAs and 6 immune cells. (A) The immune infiltration between normal group and ccRCC group; (B) correlation network between 4 lncRNAs and 6 immune cells; (C) correlation network between 10 mRNAs and 6 immune cells. The blue diamond represents mRNA or lncRNA, the Orange circle represents immune cells, the solid line represents a positive correlation, and the dotted line represents a negative correlation. *P < 0.05; ***P < 0.001; ****P < 0.0001. Abbreviation: ns, not significant.

Validation Analysis

Verified dataset was downloaded from GSE167573, which included tissue samples from 14 normal controls and 63 patients with RCC. The Wilcoxon test was adopted to analyze the expression of key 4 lncRNAs and 10 mRNAs in the RCC group and the control group. In lncRNA, it was found that PVT1 (Figure 8A) was significantly overexpressed in the RCC group. GABRD, SIGLEC8 and CDKN2A (Figure 8B–D) were significantly overexpressed, while ESRRB, ELF5 and UMOD (Figure 8E–G) were significantly down-regulated in the RCC group, which was consistent with the expression trend in our analysis. In addition, survival analysis was conducted for these genes in the validation set. According to the median expression value of the related genes in each sample, they were divided into high-expression and low-expression groups. Kaplan–Meier was used in R to analyze the impact of these genes in the validation set for survival. Among these genes, we found that the overall survival time of patients with high expressed ELF5, ESRRB, and UMOD (Figure 9A–C) was significantly higher than that those with low expression; while in SPAG4 and CDKN2A (Figure 9D and E) genes, the overall survival time of patients in high-expression group was significantly lower than that in low-expression group.

Figure 8 The expressions of key lncRNAs and key mRNAs in GSE167573 dataset. (A–G) expressions of PVT1, SIGLEC8, CDKN2A, ESRRB, ELF5 and UMOD. *P < 0.05; ****P < 0.0001. Abbreviation: ns, not significant.

Figure 9 The survival analysis of key lncRNAs and key mRNAs in GSE167573 dataset. (A–E) survival analysis of ELF5, ESRRB, UMOD, SPAG4 and CDKN2A. Expressions are divided into high expression and low expression groups according to the median expression value.

Prediction of ccRCC-Associated Drugs and Molecular Docking Between Drugs and Key Genes

The aforementioned results prompted us to focus on the therapeutic utility of the key genes. Thus, 6 key genes (GABRD, SIGLEC8, CDKN2A, ESRRB, ELF5, and UMOD) were imported into DGIdb database to obtain the candidate drugs. As demonstrated in Figure 10, four (CDKN2A, ESRRB, GABRD, and UMOD) of the 6 key genes were targeted by the promising 84 candidate drugs that were all approved by the FDA. Most candidate drugs might interact with the GABRD (61/84), followed by CDKN2A (21/84). The interaction types, sources, PMIDs of candidate drugs and key genes are listed in Supplementary Table 3.

Figure 10 Prediction of candidate drugs.

The impact of GABRD in ccRCC remains unclear, while the mutation of CDKN2A is related to poor prognosis of RCC.¹⁹ In order to explore the interaction relationship and mechanism between candidate drugs and targets, two drugs (ABEMACICLIB and RIBOCICLIB) and CDKN2A were selected for molecular docking to calculate their binding energy. The docking results are shown in Figure 11, the binding energy of the two drugs to CDKN2A was lower than −5.0 kJ/mol, and the binding was relatively stable. Among which, the binding energy of ABEMACICLIB and CDKN2A was −7.53 kcal/mol, formed four hydrogen bonds with ASP-105 and ASP-156 residues. The binding energy of RIBOCICLIB and CDKN2A was −6.25 kcal/mol, formed four hydrogen bonds with GLU-88, MET-53 and VAL-51 residues. The conformations with the lowest energy were selected as the optimal conformation by the dock binding free energy, and display the interaction mode between drug molecules and core targets (Figure 11).

Figure 11 The three-dimensional structures of the candidate drugs and key genes based on molecular docking.

Discussion

CcRCC is a common cancer of the urogenital system. Even with proper treatment, it is also prone to unpredictable metastasis. Unfortunately, the 5-year survival rate of ccRCC is less than 10%.²⁰ Accurate classification plays important roles in guiding the clinical diagnosis and treatment of tumors. Traditionally, the classification, grade, and stage of tumors are mainly judged according to histopathological characteristics, which are used to guide the formulation of treatment plans and predict the prognosis of patients. However, traditional tumor classification methods are not completely accurate in judging the prognosis and guiding the treatment of patients. Even if patients with the same classification adopt the same treatment scheme, the prognosis may vary greatly. The tumor is a complex disease with high heterogeneity at the molecular level. Tumor cells usually have abnormal changes in the genome, epigenetic group, transcriptome, and protein level. The diagnosis and even typing based on the molecular characteristics of tumor cells can be more objectively and accurately reflect the internal essential characteristics of tumor occurrence and development. The research on the internal molecular mechanism of tumors is also important for predicting the prognosis of patients and realizing personalized treatment. Various molecular markers of ccRCC play important roles in clinical practice. MiRNAs are involved in cell proliferation and death, migration, epithelial–mesenchymal transition, tumor invasion and metastasis, and play an important role in the diagnosis, treatment, and prognosis of RCC.²¹ Cochetti, G et al developed a 7p-urinary score based on urinary miR-122-5p, miR-1271-5p, and miR-15b-5p levels, which is a potential tool for the early diagnosis of ccRCC.^22,23 It has been reported that lncRNA is abnormally expressed in renal cell carcinoma, suggesting that it may act as an oncogene or tumor suppressor. In addition, lncRNA may be suitable for the diagnosis of ccRCC and even for predicting prognosis.^24,25 TCGA database facilitates the analysis of high-throughput data from different genomic changes, including non-coding RNA. Recently, researchers have developed multiple prognostic models and nomograms of ccRCC based on a combination of tumor stage, grade, subtypes, clinical features, and presentation status, which might provide reliable prognostic evaluation tools for clinician to assist in clinical treatment decisions.²⁶ The nomogram constructed by Zhang et al based on key genes can effectively predict the survival of ccRCC patients at 1, 3, and 5 years, but he only focused on large-scale mRNA sequencing data from the TCGA platform.²⁷ In this study, we analyzed the differential expression of lncRNA and mRNA in 530 ccRCC samples and 72 adjacent non-tumor samples. The DElncRNAs and DEmRNAs were detected by R toolkit, and the machine learning algorithm was used to identify the lncRNA biomarkers most suitable for diagnosis and prognosis of ccRCC. LINC02038, FAM242C, LINC01762, and PVT1 could be used for the optimal diagnostic lncRNA for ccRCC, and they demonstrated high specificity and sensitivity in 3 sets of models. In addition, lncRNAs and mRNAs co-expression analysis and drug prediction were performed, providing a new perspective for the treatment of ccRCC.

CcRCC is essentially a metabolic disease characterized by metabolic reprogramming involving multiple metabolic pathways such as glucose metabolism and pentose phosphate.^28,29 CcRCC metabolic reprogramming enables tumor cells to survive under energy deprivation and hypoxia conditions,³⁰ the metabolic flux through glycolysis is partitioned, and mitochondrial bioenergetics and oxidative phosphorylation are impaired, as well as lipid metabolism.^31–34 Notably, many lncRNAs such as PVT1 are involved in the regulation of metabolic pathways and can be used as diagnostic and prognostic biomarkers of ccRCC. PVT1 is located on human chromosome 8q24. In a meta-analysis of 9972 cancer patients, He et al³⁵ found that the expression of PVT1 in tumor tissues was higher than in paired non-cancer tissues and was related to a poorer prognosis in cancer patients. Consistent with their research, our analysis showed that the expression of PVT1 was significantly up-regulated in ccRCC and further increased with the increase of clinical stage. Its expression was negatively correlated with the OS of ccRCC patients. In addition, our study also demonstrated that the cytokine–cytokine receptor interaction pathway was one of the KEGG enrichment of DEmRNAs. Moreover, PVT1 was involved in this pathway. This is consistent with the current research that PVT1 is widely involved in the regulation of a variety of tumors.³⁶ Zhang et al³⁷ recently reported that PVT1 promoted ccRCC development by stabilizing hypoxia-inducible transcription factor 2α (HIF2α), one of the core hypoxia-inducible factors, and HIF2α further regulated the expression of downstream genes after entering the nucleus. The results of immune infiltration analysis showed that PVT1 was associated with 4 key immune cells. These results suggest that PVT1 may be involved in the disease process of ccRCC extensively. Linc02038 is located on human chromosome 3q12. At present, there are few reports on its role in tumors. In an RNA-Seq and miRNA-Seq data analysis, Chen et al³⁸ revealed LINC02038 and AC104823.1 were positively correlated with OS of patients with colorectal cancer. The function of LINC01762, FAM242C, and LINC02038 is not clear yet; however, they are related to one or two immune cells. This suggests that these lncRNAs may be involved in the process of immune regulation of ccRCC. The co-expression network showed that LINC01762 interacted more with DEmRNAs, suggesting that we can pay more attention to this lncRNA in later research to clarify its function or role in ccRCC. Our research also indicated that lncRNA PVT1 and mRNA SIGLEC8 were highly correlated with MDSC. MDSC plays an important role in tumor-mediated immune evasion. Guan et al³⁹ elucidated that the accumulation of MDSC subsets might be associated with cytokine-induced migration in RCC. Combined with our results, PVT1 may play important roles in the diagnosis and disease progression of ccRCC.

CA9 is a cell surface glycoprotein and a tumor-associated antigen overexpressed in almost all ccRCC following hypoxia.⁴⁰ CA9 catalyzes the conversion of carbon dioxide and water to bicarbonate and hydrogen ions, which promotes cell survival under hypoxic conditions.⁴¹ Our results showed that CA9 expression was significantly up-regulated, and this up regulation may be the result of transmembrane transport in response to up-regulated proteins. Expression of SPAG4 is up-regulated in RCC via the HIF and could be observed within the nucleolus.⁴² This is consistent with the reported mechanism of RCC;⁴³ hypoxia and hypoxia-related pathways play a critical role in the development and progress of RCC.⁴⁴ It is demonstrated that mutation of CDKN2A is related to poor prognosis of RCC.¹⁹ Our study showed that CDKN2A was up-regulated in ccRCC, while whether this high expression was caused by mutation needs to be further confirmed by follow-up studies. Consistent with the previous report,⁴⁵ we found that ELF5 was down-regulated in RCC, suggesting that it can be used as a tumor suppressor in the kidney. This finding deserves further research verification to explore its protective role in ccRCC. In contrast, the impact of GABRD, SIGLEC8, MFSD4A, ACP3, ESRRB, and UMOD in ccRCC remains unclear, which needs to be clarified in the follow-up study.

Abnormal blood vessels supplying and impaired blood flow are the main characters of the solid tumor. During the formation of tumors, cancer cells release various cytokines and growth factors around and reprogram many other types of cells to establish a tumor microenvironment.⁴⁶ These cytokines regulate the expression of downstream genes after entering the nucleus,⁴⁷ which further results in the occurrence of mutant cancer cells and impairs the effectiveness of the immune system in attacking cancer cells.⁴⁸ Our functional enrichment analysis showed that the DEmRNAs were mainly involved in the cytokine–cytokine receptor interaction pathway and immunodeficiency pathway. Also, the kidney is an important organ for energy metabolism and regulation in the human body. RCC will lead to abnormal regulation and metabolism of sodium in the kidney. In accordance with Qi’s report,⁴⁹ our results demonstrated that the DEmRNAs were rich in aldosterone-regulated sodium reabsorption. Therefore, we speculate that these mRNA enrichment pathways, including aldosterone-regulated sodium reabsorption, are mainly caused by ccRCC. On the other hand, our validation results did not show the survival difference of PVT1 in the high- and low-expression group. The sample size of the validation data set is relatively small, which may be the main reason for the absence of the survival difference. Thus, further studies are needed to identify possible molecular interactions between the current reported lncRNAs and mRNAs in ccRCC. In particular, DElncRNAs and DEmRNAs with unclear functions can provide more choices for the treatment of ccRCC.

In addition to the bioinformatics database analysis, drug repurposing was also investigated via structure exploration, hoping to provide a new perspective for disease diagnosis, treatment and research. Molecular docking is computational approaches that dock small molecules into the structures of macromolecular targets and scores their potential complementarity to binding sites through computer simulation.⁵⁰ In this study, DGIdb database was used to identify several potential drugs with therapeutic efficacy on ccRCC key genes, and molecular docking was used to explore the interaction and mechanism between candidate drugs and targets. A total of 84 candidate drugs were screened out by 4 key genes (CDKN2A, ESRRB, GABRD, and UMOD). Among candidate drugs, Cyclosporine A has been shown to reverse multidrug resistance in renal cell carcinoma cell lines in vitro.⁵¹ ABEMACICLIB causes decreased cellular viability, increased apoptosis, and alterations in autophagy in renal cell carcinoma cell lines.⁵² RIBOCICLIB enhanced chemotherapy and immunotherapy of RCC by inhibiting CDK4/6-Cyclin D-Rb pathway.⁵³ ABEMACICLIB and RIBOCICLIB were selected for molecular docking with CDKN2A, and the results showed that two drugs had stable binding affinity with CDKN2A. The bioinformatics database analysis and molecular docking were used to provide a novel perspective for ccRCC.

Our study has some limitations. Firstly, the sample size of the validation dataset is relatively small, which may be the main reason why there is no survival difference between the high-expression and low-expression groups of PVT1. Therefore, further studies are needed to identify possible molecular interactions between the current reported lncRNAs and mRNAs in ccRCC. Furthermore, this study identified several potential drugs with therapeutic effects on key ccRCC genes, and additional studies and clinical trials are needed to identify and explore the drugs that are effective for ccRCC treatment in future.

Conclusion

In summary, our study constructed machine learning models in ccRCC by machine learning method and found that the survival of patients with ccRCC was related to lncRNA PVT1. Furthermore, DEmRNAs were mainly involved in hypoxia and hypoxia-related pathways. The lncRNA-mRNA co-expression network demonstrated that cytokine–cytokine receptor interaction and immunodeficiency were two main KEGG enrichments of DEmRNAs in the co-expression network. A total of 84 candidate drugs were screened out by 6 key genes, of which ABEMACICLIB and RIBOCICLIB were selected for molecular docking with CDKN2A, with stable binding affinity with CDKN2A. In the future, the molecular mechanism of the DelncRNAs and mRNAs in ccRCC are needed to further explore for early detection, prognosis judgment, risk assessment and targeted therapy of ccRCC.