{'fraction_non_alphanumeric': 0.042566525260316235, 'fraction_numerical': 0.027670651754724256, 'mean_word_length': 4.901849217638691, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Novel Immune-Related ceRNA Network and Relative Potential Therapeutic Drug Prediction in ccRCC

Renal cell carcinoma (RCC) is the third common solid tumor in the urinary system with a high distant metastasis rate. The five-year survival rate of RCC has reached 75%, benefiting from the emergence and update of multiple treatments, while its pathogenesis and prognostic markers are still unclear. In this study, we committed to explore a prognostic ceRNA network that could participate in the development of RCC and had not been studied yet. We screened nine immune-related hub genes (AGER, HAMP, LAT, LTB4R, NR3C2, SEMA3D, SEMA3G, SLC11A1, and VAV3) using data of The Cancer Genome Atlas Kidney Clear Cell Carcinoma database (TCGA-KIRC) through survival analysis and the cox proportional hazard model. Next, we successfully constructed a ceRNA network of two mRNA (NR3C2 and VAV3), miRNA (hsa-miR-186-5p), and lncRNA (NNT-AS1) for ccRCC based on numerous online bioinformatics tools and Cytoscape. Finally, we predicted five potential drugs (clemizole, pentolonium, dioxybenzone, Prestwick-691, and metoprolol) based on the above results.


INTRODUCTION
Renal cell carcinoma (RCC) is one of the most common tumors all over the world, accounting for 5% of all new cancer cases in men and 3% in women (Siegel et al., 2021). With the increasing development of treatment approaches, the five-year survival rate has reached 75% for RCC in all stages and 93% for patients with localized lesions (Posadas et al., 2017). Clear cell renal cell carcinoma (ccRCC) accounts for more than 80% of RCC in different pathological subtypes and is the main concerned type (Capitanio and Montorsi, 2016). It is known that more than 30% cases are accidentally discovered during daily health checkups, instead of the typical symptoms. It is necessary to find effective biomarkers for ccRCC patients.
Nowadays, there are numerous research studies to figure out the development of ccRCC and how to treat it. In addition to the classic and abundant lipid-related treatments researched (Xiao et al., 2019), studies on immune-related treatments are also in full swing (Sharma and Allison, 2015;Barata and Rini, 2017;Wei et al., 2018). It is meaningful to seek out new immune targets for therapy of ccRCC patients.
ceRNA (competing endogenous RNA) hypothesis indicates that some molecules, like long noncoding RNA (lncRNA), can compete with the same microRNA (miRNA) response elements, thereby affecting gene expression (Salmena et al., 2011). Various studies proved that the ceRNA network could affect cancer cell proliferation and migration and serve as suitable therapeutic targets (Martens-Uzunova et al., 2014;Qu et al., 2016;Zhai et al., 2017;Zhai et al., 2019). However, the immune-associated ceRNA network was little studied in ccRCC.
In this study, we constructed a new immune-related network and validated the prognostic value of each element of the network. The potential signaling ways were analyzed by online bioinformatic tools. Several therapeutic drugs for ccRCC patients were also predicted based on the three target genes we screened.

Data Collection and Differentially Expressed Gene Screening
The mRNA-seq data, lncRNA-seq data, and clinical information of ccRCC in The Cancer Genome Atlas Kidney Clear Cell Carcinoma database (TCGA-KIRC) were downloaded from the Xena Functional Genomics Explorer (Xiao et al., 2020). All 1793 immune-related genes were collected from ImmPort (https://www.immport.org) (Bhattacharya et al., 2018). There were 2483 rows of data that were included in the file we downloaded from ImmPort; then, the duplicate genes were eliminated and 1793 genes were left. We screened the differentially expressed genes (DEGs) of KIRC using the "limma" package with parameters of adjusted p-value < 0.05 and log 2 |FC| > 1. Then, the intersection of DEGs and 1793 immune-related genes were obtained using a Venn diagram.

Functional Enrichment, Interaction Network Analysis, and Hub Gene Identification
The Database for Annotation, Visualization, and Integrated Discovery (DAVID) was used to obtain the functional enrichment results of 586 DEGs. We selected the top15 enrichment results of the biological process (BP), cellular component (CC), and molecular function (MF) to present. The top 20 signaling ways of Kyoto Encyclopedia of Genes and Genomes (KEGG) results were shown. p-value < 0.05 was considered significant. The protein-protein interaction (PPI) network of all DEGs was pictured by Cytoscape (v.3.6.1) using the "string" plugin. The confidence score was set as 0.40. Then, the PPI network was further analyzed by Cytoscape using the "MCODE" plugin to figure out the important modules. The top five modules were finally selected. Next, the "CytoHubba" plugin was conducted to screen the hub genes from all 586 DEGs with the criteria of degree >10.

Survival Analysis and Multivariate Cox Regression Analysis
Patients of ccRCC were divided into high-expression groups and low-expression groups according to the median of each gene expression. The "survival" package and R software were used to evaluate the prognostic value of the hub genes. GEPIA was utilized to verify the overall survival (OS) prognosis of the above hub genes (Tang et al., 2017). Then, multivariate cox proportional hazard regression analysis was conducted by SPSS (version 23.0). The value of p < 0.05 was significant. There were nine genes (AGER, HAMP, LAT, LTB4R, NR3C2, SEMA3D, SEMA3G, SLC11A1, and VAV3) finally screened after the above three steps. The Metascape database was used for GO enrichment analysis for hub genes (Zhou et al., 2019).

Prediction of miRNA and lncRNA and Construction of ceRNA
Starbase (Yang et al., 2011) and miRTarbase (Hsu et al., 2011) databases were used to predict miRNA reversely based on nine genes. The prediction results were merged and LAT was removed because there was no miRNA predicted. Cytoscape and plugin "CytoHubba" were conducted to obtain the top 15 miRNA according to the score of degree. The expression of miRNA was explored by UALCAN (Chandrashekar et al., 2017). The survival analysis of miRNA was identified by OncoLnc (http:// www.oncolnc.org/). Next, LncBase (Paraskevopoulou et al., 2013) and Starbase databases were used to predict lncRNA based on miRNA, both of which had different expressions and were related to patients' survival. There were 29 lncRNA in the intersection of prediction results, which were screened. Graphpad prism (v 8.2.1) was used to analyze the differential expression of lncRNA through an unpaired t-test, OncoLnc was used for survival analysis, and SPSS was used to conduct multivariate cox proportional hazard regression analysis. Finally, we constructed a new immune-related ceRNA network including NNT-AS1 (lncRNA), hsa-miR-186-5p (miRNA), NR3C2, and VAV3.

Gene Set Enrichment Analysis and Screening of Potential Drugs
LinkedOmics (Vasaikar et al., 2018) was utilized to conduct the Gene set enrichment analysis (GSEA) of two genes in the ceRNA network (NR3C2 and VAV3). Connectivity Map (Cmap) (Lamb, 2007) was a common small-molecular drug-prediction database, so we predicted potential therapeutic chemicals for ccRCC through three genes (NR3C2, VAV3, and HAMP) associated with miR-186-5p. PubChem  was used to find their structure and other information.

Immune Infiltration Analysis of VAV3 and NR3C2
In order to explore the relationship between immune-infiltration level and the two immune genes, we used CIBERSORT Kong et al., 2020;Meng et al., 2021a) (https:// cibersort.stanford.edu/) to obtain the immune infiltration level of TCGA-KIRC data. The algorithm was run using LM22 signature and 1,000 permutations. The 530 samples of TCGA-KIRC were divided into two groups according to the median of VAV3 or NR3C2. Then, R software and the "ggplot2" package were used to further present the immune cell infiltration level difference between two groups. The Wilcoxon test was used to calculate the p value. p value <0.05 was considered significantly.

Validation of Quantitative Real-Time Polymerase Chain Reaction
Paired normal and cancer tissues of ccRCC were obtained from Department of Urology, Union Hospital, Tongji Medical College, Wuhan, China, with the approval of the Ethics Committee of Huazhong University of Science and Technology. According to the instructions, we performed RNA Isolation and Real-time PCR analysis with SYBR Green mix (Thermo, Massachusetts, USA). Primers of NNT-AS1 and hsa-miR-186-5p were obtained from RiboBio (Guangzhou, China). Primers of VAV3, NR3C2, and normalized gene GAPDH were obtained from Sangon Biotech (Shanghai):

Determination of Immune-Related DEGs
We screened 2858 upregulated DEGs and 2354 downregulated DEGs in the TCGA-KIRC database, setting the parameters of p < 0.05 and log 2 |FC| >1 ( Figure 1A). Taking the intersection of DEGs and 1793 immune-related genes, 396 upregulated immune genes and 190 downregulated immune genes were finally determined and used for the following analysis ( Figure 1B).

Functional Enrichment Analysis, Construction of PPI, and Module Clustering of Immune-Associated DEGs
To determine the biological functions and potential signaling ways of DEGs, we used the DAVID database to perform GO and KEGG analysis of them. We selected the top 15 terms of GO and the top 20 terms of KEGG since there too many enrichment results. The results indicated that upregulated DEGs were enriched in biological processes (BP) of T-cell costimulation and interferon-gamma mediated signaling ways, while downregulated DEGs participated in cell proliferation and migration (Figures 2A,B). Cellular component (CC) enrichment results showed that both upregulated and downregulated DEGs were involved in the formation of various membrane and receptor complexes ( Figures 2C,D). Molecular function (MF) results suggested that upregulated DEGs were associated with the tumor necrosis factor receptor, tumor necrosis factor-activated receptor activity, and so on ( Figure 2E). Downregulated DEGs were related to transforming growth factor-beta receptor binding, fibroblast growth factor receptor binding, and so on ( Figure 2F). KEGG enrichment indicated that both groups participate in various important signaling ways in ccRCC development ( Figures 2G,H). We used Cytoscape and the "string" plugin to construct the PPI network of all 586 DEGs. Then, "MCODE" plugin was utilized to divide the DEGs into several modules. We selected the top five modules for the following enrichment analysis. Genes in module 1 were associated with neutrophils, monocytes, lymphocyte chemotaxis, and positive regulation of angiogenesis ( Figures 3A,B). Genes in module 2 were involved in cell adhesion molecules, positive regulation of B-cell proliferation, and natural killer cell-mediated cytotoxicity ( Figures 3C,D). The TNF-signaling pathway and IL8 production were enriched in module 3 ( Figures  3E,F). The NF-kappa B signaling pathway, Rap1 signaling pathway, and JAK-STAT signaling pathway were enriched in module 4 ( Figures 3G,H). It was worth mentioning that the above three pathways were also enriched in the other four modules. DEGs in module 5 were related to the PI3-Akt signaling pathway and MAPK signaling pathway ( Figures 3I,J).

Determination of Hub Genes and Survival Analysis
First, we used Cytoscape and "CytoHubba" plugin to initially screen hub genes. There were 449 genes with degree >10 selected. Then, the OS rate of these genes was evaluated through the "survival" package. GEPIA was an online website by which we verified the prognostic value of the above genes. Multivariate cox proportional hazard regression analyses were further determined by the hub genes ( Table 1). In addition, there were nine genes (AGER, HAMP, LAT, LTB4R, NR3C2, SEMA3D, SEMA3G, SLC11A1, and VAV3) finally selected. However, since later we found that there were no miRNA prediction results for LAT, we decided to remove it ( Figures 4A-H). AGER, HAMP, LTB4R, and SLC11A1 were risk factors for ccRCC, while NR3C2, SEMA3D, SEMA3G, and VAV3 were protective factors for ccRCC. We obtained the GO enrichment results of the eight hub genes using the meta scale database. Interestingly, there were two clusters in which one (AGER, HAMP, SLC11A1, and VAV3) participated in macrophage activation ( Figure 4I), and the eight hub genes were associated with cell proliferation and some immune processes ( Figure 4J).

Construction of the ceRNA Network Based on the Prediction
We predicted 673 miRNAs through Starbase and miRTarbase for all the eight hub genes ( Figure 5). Then, Cytoscape and "cytoHubbe" plugin were used again to obtain the top 15 miRNAs associated with the hub genes ( Figure 6A). Next, we explored the expression of miRNA using the UALCAN database and the prognostic value of them using the OncoLnc website. Then, hsa-miR-186-5p (alias, hsa-miR-186) was identified ( Figures 6B,C). The related genes NR3C2 and VAV3 were negatively associated with hsa-miR-186-5p, which was consistent with the known mechanism of miRNA.
As mentioned before, we predicted 1417 lncRNAs using Lncbase and 157 lncRNAs using Starbase for hsa-miR-186-5p. There were 29 lncRNAs in the intersection ( Figure 7A). We checked their expression and conducted survival analysis and multivariate cox proportional hazard regression analyses to determine the final lncRNA. In addition, according to the ceRNA hypothesis, lncRNA should be negatively related to miRNA. Finally, NNT-AS1 was the only one suitable. Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 12 | Article 755706 5 NNT-AS1 had a much lower expression level in ccRCC tissues than in normal tissues ( Figure 7B). The low expression group showed poorer survival ( Figure 7C). The results of cox regression were presented in the forest plot ( Figure 7D).
As a result, a new ceRNA network (NNT-AS1-hsa-miR-186-5p-NR3C2/VAV3) was constructed by us. What is more, GSEA results of NR3C2 and VAV3 indicated that both genes were associated with the cell cycle, cytokine-cytokine receptor interactions, the fatty acid metabolism, and the citrate cycle (TCA cycle), which made great significance in the development of ccRCC ( Figure 8A). To further verify the results from bioinformatics analysis, we performed qRT-PCR analysis on 30 pairs of clinical samples. As expected, the results showed that NNT-AS1, VAV3, and NR3C2 were remarkably downregulated in ccRCC compared to normal kidney tissues, while hsa-miR-186-5p was significantly highly expressed in ccRCC ( Figures 8B-E). Moreover, the correlation analysis results further suggested that NNT-AS1 was negatively related to miR-186-5p ( Figure 8F) but positively related to VAV3 and NR3C2 ( Figures 8G,H). In addition, MiR-186-5p was negatively associated to VAV3 and NR3C2 ( Figures 8I,J). In conclusion, the overall results were consistent with the hypothesis of the ceRNA network (NNT-AS1-hsa-miR-186-5p-NR3C2/VAV3).
To figure out the association between clinical characteristics and the two immune-related genes, we analyzed the mRNA expression levels of VAV3 and NR3C2 in different clinical subgroups. The classification standard was based on the previous literature (Meng et al., 2021b;Lou et al., 2021). The results showed that the expression of VAV3 had no significance between the elderly group (age ≤ 60 years) and young group (age > 60 years) (Supplementary Figures S3A,B). Besides, VAV3 expressed higher in females than in males. However, the expression level of VAV3 decreased with the increase of T stage, N stage, M stage, TNM stage, and G grades ( Supplementary Figures S3C-G). It may meant that VAV3 was not only downregulated in ccRCC but also decreased as the malignancy of tumors increased. Similarly, the expression level of NR3C2 was not correlated with age (Supplementary

Prediction of Potential Therapeutic Drugs
Using the Cmap database and three hub genes associated with hsa-miR-186-5p, we predicted potential drugs for ccRCC. NR3C2 and VAV3 were downregulated in ccRCC, while HAMP was upregulated. Drugs with negative connectivity scores were considered therapeutic. We finally selected five small molecular drugs in which enrichment < −0.7 and p-value < 0.01 (Figure 9; Table 2). Their structure and other information were presented through PubChem, in which Prestwick-691 could not be found.

DISCUSSION
Renal cancer is one of the top 10 cancers all over the world, and ccRCC makes up more than 75% of it. Although there had been much research studies to explore the pathogenesis and therapy for ccRCC, it still lacked suitable biomarkers and therapeutic targets. Since immunotherapy was considered a new hope for cancer therapy, we would like to figure out a novel immune-FIGURE 5 | Interaction network of eight mRNA and reverse-predict miRNA. Red cubes represented highly expressed mRNA, green cubes represent lowly expressed mRNA, and blue cubes represent miRNA.
Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 12 | Article 755706 related ceRNA network that might serve as a prognostic marker and therapeutic entry point for ccRCC ( Figure 10).
In this study, we took the intersection of DEGs of KIRC and 1793 immune genes to screen initially immune-related DEGs. Functional enrichment of all 586 DEGs was made, and we found they were associated with T-cell chemotaxis, B-cell chemotaxis, and multiple macrophage polarization-related signal pathways. The above biological processes and pathways were extensively studied and proved to be related to tumor development (Sousa et al., 2015;Chen et al., 2017;Choo et al., 2018;Locati et al., 2020). We further divided these genes into multiple modules and separately explored each modules' function. Interestingly, many classic and important pathways were enriched in various modules at the same time. Both the NF-kappa B signaling pathway (Hagemann et al., 2008;Mann et al., 2017) and JAK-STAT signaling pathway (Jiménez-García et al., 2018;Irey et al., 2019) had been proved to be key pathways in macrophage polarization.
Then, we predicted miRNA according to the hub genes and online databases. The top 15 miRNA were calculated by the "CytoHubbe" plugin and selected for the following analysis. The expression level and prognostic value of the 15 miRNA were explored again. As a result, hsa-miR-186-5p was determined. Consistent with the ceRNA hypothesis, highly expressed hsa-miR-186-5p acted as an inhibitor to silence downstream genes, which meant that low-expressed genes (NR3C2 and VAV3) in ccRCC could be the target genes in the ceRNA. Next, we predicted lncRNAs for hsa-miR-186-5p using LncBase and Starbase. There were 29 lncRNAs in the intersection, while only NNT-AS1 perfectly satisfied multiple conditions including survival analysis, cox regression analysis, and ceRNA hypothesis. Finally, we constructed a new immunerelated ceRNA network successfully, in which low-expressed NNT-AS1 downregulated NR3C2 and VAV3 to promote ccRCC through upregulating hsa-miR-186-5p. Although there had been no article proving the ceRNA network in ccRCC, the role of NNT-AS1 and the hsa-miR-186-5p axis had been explored in cervical cancer (Liu et al., 2020). The ceRNA network could also be associated with macrophage activation and polarization in the microenvironment. The Cmap database was used to predict potential drugs which might act on the ceRNA and treat ccRCC patients. Clemizole, pentolonium, dioxybenzone, and metoprolol are shown in Figure 9.
Our study had some limitations. We constructed a novel immune ceRNA network using various bioinformatics  analyses and verified the expression of NNT-AS1, has-miR-186-5p, VAV3, and NR3C2. Although the correlation analysis results of above ceRNA components meet the ceRNA hypothesis, whether NNT-AS1 could truly downregulate hsa-miR-86-5p and further upregulate VAV3 and NR3C2 needed further experiment verification. Besides, the influence of the ceRNA network to the tumor microenvironment needed further investigation. We will explore the ceRNA network based on the article.

CONCLUSION
In summary, we predicted a new ceRNA network for ccRCC, which could act as a prognostic biomarker and might contribute to the progression of ccRCC. Five potential drugs which might act on the ceRNA and treat ccRCC patients were predicted for future study.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
WL, XM, and HY performed the bioinformatics analysis and made equal contributions to the study. WX and XZ designed and supervised the study. All authors read and approved the final manuscript.