Leveraging diverse cell-death patterns to predict the prognosis, immunotherapy and drug sensitivity of clear cell renal cell carcinoma

Clear cell renal cell carcinoma (ccRCC) poses clinical challenges due to its varied prognosis, tumor microenvironment attributes, and responses to immunotherapy. We established a novel Programmed Cell Death-related Signature (PRS) for ccRCC assessment, derived through the Least Absolute Shrinkage and Selection Operator (LASSO) regression method. We validated PRS using the E-MTAB-1980 dataset and created PCD-related clusters via non-negative matrix factorization (NMF). Our investigation included an in-depth analysis of immune infiltration scores using various algorithms. Additionally, we integrated data from the Cancer Immunome Atlas (TCIA) for ccRCC immunotherapy insights and leveraged the Genomics of Drug Sensitivity in Cancer (GDSC) database to assess drug sensitivity models. We complemented our findings with single-cell sequencing data and employed the Clinical Proteomic Tumor Analysis Consortium (CPTAC) and qRT-PCR to compare gene expression profiles between cancerous and paracancerous tissues. PRS serves as a valuable tool for prognostication, immune characterization, tumor mutation burden estimation, immunotherapy response prediction, and drug sensitivity assessment in ccRCC. We identify five genes with significant roles in cancer promotion and three genes with cancer-suppressive properties, further validated by qRT-PCR and CPTAC analyses, showcasing gene expression differences in ccRCC tissues. Our study introduces an innovative PCD model that amalgamates diverse cell death patterns to provide accurate predictions for clinical outcomes, mutational profiles, and immune characteristics in ccRCC. Our findings hold promise for advancing personalized treatment strategies in ccRCC patients.


Introduction
Renal cell carcinoma (RCC) is the commonest kidney cancer, with its incidence increasing continuously [1; 2].According to traditional morphological classi cation, there are three different subtypes of RCC: clear cell, papillary, and chromophobe subtypes [3].ccRCC, which constitutes for 80-90% of RCC, has the least favorable prognosis and highest mortality rate compared to the other two subtypes [4].Despite the progress we made on managing ccRCC, it is still mainly treated with surgery, and no adjuvant therapy is found effective [4].So, the Survival outcomes of ccRCC are unsatisfactory.Among a study of 4034 patients with ccRCC in ve centers in Germany, 5-year cancer-speci c survival rates were 75.8% [5].Given the limited treatment strategies, novel targets need to be explored to improve prognosis.
Currently, a different mechanism associated with tumor formation, PCD, attracts attention gradually.PCD refers to a type of cell death with elaborate regulations and various mechanisms.Many forms of PCD have been recognized in recent decades, including apoptosis, pyroptosis, necroptosis, ferroptosis, entotic cell death, parthanatos, netotic cell death, autophagy-dependent cell death, lysosome-dependent cell death, alkaliptosis and oxeiptosis [6].Apoptosis is a manner of removing speci c cells of the body with no in ammatory responses.Speci c changes occur in apoptotic cells in order, including solidifcation, nuclear cleavage, and nucleolysis [7; 8].Originally, necroptosis is concerned as alternative form of apoptosis.Its characteristic features are necrotic cell death morphology and autophagy activation [9; 10].
As an in ammatory type of PCD, pyroptosis is triggered by speci c in ammasomes and activates cytokines like IL-18 and IL-1β [11].Ferroptosis is characterized by the accumulation of iron-dependent lipid hydroperoxides to lethal levels [12].Entotic cell death is a form of PCD that is performed non-cellautonomously by lysosomes and autophagy proteins [13].Reticular cell death is caused by the release of neutrophil extracellular traps (NETs) and depends on reactive oxygen species (ROS) produced by NADPH oxidase, making it distinct from other types of PCD [14].The execution of Parthanatos does not depend on the mediation of caspases, but on the overactivation of ribozyme, PARP-1 [15; 16].lysosomedependent cell death is hallmarked by lysosomal destabilization and dependence on lysosomal membrane permeabilization [17].Autophagy is the process of transporting endogenous or exogenous cytoplasmic substances to lysosomes for recycling the basic components of cells through degradation .Autophagy-dependent cell death is considered to be a form of PCD that is mechanistically dependent on the autophagic machinery, independent of apoptosis or necrosis.[18].Alkaliptosis is a pH-dependent form of PCD that is regulated by intracellular alkalinisation of the cell.[19].Oxeiptosis is a caspaseindependent and non-in ammatory form of PCD, induced by reactive oxygen species (ROS).Moreover, evidence show that it can induce other forms of PCD, including pyroptosis, apoptosis,ferroptosis, necroptosis and Autophagy-dependent cell death [20; 21].It has been known for decades that PCD plays a fundamental role in metazoa.However, defects in PCD is associated with development and metastasis of cancer as well as resistance to anticancer therapy [22].Therefore, more researches are needed to utilize the mechanism of PCD to ght against cancer in the future.
Despite the progress we made in the molecular mechanism of PCD, the clinical impacts of PCD in ccRCC largely remain unclear.We aimed to identify the molecular alterations and clinical relevance of PCDrelated genes in ccRCC in this study.We have developed a novel indicator, the PRS, to predict the e cacy of therapeutic interventions and prognosis in ccRCC.Using the PRS may help to select more appropriate therapeutic regimens for ccRCC patients in the future.

Data acquisition and processing
The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases were used to obtain transcriptomic and clinical data on ccRCC.From previous literature, 1078 genes related to PCD were counted.The E-MTAB-1980 was utilized as an external validation data set.From the GEO database, we downloaded three single cell sequencing datasets: GSE131685, GSE152938, and GSE171306.

Single cell sequencing data processing
There were 4 ccRCC samples and 4 normal samples with a total of 64926 cells.The R package Seurat was used to analyse the single cell data.The speci c procedure was rstly, ltering of low quality cells by the criteria that the number of expressed genes was greater than 100 and less than 6000, and the percentage of mitochondrial gene expression was less than 20.In addition, we ltered out genes with a low level of expression by the criterion that the gene must be expressed in at least 100 cells.We ended up with 17304 genes and 50201 cells.We then integrated the data to remove the batch effect using FindIntegrationAnchors (with the reduction parameter set to "rpca") and the IntegrateData function.
Next,to reduce the data dimensionality, the functions RunPCA and RunUMAP were used.The top 30 of the principal components and the top 2000 of the highly variable genes were used throughout.Finally, we clustered and grouped the cells using the functions FindNeighbors and FindClusters with a resolution parameter of 1.5.We de ned cell sub-populations by means of classical marker genes from the literature [19,20].In addition, differentially expressed genes (DEGs) between tumour cells and normal cells were screened according to the criteria of adj.p-value < 0.001 and |logFC| > 1.
Construction of the PCD-related signature First, the intersection genes of DEGs between single-cell data and TCGA data sets were screened.Univariate Cox regression analysis and LASSO analysis were used to screen for PRGs that were strongly correlated with prognosis.The R package "glmnet" was used for the construction of PRS.By linearly combining the gene expression-weighted regression coe cients, we obtained the PRS formula.The algorithm was as follows: PRS = Coef A * Gene A expression + Coef B * Gene B expression + Coef C * Gene C expression+......Coef N * Gene N expression, with Coef referring to the coe cient calculated by LASSO and gene expression referring to the expression of PRGs.The ccRCC patients were classi ed into high and low PRS groups on the basis of median PRS.In addition, patients were randomised into trainRisk and testRisk groups according to 6:4 criteria, and the E-MTAB-1980 dataset was used as an external validation set to verify the stability and accuracy of the model.To evaluate the prognostic characteristics of PRS, the time-dependent receiver operating characteristic (ROC) curve, Kaplan-Meier (KM) curve, univariate and multivariate Cox regression analyses were performed.

Construction of the PCD-related clusters and Bioinformatics analysis
On the basis of the expression pro les of the modeled PRGs, PCD-related clusters were constructed by means of consensus clustering.The KM method was used for the estimation of the difference in overall survival (OS) between clusters.DEGs were screened according to the criteria of |log 2 (fold change FC)| 2 and adj.P value 0.001 to further analyse the differences in biological pathways between clusters.Then, the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) were carried out on the DEGs.In addition, the R package GSVA and "c2.cp.kegg.v7.4.symbols" from MSigDB were used for Gene Set Variation Analysis (GSVA) enrichment analysis to evaluate the pathway enrichment for different clusters.

Immune microenvironment assessment and Mutation analysis
Multiple algorithms were used to assess immunoin ltration in each ccRCC sample.ESTIMATE was used to measure the level of immune cell in ltration (immune score) and stromal content (stromal score) for each ccRCC sample.The anticancer immune response (cancer immune cycle) in the tumour microenvironment (TME) consisted of seven steps.Tracking Tumor Immunophenotype (TIP, http://biocc.hrbmu.edu.cn/TIP/)immune activity scores for ccRCC samples were collected.Immunotherapy data for ccRCC were downloaded from The Cancer Immunome Atlas (TCIA).Differences in immunotherapy (anti-PD-1 and anti-CTLA4) between groups were further analysed.Somatic variant data were stored in Mutation Annotation Format (MAF), and the Maftools package was used to illustrate the respective mutational pro ling of the two risk levels by waterfall plot.

Statistical analysis
All analyses were performed by using R 4.2.2.All statistical tests were two-sided, and P-value < 0.05 was considered statistically signi cant unless otherwise noted.KM curve and log-rank test were utilized to analyze the correlation between PRGs and overall survival (OS) in ccRCC patients.All the adj p-value or false rate discovery (FDR) were adjusted by Benjamini-Hochberg (BH).

Results
Construction and validation of the PCD-related signature 24 differentially expressed PRGs in single cell sequencing data and TCGA were screened according to the standard | logFC | > 1 + FDR < 0.001 and the intersection was taken (Fig. 1A).Subsequently, 9 prognostic PRGs were obtained by means of univariate COX regression analysis (Fig. 1B).To construct PRS for predicting patient prognosis, 8 PRGs were selected by LASSO regression analysis (Fig. 1C and 1D).PRS was positively associated with mortality, as shown by the distribution of PRS and survival status and the KM survival curve (Fig. 1E and 1F).The heat map demonstrated the distribution of the 8 modeled gene expression pro les and clinicopathological features.(Fig. 1G).ROC curves about PRS in 1,2,3 years were 0.724, 0.663 and 0.651 (Fig. 1H).In addition, Univariate and multivariate cox analysis showed that PRS had good accuracy in predicting the prognosis of ccRCC patients (Fig. 1I and 1J).The E-MATB-1980 data set was employed as external validation set to further verify the accuracy and stability of the model.The KM curve indicated that patients in the high PRS group had a worse prognosis (Fig. 1K).ROC curves about PRS in 1,2,3 years in the E-MATB-1980 dataset were 0.851, 0.885, and 0.814, respectively (Fig. 1L).

Internal validation of the PCD-related signature
To further validate this, we randomised patients to the training and test groups according to a 6:4 ratio.The distribution of PRS and survival status and KM curves revealed that the higher PRS was related to poorer prognoses in both the training and test groups.(FigureS1A-S1D).As shown in gure S1E and S1F, the distribution of the 8 modeled gene expression pro les and clinicopathological features were summarised.The ROC curve represented optimal prognostic value in both groups (training set AUC = 0.751; test set AUC = 0.680) (Figure S1G and S1H).In both univariate and multivariate Cox regression analyses, PRS remained an independent predictive factor in both training and testing groups.(FigureS1I-S1L).The veri cation process above indicated that PRS was a relatively stable predictor for prognosis.

Identi cation of immune characteristics of the PCD-related signature
Seven algorithms were used to produce a heat map of the different immune cell components.This was shown in Fig. 2A.The correlation between PRS and immune cells according to different algorithms was shown in Fig. 2B.We further found that the expression of the three immunosuppressive cells in the high PRS group was signi cantly higher (Fig. 2C-2E).Using the ESTIMATE algorithm, immune score, stromal score and ESTIMATE score were all higher in high PRS group (Fig. 2F-2H).In order to further con rm the reliability of the PRS for immunotyping, we analysed the association between the risk scores and the immune subtypes.The expression of the PRS was higher in C1 and C6, but lower in C3, C4, and C5 (Fig. 2I).In addition, the high PRS group had markedly higher levels of all immune-related molecules than the low PRS group (Fig. 2J).The anti-tumour immune response need to be progressively initiated to effectively kill tumour cells.Therefore, in order to further investigate the role of immune cells in the progression of ccRCC, we obtained the immune activity scores at each step in the ccRCC samples from the TIP database.As shown in Fig. 2K, the frequency of anti-tumour immune cells in the high PRS group was signi cantly higher than in the low PRS group.
Drug Susceptibility Analysis of the PCD-related signature We used the "pRRophetic" package to evaluate the differences between high and low PRS group in the IC50 of ve common anti-renal cancer drugs (Sunitinib, Sorafenib, Pazopanib, Axitinib, and Temsirolimus).The results showed that the IC50 values of ve drugs were all signi cantly overexpressed in low PRS group (Fig. 3A-3E).The target genes of anticancer drugs obtained from the DrugBank database include PDGFRB, FLT3, FLT4, CSF1R, PDGFRA, RAF1, FGFR1, RET, MTOR, FGF1, SH2B3, ITK, FGF2, and FKBP1A.It was noteworthy that the expression of most target genes the high PRS group was higher than that in low PRS group (Fig. 3F).These results showed that PRS may help to select appropriate treatment options for patients and improve prognosis.

Correlation of the PCD-related signature with Clinicopathological Characteristics
In the TCGA cohort, we investigated the universality of the signature in different clinical subgroups.The PRS was signi cantly increased in advanced stage cases, advanced grade cases, advanced T stage cases, positive lymph node metastasis cases and positive distant metastasis cases (Figure S2A-S2E).Similarly, we found that patients with advanced ccRCC (grade 3 or 4, III or IV, t3 or 4, M1 and N1) were more likely to be in the high PRS group (Figure S2F-S2J).As shown in gure S2K-S2T, patients in the high-PRS group were associated with poor prognosis.This further proved that our model still had relatively accurate prediction ability in different clinicopathological stages.

Identi cation of PCD-related cluster and and Bioinformatics analysis
Based on the expression pro le of the modelled genes, the patients were divided into two clusters according to the NMF algorithm (Fig. 4A).The KM curve clearly showed that the prognosis of cluster B was worse than that of cluster A(Fig.4B).The sankey diagram further showed the correlation between classi cation, PRS and fustat.We found that most B cluster patients had high PRS, and the vast majority of deaths were in the high PRS group (Fig. 4C).This was in line with our previous analysis.The heat map showed the distribution of the 8 modeled gene expression pro les and clinicopathological characteristics (Fig. 4D).In order to further study the heterogeneity of each of the PCD-related clusters, we then performed a functional enrichment analysis for the DEGs between the clusters.KEGG and GO analysis revealed that DEGs signi cantly enriched in PI3K − Akt signaling pathway, multiple metabolic pathways, mitochondrial matrix, regulation of in ammatory response and regulation of angiogenesis (Fig. 4E and   4F).In addition, the GSVA indicated that multiple cancer-promoting pathways such as P53 signaling pathway, JAK STAT signaling pathway, CELL CYCLE, and ECM receptor interaction were signi cantly enriched in B cluster (Fig. 4G).

Mutation and Immunotherapeutic Responses of the PCDrelated signature
As TMB was signi cantly correlated with immunotherapy e cacy, we analyzed the correlation between PRS and TMB.The mutation rate was 142/184 (77.17%) in low PRS group and 126/148 (85.14%) in high PRS group (Fig. 5A and 5B).TMB was signi cantly more highly expressed in high PRS group than in low PRS group, as expected (Fig. 5C). Figure 5D revealed that PRS and TMB were signi cantly positively correlated, and cluster B had higher PRS and TMB than cluster A. The KM curve indicated that poor prognosis in high TMB group (Fig. 5E).We further analyzed the value of combining PRS and TMB in predicting the prognosis of patients with ccRCC.The KM survival curve showed that H-TMB + H-PRS had the worst prognosis and L-TMB + L-PRS had the best prognosis (Fig. 5F).Then, previous studies have shown that TMB has an association with the outcome of immunotherapy.Therefore, we further analyzed the differences in immunosuppressive checkpoint expression between the two groups.As expected, there was a signi cant over-expression of immunosuppressive checkpoints in high PRS group (Fig. 5G).Next, the relative probability of responding to CTLA4/PD-1 positive, CTLA4/PD-1 negative and CTLA4 negative/PD-1 positive treatment was higher in the low PRS group.These results showed that low PRS patients have better response to CTLA4, PD-1, and CTLA4 + PD-1 immunotherapy (Fig. 5H-5K).

Correlation of the modeled genes with Clinicopathological Characteristics
In Figure S3A, we found that the expression of 8 model genes was different between cancer and normal tissue, among which 5 were up-regulated in cancer tissue and 3 were down-regulated in cancer tissue.All 8 modeled genes had a great predictive ability according to the AUC (Figure S3B).The differential expression of the 8 genes in the different clinicopathological stages was shown in Figure S3C-S3G.We then divided patients into high-expression and low-expression groups based on median gene expression, and further analyzed survival differences between the two groups.The results showed that 5 genes were positively correlated with poor prognosis, and 3 genes were negatively correlated with poor prognosis (Figure S3-S3O).

Single-Cell Transcriptomic Context of the 8 modeled genes
For further analysis of the potential mechanism of modeled genes in ccRCC, we downloaded three singlecell sequencing datasets from the GEO database: GSE131685, GSE152938 and GSE171306.After quality control, 17,304 genes and 50,021 cells were obtained from four ccRCC and four normal samples (Figure S4A).Then, 16 different clusters and 15 groups of cells were identi ed, which included immune cells, tubular cells, endothelial cells and tumor cells (Figure S4B and S4C).In agreement with previous studies [23], renal tubular epithelial cells represented over 90% of normal renal cortical samples, whereas immune cells and tumor cells were the majority in ccRCC (Figure S4D).The majority of immune cells were identi ed in ccRCC patients, delineating the tumour immune microenvironment of ccRCC (Figure S4E).We investigated the expression pro les of eight modeled genes in different cell types (Fig. 6A-6H).We also examined the distribution of 8 genes in cancerous and paracancerous tissues.We found that 3 genes (PEBP1, NAPSA and FDX1) were highly expressed in paracancerous tissues, whereas 5 genes (SERPINE1, NOL3, CEBPB, P4HB and YBX1) were highly expressed in cancerous tissues.This was in line with our previous analysis (Fig. 6I-6P).

Validation of protein and mRNA expression of PCD-related genes in ccRCC
We obtained the protein expression data of 8 modeled genes in ccRCC from CPTAC.We found that the protein levels of 5 modeled genes (SERPINE1, P4HB, NOL3, CEBPB and YBX3) were signi cantly overexpressed in the tumor, while the remaining 3 genes (PEBP1, NAPSA and FDX1) were signi cantly underexpressed in the tumor (Fig. 7A-7H).Using qRT-PCR, differences in mRNA expression of 8 model genes were tested between ccRCC and adjacent tissue.Consistent with the trend of modeled genes protein expression, 5 modeled genes was highly expressed in ccRCC tissues and 3 genes were signi cantly underexpressed (Fig. 7I-7P).

Discussion
Traditionally, histological grade, tumour-node-metastasis (TNM) stage and clinical stage are the main factors affecting the prognosis of ccRCC.Although these indicators have been gradually improved and improved in the past few decades, they still have de ciencies in predicting patient OS and guiding anticancer treatment.New biomarkers are desperately needed for diagnosis and prognosis of ccRCC.
PCD participates in elaborate regulation, involves various mechanisms, and is crucial in tumor genesis and metastasis [24].We found that the signature we created containing eight PCD-associated genes (SERPINE1, CEBPB, NOL3, P4HB, YBX3, PEBP1, FDX1, and NAPSA) accurately predicted the overall survival of patients with ccRCC.SERPINE1, which inhibits tissue plasminogen activators and urokinase, was highly expressed in a variety of cancers, including colorectal, gastric and lung cancers, and is associated with poor prognosis, which is consistent with our ndings [25; 26; 27].CEBPB was found to be highly expressed in breast cancer, colorectal cancer and glioma,and may be involved in aerobic glycolysis to promote breast and colorectal cancer growth and metastasis[28; 29; 30].Our results also suggested that CEBPB was a pro-cancer factor in ccRCC.NOL3 encoded an anti-apoptotic protein involved in pathways such as autophagy and apoptosis regulation and signalling.He et al. found that colorectal cancer patients with high NOL3 expression had a poor prognosis [31].As a molecular chaperone, P4HB helps to improve the response of misfolded proteins to endoplasmic reticulum stress.Simultaneously, P4HB has become a new diagnostic and prognostic marker for various cancers (gastric, bladder, colorectal, etc.), and its overexpression may promote tumour progression [32; 33].YBX3 was a member of a family of transcription factors expressed primarily in endothelial cells and was essential not only in regulating cell proliferation and differentiation, but also in regulating the expression of tight junction protein, particularly in promoting the proliferation and progression of liver and breast cancer and in regulating levels of tight junction protein in melanoma and bladder cancer [34; 35].PEBP1 was a physiological endogenous inhibitor of RAF1/MEK/ERK pathway, which may inhibit cancer cell migration, proliferation and invasion of cancer cells.PEBP1 was downregulated in liver and pancreatic cancer, which may promote aggressive tumor behavior and poor prognosis [36; 37].FDX1 may not only positively regulate the speci c cuproptosis pathway, but may also be the key regulator of copper ion carrier-induced cell death [38].FDX1 down-regulated in most types of cancer and associated with better outcome [39].
Compared to the above genes, NAPSA has fewer reports, but was considered to be one of the surfaceactive metabolism-related genes related to the prognosis of primary LUAD and LUAD brain metastases.[40].
Tumor-in ltrating immune cells played a vital role in regulating cancer cell function in the TME with a high degree of consistency and plasticity, exhibiting distinct anti-tumour or pro-tumour functions [41].Previous studies and our analysis have shown that PCD was a critical regulator of the immunosuppressive tumour microenvironment (TME) [42].The appearance of PCD in the TME was paralleled by the release of intracellular components, including cytokines, mtDNA and exosomes, which were engaged in shaping the immune landscape of the TME [43; 44].These components may enrich antitumour immune cells (such as cytotoxic T and NK cells) or regulate immunosuppressive cells (such as Tregs, MDSC and macrophages), ultimately causing tumour regression/progression [45; 46] IL-1 caused chronic in ammation and promoted tumour progression and development through stimulation of epithelial-to-mesenchymal transformation, cancer cell proliferation and enrichment of immunosuppressive cell populations in TME [47].Extracellular ATP can be further converted via the exonucleotides of CD39 and CD73 on the cell membrane to the immunosuppressive metabolite adenosine, which can induce the proliferation of tumour-in ltrating macrophages [48].In our analysis, there was a signi cant over-expression of immunosuppressive cells and immunosuppressive checkpoints in the high PRS group.PRS was highly expressed in immunoassay C6 and low in immunoassay C3.Previous studies have shown that the C3 subtype has a better prognosis, while the C6 subtype has a worse prognosis.In conclusion, it was reasonable to speculate that PCD, tumour immunity and ccRCC were inextricably linked.However, the mechanism of this interaction need to be further investigated.
Our study's strength was the statistical analysis of PCD-related genes using high-throughput data and large databases, addressing the urgent need for ccRCC validation indicators.In addition, our study helped to better understand how PCD plays a role in ccRCC.Inevitably, there were some limitations to our study.First, we used traditional univariate and Lasso regression risk analyses to construct and estimate a PCDrelated risk prognostic model.Although this method has been recognised and used in many studies, more advanced methods and techniques were needed to further re ne our research in the future.Second, the clinical information in the TCGA database was not comprehensive, and we did not have access to additional parameters to validate our model, such as CT images.

Conclusion
In summary, we constructed and validated an 8-genes PCD related signature to predict the prognosis of ccRCC.This signature was signi cant for the assessment of clinical prognosis and was associated with immune in ltration.

Figure 3 Drug
Figure 3