Identification of prognostic immune-related genes in rhabdoid tumor of kidney based on TARGET database analysis

Background: Malignant rhabdoid tumor of the kidney (RTK) is a rare and highly aggressive pediatric malignancy. Immune system dysfunction is significantly correlated with tumor initiation and progression. Methods: We integrated and analyzed the expression profiles of immune-related genes (IRGs) in 65 RTK patients based on the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. Prognostic related IRGs in RTK patients were analyzed using univariate and multivariate analysis, based on which a prognostic model with IRGs was constructed. Correlation analysis between the risk score of our model and tumor-infiltrating cell were also investigated. Results: Twenty two IRGs were significantly associated with the clinical outcomes of RTK patients. Gene ontology (GO) analysis revealed that inflammatory pathways were most frequently implicated in RTK. A prognostic model was constructed using 7 IRGs (MMP9, SERPINA3, FAM19A5, CCR9, PLAUR, IL1R2, PRKCG), which were independent prognostic indices that could differentiate patients based on their survival outcomes. Furthermore, the risk scores from our prognostic model was positively associated with cancer-associated fibroblasts (CAFs). Conclusions: We screened seven IRGs of clinical significance to distinguish patients with different survival outcomes. This may enhance our understanding of the immune microenvironment of RTK, and could use to design individualized treatments for RTK patients. Background: Malignant rhabdoid tumor of the kidney (RTK) is a rare and highly aggressive pediatric malignancy. Immune system dysfunction is significantly correlated with tumor initiation and progression. Methods: We integrated and analyzed the expression profiles of immune-related genes (IRGs) in 65 RTK patients based on the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. Prognostic related IRGs in RTK patients were analyzed using univariate and multivariate analysis, based on which a prognostic model with IRGs was constructed. Correlation analysis between the risk score of our model and tumor-infiltrating cell were also investigated. Results: Twenty two IRGs were significantly associated with the clinical outcomes of RTK patients. Gene ontology (GO) analysis revealed that inflammatory pathways were most frequently implicated in RTK. A prognostic model was constructed using 7 IRGs (MMP9, SERPINA3, FAM19A5, CCR9, PLAUR, IL1R2, PRKCG), which were independent prognostic indices that could differentiate patients based on their survival outcomes. Furthermore, the risk scores from our prognostic model was positively associated with cancer-associated fibroblasts (CAFs). Conclusions: We screened seven IRGs of clinical significance to distinguish patients with different survival outcomes. This may enhance our understanding of the immune microenvironment of RTK and could use to design individualized treatments for RTK patients.


INTRODUCTION
Malignant rhabdoid tumor of the kidney (RTK) is a rare and highly lethal malignancy that primarily affected infants and young children [1]. Previous studies have reported that 10-15% of patients with RTKs had primary central nervous system (CNS) disease, which is currently designated as atypical teratoid-rhabdoid tumors [2]. The morbidity associated with RTK is extremely rare, accounting for about 2% of renal tumors in children. The prognosis of RTK patients remains to be extremely poor, with an overall 5-year survival rate of no more than 20% to 25% [3][4][5]. The main obstacle to improve the survival rate of patients with RTK might be due to poor understanding of the biological characteristic and regulatory mechanisms that underlie this fatal disease.
AGING programmed death-1 (PD-1) or cytotoxic T-lymphocyte associated antigen 4 (CTLA-4) polymorphisms and the immune escape of cancer. Immune checkpoint inhibitors (ICIs), such as pembrolizumab and nivolumab, have shown apparent efficacy in reducing tumor growth. This is mainly done by interfering with PD-1 or programmed death ligand-1 (PD-L1) interaction, thereby limiting the immune escape of cancer cells [7]. Also, certain ICIs have been shown to greatly improve the prognosis of patients [8][9][10][11]. However, only a small subset of patients responds to these therapies. Therefore, it is urgent to screen patients who are sensitive to immunotherapy [9].
Accumulating studies have suggested that immunerelated genes (IRGs) showed association with the efficacy of immunotherapy and the prognosis of patients [10]. According to a previous study, high OX-40 expression in tumor immune infiltrates indicated a favorable prognosis in patients with non-small cell lung cancer [11]. Zhang et al. have reported that PKD1 is significantly overexpressed in the tissues and remarkably associated with dismal prognosis in osteosarcoma [12]. Ryu et al. have found that the expression of p16 showed association with tumor immune microenvironment and had favorable prognosis of head and neck squamous cell carcinoma [13]. Besides, Bai et al. have revealed that the expression of BRAF V600E showed positive association with PD-L1/PD-1 in papillary thyroid carcinoma (PTC) samples, suggesting that immunotherapies that target PD-L1/PD-1 might be effective in PTC patients with BRAF V600E mutation [14]. However, there are only few studies that focused on the roles of IRGs in RTK [15]. Therefore, in this study, the Therapeutically Applicable Research To Generate Effective Treatments (TARGET) database were searched to examine the relationships of multiple immune genes with the prognosis of RTK, and aimed to construct a new prognostic model of RTK with IRGs, verifying the validity of this model in RTK patients based on TARGET-RTK cohort.

DEGs from RTK and normal controls
In this study, 7002 DEGs (3927 up-regulated genes and 3075 down-regulated genes) were found in the RTK samples when compared with normal control samples. They had a | log FC ≥ 1 and a false discovery rate (FDR) < 0.05. The heatmap and volcano diagrams of all DEGs are shown in Figure 1A, 1C.

Identification of DEIRGs
A total of 478 DEIRGs, including 254 up-regulated IRGs and 224 down-regulated IRGs, were extracted from the DEGs. The heatmap and volcano diagrams of all DEIRGs were shown in Figure 1B, 1D. As expected, the GO analysis revealed that inflammatory pathways were the most frequently implicated ones, and the receptorligand activity, cytokine activity, and cytokine receptor binding, growth factor activity, and cytokine receptor activity were the most frequent biological processes in RTK (Figure 2A). Based on the KEGG pathways analysis, our results demonstrated that DEIRGs were particularly enriched in cytokine receptor interactions of signaling pathways (p< 0.001), ( Figure 2B).

Identification of prognostic DEIRGs
As RTK is a lethal disease, monitoring the disease outcomes is important for clinical management. So, the molecular biomarkers that could serve as viable prognostic indicators were identified. After univariate Cox regression analysis, 22 DEIRGs were shown to be significantly associated with the OS of RTK patients in the TARGERT cohort (p< 0.05), (Table 1).

TF (transcription factor) regulatory network
To explore the potential regulatory mechanisms that correspond to the clinical significance of DEIRGs, the correlation between TFs and PIGRs (prognostic related immune genes) was analyzed. Firstly, the mRNA levels of TFs in RTK samples and normal controls were analyzed. A total of 131 differentially expressed TFs were identified between the two tissue types (FDR < 0.05, |log2 FC| > 1). Next, the correlations between the mRNA levels of the 131 TFs and the 22 PDEIRGs, with a correlation coefficient > 0.4 and a p-value < 0.01 as the cut-off values were identified ( Figure 3A, 3B). A total of 73 TFs were shown to be significantly associated with the abnormal expression of PIRGs (P< 0.01). To better visualize the regulatory relationships, a TF-based regulatory network was constructed using the Cytoscape software Figure 3C. Among these, the three transcription factors including HNF4A, FOSL2, and GATA4 interacted with multiple immune-related genes.

Construction of immune-related genes prognostic model for RTK
Multivariate Cox regression analysis identified 7 genes that showed significant association with the clinical outcomes (Table 2), which were used to construct a prognostic model to separate the RTK patients into two groups with discrete clinical outcomes. The formula was as follows: The median value of the risk scores was used as the cutoff value to divide patients in the TARGET-RTK cohort into high-risk group (n =28) and low-risk group (n = 30). As shown in Figure 4A, the number of deaths was significantly greater while the OS was shorter in high-risk cases when compared to the low-risk group (p< 0.001). The 1-year OS rate in the high-risk group of TARGET-RTK cohort was 16.5%, while the corresponding rate for that of the low-risk group was 79.2%. The area under the ROC (AUC) value for the prognostic model was 0.915, which suggested an excellent potential for the prognostic model based on PIRGs in monitoring the survival ( Figure 4B). The risk scores of the patients in the TARGET-RTK cohort were ranked and then their distribution was analyzed ( Figure  4C). The survival status of each patient in the TARGET-RTK cohort was marked on the dot plot in Figure 4D, showing the heatmap of the expression patterns of the risk genes in the two prognostic groups ( Figure 4E).

Independent prognostic value of the risk model in the TARGET-RTK cohort
To explore whether the prognostic model acts as an independent prognostic index in RTK patients, univariate and multivariate Cox regression analyses were performed to assess the effectiveness of this model. As expected, univariate analyses indicated that the variables of the prognostic model and clinical stage showed significant association with the prognosis of RTK patients ( Figure 5). Multivariate analysis revealed that the prognostic model was independently associated with the OS in the TARGET-RTK cohort (p< 0.001), ( Figure 6). These results indicated that the prognostic model could be used independently for predicting the prognosis of RTK patients. Consistent with the clinical observation, the clinical-stage also acts as an independent prognostic indicator for RTK in multivariate analysis (p< 0.05).

Clinical utility of the prognostic model
To validate the clinical utility of our model in predicting the progression of RTK, the relationships were analyzed between our model (the risk score and risk genes) and clinical and demographic characteristics (age, gender, stage). The expression of MMP9 was lower in the advanced stage cases, while the risk scores were significantly higher in the advanced stage cases and female patients (Figure 7), and no difference was observed in other variables.

Tumor-infiltrating cell
To verify whether our prognostic index model could reflect the status of the tumor immune microenvironment in RTK patients, the relationships between the risk score and tumor-infiltrating cells in the TARGET-RTK cohort were analyzed. There were fewer B cells, CD4 + T cells, CD8 + T cells, and endothelial cells in the tumor samples, while cancer-associated fibroblasts (CAFs) and macrophages were more in the tumor samples than in the normal controls ( Figure 8). Furthermore, correlation analysis showed that the risk score was positively associated with CAFs, and no difference was observed between the risk score and other infiltrating cells ( Figure 9).

DISCUSSION
In the present study, integrated analysis of the TARGET-RTK cohort was conducted and the results revealed that 478 IRGs showed aberrant expression in RTK patients when compared to normal controls. Next, univariate and multivariate COX regression analyses were performed to explore the PIRGs, and the results showed that a total of seven PIRGs were used to construct an independent prognostic model. As expected, compared with the low-risk group, the number of deaths in the high-risk group was significantly greater, and the OS time was shorter. Furthermore, the expression profiles of tumorinfiltrating cells in RTK were analyzed. The results showed that there were more CAFs and macrophages in the tumor tissues than in the normal controls, while B cells, CD4 + T cells, CD8 + T, and endothelial cells were fewer than normal controls. Also, the relationship between the risk score and tumor-infiltrating cells was explored, and our results showed that the risk score was positively correlated with infiltrating CAFs.
Cancer research has uncovered it as a disease that involves a succession of alterations to the genome [16].
Our study explored the alterations to reveal the relationship between immunogenomic profiles and immune microenvironment and to uncover the potential clinical implications. GO analysis and KEGG indicated that these altered IRGs were mainly involved in cytokinecytokine receptor interactions, MAPK, and PI3K-Akt signaling pathways. Previous studies have revealed that cytokine and their receptors actively participate in various types of cancers, including RTK [17,18]. Some of these cytokines have been demonstrated to exert antitumor or tumor-promoting effects on rhabdoid tumor cell lines [18]. As such, these altered IRGs might serve as diagnostic markers or therapeutic targets for RTK. However, comprehensive and an in-depth exploration with regard to the role of these IRGs in RTK is needed.
The results of KEGG analysis indicated that PI3K-Akt and MAPK signaling pathways were active in RTK, which were consistent with that of the previous studies [19,20]. We hypothesized that inhibitors of PI3K-Akt and MAPK signaling pathways might be effective in RTK. Experiments and clinical trials are needed to further verify this conjecture.
Previous studies have revealed that TFs play an important role in the development and differentiation of AGING  immune cells and immune response by regulating IRGs [21,22]. In our study, 73 TFs in the transcription factor network interacted with PIRG, indicating that these PIRGs might be regulated by these 73 TFs in RTK. Among these, three TFs including HNF4A, FOSL2, and GATA4 interacted with multiple IRGs. Increasing studies have demonstrated that HNF4A, FOSL2, and GATA4 showed aberrant expression in various malignant tumors, playing a role in promoting or inhibiting malignancies [8,[23][24][25], and our study results were consistent with these studies. We hypothesized that the abnormal expression of these three TFs might be crucial for the development and progression of RTK. However, there might exist other underlying regulatory mechanisms that are independent of these TFs. The molecular regulatory mechanism of these IRGs requires further research.

AGING
Previous studies analyzing the prognostic factors are limited because of the low incidence of RTK. Tomlinson G et al. have found that young age at diagnosis showed association with dismal survival of RTK patients [5]. However, little is known about the clinical importance of IRGs in RTK. In the present study, 7 IRGs that were closely associated with the clinical outcomes of children with RTK were found, and a prognostic model was constructed by using these to accurately discriminate between patients with different survival outcomes. Moreover, the risk scores from our prognostic model were significantly higher in advanced stage cases and female patients. It is well known that the prognosis of RTK patients with advanced-stage remained very poor, and our results were consistent with the clinical observation [26]. However, multivariate Cox regression analysis indicated that gender was not an independent prognostic factor. This inconsistency might be attributed to the small sample size. Based on our research, we hypothesized that this model could be applied to identify high-risk RTK patients, enabling early, intensive interventions and regular monitoring of disease recurrence to improve the prognosis of patients. Moreover, as revealed in the previous studies that tumor-infiltrating immune cells are important determinants of the prognosis and response to therapy [27,28]. Our results indicated that ICIs might be effective for low-risk RTK patients based on our model. The comprehensive and integrated analyses of IRGs in RTK deepens our understanding with regard to their clinical significance and illuminates the underlying regulatory mechanism of these IRGs. There were several advantages in our study over other reports. Firstly, the ImmPort database used in this study was a specialized immunology database, allowing us to analyze as many IRGs as possible. Secondly, although previous studies have demonstrated that IRGs could be applied to predict the prognosis of patients with cancer, this is the first study to our knowledge that used such a great number of samples to comprehensively explore the clinical importance of IRGs in RTK patients. Thirdly, the novel immune-related prognostic model exhibited a prominent performance for predicting the OS based on the TARGET database. To our knowledge, this is the first prognostic model developed for RTK patients.
Studies have revealed that regulation of tumor immune microenvironment plays a crucial role in tumor progression and metastasis [29,30]. To characterize the tumor immune microenvironment status, tumorinfiltrating immune cells were investigated. There were more CAFs and macrophages, but fewer B cells, CD4 + T cells, CD8 + T cells, and endothelial cells in the tumor tissues than normal controls. In addition to the difference in the number of tumor-infiltrating cells, the relationships between tumor-infiltrating cells and risk scores were further evaluated. Also, CAFs in the tumor stroma were found to be positively correlated with risk score, but no difference was observed between the risk score and other infiltrating cells. This might be due to relatively high content of CAFs and fewer tumorinfiltrating cells in RTK tissues. Many studies have AGING     www.aging-us.com 5471 AGING suggested a prominent functional role of CAFs in cancer progression and metastasis over the past few years. CAFs in tumor microenvironment induce angiogenesis by secreting CAF derived stromal cell-derived factor 1 and recruiting bone marrow-derived endothelial cells to promote tumor growth [31,32]. CAFs prompted epithelial to mesenchymal transition (EMT) and invasiveness of adjacent cancer cells by producing ECM-degrading proteases such as the MMPs [33]. CAFs promoted immunosuppressive TME, which is conducive to the immune escape of tumor cells [34]. Our result revealed that CAFs might play a tumor-promoting role in RTK, and this was consistent with the function of CAFs in other tumors [35,36]. However, the specific role and regulatory mechanism of CAFs in RTK should be further confirmed by experiments.
Admittedly, our study also had some limitations. Firstly, the number of RTK tissue samples in the TARGET cohort was relatively small, which might in turn cause selection bias. Secondly, transcriptomic analysis could not reflect the global alterations of immune status.
Thirdly, the prognostic model we constructed was not validated in the prospective clinical trials. Additionally, the reliability of IRGs identified requires further validation in both in vivo and in vitro experiments.
In conclusion, a prognostic model using 7 IRGs that showed significant correlation with clinical outcomes was constructed. The risk score acts as an independent prognostic index that could distinguish patients with different survival outcomes, and this could be applied to predict patients who were sensitive to immunotherapy.

Database
The required original data were obtained from TARGET (https://ocg.cancer.gov/programs/target) database. A total of 77 samples were screened, which included 65 tissues from RTK patient and 12 samples of normal kidney tissues. All data were processed using the R software (https://www.r-project.org/) and normalized using the limma package.

Analysis of differentially expressed genes (DEGs) by R software
DEGs were identified using the limma package [37]. The absolute value of | log2 fold change (FC) | was set to >1, and the cutoff value was adjusted to a P-value of < 0.05, which was considered to be statistically significant. The FDR and P values were screened to obtain all DEGs. The ggplot2 and the heatmap package were used to plot volcano diagrams and heatmap of DEGs, respectively.

Identification of differential expression of immunerelated genes (DEIRGs)
A total of 2498 IRGs were retrieved from the ImmPort database (https://www.immport.org/home) [38]. The DEIRGs were obtained from the DEGs using the R software. The absolute value of | log2 fold change (FC) | was set to >1, and the P-value was adjusted to <0.05 as cut-offs for filtering the DEIRGs. The ggplot2 and the heatmap package were used to plot the heatmap and volcano diagrams of the DEIRGs, respectively.

Univariate and multivariate analysis of the DEIRGs
To explore the prognostic significance of each DEIRG, univariate Cox analysis was conducted to screen the IRGs with a significant prognostic value (P< 0.05) as the candidates of individual risk score. A multivariate Cox analysis was performed to identify robust IRGs (P< 0.05), showing association with the overall survival (OS) to build the IRGs prognostic model. The forest map was drawn with ggplot2 package.

Molecular characteristics of prognostic immunerelated genes (PIRGs)
DEIRGs that showed significantly relation to the OS of RTK patients were regarded as PIRGs. It is well known that the transcription factors (TFs) act as important biological regulators of gene function in response to various internal and external stimuli. To better understand the regulatory mechanism of these immune genes, it is necessary to identify some TFs in patients with RTK that have the potential ability to regulate these PIRGs. Cistrome Cancer is a comprehensive resource for predicting the targets and enhancer profiles of TFs in cancers. This database contains a total of 318 TFs, and these are considered as a precious resource for experimental and computational cancer biology research. After that, clinically relevant TFs were extracted to construct the regulatory network of clinically relevant IRGs and potential TFs.

Generation of individual risk score
To assess the risk in each patient, a new prognostic model based on the expression data by multiplying with the Cox regression coefficient was constructed. The prognostic model was used to measure the prognostic risk of each patient with RTK, and the median risk score of the cohort was used as the cut-off value to divide all RTK patients into two groups: the high and low-risk groups. High-risk score indicates a poor prognosis in RTK patients.

Analysis of tumor-infiltrating cell
EPIC is an online database that is designed to estimate the proportion of immune and cancer cells from bulk tumor gene expression data (https://gfellerlab. shinyapps.io/EPIC_1-1/) [39]. After analyzing the immune infiltrating cells, including B cells, CAFs, CD4 + T cells, CD8 + T cells, endothelial, macrophages, and natural killer (NK) cells, the correlation between the risk score and tumor-infiltrating cells was analyzed by the R software.

Statistical analysis
The R software cluster Profiler package was used to perform GO analysis for identifying the biological themes among the gene clusters. The Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to screen the pathway enrichment of IRGs using the cluster Profiler package. The survival ROC R software package was used to calculate the AUC of the survival ROC curve to measure the performance of the prognostic model. The differences among the clinical parameters were tested using independent t-tests. R software was used to perform all statistical analyses, and p< 0.05 was considered to be statistically significant.