Identification of a Novel Myc-Regulated Gene Signature for Patients with Kidney Renal Clear Cell Carcinoma

Given that myc was known to be a cancer-causing gene in several cancers including kidney renal clear cell carcinoma (KIRC). We aimed to construct myc-regulated genes (MRGs)-based prognostic signature. We obtained the mRNA expression and clinical data of KIRC from The Cancer Genome Atlas (TCGA) database and MRGs from the Molecular Signature Database (MSigDB). Then, a prognostic signature consisting of 8 MRGs (IRF9, UBE2C, YBX3, CDKN2B, CKAP2L, CYFIP2, FBLN5, and PDLIM7) was developed by differential expression analysis, cox regression analysis, and least absolute shrinkage and selection operator (lasso) analysis. Patients with KIRC were divided into high- and low-risk groups based on risk scores of MRGs-based signatures. Patients in the high-risk group showed inferior clinical characteristics and survival. In addition, the risk score was an independent prognostic factor for KIRC, and the risk score=based nomogram displayed satisfactory performance to predict the survival of KIRC. The MRGs-based signature is also correlated with immune cell infiltration and the mRNA expression of important immune checkpoints (IDO2, PDCD1, LAG3, FOXP3, and TIGIT). The tumor mutation burden (TMB) landscape between the high- and low-risk groups showed higher levels of TMB in the high-risk group than in the low-risk group and that higher levels of TMB predicted a poorer prognosis in KIRC. Furthermore, patients with KIRC in the high-risk group are more likely to experience immune escape. At last, we found patients with KIRC in the high-risk group were more sensitive to several chemotherapy drugs such as sunitinib, gefitinib, nilotinib, and rapamycin than patients with KIRC in the low-risk group. Our study successfully constructed and validated an MRGs-based signature that can predict clinical characteristics, prognosis, level of immune infiltration, and responsiveness to immunotherapy and chemotherapy drugs in patients with KIRC.


Introduction
Renal cell carcinoma (RCC) is global cancer that afects more than 300,000 people worldwide each year [1]. Among all RCCs, kidney renal clear cell carcinoma (KIRC) is the most common pathological subtype, accounting for more than 70% of all RCCs [2]. Although KIRC is common, its onset is insidious and the early clinical symptoms are mild or nonexistent, and once patients developed typical clinical symptoms such as hematuria, abdominal mass, and back pain, it often indicates that cancer has entered an advanced stage and the best time for treatment has been missed [3,4]. In addition, KIRC is an invasive tumor that can often develop distant organ metastases [5][6][7]. Delay in the diagnosis of KIRC and susceptibility of KIRC to metastasis contribute to the unsatisfactory prognosis of patients. Terefore, the identifcation of prognostic signatures and customization of treatment options are necessary to allow a satisfactory prognosis for patients with KIRC.
Myc is one of the most extensively studied oncogenes and is closely correlated with the initiation, maintenance, and progression of several cancers [8]. It encodes a protein that functions primarily as a transcriptional regulator of genes involved in the regulation of several cellular processes such as cell growth, cell cycle, cell diferentiation, apoptosis, angiogenesis, metabolism, and immune response [4,9].
Dysregulation of the expression of myc has been detected in several types of cancer including KIRC, and studies have demonstrated that myc played an important role in the progression of KIRC [10][11][12]. In view of the important role played by myc in KIRC progression, a comprehensive analysis of genes regulated by myc was performed in this study.
In the present study, we identifed diferentially expressed myc-regulated genes (MRGs) and explored their potential functions. Tese MRGs were further screened, and we fnally obtained 8 MRGs (IRF9, UBE2C, YBX3, CDKN2B, CKAP2L, CYFIP2, FBLN5, and PDLIM7) based on which we successfully constructed and validated a prognostic signature of KIRC. We uncovered the correlation of MRGs-based signature with the clinical characteristics and prognosis of KIRC. Based on the evidence that myc was associated with the regulation of immune responses in cancer [13][14][15], we also explored and unveiled the correlation of MRGs-based signature with immune cell infltration, immune checkpoint expression, tumor mutation burden (TMB), and immune treatment response in KIRC. In the end, we also found that the MRGs-based signature was signifcantly correlated with the sensitivity of KIRC to chemotherapeutic agents.

Identifcation of Diferentially Expressed MRGs.
Analysis of diferential mRNA expression of MRGs between KIRC samples and normal control samples was performed by running the limma R package. Diferentially expressed MRGs were defned as the diference in mRNA expression of MRGs in KIRC samples and normal control samples that met both |log 2 fold change (FC)| > 1 and false discovery rate (FDR) < 0.05.

Enrichment Analysis of MRGs.
Gene Ontology (GO) analysis focuses on the molecular function, biological process, and cellular component of the gene product; Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis focuses on the metabolic pathways in which the gene product is involved. In order to explore which functions and metabolic pathways MRGs are mainly involved in, GO and KEGG analyses were performed by running cluster Profler R package, org.Hs.eg.db R package, and enrichplot R package.

Establishment and Validation of MRGs-Based Prognostic
Signature. Samples included in the construction and validation of the risk model were subject to the two criteria of [1] having complete mRNA expression of MRGs and [2] having a minimum overall survival (OS) time of no less than 30 days. Te samples eligible for inclusion in the model construction and validation were randomly divided into training and validation sets in a 7 : 3 ratio by running the caret R package.
First, to identify MRGs associated with the prognosis of KIRC in the training cohorts, univariate cox regression was performed. To avoid overftting of MRGs with the model, least absolute shrinkage and selection operator (lasso) regression was then performed to remove MRGs that overftted with the model. Later, multivariate cox regression analysis was performed to select risk MRGs independently from other factors. Te risk score for each KIRC sample can be estimated by the following formula: Risk score � (Coef1 * mRNA expression of gene1) (1) Te "Coef" is the coefcient in the lasso regression model. Te KIRC samples in the training sets can be divided into high-and low-risk groups based on the median value of risk scores for all KIRC samples in the training sets. Principal component analysis (PCA) was performed to downscale this multigene signature model. Te Kaplan-Meier survival analysis and time-dependent receiver operating characteristic (ROC) analysis were then performed to assess the ability of this signature in predicting the prognosis of KIRC. To evaluate the reliability and stability of the signature, the same analysis was performed on the test set and the total set. Te "glmnet," "survival," "survminer," and "time ROC" R packages were applied to perform these analyses.

Establishment of Risk Score-Based Nomogram.
Te correlation of MRGs-based signatures with the clinical characteristics of KIRC samples was frst explored. Next, Cox regression analysis was performed to analyze whether the risk score of MRGs-based signature was an independent prognostic factor for the KIRC sample. A nomogram predicting the overall survival (OS) of the KIRC sample at 1, 2, and 3 years was constructed based on the independent factors of KIRC. Ten, calibration curves were plotted to evaluate the performance of the nomogram in predicting KIRC for 1, 2, and 3 years of OS. Te "survival" and "rms" R packages were used to perform these analyses.
2.6. Immune Analysis of MRGs-Based Signature. Given the importance of the tumor microenvironment (TME) in tumorigenesis and progression, the ESTIMATE algorithm was applied to analyze the level of infltration of two important components of the TME, stromal cells, and immune cells.

Analysis of TMB and TIDE Based on MRGs Signature.
Te TMB in KIRC was evaluated based on this formula: TMB � (total count of variants)/(the whole length of exons). We analyzed the respective mutational profles of the KIRC samples in the high-and low-risk groups by applying the Maftools R package. Diferences in TMB between KIRC patients in high-and low-risk groups, diferences in Kaplan-Meier survival for KIRC patients with high and low levels of TMB, and diferences in Kaplan-Meier survival across subgroups stratifed by H−TMB + high risk, H−TMB + low risk, L−TMB + high risk and L−TMB + low risk were further analyzed by applying the limma, ggpubr, survival, and survminer R packages.
In addition, the diference in the therapeutic efcacy of immune checkpoint inhibitors between the high-and lowrisk groups was obtained by calculating the diference in tumor immune dysfunction and exclusion (TIDE) scores in the TIDE database (http://tide.dfci.harvard.edu/) between the high-and low-risk groups.

Drug Sensitivity Analysis.
To further investigate the sensitivity of KIRC to common chemotherapeutic drugs in the high-and low-risk groups, we searched the Genomics of Drug Sensitivity in Cancer (GDSC, http://www. cancerRxgene.org) database combined with the application of the R package to compare the diferences in halfmaximal inhibitory concentration (IC50) values.

Statistical Analyses.
All statistical analyses in the present study were performed by R (version 4.0.0). Wilcoxon ranksum test was used to compare gene expression diferences between KIRC samples and normal control samples. As for the correlation analysis, the Pearson correlation coefcient was adopted when the data met the three conditions of continuous variables, normal distribution, and linear relationship, otherwise, the Spearman correlation coefcient was used. P values less than 0.05 were considered statistically diferent.

Identifcation of the Diferentially Expressed MRGs.
We obtained 480 MRGs from the MSigDB database, and the mRNA expression levels of MRGs were obtained by comparing with RNA-seq data downloaded from the TCGA database including 539 KIRC samples and 72 normal control samples. By diferential expression analysis, we fnally successfully identifed 94 diferentially expressed MRGs including 18 downregulated MRGs and 76 upregulated MRGs based on the screening criteria of |log 2 FC| > 1 and FDR < 0.05. Te results of this diferential expression analysis are presented in the form of a heat map and a volcano map, respectively (Figures 1(a), 1(b)).

Functional Enrichment Analysis of Diferentially
Expressed MRGs. To uncover the potential molecular functions and mechanisms of diferentially expressed MRGs, we performed GO and KEGG pathway enrichment analysis on diferentially expressed MRGs. According to the low to high p values and high to low correlation coefcients, we found that the top three GO terms in biological process (BP) were response to oxygen levels, wound healing, and response to decreased oxygen levels. In the cellular component (CC)  category, collagen−containing extracellular matrix, apical part of the cell, and secretory granule membrane occupied the top three GO terms. Extracellular matrix structural constituent, growth factor binding, and G protein−coupled receptor binding were the top three GO terms in the molecular function (MF) category. Te results of the KEGG pathway enrichment analysis showed that diferentially expressed MRGs were closely associated with Epstein−Barr virus infection, human papillomavirus infection, and PI3K−Akt signaling pathway. Te results of GO and KEGG pathway enrichment analysis are shown together in Figure 2.    In this way, we can calculate the median value of risk scores for all KIRC samples in the model, and based on the median value, the KIRC samples can be divided into high-risk and low-risk groups. Te results of PCA analysis showed a signifcant discrete trend in the three-dimensional plane for the high-and low-risk groups ( Figure 4).

Construction and Validation of a Prognostic MRGs-Based
Te mRNA expression of the eight MRGs in the high-and low-risk groups was presented as a heatmap in Figure 5(a). Tere was a signifcant diference in OS between the high-and low-risk groups, as shown in Figure 5(b), where the Kaplan-Meier survival analysis between the high-and lowrisk groups showed that OS was signifcantly lower in the high-risk group than in the low-risk group (p � 2.331e − 15). Consistent with this, the survival status of the high-risk group (percentage of death status 54%) was also signifcantly worse than that of the low-risk group (percentage of death status 17%) ( Figure 5(c)). Te time-dependent ROC showed the accuracy of the risk score in predicting OS at 1, 2, and 3 years in the KIRC sample to be 0.743, 0.726, and 0.739, respectively ( Figure 5(d)). Te test set and the overall set were performed with the same analysis to verify the stability and reliability of the model constructed from the training set. MRGs incorporated into the constructed risk model exhibited

Correlation between MRGs-Based Signature and Clinical
Characteristics. We also explored the correlation between MRGS-based signatures and clinical characteristics of the KIRC samples. Te heatmap (Figure 8) showed that the mRNA expression levels of IRF9, UBE2C, YBX3, CKAP2L, and PDLIM7 were higher in the high-risk group than in the low-risk group, while the mRNA expression levels of CDKN2B, CYFIP2, and FBLN5 were lower in the high-risk group than in the low-risk group. In addition, an important fnding was that samples in the high-risk group possessed inferior clinical characteristics including higher grade, stage, T-staging, and M-staging in the TNM staging system than those in the low-risk group. Furthermore, the relationship between MRGs-based signature and the prognosis of KIRC in each clinical subgroup stratifed by age (≤65 years or >65 years), gender (female or male), grade (G1 + 2 or G3 + 4), stage (stage I + II or III + IV), T stage (T1 + 2 or T3 + 4), M stage (M0 or M1), and N stage (N0 or N1) was analyzed. Te Kaplan-Meier analysis showed that higher risk scores were correlated with worse prognosis in multiple subgroups including age ≤65 years or >65 years, female or male, G1 + 2 or G3 + 4, Stage I + II or III + IV, T1 + 2 or T3 + 4 stage, M0 or M1 stage, and N0 stage compared to lower risk scores (Figure 9). It should be noted that in the N1 subgroup (Figure 9), the risk score was not signifcantly associated with the prognosis of the KIRC sample, which may be due to the small number of KIRC samples in the N1 subgroup.   Figure 10(b)). In addition, we compared the accuracy of risk and traditional clinical characteristics in predicting prognosis in the KIRC. Te results showed that the accuracy of risk in predicting one-year OS in the KIRC sample was 0.743, which was only less accurate than that of To further predict the prognosis of KIRC, we constructed an MRGs-based nomogram. As shown in Figure 11(a), the nomogram predicted the OS of the KIRC sample at 1, 2, and 3 years, respectively. Te C-index of 0.747 implied that the nomogram had moderate accuracy in predicting KIRC 1year, 2-year, and 3-year OS. In addition, the calibration curves also showed that the 1-year, 2-year, and 3-year survival probabilities of the patients with KIRC predicted by the nomogram were consistent with the actual survival probabilities of the patients (Figures 11(b)-11(c)).

Diference of Immune Characteristics of MRGs-Based
Signature. To explore the correlations between MRGs-based signature and immune characteristics in the TME, we performed a series of analyses including the correlation analyses between MRGs-based signature and immune cells, immune checkpoints, and immune escape. Te results of stromal cell and immune cell infltration analysis in TME  Journal of Oncology showed that higher levels of immune cell infltration in the high-risk group than in the low-risk group in TME (p � 2.8e − 07) (Figures 12(a)-12(c)). Te results of immune cell infltration analysis based on TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC algorithms showed that the risk score of MRGs-based signature was closely correlated with the infltration of immune cells (Figure 12(d)). Enrichment analysis of immune-related functions in the high-and lowrisk groups showed that APC_co_stimulation, checkpoint, cytolytic_activity, HLA, Infammation−promoting, T_cell_coinhibition, T_cell_costimulation, and Type_-I_IFN_response were more active in the high-risk group than in the low-risk group (Figure 12(e)).

Diference between TMB and TIDE Based on MRGs
Signatures. We further explored mutational characteristics between high-and low-risk groups. Te TMB landscape in the high-and low-risk groups showed that the frequency of mutations in SETD2, BAP1, and MTOR was higher in the high-risk group than in the low-risk group (Figure 14(a)). In addition, the results of the correlation analysis between risk and TMB levels showed that the levels of TMB in the highand low-risk samples were close to being statistically different (p � 0.069) (Figure 14(b)). Te Kaplan-Meier survival analysis showed higher levels of TMB were associated with a worse prognosis for KIRC patients (p � 0.001) (Figure 14(c)). Results of survival analysis for subgroups stratifed by H−TMB + high risk, H−TMB + low risk, L−TMB + high risk, and L−TMB + low risk showed that patients with high levels of TMB and high risk have the worst prognosis, while those with low levels of TMB and low risk have the best prognosis (p < 0.001) (Figure 14(d)). Immune checkpoint inhibition therapy has become a hot treatment modality for cancer [16]. Te TIDE score is gaining popularity because it is more accurate than a single biomarker for predicting the efcacy of immune checkpoint blockade for cancer [17]. We analyzed the diference in TIDE scores between the high-and low-risk groups, and the results showed that the TIDE scores were signifcantly higher in the high-risk group than in the low-risk group (p < 0.01) ( Figure 15).

Analysis of MRGs-Based Signature and Chemotherapy
Drug Sensitivity. Te results of the correlation analysis between MRGs-based signature and chemotherapeutic drug sensitivity showed that the patients in the high-risk group had lower IC50 values compared to patients in the low-risk group for sunitinib, geftinib, nilotinib, rapamycin, mitomycin.C, paclitaxel, vinblastine, salubrinal, parthenolide, and metformin. However, the IC50 values for embelin and thapsigargin were lower in the low-risk group than in the high-risk group (Figure 16).

Discussion
In this study, we identifed 94 MRGs and analyzed the potential functions of these MRGs. Te results of the GO and KEGG pathway enrichment analysis displayed the complexity of the function of the MRGs signature, suggesting the application potential of the MRGs signature. Combining cox regression analysis and lasso regression analysis we successfully constructed a prognostic signature consisting of eight MRGs (IRF9, UBE2C, YBX3, CDKN2B, CKAP2L, CYFIP2, FBLN5, and PDLIM7). Te results of PCA and t-SNE analysis displayed that the signature received good dimensionality reduction, which implied the reliability of the risk model. Te results of the prognostic analysis of the samples in the test and total sets further increase the credibility of our risk model constructed using the training set. Te results of the correlation analysis between MRGsbased signature and clinical characteristics revealed that patients in the high-risk group were signifcantly associated with malignant clinical characteristics such as high grade, high stage, high T stage, and high M stage, which implied a poor prognosis for patients. In addition, the results of the Kaplan-Meier survival analysis for subgroups showed signifcant diferences in OS between patients in high-and lowrisk groups across multiple subgroups. Te results of these analyses implied the satisfactory performance of MRGsbased signatures to predict the prognosis of KIRC. Te results of univariate and multivariate Cox regression     IDO2  ICOS  PDCD1  CD70  LAIR1  CD28  CD40  VTCN1  CD160  TNFSF9  LAG3  HHLA2  BTLA  CD48  CD44  CD200R1  TIGIT  TNFSF4  KIR3DL1  TMIGD2  TNFRSF14  TNFSF14  LGALS9  TNFSF18  TNFRSF9  CD86  TNFSF15  CD244  TNFRSF25 Gene expression Risk low high analyses demonstrated that risk score could independently predict the prognosis of patients with KIRC, and risk scorebased nomogram and calibration curves suggested satisfactory accuracy of MRGs-based signature in predicting the prognosis of KIRC patients. Te TME is highly complex, and the study of TME can help provide new ideas for the treatment of tumors [18]. Stromal cells and immune cells are two important cell members of the TME [19]. Stromal cells have been demonstrated to be closely associated with the growth, metastasis, drug resistance of several cancers [20][21][22], and immune cells, depending on their type, can play an important role in fghting tumors and in promoting tumor progression or immune escape, respectively [23,24]. Analytical results of TME based on MRGs signature showed that higher risk scores were associated with higher levels of immune cell infltration, with no signifcant correlation between stromal cells and risk scores. Tis suggested that MRGs-based signature may afect the prognosis of KIRC by infuencing the immune cell landscape in TME. Te higher level of immune cell infltration and immune function enrichment scores in the high-risk group further suggested that the infltration of immune cells in TME may have contributed to tumor progression. Te results of the analysis of the diference in the expression of immune checkpoints between high-and low-risk groups implied an immunosuppressed status in the KIRC samples from the high-risk group. Studies have shown that IDO2 was associated with Bcell immunity and can regulate T-cell-related immunity by afecting B-cell intrinsic mechanisms [25,26]. PDCD1, also known as PD1, were known as programmed death ligands and receptors, respectively, and its combination with PD-L1 allows tumor cells to evade the body's immune surveillance [27]. FOXP3 is an important marker molecule of Tregs that   Figure 16: Diferential analysis of IC50 value of sunitinib, geftinib, nilotinib, rapamycin, mitomycin.C, paclitaxel, vinblastine, salubrinal, parthenolide, metformin, embelin, and thapsigargin between high-and low-risk groups. IC50: half-maximal inhibitory concentration. directly or indirectly regulates the activity and function of Tregs, and changes in its protein levels have been shown to be closely associated with a variety of human diseases including tumorigenesis and metastasis [28][29][30][31]. LAG-3 was demonstrated to negatively regulate T-cell function, and antibodies to LAG-3 could relieve the inhibition of T-cell function by Tregs [32]. TIGIT is expressed in various T-cell subsets (CD4+ T, CD8+ T, Tregs) and can suppress the body's innate and adaptive immunity through various mechanisms and is considered a promising target for immunotherapy [33]. Takamatsu et al. demonstrated that profles of LAG-3, TIM-3, and TIGIT were valuable for predicting the prognosis and TME of KIRC [34]. Te signifcant positive correlation of risk scores with these key immune checkpoints further confrmed that upregulation of the expression of immune checkpoint in the KIRC patients from a high-risk group may be responsible for their poor prognosis.
With the widespread use of high-throughput sequencing technology, TMB has become a marker for predicting the response to immune checkpoint blockade in several types of cancer [35]. Although, in general, higher levels of TMB in cancer can lead to more infltration of CD8+ T cells and thus contribute to a better prognosis by exerting an antitumor efect, KIRC has been shown to be cancer that challenged conventional thinking about cancer immunology based on this evidence that KIRC has a modest mutation burden but is responsive to immunotherapy and higher CD8+ T-cell infltration is usually associated with a poorer prognosis [36,37]. Our study displayed that patients with KIRC in the high-risk group had higher levels of TMB than those with KIRC in the low-risk group and that the higher the TMB the worse the prognosis for KIRC, which further validated the unusual relationship between the TMB and immunology in KIRC. SETD2, BAP1, and MTOR possessed higher mutation frequencies in the high-risk group than in the low-risk group implying that these three cancer-driven mutations may promote the progression of KIRC. Te diference in TIDE scores between KIRC in the high-and low-risk groups refected the efectiveness of immune checkpoint blockade, with higher TIDE scores indicating poorer immune checkpoint blockade, and higher TIDE scores for KIRC patients in the high-risk group than those in the low-risk group suggesting that patients in the high-risk group were more likely to experience immune escape, which was consistent with our analysis results that KIRC patients in the high-risk group had a worse prognosis than those in the lowrisk group.
Of these MRGs, some have been shown to be associated with KIRC progression. IRF9 was thought could predict the prognosis and immune characteristics of KIRC as a member of a multigene signature [38,39]. UBE2C was shown to be closely associated with the proliferation and invasion of KIRC and could predict the immune characteristics and prognosis of KIRC [40][41][42]. Jafri et al. identifed germline CDKN2B mutations as a novel causative agent of familial KIRC, suggesting an important role for CDKN2B in the initiation of KIRC [43]. CYFIP2 was predicted to be downregulated in KIRC and downregulation of expression of CYFIP2 was associated with poor prognosis of KIRC [44]. In a rat model, FBLN5 was found to be associated with the metastasis of KIRC [45]. Other MRGs were not reported in KIRC, so our future studies will focus on these MRGs.
Te results of the analysis of diferences in sensitivity to chemotherapeutic agents between high-and low-risk groups revealed that patients in the high-risk group benefted more from treatment with sunitinib, geftinib, nilotinib, rapamycin, mitomycin.C, paclitaxel, vinblastine, salubrinal, parthenolide, and metformin, while patients in the low-risk group benefted more from treatment with embelin and thapsigargin. Tese fndings provide a theoretical basis for individualized pharmacological treatment of patients with KIRC.
Te lack of prognostic information in all KIRC-related datasets in the GEO database resulted in our inability to validate our results with an independent external dataset, which is the main limitation of this study. Further studies need to be performed to explore the molecular mechanisms and biological functions of MRGs-based signatures.

Conclusion
In this study, a prognostic signature consisting of 8 MRGs (IRF9, UBE2C, YBX3, CDKN2B, CKAP2L, CYFIP2, FBLN5, and PDLIM7) in KIRC was successfully constructed and validated. Te MRGs-based signature can predict the clinical characteristics, prognosis, immune characteristics, and sensitivity to chemotherapeutic agents of patients with KIRC and has the potential to be applied in the clinical setting.