Nonsense-Mediated Decay Targeted RNA (ntRNA): Proposal of a ntRNA–miRNA–lncRNA Triple Regulatory Network Usable as Biomarker of Prognostic Risk in Patients with Kidney Cancer

The most prevalent subtype of renal cell carcinoma (RCC), kidney renal clear cell carcinoma (KIRC) may be associated with a poor prognosis in a high number of cases, with a stage-specific prognostic stratification currently in use. No reliable biomarkers have been utilized so far in clinical practice despite the efforts in biomarker research in the last years. Nonsense-mediated mRNA decay (NMD) is a critical safeguard against erroneous transcripts, particularly mRNA transcripts containing premature termination codons (called nonsense-mediated decay targeted RNA, ntRNA). In this study, we first characterized 296 differentially expressed ntRNAs that were independent of the corresponding gene, 261 differentially expressed miRNAs, and 4653 differentially expressed lncRNAs. Then, we constructed a hub ntRNA–miRNA–lncRNA triple regulatory network associated with the prognosis of KIRC. Moreover, the results of immune infiltration analysis indicated that this network may influence the changes of the tumor immune microenvironment. A prognostic model derived from the genes and immune cells associated with the network was developed to distinguish between high- and low-risk patients, which was a better prognostic than other models, constructed using different biomarkers. Additionally, correlation of methylation and ntRNAs in the network suggested that some ntRNAs were regulated by methylation, which is helpful to further study the causes of abnormal expression of ntRNAs. In conclusion, this study highlighted the possible clinical implications of ntRNA functions in KIRC, proposing potential significant biomarkers that could be utilized to define the prognosis and design personalized treatment plans in kidney cancer management in the next future.


Introduction
Kidney cancer is a potentially lethal disease. According to the International Agency for Research on Cancer's release of global cancer statistics for 2020, there were 179,368 fatalities and 431,288 new cases of kidney cancer worldwide [1]. The most prevalent histological type of kidney cancer is clear cell carcinoma, present in about 75% of cases, with an estimated 20-30% of patients presenting at diagnosis with distant metastasis and 20-40% of cases developing metastatic spread after primary surgical treatment [2,3]. Different

Definition and Identification of the Differentially Expressed ntRNAs That Are Independent of Genes and Detection of Differentially Expressed Genes (DEGs)
Based on the two versions of human genome annotation files that were downloaded from the NCBI (GRCh37.p13) and Ensembl (GRCh37.73), we defined transcripts as the ntRNAs if they met one of the following criteria: for NCBI transcripts, the one with the prefix "NR_" was considered as ntRNA; the one with the prefix "NM_" was considered as a normal transcript, which will produce protein. For Ensembl transcripts, one of three criteria must be met before the transcript was defined as a

Definition and Identification of the Differentially Expressed ntRNAs That Are Independent of Genes and Detection of Differentially Expressed Genes (DEGs)
Based on the two versions of human genome annotation files that were downloaded from the NCBI (GRCh37.p13) and Ensembl (GRCh37.73), we defined transcripts as the ntRNAs if they met one of the following criteria: for NCBI transcripts, the one with the prefix "NR_" was considered as ntRNA; the one with the prefix "NM_" was considered as a normal transcript, which will produce protein. For Ensembl transcripts, one of three criteria must be met before the transcript was defined as a ntRNA: (1) It is labeled with "not protein coding" by Ensembl; (2) it is labeled with "protein coding" by Ensembl, but it does not have a CCDS number; or (3) it is labeled with "protein coding" by Ensembl, but its length is less than 75% of the longest NM.
Due to AS, ntRNA may have ntRNA tags compared with NM ( Figure 2A). After ntRNAs were defined, we compared them with all the NM transcripts of the gene to define the ntRNA tags that are unique in the ntRNA. Then, we searched for ntRNA unique tags in the ASE dataset. In subsequent experiments and analysis, ASEs corresponding to ntRNA unique tags were used as ntRNA. Our study focused on the ntRNAs that were differentially expressed between tumor and normal samples, but the normal transcripts of the corresponding genes were not differentially expressed. We selected genes that were not differentially expressed using the 50th percentile of the false discovery rate (FDR) as the cutoff value. The differentially expressed ntRNAs were determined based on the 10th percentile of the FDR and |log 2 | > 1 as criteria. Similarly, we determined the differentially expressed miRNAs (DEmiRNAs) and differentially expressed lncRNAs (DElncRNAs) with thresholds of FDR < 0.05 and and |log 2 | > 1. Additionally, to investigate the potential biological pathways and mechanisms of ntRNAs, we performed a functional enrichment analysis of the genes containing ntRNAs in Metascape (https://metascape.org/gp/index.html#/main/step1 (accessed on 5 January 2022)) [27].

Establish a Triple Regulatory Network in KIRC
Following the theory that by competing with ntRNA for miRNA binding targets, lncRNA could indirectly regulate the expression of ntRNA [21], a triple regulatory network was created by performing the following procedures: (1) TargetScan is a widely used database that searches for conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA to anticipate biological targets for miRNAs [28]. In this paper, the targeted DEmiRNAs of the ntRNA were predicted Our study focused on the ntRNAs that were differentially expressed between tumor and normal samples, but the normal transcripts of the corresponding genes were not differentially expressed. We selected genes that were not differentially expressed using the 50th percentile of the false discovery rate (FDR) as the cutoff value. The differentially expressed ntRNAs were determined based on the 10th percentile of the FDR and |log 2FC| > 1 as criteria. Similarly, we determined the differentially expressed miRNAs (DEmiRNAs) and differentially expressed lncRNAs (DElncRNAs) with thresholds of FDR < 0.05 and and |log 2FC| > 1. Additionally, to investigate the potential biological pathways and mechanisms of ntRNAs, we performed a functional enrichment analysis of the genes containing ntRNAs in Metascape (https://metascape.org/gp/index.html#/main/step1 (accessed on 5 January 2022)) [27].

Establish a Triple Regulatory Network in KIRC
Following the theory that by competing with ntRNA for miRNA binding targets, lncRNA could indirectly regulate the expression of ntRNA [21], a triple regulatory network was created by performing the following procedures: (1) TargetScan is a widely used database that searches for conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA to anticipate biological targets for miRNAs [28]. In this paper, the targeted DEmiRNAs of the ntRNA were predicted using TargetScan8.0 (http://www.targetscan.org/vert_80/ (accessed on 20 September 2021)).
As determined by the expression profiles of the selected ntRNAs, DEmiRNAs, and DElncRNAs, the R package WGCNA (short for weighted gene correlation network analysis), a scientific methodology in biology that compares the gene association patterns between different samples, was employed to create the weighted gene coexpression network and examine the modules that are associated to the prognosis [30]. The module that significantly correlated with sample traits (patients' overall survival and stage, for example) was used for core gene selection. Using the Cytoscape (https://cytoscape.org (accessed on 10 October 2021)) plug-in cytoHubba [31], the hub triple regulatory network was then further identified. Only genes that ranked in the top 15 for all five methods, namely BC (Betweenness), Clo (Closeness), Degree, EPC, and RAD (Radiality) [32], were considered core genes.

Establishment of the Risk Prediction Model for the Hub Triple Regulatory Network
We combined the patient survival data for core genes with their expression, and then all samples were randomly split into training and test groups in a ratio 4:1. Following this, we used LASSO Cox regression to select the most appropriate variables from the hub triple regulatory network and finally obtained the risk score calculation formula: where Expression of Gene i denotes the patient's ith gene expression value, and Coefficient of Gene i denotes Gene i 's regression coefficient in the LASSO Cox regression model. Using the median risk score as the cutoff value, the samples were then split into high-and low-Riskscore gene groups. To compare the survival rates of high-and low-Riskscore gene groups, the Kaplan-Meier (KM) survival curve was plotted.

Establishment of Risk Prediction Model for the Immune Cells Associated with the Network
Kaplan-Meier survival curves and the log-rank test were performed to screen immune cells associated with the KIRC patients' overall survival. Furthermore, the differences between high-and low-Riskscore gene groups in each immune cell were examined using the Wilcoxon rank-sum test. Finally, we screened out immune cells that were associated with both the hub triple regulatory network and overall survival and constructed a LASSO Cox regression model to quantitatively describe the association between immune cells and patients' prognosis: where Expression of immune cell i denotes the patient's ith immune cell expression value, and Coefficient of immune cell i denotes immune cell i 's regression coefficient in the LASSO Cox regression model. Using the median risk score as the cutoff value, the samples were then split into high-and low-Riskscore immune groups. To compare the survival rates of high-and low-Riskscore immune groups, the Kaplan-Meier (KM) survival curve was plotted.
2.6. Using Core Gene and Immune Groups to Construct and Verify Cox Model Multivariate Cox regression analysis was conducted including two risk groups (both gene and immune) and clinical information (grade, age, and TMN). Cox regression prognostic factors were then incorporated into a calibration chart and nomogram design using the R computer package "rms", which would facilitate the application of our model to the clinical prognosis. To further assess the regression model's capacity to forecast patient survival, the C-index of the regression analysis model was established. The consistency between the anticipated survival and observed survival was examined using calibration curves. Finally, five other KIRC studies in the last 3 years were replicated on the same dataset, and the results were compared using receiver operating characteristic (ROC) curves [33][34][35][36][37].

Methylation Analysis of ntRNAs in Hub Triple Regulatory Network
DNA methylation affects the behavior of cancer cells through three DNA methyltransferases (DNMT1, DNMT3A, and DNMT3B). The Wilcoxon rank-sum test was conducted to examine the significant difference of the expression of each DNA methyltransferases between the high-and low-ntRNA groups. After methylation probe data were downloaded from UCSC Xena (https://xenabrowser.net/datapages/ (accessed on 15 November 2021)), the Pearson correlation coefficient was used to analyze the relationship between the methylation probe and ntRNAs, as well as the connections between the methylation probe and patient overall survival to find the key methylation sites.

Summary of Datasets
A total of 60,484 genes including 20,501 mRNAs, 14,353 lncRNAs, and 4077 miRNAs were found in The Cancer Genome Atlas (TCGA) database for 533 KIRC samples and 72 normal samples (Tables S1 and S2). We excluded genes that were not expressed in at least 60% of the samples, leaving 27,146 genes, namely 18,101 mRNAs, 8668 lncRNAs, and 377 miRNAs, for the study. The spectrum of ASEs was composed of 46,415 ASEs and involved 10,600 genes.
There were 530 cases with both clinical and immune information, and the samples were proportionately split into the training and test groups (4:1) at random. Table 1 displays the features of the samples in the training and test groups. In addition, 341 samples had methylation probe information.

Identification of Differentially Expressed ntRNAs and DEGs
A total of 102,309 ntRNA unique tags involving 15,045 genes were identified by our Python script, and 6720 ntRNAs involving 4598 genes were identified by the criteria (see method Section 2.2). log2FC and FDR techniques were utilized to detect differentially expressed ntRNAs independent of the corresponding genes, DEmiRNAs, and DElncRNAs; and, in total, 296 ntRNAs involving 263 genes, 261 DEmiRNAs, and 4653 DElncRNAs were identified in KIRC samples. The list of ntRNAs is shown in Table S3. In addition, volcano plots were drawn to display the distribution of ntRNAs, DEmiRNAs, and DElncRNAs ( Figure 2B-D).
We performed a functional enrichment analysis utilizing Metascape for 263 genes containing ntRNAs. They were highly enriched in proteolysis, which is involved in the cellular protein catabolic process, small molecule catabolic process, carbohydrate metabolic process, and so on ( Figure 2E). Therefore, ntRNAs play a vital role in a number of crucial biological pathways and mechanisms related to cancer.

Construction of the Hub ntRNA-miRNA-lncRNA Triple Regulatory Network in KIRC
The TargetScan database was firstly used to identify potential DEmiRNAs targeting transcripts associated with ntRNAs. Forty-six ntRNA-miRNA interaction pairs were identified, including 18 DEmiRNAs and 20 ntRNAs. These 20 ntRNAs had significantly higher expression in KIRC samples compared with normal samples. Secondly, the LncBase2 database was used to predict 18 DEmiRNAs and targets of these DEmiRNAs. A total of 2759 miRNA-lncRNA interaction pairs were discovered, containing 18 DEmiRNAs and 1624 DElncRNAs. Finally, we constructed a ntRNA-miRNA-lncRNA triple regulatory network containing 20 ntRNAs, 18 DEmiRNAs, and 1624 DElncRNAs with 2805 interactions (Tables S4 and S5 and Figure S1).
A WGCNA algorithm was used to select the most critical genes. First, we chose the soft threshold power 5 using the approximation scale-free topology criterion ( Figures 3A and S2). Then, a total of 15 non-gray modules were found by setting the minimum number of module genes to 10 and merging the modules with a similarity lower than 0.85. Next, we determined the Pearson correlation coefficient (Cor) between each module and clinical data, and we discovered that the turquoise module was the most closely associated with a patient's overall survival (OS) (Cor = 0.25, p = 1 × 10 −8 ) and OS time (Cor = −0.25, p = 5 × 10 −5 ). In addition, there was a significant positive connection between the stage and turquoise module (Cor = 0.18, p = 2 × 10 −5 ) ( Figure 3B). To sum up, we regarded the turquoise module as the most important module. Based on the turquoise module, a ntRNA-miRNA-lncRNA triple regulatory network including 13 ntRNAs, 5 DEmiRNAs, and 392 DElncRNAs was constructed ( Figure 3C).

Construction of a Prognostic Risk Score Model Based on the Hub Triple Regulatory Network
To prove that the genes in the hub triple regulatory network can be employed as possible biomarkers for the KIRC patients' prognosis, we performed LASSO Cox regression on 11 genes in the hub triple regulatory network in the training set and selected six genes that were strongly related to the patient prognosis ( Figure S3A). Finally, a risk score model was developed using the selected six genes to forecast the prognosis of the KIRC patients by the following formula:  In addition, the hub triple regulatory network was further identified using the Cytoscape plug-in cytoHubba. By comparing the intersection of the top 15 scores for each of the five methods (see method Section 2.3) (Table S6), three ntRNAs (KLC2_16990_AP, KLC2_16992_RI, and ABAT_33905_AP), four DEmiRNAs (hsa-miR-25-3p, hsa-miR-125b-5p, hsa-miR-148a-3p, and hsa-miR-499a-5p), and four DElncRNAs (RP6-159A1.4, LL0XNC01-7P3.1, CTD-3162L10.1, and AC004448.5) were identified to establish a hub triple regulatory network ( Figure 3D). The three ntRNA structures are shown in Figure 3E. The target sites of hsa-miR-125b-5p were predicted to pair with KLC2_16990_AP, KLC2_16992_RI, ABAT_33905_AP, AC004448.5, RP6-159A1.4, LL0XNC01-7P3.1, and CTD-3162L10.1 by TargetScan and LncBase2, and Figure 3F shows the binding sites between these genes.

Construction of a Prognostic Risk Score Model Based on the Hub Triple Regulatory Network
To prove that the genes in the hub triple regulatory network can be employed as possible biomarkers for the KIRC patients' prognosis, we performed LASSO Cox regression on 11 genes in the hub triple regulatory network in the training set and selected six genes that were strongly related to the patient prognosis ( Figure S3A). Finally, a risk score model was developed using the selected six genes to forecast the prognosis of the KIRC patients by the following formula: where the expression of each gene is indicated by the gene symbol, and the number before each gene symbol indicates the gene's regression coefficient in the optimal LASSO Cox regression model. After determining each patient's risk score using Formula (3), we divided the patients into high-and low-Riskscore gene groups with a cutoff value of the risk score's median. A KM plot revealed significant differences in survival between the high-and low-Riskscore gene groups ( Figure S4), confirming that the genes in the hub triple regulatory network we found are potential biomarkers for KIRC, which can be used to classify the prognostic risk of KIRC.

Correlation between Immune Infiltration and the Hub Triple Regulatory Network in KIRC
Nine immune cells were significantly connected with the OS of the KIRC patients, according to the results of the survival analysis ( Figure 4A). Then, we analyzed the immune infiltration levels between the high-and low-Riskscore gene groups by Wilcoxon rank-sum test, and 15 immune cells were differentially expressed under the standard p < 0.05 ( Figure 4B). Finally, eight immune cells (i.e., Tex, iTreg, Th1, Tem, CD8_T, Th2, Tcm, and Tgd) were both associated with OS and differently expressed in between the high-and low-Riskscore gene groups.

Establishment and Validation of a Nomogram for OS Prediction in KIRC
We performed a multivariate Cox regression that included our recorded risk groups (both Riskscoregene and Riskscoreimmune) and clinical information (grade, age, and TMN) from the samples in the training set. In the multivariate Cox model, Riskscoregene and Riskscoreimmune were associated with survival, proving that these two factors can be Furthermore, LASSO Cox regression was performed on the eight immune cells to investigate the most significant immune cells for the prognosis. Ultimately, a risk model was created using four immune cells that were highly related with the prognosis ( Figure S3B). Based on the immune infiltration of the four immune cells in the KIRC patients, we constructed the following formula to determine the risk scores: where the expression of each immune cell is indicated by the immune cell symbol, and the number before each immune cell symbol indicates the immune cell's regression coefficient in the optimal LASSO Cox regression model. After determining each patient's risk score using Formula (4), we divided the patients into high-and low-Riskscore immune groups with a cutoff value of the risk score's median. A KM plot revealed significant differences in survival between the high-and low-Riskscore immune groups ( Figure 4C), confirming that the four immune cells can be used as a prognostic biomarker for KIRC.

Establishment and Validation of a Nomogram for OS Prediction in KIRC
We performed a multivariate Cox regression that included our recorded risk groups (both Riskscore gene and Riskscore immune ) and clinical information (grade, age, and TMN) from the samples in the training set. In the multivariate Cox model, Riskscore gene and Riskscore immune were associated with survival, proving that these two factors can be used as prognostic biomarkers of KIRC ( Figure 5A). Then, a prognostic nomogram was developed to predict the survival rates for 1, 3, and 5 years of the KIRC patients according to the Cox model ( Figure 5B), which may provide a new clinical prognostic model reference for doctors. For 1-, 3-, and 5-year OS rates, the calibration curve of this nomogram showed a fair agreement between the forecast and actual observations ( Figure 5C), indicating that the generated nomogram had high accuracy for predicting OS in KIRC.
To further demonstrate the effectiveness of our multivariate Cox prognostic model, we compared our result with those of five other researchers' works on the same training and test sets (Table 2 and Figure 5D). Notably, our prognostic model had the highest area under the curve (AUC) for both the training set (AUC = 80.8%) and the test set (AUC = 86.3%) as compared with other prognostic models, including the models with a single mRNA [34,37], several DEGs [36], and a class of genes [33,35], which demonstrated the excellent predictive OS ability of our prognostic model; moreover, the novel findings provide biomarkers with a potential value for prognostic stratification of patients and may improve personalized approaches to treatment.

DNA Methylation Involved in Regulating the Expression of KLC2_16990_AP and KLC2_16992_RI
We explored the correlation between the expression level of three ntRNAs (i.e., KLC2_16990_AP, KLC2_16992_RI, and ABAT_33905_AP) contained in our hub triple regulatory network and the methylation status. First, a comparison of the expression of

DNA Methylation Involved in Regulating the Expression of KLC2_16990_AP and KLC2_16992_RI
We explored the correlation between the expression level of three ntRNAs (i.e., KLC2_16990_AP, KLC2_16992_RI, and ABAT_33905_AP) contained in our hub triple regulatory network and the methylation status. First, a comparison of the expression of three DNA methyltransferases (DNMT1, DNMT3A, and DNMT3B) in groups of ntRNAs with high and low expression revealed that the expressions of DNMT1, DNMT3A, and DNMT3B were significantly higher in KLC2_16990_AP, KLC2_16992_RI, and ABAT_33905_AP highexpression groups than the low-expression groups ( Figure 6A). Furthermore, three methylation sites (cg25736830, cg1529835, and cg07264522) in KLC2 DNA sequences were negatively correlated with two ntRNA expression levels and patients' survival state ( Figure 6B). The results suggested that some ntRNA expressions in KIRC may be dysregulated as a result of DNA methylation. The expression level of ABAT 33905 AP, however, did not significantly correlate with the methylation sites in ABAT (Table S7). The reason for the abnormal expression of ntRNAs in KIRC needs to be further studied.
patients' survival state ( Figure 6B). The results suggested that some ntRNA expressions in KIRC may be dysregulated as a result of DNA methylation. The expression level of ABAT 33905 AP, however, did not significantly correlate with the methylation sites in ABAT (Table S7). The reason for the abnormal expression of ntRNAs in KIRC needs to be further studied.

Discussion
In this study, we constructed a hub ntRNA-miRNA-lncRNA triple regulatory network to investigate the role of ntRNAs in KIRC, which could provide new perspective for the application of ntRNA data in clinical practice in urology. First, an in silico analysis was used to create a hub ntRNA-miRNA-lncRNA triple regulatory network with three ntRNAs, four DEmiRNAs, and four DElncRNAs. Then we analyzed the relationship between the network with a patient's prognosis and immune infiltration level through LASSO Cox regression, multivariate Cox regression, and correlation analysis. The results suggested that the genes in our hub triple regulatory network were independent prognostic factors of KIRC and significantly affected the immune infiltration expression. A prognostic nomogram was developed for a multivariate model that includes Riskscore gene and Riskscore immune groups and clinical information (grade, age, and TMN). The results demonstrated its ability to predict KIRC patients' prognosis with better performance as compared with other models [33][34][35][36][37]. These proved that our study identified novel biomarkers that could improve the prognostic stratification for KIRC patients. Finally, the results of methylation analysis showed that methylation influenced the expression of two ntRNAs in the hub triple regulatory network, suggesting that some ntRNAs' expression in KIRC may be dysregulated due to DNA methylation.
This study focused on the independently high expression of ntRNA to find some interesting results that cannot be found by analyzing gene expression. Many of the enriched terms of the 263 genes containing ntRNAs have important roles in cancer initiation and development ( Figure S5). Proteolysis involved in cellular protein catabolic process (GO: 0051603) has been reported in conjunctival melanoma and acute promyelocytic leukemia [38,39]. In another paper exploring hepatocellular carcinoma, small molecule catabolic process (GO: 0044282) was identified as the first significantly enriched gene ontology terms of DEGs [40]. It has been proven that PBRM1, the second-most often mutated gene in KIRC, affects the functioning of the carbohydrate metabolic pathway (GO: 0005975), a crucial biological process [41]. The process has also been reported in prostate cancer, lung cancer, and pancreatic ductal adenocarcinoma [42][43][44]. The level of glycerol in control renal tissue samples was 2.2 times higher than that in RCC, suggesting that glycerolipid metabolism (hsa00561) may be a key pathway in RCC [45], which was also confirmed by Lucarelli, G. et al. [46]. A recently identified form of cell death called ferroptosis (hsa04216) is caused by a buildup of iron-dependent lipid peroxide [47], and it has been observed that inhibiting ferroptosis in KIRC was found to inhibit tumor growth [48].
We constructed a hub ntRNA-miRNA-lncRNA triple regulatory network that statistically correlates with KIRC prognosis. Many genes in the network have been previously reported to be linked to kidney cancer growth and clinical behavior. For example, miR-125b-5p is in charge of maintaining the smooth muscle phenotype of juxtaglomerular cells [49]. In addition, miR-125b-5p was found elevated in RCC tissues and regulated RCC cell motility, invasion, and apoptosis, demonstrating that miR-125b-5p played a crucial role in RCC [50]. The expression of miR-148b-3p is upregulated in RCC, and miR-148b-3p is related with the DNA damage response, cell cycle progression, and apoptosis [51]. Meanwhile, miR-148b-3p directly binds to the 3 UTR of FOXK2 to inhibit FOXK2 from being expressed, and a poor prognosis is indicated by the low expression level of FOXK2 [52]. Similarly, the circPSD3/miR-25-3p/FBXW7 and miR-25-3p/IMPA2 axes increase KIRC cell invasion, motility, and the epithelial-mesenchymal transition (EMT) [53,54]. Heineman et al. noted that the expression value of miR-499a-5p was upregulated in RCC compared with benign renal tumors (BRTs) [55]. Interestingly, the gene ABAT acts as a tumor suppressor of KIRC [56], and ABAT_33905_AA expression in KIRC patients is significantly higher than that in normal samples, which deserves further study and exploration. Furthermore, based on the CTD (Comparative Toxicogenomics Database, https://ctdbase.org/ (accessed on 3 September 2022)) [57], we predicted 11 potential drugs, such as acetaminophen, doxorubicin, folic acid, celecoxib, and resveratrol, that can target five genes, namely ABAT, KLC2, hsa-miR-148a-3p, hsa-miR-25-3p, and hsa-miR-499a-5p, in KIRC, forming a gene-drug network ( Figure S6). Furthermore, we extended the ntRNA-miRNA-lncRNA network to show the genes that can be targeted by drugs, which can more clearly demonstrate the relationship between drugs and the network ( Figure S7). These drugs have been previously reported to be associated with KIRC. For instance, Cho et al. reported that taking acetaminophen raised the risk of significantly developing RCC [58]. In a xenograft model, doxorubicin combined with the TGase 2 inhibitor reduced the RCC growth [59]. Kabel et al. experimentally confirmed that administration of resveratrol can ameliorate KIRC [60]. Our study may offer helpful insight into the targeted treatment in KIRC. Our results provide advice for clinicians on drug choices in the treatment of KIRC patients.
Prior research has recognized the importance of immune infiltration for prognosis [61,62]. Moreover, it has recently been hypothesized that therapy, or combination of therapies, in cancer has to be adjusted to the tumor immunobiology; the potential use of immuno-based biomarkers, such as the ones we discovered, has been recently proposed in other urological tumors, such as bladder cancer [63]. In our study, eight of the 24 different types of immune cells that were associated with our hub triple regulatory network and overall survival of KIRC patients were explored. Vesicles with a diameter of 30 to 150 nm are known as tumor-derived exosomes (Texs). These vesicles are important mediums for communication between tumor immune cells [64], and Tex suppresses KIRC by inducing NK cell dysfunction [65]. The Th1/Th2 cytokines control the immune response to malignancy [66]. Patients with RCC showed a shift from Th1 to Th2 cells, indicating unbalanced Th1 and Th2 profiles [67]. CD8_T cells from normal donors recognized more antigens than those from patient sources and showed superior antitumor activity [68]. Therefore, these immune cells, which are linked to our hub triple regulatory network, have a significant impact on the development of KIRC.
To further explore the impact of our genes and immune cells on the prognosis of KIRC patients, we performed a multivariate Cox regression analysis. Our model included gene expression, patient characteristics, and immune infiltration at the same time, and the ROC curve of the model in the training set (AUC = 0.808) and test set (AUC = 0.863) achieves the best performance in comparison with other studies [33][34][35][36][37]. Moreover, few other models take into account the impact of gene expression, patient characteristics, and immune infection on patient prognosis at the same time, which may help to increase the model's robustness.
The multiclass feature data improve the ability of our prognostic model to predict OS, but they also limit our selection of an external validation set. We failed to find any additional KIRC datasets that contain ntRNAs, gene expression, immune infection, and clinical data. To address this problem, we used the KIRP (kidney renal papillary cell carcinoma) data from TCGA as an external validation dataset. We substituted patient stage data for grade data (e.g., stage I is regarded as grade I, and stage II is regarded as grade II) since the KIRP dataset lacks patient-grade information. The AUC of our prediction models for KIRP was 0.703, the highest AUC among other prediction models [33][34][35][36][37] (Figure S8), which further demonstrated its ability to predict the prognosis of kidney cancer. There were probably two reasons why AUC in KIRP was lower than that in KIRC (AUC = 0.808): (1) we used stage instead of grade for verification, and (2) KIRP is still not a perfect external validation set for KIRC.
Our study has a number of limitations. First, additional experiments on the role of core genes are required. Second, the function and method by which lncRNA affects miRNA regulation on ntRNA need to be further elucidated.

Conclusions
This study focused on ntRNA, a novel potential useful biomarker for kidney cancer management, which was previously defined as the nonsense-mediated decay targeted RNA. The results of pathway enrichment analysis suggested that ntRNAs were strongly associated with the clinical evolution of KIRC. Therefore, we constructed a ntRNA-miRNA-lncRNA triple regulatory network that showed to be linked to the degree of immune infiltration and prognosis of KIRC. The prognostic model based on the Riskscore gene and Riskscore immune groups and clinical information achieved significant advantages, compared with other previous models, proving that our data can be utilized to design reliable prognostic biomarkers and promising therapeutical targets. Furthermore, we performed methylation analysis of ntRNAs in the network, which proved that some ntRNAs were regulated by methylation, providing a direction for future studies on the causes of the abnormal expression of ntRNAs.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13091656/s1, Figure S1: Triple regulatory network in KIRC. Red squares indicate HEntRNAs, yellow V-shapes indicate DEmiRNAs, and blue circles indicate DElncRNAs; Figure S2: Analysis of the scale-free index for various soft-threshold powers β in WGCNA; Figure S3: (A) Cross-validation plot for gene LASSO Cox model tuning parameter selection and profiles of LASSO coefficients for 11 genes in the network hub triple regulatory network. (B) Cross-validation plot for immune LASSO Cox model tuning parameter selection and profiles of LASSO coefficients for eight immune cells; Figure S4: KM plot (high-Riskscore gene vs. low-Riskscore gene groups) for KIRC patients; Figure S5: 55 enrichment Gene Ontology (GO) and pathways of 263 genes containing HEntRNA, including 33 GO biological processes, nine reactome pathways, seven KEGC pathways, five WikiPathways, and one CORUM pathway; Figure S6: Gene-drug network for genes in the hub triple regulatory network; Figure S7: Extended hub triple regulatory network with drug information; Figure S8: ROC curve comparison with other studies in KIRP dataset; Table S1: List of 533 KIRC samples; Table S2: List of 72 normal samples; Table S3: List of 296 HEntRNAs; Table S4: ntRNA-miRNA interaction pairs from TargetScan database; Table S5: miRNA-lncRNA interaction pairs from LncBase2 database; Table S6: Core genes obtained by cytoHubba's five methods; Table S7: Correlation between ABAT methylation probes and ABAT_33905_AP expression in KIRC.  Conflicts of Interest: Nelson L. S. Tang is a director and hold shares of the company Cytomics Ltd. Cytomics Ltd. holds a licence to use the patent related to single cell type specific gene expression analysis.