Pyroptosis-Related Gene Signature Predicts the Prognosis of ccRCC Using TCGA and Single-Cell RNA Seq Database

Clear cell renal cell carcinoma (ccRCC) is the most prevalent type of renal carcinoma, which is not sensitive to both radiotherapy and chemotherapy. The objective response rate of metastatic renal cancer to targeted drugs and immunotherapy is unsatisfactory. Pyroptosis, proven as an inflammatory form of programmed cell death, could be activated by some inflammasomes, while could create a tumor-suppressing environment by releasing inflammatory factors in the tumor. To explore indicators predicting the prognosis of ccRCC and the effect of antitumor therapy, we constructed a pyroptosis risk model containing 4 genes after 11 pyroptosis-related genes of 516 ccRCC cases in the TCGA database were scanned. Based on the risk score, 516 ccRCC cases were divided into two groups for functional enrichment analysis and immune profile to seek functional pathways and potential therapeutic targets. Besides, those results were verified in GSE29609 and single-cell transcriptomic data. The study suggests that the conducted pyroptosis model could predict the prognosis of ccRCC and reflect the immune microenvironment, which may help in immune checkpoint inhibitor treatment.


Introduction
Renal cell carcinoma (RCC), which originates from renal proximal convoluted tubule epithelial cells, accounts for about 90% of all primary renal malignancies. RCC is one of the most prevalent malignant tumors of the genitourinary system. In recent decades, the incidence of RCC continues to rise [1]. According to histological classification, clear cell renal cell carcinoma (ccRCC) is the most common type of renal carcinoma (about 75-80%), followed by papillary carcinoma (15%) and chromophobe cell carcinoma (5%) [2]. e onset of RCC is insidious and lacks specific clinical manifestations and features at the early stage. About 20-30% of RCC has been found with metastasis at the time of initial diagnosis [3]. RCC is not sensitive to both radiotherapy and chemotherapy, which currently mainly relied on surgical resection. Despite the development of targeted drugs and immunotherapy in recent years, the objective response rate of metastatic renal cancer is only about 30% [4]. Although some clinical indicators and pathological results have been used to predict the treatment and prognosis of ccRCC, their pyroptosis could create a tumor-suppressing environment with released inflammatory factors in different tumor types. However, pyroptosis can also debilitate the own immune response to cancer cells and quicken tumor growth [12][13][14][15].
ccRCC with metastasis is usually incurable by surgical resection and requires systemic treatment [16]. However, metastatic RCC shows insensibility to radiotherapy and systematic treatment in the later stages of treatment, including hormone therapy, chemotherapy, and interleukin-2-based immunotherapy [17].
e study of the Cancer-Genome Atlas has significantly advanced the molecular classification of renal cell carcinoma to guide the treatment and prognosis. Among them, the activation of protein kinase B (PKB/Akt), the mammalian target of the rapamycin (mTOR) pathway is a key driver of RCC. e expression and activity of mTOR downstream effectors in RCC are unbalanced, which lays a theoretical foundation for the clinical application of ccRCC-targeted therapy [18]. In the last few years, with the indepth study of cytotoxic T cell inhibitory molecules such as cytotoxic T lymphocyte-associated protein 4, programmed cell death receptor 1 (PD-1), and programmed cell death ligand-1 (PD-L1), the immunomodulators have been applied in clinical practice. Although these treatments improved the prognosis of ccRCC, drug resistance and recurrence still occurred [19,20]. Due to the lack of a well-established subgroup classification of ccRCC, there is still a lack of molecular subtypes to guide clinical practice. Consequently, it is urgent to construct an effective genetic signature to guide subgroup classification and predict prognosis.
In view of the important part of pyroptosis in the development and treatment of ccRCC, we constructed a pyroptosis risk model to classify ccRCC in the present study to predict prognosis and treatment. 11 pyroptosis-related genes were selected, and survival analysis and GSVA analysis were performed on 516 ccRCC cases in TCGA database. e pyroptosis risk model containing 4 genes was constructed by the multivariate COX regression analysis. Using the risk score, 516 ccRCC cases were grouped into two groups for functional enrichment analysis and immune profile. ose results were verified in GSE29609 and single-cell transcriptomic data. Our findings suggest that the pyroptosis model could predict the prognosis of ccRCC and reflect the immune microenvironment. Figure 1 shows the flowchart of the study.

Database.
e patients' characteristics and renal cell carcinoma patients with level 3 gene expression profiles were downloaded from the TCGA database (June 2020) (https:// cancergenome.nih.gov) [21]. We selected the 516 cancer cases whose pathological diagnosis is clear cell renal cell carcinoma (ccRCC). Cases without pathological or clinical information were excluded. GSE29609 dataset [22], as a validation cohort, was obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo). e limma package in R was used to normalize the gene expression.
Single-cell transcriptome profiling data for analyses were downloaded from the supplemental data in the published article [23]. e Seurat package in R (version 4.0.4) was applied to process the single-cell RNA-seq data. Cell clusters were recognized by Uniform Manifold Approximation and Projection (UMAP) with a resolution of 0.5 [24]. e function FeaturePlot and VlnPlot of the Seurat package were used for visualization of the expression profiling of the genes.

Construction of a Risk Model.
e heatmaps of pyroptosis-related genes were generated by the pheatmap package in R. Pyroptosis pathway enrichment was performed based on the pyroptosis-related signatures and Gene Set Variation Analysis (GSVA) [33]. R package survival was used for the Kaplan-Meier survival analysis of these 11 pyroptosis-related genes. e function coxPH of the survival package was utilized to generate the Cox proportional hazards regression model. We selected the prognosisrelated genes and formed the formula: risk score � β1gene1 × expression of gene1 + β2gene2 × expression of gene2 +. . .+ βngenen × expression of genen. e R package survival and ROCR were applied to form the Kaplan-Meier analysis and the receiver operating characteristic (ROC) curves. e predictive value of the new risk model was validated using the GSE29609 dataset downloaded from the GEO database.

Identification of Differentially Expressed Genes and
Functional Enrichment Analysis. Cancer cases were grouped into two groups, the high-risk group, and the low-risk group, following the median value of the risk score calculated. e differentially expressed genes (DEGs) between the two groups were identified using the limma package in R with the fold change (|fold change| ≥ 1.5) and adj. P < 0.05. e functional enrichment analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis were conducted using the clusterProfiler package in R.

Assessment of Immune Cell Type Fractions.
e analytical web server tool CIBERSORT (https://cibersort.stanford.edu/) was applied to estimate the immunologic cell abundances in the cancer immune microenvironment [34]. e leukocyte gene signature matrix termed LM22 was used to distinguish the 22 immune cell types between the high-and low-risk score groups. e 22 immune cell types including CD8 T cells, naive CD4 T cells, resting memory CD4 T cells, activated memory CD4 T cells, naive B cells, memory B cells, plasma cells, follicular helper T cells, T cells regulatory (Tregs), gamma delta T cells, resting NK cells, activated NK cells, monocytes, macrophages M0, macrophages M1, macrophages M2, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils.

Assessment of Immunomodulators and Immunosuppressive Cytokines' Expression Profile.
We quantified a group of key immunomodulators and tumor immunosuppressive cytokines.

Overview of Pyroptosis-Related Genes in ccRCC.
e selected 516 cancer cases from the TCGA database were pathologically diagnosed as ccRCC. e basic patient information and characteristics were shown in Table 1. Based on the Kaplan-Meier survival curves, the 11 pyroptosis-related genes were all significantly related to the overall survival (OS) outcome of the cancer cases with the log-rank test P < 0.05 ( Figure 2). e distinct gene expression patterns of the 11pyroptosis-related genes in these cancer cases were presented in the heatmap (Figure 3(a)). Pyroptosis activity was calculated based on the pyroptosis signatures and GSVA. As shown in Figure 3(b), cancer cases with T staging III/IV, according to tumor node metastasis (TNM) classification, had higher GSVA scores than cancer cases with T staging I/II (P < 0.001). Figure 3(c) exhibited that no significant difference was observed in GSVA scores when the cancer cases were grouped according to with or without regional lymph node metastasis. When compared to cancer cases without metastasis (M0), cancer cases with distant metastasis (M1) had significantly higher GSVA  ). e cancer cases were divided into low-and high-GSVA score groups, based on the optimal cut-off value calculated by the survminer package in R. As shown in Figure 3(e), the cancer cases in the high-GSVA score group had a poorer prognosis, while cancer cased in low-GSVA score group had better OS (log-rank test P < 0.001).

Construction and Evaluation of the Pyroptosis Risk Model.
Multivariate Cox regression analyses were used for the pyroptosis risk model establishment (Supplementary  Table S1). Using the genes with a P value less than 0.1 within the supplementary table, the pyroptosis risk model was Following the median value of risk score (P < 0.05), all 516 ccRCC cases were divided into low-and high-risk groups. e whole of the 11 pyroptosis-related genes was upregulated in the high-risk score group (Figure 4(a)). ROC curves and Kaplan-Meier analysis were applied to assess the pyroptosis risk model. e risk model had an accuracy of 0.674 (95% CI: 0.624-0.724) in the TCGA cohort (Figure 4(b)). e cancer cases in high-risk score group had significantly poor OS (P < 0.001) (Figure 4(c)). To reveal the independent predictability of the risk model   Table 2. e hazard ratio (HR) was 2.444 (95% CI 1.863-3.205) (P < 0.001).
Using the dataset GSE29609, external validation was performed. As shown in Figure 4(d), the area under the ROC curve (AUC) was 0.679 (95% CI: 0.506-0.852). Figure 4(e) presented the same results that the cancer cases in high-risk score group had a significantly poor OS (P � 0.045).

Functional Enrichment Analyses.
DEGs between lowand high-risk score groups were identified. e profiles of DEGs expression for each group were exhibited in the heatmap (Figure 5(a)). GO and KEGG analyses were taken to appraise the biological involvement of the DEGs. As highlighted in Figure 5(b), the top GO terms comprised acute-phase response, carboxylic acid transport, humoral immune response, etc. Furthermore, KEGG analysis exposed that the DEGs were chiefly involved in carbohydrate digestion and absorption, complement and coagulation cascades, glycolysis/gluconeogenesis, PPAR signaling pathway, etc. (Figure 5(c)).

Immune Microenvironment of Low-and High-Risk Score
Groups. In consideration of the established pyroptosis risk model that could also reflect the immune microenvironment of ccRCC, the disparate immune cell fraction between low-and high-risk score groups was studied. e diverse immune cell fraction upshot of the 516 ccRCC cancer cases grouped into different risk score groups was depicted in Figure 6(a). In spite of the higher multiple effector immune cells (e.g. plasma cells, CD8+ T cells) in high-risk score groups, the immunosuppressive cells (e.g. regulatory T cells) were significantly higher in the same group (Figure 6(b)) (P < 0.05). is status may imply the immunosuppressive microenvironment in high-risk-score cancer cases. In addition, we ferreted out the dissimilar expression of immunomodulators and immunosuppressive cytokines between low-and high-risk score groups. e expression of ICOS, KLRC1, PDCD1, TIGIT, ICAM1, IFNB1, CTLA4, and LAG3 were all significantly upregulated in the high-risk score group (Figure 6(c)) (P < 0.05). We discovered that several immunomodulators (e.g. TGFβ1, IL10) were upregulated in the high-risk score group, while NOS2 and NOS3 TMEM27  SLC16A9  HAO2  GLYATL1  CXCL5  PITX2  SLITRK2  SLN  SLC38A5  CPA4  IGDCC4  IGF2BP2  GPRC5A  CLMP  HMGA2  NIPAL4  MMP12  PPP2R2C  PANX2  PI3  RARRES1  SLPI  MFSD2A  MAT1A  PPP1R1A  MOCOS  ITPKA  MELTF  GFPT2  PKP3  TNNT1  PTPRH  IGF2BP3  PCBP3  ADGRF4  PADI3  KCNS1  WFDC5  UGT1A10  PITX1  FGG  FGB  FGA  PRAME  AQP9  SAA2.SAA4  SAA2  SAA1  HP  F2  LBP  CIDEC  IL20RB  IGFBP1  PAEP  PLA2G2A  IL6      usly, the cancer cases with a higher pyroptosis risk score may exist in an immunosuppressive microenvironment.

Single-Cell Transcriptomic Analysis of the Four Modeling
Genes. In order to delve into the extra interrelation among the four modeling genes in ccRCC, single-cell transcriptomic data were exploited for further analysis. We identified 18 different cell clusters, including cancer cells, renal tubule cells 1, renal tubule cells 2, renal tubule cells 3, CD8+ T cells, CD4+ T cells, regulatory T cells (Treg cells), natural killer cells (NK cells), macrophages/dendritic cells (MACDC) 1, MACDC 2, B cells, neutrophils, fibroblasts (FIB), endothelial cells (EC) 1, EC 2, EC 3, collecting duct cells 1, and collecting duct cells 2 (Figure 7(a)). e different expression profiles of the four modeling genes in different types of cells were scrutinized. FeaturePlot revealed that the expression of the four modeling genes was higher in cancer samples than that in adjacent non-neoplastic samples (Figures 7(b)-7(e)).
Due to the important role of pyroptosis-related genes in the occurrence, development, and prognosis of ccRCC, we constructed a pyroptosis risk model by multivariate COX regression analysis. e pyroptosis risk model, containing 4 genes namely CASP3, CASP4, GSDMB, and GZM, had the independent predictability of ccRCC prognosis in the TCGA cohort as well as the dataset GSE29609. In the functional enrichment analyses, the glycolysis signaling pathway was exhibited. It had been reported that glycolysis could play a key role in the process of proinflammatory activation during cell pyroptosis, in which interleukin (IL)-1β and IL-18 are released from plasma membranes. e metabolism of macrophages could switch from oxidative phosphorylation to glycolysis following proinflammatory activation [45]. Besides, the insulin resistance pathway was enriched in our study. Pyroptosis occurs not only in monocytes and dendritic cells but also in nonmacrophage cells [46][47][48]. It had been reported that adipose tissue also experiences pyroptosis [49,50]. e intracellular concentration of LPS-inducing pyroptosis determined the adipocyte death size [51]. e adipocyte overexpansion induces a stress state, leading to obese adipocyte pyroptosis, which in turn recruits macrophages into adipose tissue and induces inflammation and insulin resistance in obese mice [50].
In addition, we explored the immune microenvironment of low pyroptosis and high pyroptosis risk score groups. In high pyroptosis risk score groups, the multiple effector immune cells such as plasma cells and CD8+ T cells were higher, possibly due to stimulation by inflammatory factors released during pyroptosis. However, the immunosuppressive cells such as regulatory T cells were also higher in the same groups, implying the immunosuppressive microenvironment induced by pyroptosis. In addition, we ferreted out the differentially expressed immune checkpoints. e expression of ICOS, KLRC1, PDCD1, TIGIT, ICAM1, IFNB1, CTLA4, and LAG3 were all significantly upregulated in the pyroptosis high-risk score group, indicating that pyroptosis could stimulate the activation of immune systems. In terms of immunomodulators, TGFβ1 and IL10 were upregulated, while NOS2 and NOS3 were reduced in the pyroptosis high-risk score group. e discovered phenomenon of tumor immune microenvironment characterized the multiple impacts of pyroptosis. e pyroptosis risk score model may be useful in immunotherapy response grading and classification. e ccRCC with higher pyroptosis risk may benefit from the renovation of the immunosuppressive condition of the tumor microenvironment. ese may provide ideas for preventing drug resistance and increasing the efficiency of immunotherapy.
However, there are still some defects in the study. We cannot avoid the potential for selection bias, since we drew the results based on the data downloaded from the TCGA database and GEO database. We cannot obtain additional detailed clinical information for further analysis. In this study, we exploited single-cell transcriptomic data to identify 18 different cell clusters, validating the extra interrelation among the four modeling genes in ccRCC. Further clinical trials and single-cell transcriptomic-based analysis should be performed to validate the pyroptosis risk model. Despite these defects listed above and the lack of further validation, the presented findings still proved the predicting ability of the pyroptosis risk model statistically and its potential application.

Conclusion
e study developed a pyroptosis-related risk model based on 4 identified pyroptosis-related genes. e conducted pyroptosis model could predict the prognosis of ccRCC and reflect the immune microenvironment, which may help in prognostic biomarker discovery in ccRCC patients and immune checkpoint inhibitor treatment in the future.
Data Availability e datasets used and/or analyzed during the current study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. .

Authors' Contributions
Ying Gan and Zhenan Zhang obtained the data and performed the analysis. Xiaofei Wang and Aolin Li checked the statistical method. Ying Gan, Zhenan Zhang, Xiaofei Wang, and Aolin Li performed the figures. Ying Gan and Zhenan Zhang wrote the manuscript. Qian Zhang and Yu Fan designed and supervised the study. All authors contributed to the article and approved the submitted version. Ying Gan and Zhenan Zhang contributed equally to this work.