Identification and Validation of Cytotoxicity-Related Features to Predict Prognostic and Immunotherapy Response in Patients with Clear Cell Renal Cell Carcinoma

Background Clear cell renal cell carcinoma (ccRCC) is a renal cortical malignancy with a complex pathogenesis. Identifying ideal biomarkers to establish more accurate promising prognostic models is crucial for the survival of kidney cancer patients. Methods Seurat R package was used for single-cell RNA-sequencing (scRNA-seq) data filtering, dimensionality reduction, clustering, and differentially expressed genes analysis. Gene coexpression network analysis (WGCNA) was performed to identify the cytotoxicity-related module. The independent cytotoxicity-related risk model was established by the survival R package, and Kaplan–Meier (KM) survival analysis and timeROC with area under the curve (AUC) were employed to confirm the prognosis and effectiveness of the risk model. The risk and prognosis in patients suffering from ccRCC were predicted by establishing a nomogram. A comparison of the level of immune infiltration in different risk groups and subtypes using the CIBERSORT, MCP-counter, and TIMER methods, as well as assessment of drug sensitivity to conventional chemotherapeutic agents in risk groups using the pRRophetic package, was made. Results Eleven ccRCC subpopulations were identified by single-cell sequencing data from the GSE224630 dataset. The identified cytotoxicity-related T-cell cluster and module genes defined three cytotoxicity-related molecular subtypes. Six key genes (SOWAHB, SLC16A12, IL20RB, SLC12A8, PLG, and HHLA2) affecting prognosis risk genes were selected for developing a risk model. A nomogram containing the RiskScore and stage revealed that the RiskScore contributed the most and exhibited excellent predicted performance for prognosis in the calibration plots and decision curve analysis (DCA). Notably, high-risk patients with ccRCC demonstrate a poorer prognosis with higher immune infiltration characteristics and TIDE scores, whereas low-risk patients are more likely to benefit from immunotherapy. Conclusions A ccRCC survival prognostic model was produced based on the cytotoxicity-related signature, which had important clinical significance and may provide guidance for ccRCC treatment.


Introduction
Renal cancer is a frequent cancer of the urinary system accounting for about 2-3% of adult malignancies [1].Clear cell renal cell carcinoma (ccRCC) is a rather prevalent and aggressive histological subtype representing about 80-90% of metastatic renal cancer cases [2].Patients with ccRCC typically have symptoms such as fank pain, abdominal masses, and hematuria, but most ccRCC patients are asymptomatic at the early stage, so about 1/3 of patients have distant metastases at diagnosis and about 1/4 will experience relapse and metastases after surgery, making to a low overall survival and poor prognosis for ccRCC patients [3,4].Currently, a substantial proportion of patients are insensitive to chemotherapy and conventional radiotherapy and surgical intervention is still the main treatment for ccRCC patients in the early stage [5].However, about 25% of patients will recur or develop tumor metastasis after surgery.In other words, surgery cannot completely address the treatment needs of these ccRCC patients [6,7].Currently, there are studies confrming that immunotherapeutic approaches with diferent targets are efective in improving the survival rate of cancer patients [8,9].Among these, immune checkpoint blockade is a novel immunotherapy that reduces inhibitory signaling and restores tumor-specifc T-cellmediated immune responses [10].Several anti-PD-1/PD-L1 drugs have been approved for the treatment of advanced renal cell carcinoma and have shown acceptable efcacy [11].However, the efectiveness of the treatment is limited due to the complex tumor microenvironment (TME), and investigations on pathogenesis and prognostic predictors of ccRCC including new efective immunotherapeutic targets for improved clinical outcomes are still imperative.
Immune cells in TME exert a vital role in tumor progression regulation and could be used as a therapeutic target [12].Te characteristics of TME strongly infuence the response of ccRCC patients to immunotherapy because the tumor is prone to immune infltration [13].Te clinical success of antitumor immune response involved the activation and synergistic action of multiple tumor-infltrating lymphocytes [14].T cells as a type of tumor-infltrating cells are closely related to the immunosuppressive properties of ccRCC [15].Recent studies based on various biological aspects of ccRCC have revealed features associated with immune infltration.Specifc lymphocyte-related characteristics such as CD8+ Tcells [16], TNFRSF9 + CD8+ Tcells [17], and the CXCL13+ CD8+ T cells [18] have been found in ccRCC.Te CD8+ T cells are cytotoxic T lymphocytes, which are key players in eliminating the pathogen-infected cells or tumorigenic via the secreted cytotoxic proteins (perforin and granzymes) [19,20].In general, T-cell exhaustion represents an immune-functionally impaired status converted formation antitumor status of CD8+ in the TME, which is considered to be one of the main factors contributing to the low response rate to immunotherapy [21].Terefore, immunotherapy targeting the conversion of exhausted T cells recovering to an activated state in ccRCC has recently received great research attention.Te infltration of CD8+ T cells in the TME of ccRCC presents highly heterogeneous phenotypes, which are closely associated with the immunotherapy response [22].High levels of immune-evasive biomarkers and enhanced immunosuppressive infltrations are often associated with a poor prognosis in patients with ccRCC [23,24].However, interrelated association between the efectiveness of immunotherapy and the degree of CD8+ T-cell infltration in ccRCC is still not clear [25,26].Terefore, additional analysis of biomarkers linked to CD8+ T cells is urgently required in order to fnd novel prognostic markers that will guide immunotherapy for ccRCC.
Tis study characterized CD8+ T-cell-associated molecule clusters by the cytotoxic score in ccRCC; the WGCNA was used to distinguish the modules related to the cytotoxic score in the TCGA dataset.Subsequently, three distinct cytotoxicity-related molecular clusters with diferent prognoses and clinical characteristics were categorized by unsupervised cluster analysis.Finally, a risk model consisting of six cytotoxicity-related genes that afect prognosis was developed to guide prognosis and provide new insights for personalized immunotherapy.

Data Collection and
Processing.Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) is a public database that provides the ccRCC single-cell dataset GSE224630, which included six tumor samples.We used the tool of Te Cancer Genome Atlas (TCGA) GDC API to download the RNA-sequencing data and clinical grades on Kidney Renal Clear Cell Carcinoma (KIRC, https://cancergenome.nih.gov/) for the training set, which contained a total of 72 normal samples and 530 primary tumor samples (513 tumor samples had complete survival time).Te gene expression profles of RECA-EU/Renal cell carcinoma as a validation set were acquired from the database of the International Cancer Genome Consortium (ICGC), including 91 primary tumor samples.

Te Analysis of scRNA Data and Identifcation of
Cytotoxicity-Related Cluster.Te single-cell data fltering of the dataset GSE224630 was performed under the criteria that each cell expressed at least 200 genes and each gene is expressed in at least three cells.Te proportion of mitochondria and rRNA was calculated by the PercentageFea-tureSet function, ensuring that of the 100 <mRNA of each cell <5000, the mitochondrial content is <15%.Te samples were normalized by log-normalization, and the highly variable genes were identifed by the FindVariableFeatures function based on the variance-stabilization transformations (selection.method� "vst").Te FindIntegrationAnchors function (reduction � "cca") was used to remove the six samples' batch efect, and the data were consolidated using the IntegrateData method.Te FindNeighbors and FindClusters functions were applied for cell unsupervised clustering.In addition, we collected the cell markers from the CellMarker2.0website, and the cells were reclassifed according to the expression of these marker genes.tSNE dimensionality reduction of cells was performed by using the RunTSNE method.Ten, we screened these marker genes in diferent species cells using FindAllMarkers (setting by |log(fold change, FC)| � 0.5, p value < 0.05, and min.pct� 0.25).Te clusterProfler was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis.Te above-used functions are in the Seurat package [27].A signature gene set of T cells was obtained from the previous article [28], and the AUCell R package was used to calculate the cytotoxic score of each T cell [29].Moreover, the cytotoxic scores of TCGA samples were evaluated by the method of single-sample Gene Set Enrichment Analysis (ssGSEA) through the GSVA package [30], and the samples were grouped based on the median cytotoxic score.

Te WGCNA for the Cytotoxicity-Related Module Genes.
Te limma package in R (|log2(FC)| > 1 and FRD < 0.05) [31] was employed to screen diferentially expressed genes (DEGs) in two types of TCGA samples.Te gene modules the most associated with the cytotoxic score were identifed by WGCNA.After excluding the top 50% of genes with the smallest medium absolute deviation (MAD) in the gene 2 Genetics Research expression profle, the function of pickSoftTreshold in the WGCNA package was used to determine the soft threshold β [32], and the gene modules were identifed by hierarchical clustering (the minimum module size was at least 50 genes (minModuleSize � 50) and modules were merged if the distance is <0.2).Te correlation between modules and clinical features was evaluated by the Spearman method.

Defnition of the Cytotoxicity-Related Molecular Subtypes.
Te intersection between the C0 cluster of T cells, magenta module, and DEGs of ccRCC resulted in 18 overlapping genes.According to the expression of these overlapping genes, a consensus matrix of the expression profle data of TCGA-KIRC was created using the ConsensusClusterPlus R package for sample classifcation [33].Te "pam" acts as the clustering algorithm, "pearson" acts as the metric distance, and 80% of patients of TCGA-KIRC were incorporated into each bootstrap within 500 bootstraps.Cumulative distribution function (CDF) curves between the number of clustering from 2 to 10 were used to determine the optimal number of clustering and cytotoxicity-related molecular subtypes.

Immune Landscape Analysis of Cytotoxicity-Related
Molecular Subtypes.Based on the optimal number of clustering, the TCGA samples were divided into diferent subtypes.Te KM survival analysis [34] was conducted, and subsequently, the survival state of diferent subtypes was determined.Te CIBERSORT algorithm in the estimate R package [35] was used to calculate the relative abundance of 22 types of immune cells.Te MCP-counter function of the MCP-counter R package [36] and the TIMER method [37] were used to evaluate the relative proportion of 10 immune cells and the six immune cell scores, respectively.Te ES-TIMATE algorithm was used to evaluate the immune cell infltration, including the immune score, stromal score, and ESRIMATEScore, which are positively correlated with immunity, stroma, and the sum of immunity + stroma, respectively [38].To study the pathways of biological processes, a GSEA (https://software.broadinstitute.org/gsea/index.jsp)was performed using the gene set (h.all.v7.5.1.symbols.gmt,a false discovery rate (FDR) < 0.05) in the MSigDB database through the clusterProfler R package.At the same time, this gene set was used to calculate the pathway ssGSEA score of diferent risk groups using the GSVA R package.
2.6.Screening of Cytotoxicity-Related Hub Genes and Establishment of RiskScore Model.In order to establish a Risk-Score model, the genes with signifcant prognosis were selected by comparing the diferently expressed genes (C1 vs (C2, C3), C2 vs (C1, C3), C3 vs (C1, C2)) among three molecular subtypes using the limma package (setting FDR <0.05 and |log2(FC)| > 1).Te coxph function of the survival R package [39] was used in the univariate Cox regression analysis to flter key prognostic genes with p < 0.05, followed by applying the least absolute shrinkage and selection operator (Lasso) Cox regression analysis in the glmnet function of the R package to reduce the total gene number [39].Te multivariate Cox regression analysis with the stepwise regression method was performed to determine the fnal risk factor.After that, a cytotoxicity-related scoring system for ccRCC patients was established by the multivariate Cox results and the expression of genes: cytotoxicityrelated RiskScore � Σβi * Expi (i represents the expression of a risk gene, and β is the Cox regression coefcient of the gene).Te RiskScore of each patient and the optimal cutof were calculated by the survminer package [40] for patient's risk classifcation.[41].Te receiver operating characteristic analysis of 1, 3, 5 years with the AUC was applied to identify the classifcation efciency of the model using the timeROC R package [42].Te calibration curve was used to assess the predictive accuracy of nomogram model at 1, 3, and 5 years.Te decision curve analysis (DCA) was performed to evaluate the reliability of the nomogram.

Immunotherapy Evaluation of Risk Groups.
To elucidate the immunotherapy diference in diferent risk groups, the gene expression of the immune checkpoint and Tumor Immune Dysfunction and Exclusion (TIDE, https://tide.dfci.harvard.edu/)score were analyzed [1].Te TIDE score refected that the tumor inhibits the function and infltration of cytotoxic T lymphocytes to achieve immune escape [43], and a high TIME score is not benefcial for immunotherapy.
2.9.Drug Sensitivity of Risk Groups.Considering that there exists a prognostic diference in diferent risk groups, we performed further analysis of conventional drug sensitivity.Te lower 50% inhibiting concentration (IC50) value represented enhanced susceptibility to drugs.Utilizing the pRRophetic R package [44] and the pharmacogenomic data of the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/)[45], we calculated the drug sensitivity in risk groups.

Pathway and Mutation Characteristic Analysis of Risk
Groups.Te ssGSEA score of each sample in the HALL-MARK pathway and the correlation between the RiskScore and ssGSEA score were calculated.Moreover, the molecular characteristics of TCGA-KIRC in the previous pan-cancer study [46] and the mutect2-processed TCGA mutation dataset were used for mutation characteristic analysis.Te Fisher test was used to screen genes showing signifcant high-frequency mutations in diferent groups.

Statistical Analyses.
All the statistical analyses and fgures were produced in the R environment (version 3.6.3).Te Two-tailed Wilcoxon rank sum test was applied to calculate the diferences between two sets of continuous variables.Te Pearson or Spearman correlation was used for calculating the correlation matrices.Te survival diferences were depicted using KM curves with a log-rank test.Te Fisher test was used to screen the signifcant high-frequency mutation genes.Sangerbox (https://sangerbox.com/home.html), which is an interaction-friendly bioinformatics analysis platform, ofered analysis support in this paper.A p value <0.05 was considered statistically signifcant.

Single-Cell Dimension Reduction and Identifcation of
Cytotoxicity-Related Cluster.We conducted a single-cell analysis of six tumor samples from GSE224630.A total of 26906 cells were identifed (Figures S1A-S1C), and these cell samples were further divided into 11 subgroups (clusters 0-10) by dimensionality reduction cluster analysis (Figure 1(a)).According to the expression of 14 marker genes taken from the Cellmarker 2.0, the cells were reclassifed into six cell types, including the fbroblasts, T cells, epithelial cells, endothelial cells, smooth muscle cells, and B cells (Figure 1(b)).Among them, epithelial cells, smooth muscle cells, and endothelial cells accounted for the largest proportion, while the T cells had the smallest proportion in the six samples (Figure 1(c)).Te expressions of the top fve diferential marker genes with the most outstanding contributions in T cells are NKG7, CCL4, CCL5, GNLY, and KLRB1 (Figure 1(d)).Te KEGG enrichment analysis of these marker genes revealed that the pathway of regulation of actin cytoskeleton, Yersinia infection, Fc gamma R-mediated phagocytosis, endocytosis, and natural killer cell-mediated cytotoxicity were enriched in T cells (Figure 1(e)).Te T cells were further divided into 2 T cell subsets by FindClusters (resolution � 0.2), and the C0 cluster had a higher cytotoxic score than the C1 cluster (Figures 1(f ) and 1(g)).Te classifer based on the cytotoxic score exhibited favorable classifcation performance (Figure 1(h)); thus, the C0 cluster was regarded as the cytotoxicity-related cluster.

Te Cytotoxic Classifcation in TCGA Cohort and KEGG
Enrichment Analysis of T Cells.In addition, we calculated the cytotoxic score in the TCGA cohort and found that the tumor had a higher cytotoxic score via the Wilcoxon rank sum test (Figure 2(a)).Te samples were further divided into high-and low-score groups through the median cytotoxic score, and the KM survival analysis showed that the lowscore groups had better prognosis (Figure 2(b)).Te bubble plot presented the top 10 diferential marker genes in two Tcell clusters, and the antibody proteins such as FGFBP2, FCGR3A, SPON2, GZMH/B, GNLY, PLAC8, AKR1C3, PRF1, and ITGB2 were highly expressed in the C0 cluster (Figure 2(c)).Te KEGG analysis suggested that the cytotoxicity is closely related to antitumor because from the results, we observed that the C0 cluster was closely associated with the Fc gamma R-mediated phagocytosis pathway, regulation of actin cytoskeleton, Yersinia infection, and natural killer cell-mediated cytotoxicity (Figure 2(d)).Te intersection of DEGs and magenta module genes and cluster 0 resulted in 18 overlapping genes (Figure 3(c)).To further identify the subtypes, the consensus clustering analysis was performed on ccRCC samples from TCGA based on the expression profles of 18 overlapping genes.

Identifcation of
From the results of the CDF Delta area, it has a relatively stable clustering efect when the clustering is 3. Terefore, considering that the optimal clustering number (k value) of 3 is a better choice, we categorized the cohort into three (C1, C2, and C3) clusters (Figure 3(d)).Te KM survival analysis of the three subtypes demonstrated that the C1 cluster had the best prognosis, while C3 had the worst prognosis (Figure 3(e)), and the survival state in the C1 cluster was signifcantly higher than in other clusters (Figure 3(f )).

Characterization of the TME in Diferent Subtypes.
To explore the tumor microenvironment (TME) of the three molecular subtypes, we evaluated the relative abundance of 22 types of immune cells using the CIBERSORT, and the MCP-counter and the TIMER method were used to evaluate immune cell infltration.We observed that the antitumor immune cells including T-cell gamma delta, macrophage M1, and T-cell CD8 are mainly contributing to the TME score and better prognosis in the C1 cluster (Figure 4(a)), which also had a higher immune cell infltration score than patients in C2 and C3 groups (Figures 4(b) and 4(c)).Te Pathway enrichment analysis uncovered that the metabolism-related pathways in the C1 group, and the immune, cell cycle, and some tumor-related pathways in the C3 group were activated, and most of the pathways in the C2 group were suppressed (Figure 4(d)).Taken together, the C1 group was characterized by the best prognosis, higher survival rate, and higher immune cell infltration, and the C3 group was characterized by the worst prognosis.).Among these genes, SOWAHB, SLC16A12, PLG, and HHLA2 had negative coefcients in the model, indicating that upregulating of their expression levels can improve the survival time of ccRCC patients.

Validation of Model Prediction Performance.
According to the above RiskScore, high-risk and low-risk patients in TCGA-KIRC were grouped by the optimal cutof point.Te receiver operating characteristic (ROC) analysis revealed a high accuracy of the RiskScore to predict the longterm prognosis of ccRCC (AUC � 0.79, 0.73, and 0.73 at 1, 3, and 5 years, respectively) (Figure 5(d)).Patients with a high RiskScore tended to show the worst survival rate, as shown by the KM curves, with a 5-year survival of 21% (Figures 5(e) and 5(f )).Te prognosis accuracy of the RiskScore was also excellent (AUC � 0.76, 0.68, and 0.62 at 1, 3, and 5 years, respectively) (Figure 5(g)), with the high-risk group showing a signifcantly lower 5-year survival of 47% (Figures 5(h) and 5(i)) in the validation set of ICGC.A comparison of the diferences between the RiskScore and clinical grades shows that a higher RiskScore was associated with a higher clinical grade (Figure 5(j)).As a result, an efective six-gene signature to assess ccRCC prognosis was successfully created.A nomogram combining RiskScore, staging, and age was constructed, which had a wide range of RiskScore scores and contributed the most to the total score and thus had the greatest impact on survival prediction (Figure 6(c)) (p < 0.01).

Screening Independent
Te calibration curve of 1, 3, and 5 years was close to the standard curve (Figure 6(d)), indicating that the nomogram can efectively predict the actual survival outcomes.Te DCA analysis showed that the net beneft of the nomogram and RiskScore was signifcantly higher than the extreme curves, suggesting that the model had a better reliability (Figure 6(e)).Te timeROC analysis revealed a higher AUC of the nomogram and RiskScore than other clinical indicators (Figure 6(f)), which suggested that the nomogram and RiskScore had the strongest prediction ability.

Immune Characteristics between Risk Groups. Te
RiskScore was signifcantly positively (p < 0.05) correlated with the macrophage M0, T-cell CD4 memory activated, plasma cells, and T-cell regulatory (Tregs) and signifcantly negatively (p < 0.05) correlated with the dendritic and mast cell resting, monocytes, T-cell CD4 memory, NK cell resting, and macrophage M1, suggesting that the higher the Risk-Score score, the weaker the immune-killing ability of the body.Tis applies to the IL20RB and PLG as well, and downregulation of these genes can lead to a low RiskScore (Figure 7(a)).Te Analysis of the infltration of immune cells using ESTIMATE showed higher stromal scores, immune scores, and ESTIMATEScore of the high-risk group (Figure 7(b)).Te correlation between the RiskScore and immune infltration of MCP-counter also confrmed that the RiskScore was signifcantly negatively linked with the infltration of most immune cells (Figure 7(c)), suggesting that patients with a high RiskScore were less immune to tumors.Subsequently, we found that seven immune checkpoint genes showed a high expression in the high-risk group (Figure 7(d)).However, the potential clinical immunotherapy response analysis showed a higher TIDE score in high-risk patients (Figure 7(e)), suggesting a greater possibility of immune escape and relatively limited beneft from immunotherapy.In addition, the analysis of drug sensitivity showed a higher sensitivity in the low-risk group to BMS-509744, erlotinib, rapamycin, and sorafenib, while patients in the high-risk group were more sensitive to dasatinib, cisplatin, paclitaxel, and GNF-2 (Figure 7(f )).Tese results may provide guidance for the drug selection for the treatment of ccRCC.

Analysis of Pathway and Mutation Feature between Risk
Groups.Te correlation analysis between the pathway of the ssGSEA score and the RiskScore demonstrated a positive relation of the RiskScore to cell cycle-related pathways and a negative relation to metabolism-related pathways (Figure 8(a)).In addition, based on a previous pan-cancer study to characterize the mutation, the results showed that the high-risk groups were closely associated with the fraction altered, aneuploidy score, number of segments, and homologous recombination defects (Figure 8(b)).Te Fisher analysis of gene mutation characteristics in the mutect2treated dataset revealed that the PBRM1 is a highly mutated gene that afected tumorigenesis in the two groups, the BAP1 was a high-frequency mutation gene in the high-risk group (Figures 8(c) and 8(d)), and their function should be further studied.

Discussion
Renal cancer remains a serious urinary system problem with high mortality and morbidity [47].Te tumor microenvironment of ccRCC is usually accompanied by a high level of CD8+ T-cell infltration [25], which is closely associated with patients' prognosis and immunotherapy efcacy in ccRCC.To enhance the understanding of biological functions of CD8+ T cell in TME, this study developed a cytotoxicity-related signature to provide therapeutic guidance and prognostic prediction for ccRCC patients.In this study, we identifed T-cell subtypes based on single-cell sequencing data and used WGCNA to identify the gene modules most associated with cytotoxicity.In particular, in the TCGA dataset, the tumor tissues had a high-cytotoxic score and inferior prognosis, suggesting that the cytotoxic characteristics of TME are closely correlated with tumorigenesis.Omics analysis has facilitated the discovery of disease markers [48].In addition paper, using transcriptomics, we constructed a cytotoxicity-related prognostic prediction risk model consisting of six key genes (SOWAHB, SLC16A12, IL20RB, SLC12A8, PLG, and HHLA2), precise and concise, showing a preferable application prospect that benefts the prognostic evaluation personalized treatment for ccRCC patients.Te SOWAHB, SLC16A12, PLG, and HHLA2 are considered protective factors, and the IL20RB, SLC12A8 are risk factors for ccRCC patients.Te SOWAHB was reported as a potential regulator in ccRCC progression, but its functionality needs further verifcation.SLC16A12 is an identifed creatine transporter [49], the lack of SLC16A12 usually leads to the low levels of creatine in chronic renal failure [50], and the comparatively lower expression of SLC16A12 in tumors often correlates with the worst prognosis [51].Tis implies that SLC16A12 may exert its antitumor efects in ccRCC by maintaining the stability of energy metabolism in renal cells.In addition, the PLG is an antitumorigenic factor, and its hydrolysate contains angiotensin that will function against cancer progression [52].
Te expression of PLG in ccRCC was lower than that in the adjacent normal tissue, indicating low levels of PLG in favor of ccRCC progression, leading to the worst overall survival (OS) [53].Te HHLA2 protein is expressed on antigenpresenting cells, and elevated HHLA2 is thought to be associated with a more severe pathology and poor prognosis in cancer patients [54][55][56][57], whereas several studies showed that patients with higher HHLA2 expression had better survival rate [58,59]; this paradox can be explained, in part, by the dual role of HHLA2 in immunity and by the fact that the HHLA2 acts as a protective factor in this article.Te high expressions of IL20RB and SLC12A8 usually predict the worst survival and were unfavorable prognostic biomarkers for ccRCC [60,61].Tis implies that all of these key genes have important roles in promoting the development of ccRCC, and their associated mechanisms of action still need to be explored in depth in future studies.
To evaluate the diference in the TME in the high-and low-risk groups, we performed the immune-infltration landscape analyses.Te high-risk group exhibited elevated immune infltrations, including the stroma and immune score, which are associated with the worst prognosis [62,63].Meanwhile, the RiskScore was positively correlated with the immunosuppressive cells (M0 dormant macrophage and Tcell regulatory) and negatively correlated with the immuneactivated cells (macrophage M1 and neutrophils), which explained an unfavorable ccRCC prognosis in patients showing high immune cell infltration [25,64].Te immunoediting theory suggests that the lack of immune cells, the presence of immunosuppressive cells, and high levels of immunosuppressive cytokines and fbrosis in TME might contribute to the immune escape of tumors [65][66][67].Te high-risk group had more infltrations of immune cell and stroma cell as well as the high expression of immune checkpoint genes and a high TIME score, indicating the complex immune escape mechanism in the high-risk group.In the HALLMARK pathway analysis, the RiskScore was positively correlated with cell cycle-related pathways that can help to fnd more cell cycle check points.In addition, the high-risk group exhibited more variety of chromosome mutations and high mutant frequency of tumor suppressor BAP1 [68], while knockdown of BAP1 inhibited tumorigenicity and lung metastasis [69]; thus, the role of BAP mutations in the development of renal cancer needs to be further elucidated.PBMR1 is a hub gene that had higher mutation frequency in high-and low-risk groups.Usually, unbalanced PBMR1 leads to a shift from immune activation to immune suppression and is deregulated in tumors [70]; this could be a target gene for ccRCC patients.
Finally, there are some problems and limitations in this paper.First, our study was based only on samples from public databases due to the limited sample size; for this reason, future studies will incorporate more samples, including patients of diferent ethnicities and geographic regions, in order to improve the model's generalization ability and applicability.In addition, the specifc functions and substrates of our screened key genes in ccRCC have not yet been thoroughly validated.Terefore, further in vivo and in vitro experiments, including cellular and animal models, are necessary to investigate the mechanism of action of the key genes.Finally, inferring partial immune characteristics from gene expression data alone may fail to fully resolve the full complexity and dynamics of the immune microenvironment.In the future, we will combine single-cell multiomics technology and spatial transcriptomics to deeply analyze the complex changes in the immune microenvironment of ccRCC.

Conclusion
Cytotoxicity-related prognosis genes were selected by performing univariate/multivariate and Lasso-Cox regression analysis and further used to build a risk model for accurately predicting the immunotherapy response and clinical outcomes of patients sufering from ccRCC.High-and low-risk ccRCC patients with diferent clinical features and immunogenomic landscapes were grouped by the gene classifer (RiskScore model), which can provide therapeutic guidance for ccRCC patients and improve the current individuated treatment options.

Figure 1 :Figure 2 Figure 2 :Figure 3 Figure 3 :Figure 4 Figure 4 :
Figure 1: Single-cell analysis and cytotoxicity-related T-cell clustering.(a) Te tSNE plot of 11 subgroups' distribution after clustering.(b) tSNE diagram of six species cells' distribution after annotation.(c) Te proportion of six species cells in diferent samples.(d) Te heatmap of top fve marker genes of six species cells.(e) Te KEGG enrichment results of the marker gene in six subgroups.(f ) Te diference of cytotoxic scores in T-cell subsets.(g) Te tSNE plot of two T-cell subsets.(h) Te AUC values of T cells with diferent cytotoxic scores.

Figure 5 :
Figure 5: Establishment and validation of cytotoxicity-related risk model.(a) Te trajectory of each independent variable as lambda changes.(b) Te lambda confdence interval analysis.(c) Multifactor forest map of model genes.(d and e) Te ROC curve and Kaplan-Meier curve of risk model in TCGA datasets.(f ) Survival state diferences between risk groups in TCGA datasets.(g and h) Te ROC curve and Kaplan-Meier curve of risk model in ICGC datasets.(i) Survival state diferences between risk groups in ICGC datasets.(j) RiskScore diference among diferent clinical grades in TCGA datasets.

Figure 6 :Figure 7 :Figure 7 :Figure 8 :
Figure 6: Identifcation of the independent risk factor and a nomogram development.(a) Te forest plot of univariate Cox regression analysis between the RiskScore and clinical features.(b) Te forest plot of the multivariate Cox analysis between the RiskScore and clinical features.(c) Nomogram model integrating the RiskScore, age, and stage.(d) Calibration plots for 1, 3, and 5 years of the nomogram.(e) Te decision curve of a nomogram.(f ) Te timeROC curves with AUC of a variety of clinical features for overall survival (OS) at 1-5 years.

Figure 8 :
Figure 8: Analysis of pathway diference and mutations between risk groups.(a) Correlation between the RiskScore and HALLMARK pathway of the ssGSEA score.(b) Molecular characterization scores between diferent risk groups.(c) Diferentially mutated genes between diferent risk groups.(d) Gene waterfall map of diferential mutations between diferent risk groups.* * * * p < 0.0001.