In this study, clinical data and bulk transcriptomic data of 602 samples from the TCGA-KIRC cohort were downloaded from UCSC Xena (http://xena.ucsc.edu/, accessed on 2 January 2023) (18). The bulk transcriptomic data comprised the raw gene counts matrix and the corresponding fragments per kilobase of exon model per million mapped fragments (FPKM) matrix. The clinical data included sample pathology status, survival time, survival outcome, age, gender, and other relevant information. Additionally, microarray data from the GSE53757 (n=144) (19), GSE36895 (n=52) (20), GSE22541 (n=24) (21), and GSE67501 (n=11) (22) cohorts were obtained from the gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/, accessed on 20 January 2023). The microarray data from E-MTAB-1980 (n=101) (23) were sourced from the EMBL-EBI database (https://www.ebi.ac.uk/, accessed on 12 February 2023). Sample information provided by data contributors was used to distinguish ccRCC samples and normal tissue samples in the GSE53757 and GSE36895 cohorts. Survival outcome information for ccRCC samples in the E-MTAB-1980 and GSE22541 cohorts was extracted from the corresponding supplementary files of the respective studies. The response data of 11 ccRCC patients to nivolumab in the immunotherapy cohort GSE67501 was provided by the authors. We retained only those samples in all cohorts that were explicitly labeled as ccRCC tumors or normal kidney tissues. Detailed baseline clinical data for all ccRCC cohorts are summarized in Table 1 .

The FPKM matrix from the TCGA-KIRC cohort was transformed into a transcripts per kilobase of exon model per million mapped reads (TPM) matrix for subsequent analyses, such as the evaluation of the level of immune cell infiltration. Specifically, for each gene within a sample, the FPKM value of that gene was divided by the sum of FPKM values for all genes within that sample. Subsequently, the library size of each sample was uniformly scaled to 1 000 000.

The transcriptomic data provided by the E-MTAB-1980 cohort was directly utilized for our analysis. Additionally, we downloaded microarray data from the GSE53757, GSE36895, GSE22541, and GSE67501 cohorts. We computed the signal intensities of each probe using the affy v1.74.0 package in R with the robust multichip average (RMA) algorithm. The GSE53757, GSE36895, and GSE22541 cohorts were all based on the GPL570 platform. We utilized the AnnoProbe v0.1.7 package in R to obtain annotations for the probes of the GPL570 platform and matched probe IDs with gene names. In cases where the same gene corresponded to multiple probes, the data from the probe with the highest signal intensity were retained as the final gene expression level. The GSE67501 cohort was based on the GPL14951 platform. We manually downloaded annotation information of the GPL14951 platform from GEO and employed the same method as mentioned above to convert the probe signal intensity matrix to a gene expression matrix. In addition, the microarray data from the Caki-1 and RPTEC/TERT1 cell lines in the GSE232951 dataset (GPL17692 platform) were also processed in the same manner as described above.

We extracted normal control samples and matched KIRC samples from 71 patients in the TCGA-KIRC cohort. Using the raw gene counts matrix of the samples, we utilized the DESeq2 v1.36.0 (24) package in R to identify differentially expressed genes. Genes with | log2(FoldChange) | > 1 and an adjusted p-value< 0.05 were considered to have significantly changed expression levels. Pathway information was obtained from the kyoto encyclopedia of genes and genomes (KEGG) (https://www.kegg.jp, accessed on 5 March 2023) (25). We used the clusterProfiler v4.4.1 (26) package in R to perform gene functional enrichment analysis and gene set enrichment analysis (GSEA). Pathways with an adjusted p-value< 0.05 were considered significantly enriched in gene functional enrichment analysis. In the GSEA, we first sorted the genes according to their log2(FoldChange) from largest to smallest, and then carried out enrichment analysis. Visualization of the results from GSEA was completed using the enrichplot v1.16.1 package in R.

The differential gene expression analysis between Caki-1 and RPTEC/TERT1 cell lines was performed using limma v3.56.2.

We obtained single-cell sequencing data of 7 ccRCC tumor samples and 6 normal samples in GSE159115 (accessed on 25 February 2023) (27) and incorporated them into the analysis pipeline using the Seurat v4.1.1 package in R (28). We retained only high-quality cells with more than 200 expressed genes and mitochondrial gene expression accounting for less than 10% of the total expression. The raw gene counts matrix was normalized using the NormalizeData() function with LogNormalize method. We then selected the top 2000 highly variable genes for principal component analysis (PCA) using the FindVariableFeatures(), ScaleData(), and RunPCA() functions. Batch differences between samples were eliminated using the harmony v0.1.1 package (29). The FindNeighbors(x, reduction = “harmony”) and FindClusters() functions were employed for cell clustering. Subsequently, the FindMarkers() and FindAllMarkers(x, min.pct = 0.3) function was used to identify differentially expressed genes in each cell cluster. Genes with | logfc.threshold | > 0.1 and an adjusted p-value< 0.05 were considered as differentially expressed genes. We obtained a list of marker genes for multiple cell types from Cell Taxonomy (https://ngdc.cncb.ac.cn/celltaxonomy) (30) and annotated cell types for each cell cluster based on the differentially expressed genes in each cell cluster. The uniform manifold approximation and projection method (UMAP) algorithm was then employed to generate two-dimensional mappings of all cells.

To study the activation levels of metabolic pathways within each cell, we utilized the scMetabolism v0.2.1 package (31) for evaluation. scMetabolism encapsulates multiple algorithms for quantifying metabolic activity at the single-cell resolution and is compatible with Seurat. We used the AUCell algorithm to assess the activation levels of KEGG metabolic pathways integrated within scMetabolism.

Caki-1 and RPTEC/TERT1 cell lines were purchased from the Cell Bank of Chinese Academy of Sciences (China). Caki-1 cells were cultured in McCOY’s 5A media (Sigma-Aldrich, USA) supplemented with 10% fetal bovine serum (HyClone, Logan, UT), 1% L-glutamine (Invitrogen, Carlsbad, CA), and 1% penicillin-streptomycin (Invitrogen, Carlsbad, CA), and incubated in 5% CO₂ with a balance of air at 37°C. The RPTEC/TERT1 cells were cultured in DMEM/F12 medium (Invitrogen, Carlsbad, CA) containing 10% fetal bovine serum (HyClone, Logan, UT) and 1% L-glutamine (Invitrogen, Carlsbad, CA), and 1% penicillin-streptomycin (Invitrogen, Carlsbad, CA). Culture media were replaced every 2–3 days until the cells reached 80% confluence, and 0.25% trypsin was then added for cell passaging.

Total RNA was isolated from Caki-1 and RPTEC/TERT1 cells using RNAiso plus (Takara, 108–95-2). The RNA concentration and purity were determined through a NanoDrop One Spectrophotometer and 2100 Bioanalyzer (Agilent Technologies, Wilmington, DE, USA). A total of 500 ng RNA was reverse transcribed to cDNA and purified using the SuperScript First-Strand Synthesis System for RT-PCR (Invitrogen, 11904–018). Real-time RT-PCR was performed on the BioRAD CFX96™ Touch (Bio-Rad Laboratories Inc.), using SYBR Green Supermix (Vazyme, Q312–02). qPCR reactions using primers targeting ACAT1 (forward 5′-GGAGGTGAAGGACAAGCTCC-3′ and reverse 5′-TCTACAGCAGCGTCAGCAAA-3′), ECHS1(forward 5′- TCCTGACTGGAGCACCTTCT3′ and reverse 5′-GCATCTGTATGAAGGCAGCA3′) and ACADM (forward 5′- ACAACGTGAACCAGGATTAG-3′ and reverse 5′-TGGCAAATTTACGAGCAGTA-3′), with GAPDH (forward 5′-GCACCACCAACTGCTTA-3′ and reverse 5′-AGTAGAGGCAGGGATGAT-3′) as the loading control.

The survival v3.3.1 package in R was utilized for conducting kaplan-meier analysis, and constructing univariate Cox regression models (Cox models) and multivariate Cox models. Based on the prognostic biomarkers we found, we constructed a lipid metabolic risk score (LMRS):

β_i is the regression coefficient of gene i in the univariate Cox regression model and is used as the weight of gene i in the calculation formula. EXP_i indicates the expression level of gene i.

Univariate Cox models were built based on survival time, survival outcome, gene expression levels or the LMRS. Multivariate Cox models were constructed based on survival time, survival outcome, gene expression levels or LMRS, along with age and gender. For Cox models, a model was considered reliable when the p-value from the model’s wald test was less than 0.05. Factors with wald test p-values less than 0.05 were considered to significantly influence survival in a reliable Cox model. Kaplan-meier analysis was employed to compare the survival rates of different sample groups and to plot survival curves. The log-rank test was utilized as a non-parametric method for comparing survival rates, with a p-value less than 0.05 indicating significant differences in survival rates between two groups. Additionally, the time-dependent prognostic significance of LMRS was investigated using the timeROC v0.4.0 package in R.

We used the glmnet v4.1 package in R to build a lasso classification model using expression data from all ccRCC samples obtained from the TCGA-KIRC cohort, aiming to select important variables significantly associated with patient prognosis. Variables with non-zero coefficients in the lasso model are considered important for model.

In addition to this, the samples from the TCGA-KIRC cohort were divided into a training set and an internal validation set in a 7:3 ratio. We used the caret v6.0.92 package in R to build machine learning models. To construct the best-performing models, we included ten of the most common machine learning algorithms: bayesian generalized linear model (bayesglm), random forest (rf), neural network (nnet), k-nearest neighbors (knn), support vector machines with linear kernel (svmLinear), support vector machines with radial basis function kernel (svmRadial), classification and regression tree (rpart), linear discriminant analysis (lda), linear discriminant analysis with stepwise feature selection (stepLDA), and bagged AdaBoost (AdaBag). The caret package automatically selects the optimal values for each important parameter affecting model performance during training, thus constructing the best-performing machine learning models. Subsequently, we validated the ten trained models using the internal validation set. By plotting the receiver operating characteristic (ROC) curves, we selected the model with the highest area under curve (AUC) value as the final best machine learning model. We then applied the best model to the external validation set, and the AUC value reflected the predictive performance of the model.

We used the immunedeconv v2.0.3 (32) package in R to evaluate the immune cell infiltration levels in the TCGA-KIRC cohort. Immunedeconv provides a unified interface for multiple immune cell infiltration algorithms, such as quantiseq and cibersort, and accepts the TPM matrix as input for evaluation. We used the quantiseq (33) algorithm to assess the relative infiltration levels of 11 cell types in the TME of each sample. Additionally, we employed the cibersort (34) algorithm to assess the relative infiltration levels of 22 cell types in the TME of each sample. Furthermore, the estimate v1.0.13 (35) package in R was used to evaluate the StromalScore, ImmuneScore, and ESTIMATEScore of the TME for each sample in the TCGA-KIRC cohort. The tumor immune dysfunction and exclusion (TIDE) scores for each sample in the TCGA-KIRC cohort were obtained from TIDE (http://tide.dfci.harvard.edu/, accessed on 17 October 2023) (36).

Previously published work provided a high-quality assessment of the somatic mutation landscape in the TCGA cohort and made the results available for download in the form of MAF files (accessed on 23 June 2023) (37). We used the maftools v2.12.0 package in R to read the downloaded MAF files. Oncoplot was used to display the mutation frequencies of genes. Additionally, we evaluated the tumor mutation burden (TMB) for each sample.

Transcriptomic alterations in cancer patients strongly influence their response to anticancer drugs (38). The development of the oncoPredict v0.1.0 (39) package in R helped us predict patient responses to anticancer drugs based on the bulk transcriptomic data. Using gene expression matrices and drug response matrices from GDSC (accessed on 28 May 2023) (40) as references, we evaluated the drug sensitivity (IC50 values) of the TCGA-KIRC cohort’s ccRCC samples to 449 drugs. The assessed drug sensitivity values represented the IC50 values of the samples’ response to the drugs. A smaller IC50 value indicates a stronger response to the drug. Additionally, we obtained detailed information on the drug targets from GDSC.

All statistical analyses were conducted using R v4.2.2. The psych v2.2.5 package was used for correlation analysis. Unless otherwise specified, all correlation analyses were pearson correlations. The stats v4.2.2 package was used for related statistical tests, such as the T-test. Plotting was primarily done using packages such as pheatmap v1.0.12, ggplot2 v 3.3.6, and ggpubr v 0.4.0 in R.

In our study, we incorporated 14 lipid metabolism-related pathways from KEGG, encompassing a total of 372 lipid metabolism-related genes ( Supplementary Table S1 ). To investigate the dysregulation of lipid metabolism-related genes in clear cell renal cell carcinoma (ccRCC), we collected bulk transcriptome data from 602 samples within the TCGA-KIRC cohort (531 tumor samples; 71 normal samples) and conducted a comprehensive analysis. The analysis of differentially expressed genes (DEGs) between tumor and normal samples revealed significant upregulation of expression in 6336 genes and significant downregulation in 5293 genes within the ccRCC samples ( Figure 1A ). Functional analysis of the DEGs indicated that the upregulated genes in ccRCC were primarily enriched in pathways such as HIF-1 signaling pathway and FoxO signaling pathway, whereas downregulated genes were predominantly enriched in pathways including retinol metabolism, carbon metabolism, and fatty acid degradation ( Figures 1B, C ). Gene set enrichment analysis (GSEA) corroborated these findings, revealing significant upregulation of pathways such as the NOD-like receptor signaling pathway and FoxO signaling pathway in ccRCC ( Figure 1D ). Conversely, pathways like carbohydrate digestion and absorption, fatty acid degradation, and pyruvate metabolism were significantly associated with the phenotypes of normal samples and downregulated in ccRCC samples ( Figure 1E ). Moreover, we identified a total of 134 dysregulated lipid metabolism genes, of which 50 genes exhibited significantly increased expression and 84 genes displayed significant downregulation in ccRCC ( Supplementary Table S2 ).

We constructed univariate and multivariate Cox regression models to explore the association between dysregulated lipid metabolism genes and the prognosis of ccRCC patients ( Supplementary Table S3 ). Under the condition that the significance level of the univariate Cox model and multivariate Cox model were met at the same time, we discovered that the expression of 30 genes among 134 dysregulated lipid metabolism genes identified by differential expression analysis was significantly associated with the prognosis of ccRCC patients. Among the 30 lipid metabolism dysregulation genes associated with prognosis (LMDPs), 9 genes were implicated in the steroid hormone biosynthesis pathway, followed by 8 genes in the fatty acid degradation pathway ( Figure 1F ). We further constructed a lasso classification machine learning model to select important variables associated with prognosis ( Supplementary Table S4 ). Under the optimal lambda parameter, a total of 23 significant LMDPs were identified ( Figures 1G, H ). These LMDPs serve as our primary focus for subsequent analysis.

Bulk RNA-Seq reflects a complex transcriptional landscape of all cell populations within the tumor microenvironment (TME). To finely explore the landscape of lipid metabolism reprogramming in various cells within the tumor microenvironment (TME) of ccRCC patients, we conducted a study on single-cell transcriptomes of 7 ccRCC samples and 6 normal samples from the GSE159115 dataset. After quality control, 25,558 high-quality cells were retained and clustered into 13 major cell types: Tumor/Epithelial cell, Plasma cell, Monocyte, Memory B cell, Mast cell, Macrophage, Gamma-delta T cell, Alpha-beta T cell, Endothelial cell, Smooth muscle cell, Intercalated cell-tran-principal cell, Mesangial cell, Fibroblast ( Figures 2A, B ). High lipid absorption, storage, and fat generation occur in various cancers, contributing to the rapid growth of tumors (41). We evaluated the activation levels of 14 lipid metabolism-related pathways in various cell subtypes. Subsequently, we used the U-test to compare the differences in the activation levels of metabolic pathways between tumor samples and normal samples ( Supplementary Table S5 ; Figure 2C ). For Tumor/Epithelial cells, we found that the activation levels of 9 lipid metabolism pathways were significantly higher in tumor samples than in normal samples. In addition, lipid metabolism undergoes varying degrees of reprogramming in seven other immune cell subtypes among cancer patients. Thus, within the nutritionally deficient tumor microenvironment, as nutrient levels fluctuate, tumor cells utilize lipid metabolism reprogramming to enhance their survival environment, leading to metabolic competition that may affect the normal function of immune cells.

Based on single-cell transcriptome data, we analyzed the transcriptional variation landscape of various cell subtypes in cancer samples compared to normal samples ( Supplementary Table S6 ). We then overlapped all identified differentially expressed genes in each cell subtype with the 23 LMDPs identified through bulk transcriptome analysis, revealing 5 overlaps. These included 4 protective genes (ACADM, ACAT1, ECHS1, HPGD) and 1 risk gene (DGKZ) ( Figure 2D , Supplementary Figure S1A ). The expression levels of three protective genes (ACAT1, ACADM, and ECHS1) were significantly higher in Tumor/Epithelial cells of tumor patients than in normal samples at single-cell level. To verify this result, we collected microarray data from the ccRCC cell line (Caki-1) and the human renal epithelial cell line (RPTEC/TERT1), and found that the expression levels of three genes (ACAT1, ACADM, and ECHS1) were significantly higher in Caki-1 ( Supplementary Table S7 ; Figure 2E ). Next, we also found that the expression levels of three genes (ACAT1, ACADM, and ECHS1) were significantly higher in Caki-1 cell line than RPTEC/TERT1 through qRT-PCR experiments ( Figure 2F ). This was consistent with our findings at the single-cell level. Conversely, in adaptive immune cells of tumor samples (Alpha-beta T cell, Plasma cell, Memory B cell), the expression level of ACAT1 was significantly lower than in normal samples. Bulk transcriptome analysis based on the TME indicated that relative to normal samples, the expression levels of the 4 protective genes were significantly downregulated in ccRCC samples, while the expression level of the risk gene DGKZ was significantly upregulated ( Supplementary Figure S1B ), consistent with the differentially expressed genes results based on ANOVA analysis from GEPIA2 ( Supplementary Figure S1C ). In the breakdown of fatty acids via the beta-oxidation pathway into acetyl-CoA, the expression product of ACAT1 catalyzes the final step of this reaction. ACAT1, ACADM, and ECHS1 are all involved in the fatty acid degradation pathway. ACADM and ECHS1 also play critical roles in this beta-oxidation process. The expression product of HPGD catalyzes the dehydrogenation reaction of a series of hydroxylated polyunsaturated fatty acids, participating in the arachidonic acid metabolism pathway. The risk gene DGKZ is a diacylglycerol kinase involved in the glycerolipid metabolism and glycerophospholipid metabolism metabolic pathways. In summary, these 5 key LMDPs identified through bulk and single-cell transcriptome analysis may play important roles in lipid metabolism reprogramming within the TME and could serve as prognostic signatures for ccRCC.

Renal cell carcinoma has been described as a metabolic disease, with alterations in metabolism leading to diverse cancer etiologies, particularly changes in lipid metabolism (42). Therefore, lipid metabolism may hold potential value in the diagnosis of ccRCC occurrence and progression. To explore the applicability of lipid metabolism-related genes in predicting the occurrence and prognosis of ccRCC, we constructed machine learning models to assess their diagnostic performance. Leveraging the 5 key prognostic signatures identified through combined bulk and single-cell transcriptome analysis, we initially developed a machine learning model for predicting the occurrence of ccRCC, distinguishing ccRCC samples and normal samples. We partitioned the entire TCGA-KIRC cohort into a training set and an internal validation set in a 7:3 ratio, and found no significant differences in the distribution of various clinical features between the two datasets ( Figure 3A ). To build the optimal machine learning model, we included 10 of the most common machine learning algorithms. Under the optimal parameters for each model, the bayesian generalized linear model (bayesglm) exhibited the best performance, with an AUC value of 0.946 for predictions on the internal validation set ( Figure 3B ). For the bayesglm model, the risk gene DGKZ is the most important factor ( Figure 3C ). Furthermore, to further evaluate the performance of the model we constructed, we included two external validation datasets, GSE53757 (72 ccRCC tumor samples; 72 normal kidney samples) and GSE36895 (29 ccRCC tumor samples; 23 normal kidney cortex samples). In both datasets, our bayesglm model effectively distinguished ccRCC samples and normal samples, with AUC values of 0.938 ( Figures 3D ) and 0.879 ( Figures 3E ), respectively, indicating the robustness and accuracy of the model.

The results of survival analysis from HPA (www.proteinatlas.org/) confirmed the association between the expression of five key prognostic signatures and the survival of ccRCC patients, which was consistent with our findings ( Supplementary Figure S2 ). To evaluate the potential of prognostic signatures we identified in predicting ccRCC patient prognosis, we employed a similar approach to construct a machine learning model for predicting the survival outcomes (death or survival) of ccRCC patients. We retained the cancer samples from the TCGA-KIRC cohort and split the dataset into a training set and an internal validation set in a 7:3 ratio, with no significant differences in the distribution of various clinical features between the two datasets ( Figure 3F ). Under the optimal parameters for each of the 10 machine learning algorithms, the neural network (nnet) model exhibited the best performance, with an AUC value of 0.831 for predictions on the internal validation set ( Figure 3G ). Regarding the constructed nnet model, the protective gene ACADM demonstrated the highest importance ( Figure 3H ). To validate the model’s performance, we included two external validation datasets, GSE22541 (6 dead samples; 18 alive samples) and E-MTAB-1980 (23 dead samples; 78 alive samples). Both datasets investigated the transcriptomes of multiple ccRCC patients and provided clear survival outcome information for each sample. Based on the transcriptome data from the two cohorts, we utilized the constructed nnet model to predict the survival outcomes of each sample. We found that the nnet model we developed performed well in both datasets, with AUC values of 0.806 ( Figures 3I ) and 0.762 ( Figures 3J ), respectively. These findings suggested that prognostic signatures we found had great application value in indicating the occurrence of ccRCC and patient prognosis.

To further explore the clinical application value of lipid metabolism-related prognostic signature, we constructed the lipid metabolism risk score (LMRS) as a risk index in ccRCC. LMRS is a simple linear model based on the gene expression levels of key prognostic signatures and their univariate Cox regression coefficients. The calculation formula for LMRS is as follows: Risk score = 0.0272 * Exp (DGKZ) + (−0.0059) * Exp (ACADM) + (−0.0033) * Exp (ACAT1) + (−0.0012) * Exp (ECHS1) + (−0.0186) * Exp (HPGD). Using the median LMRS of the patients, the ccRCC samples from the TCGA-KIRC cohort were divided into two groups: High-LMRS and Low-LMRS ( Supplementary Table S8 ). Kaplan-meier analysis revealed a significantly poorer prognosis in the High-LMRS group (log-rank test, p value< 0.0001) ( Figure 4A ). Subsequently, we constructed a univariate Cox regression model and found that LMRS was significantly associated with poor prognosis in the KIRC samples (HR = 3.56, p value< 0.0001). In the multivariate Cox regression model considering the LMRS, patient age, and patient gender, the relationship between LMRS and poor prognosis remained significant ( Figure 4B ). Furthermore, there were no significant differences in LMRS among different age and gender groups ( Supplementary Figures S3A, B ). These findings suggest that LMRS is an independent prognostic factor and is significantly associated with poor prognosis in ccRCC patients. Time-dependent ROC curve analysis revealed that the LMRS predicted the one-year, three-year, five-year, and ten-year survival outcomes of ccRCC patients with AUC values of 0.693, 0.669, 0.719, and 0.710, respectively ( Figure 4C ). This indicates that LMRS has good sensitivity and specificity for predicting patient prognosis. Then, we also observed that the high LMRS group had a lower survival probability in the E-MTAB-1980 (log-rank test, p value = 0.018) and GSE22541 cohorts (log-rank test, p value = 0.092) ( Supplementary Table S9 ; Figure 4D ).

Genomic variations are not only the driving forces behind the development of tumors but also major characteristics of cancer (43). The proportion of mutation events was found to be higher in the High-LMRS group (92.98%) compared to the Low-LMRS group (87.96%) ( Figure 4E ). The most common mutations in both the High-LMRS and Low-LMRS groups were VHL, PBRM1, and TTN. The expression product of the von Hippel-Lindau tumor suppressor (VHL) plays a critical role in cellular oxygen sensing by targeting hypoxia-inducible factors for ubiquitination and proteasomal degradation, and it has emerged as a potential therapeutic target for the treatment of late-stage ccRCC (44). Polybromo 1 (PBRM1) encodes a subunit of the ATP-dependent chromatin remodeling complex, and its mutations are associated with primary ccRCC (45, 46). Titin (TTN) encodes a highly abundant sarcomeric protein, and its mutations are implicated in the development of various cancers (47–49). Tumor mutational burden (TMB) refers to the number of variant bases per million bases in tumor tissue and is an emerging biomarker increasingly used for predicting patient prognosis. We found a significant positive correlation between TMB and LMRS (Coef = 0.21, p value< 0.001) ( Supplementary Table S10 ; Figure 4F ). Kaplan-meier analysis revealed that ccRCC patients in the High-TMB group had a poorer prognosis (log-rank test, p value< 0.0001) ( Figure 4G ). Furthermore, when combining TMB and LMRS, a better prediction of prognosis was achieved, as the ccRCC samples in the High-TMB with High-LMRS group exhibited significantly poorer prognosis (log-rank test, p value< 0.0001) ( Figure 4H ).

Different types of infiltrating immune cells in the tumor microenvironment (TME) have varying impacts on tumor progression, which may differ depending on the cancer type (50). Clinical studies on immune cells infiltrating the TME have confirmed the roles of cytotoxic T cells and tumor-associated macrophages in cancer progression (51, 52). To further investigate the immune characteristics of different LMRS groups, we employed the quantiseq algorithm to evaluate the relative infiltration levels of 11 immune cell types (B cell, macrophage M1, macrophage M2, monocyte, neutrophil, natural killer (NK) cell, T cell CD4+ (non-regulatory), T cell CD8+, T cell regulatory (Tregs), myeloid dendritic cell, uncharacterized cell) in the TCGA-KIRC cohort of ccRCC samples ( Supplementary Table S11 ). The uncharacterized cell was not included in our subsequent analysis. We observed a significant positive correlation between LMRS and the infiltration levels of macrophage M1, T cell CD8+, T cell regulatory (Tregs), myeloid dendritic cell, and macrophage M2 ( Figure 5A ), while a significant negative correlation was found with the infiltration levels of NK cells, B cells, and neutrophils. Furthermore, compared to the Low-LMRS group, the High-LMRS group exhibited significantly higher infiltration levels of macrophage M1, T cell CD8+, and T cell regulatory (Tregs), and significantly lower infiltration levels of neutrophils and NK cells ( Figure 5B ).

In numerous human malignancies, the presence of T cells has been associated with improved patient prognosis (53). However, we found a significant positive correlation between LMRS and the infiltration levels of T cell CD8+ (Coef = 0.24, p value< 0.001) ( Figure 5C ), with the High-LMRS group exhibiting a poorer prognosis. To elucidate this phenomenon, we re-evaluated the immune infiltration levels of the TCGA-KIRC cohort of ccRCC samples using the cibersort algorithm and similarly observed a significant positive correlation between LMRS and T cell CD8+ infiltration levels (Coef = 0.23, p value< 0.001) ( Figure 5D ). Indeed, research evidence suggests elevated levels of T cells CD8+ in late-stage ccRCC, impacting patient response to immunotherapy (54). The TIDE algorithm evaluated the tumor-infiltrating cytotoxic T cell exclusion score and dysfunction score. We found that the High-LMRS group exhibited significantly higher T cell dysfunction scores (T test, p value< 0.0001). Although the T cell exclusion score was higher in the High-LMRS group, there was no statistically significant difference between the two groups ( Figure 5E ). Additionally, studies have indicated a link between lipid accumulation and T cell dysfunction (55). Therefore, we speculate that the tumor-infiltrating T cell CD8+ in ccRCC samples may primarily be in an immunosuppressive state, unable to exert normal cytotoxic functions, thus affecting patient prognosis. Furthermore, we utilized the ESTIMATE algorithm to assess the stromal score, immune score, and overall ESTIMATE score of the TCGA-KIRC cohort of ccRCC samples, finding all three scores significantly higher in the High-LMRS group ( Figure 5F ) (T test, p value< 0.0001). Overall, the High-LMRS group exhibited higher levels of immune cell infiltration in the TME, potentially benefiting more from immunotherapy.

Therapeutic interventions targeting lipid metabolism can restrain tumor cell growth, and alleviate immune suppression within the tumor microenvironment, thereby enhancing responsiveness to immune checkpoint blockade therapy (56). To evaluate the predictive potential of LMRS for the efficacy of immune checkpoint blockade therapy, we collected 36 widely utilized immune checkpoint genes (57). We found a widespread and significant positive correlation between LMRS and immune checkpoint genes in the TCGA-KIRC cohort ( Figure 6A ). Additionally, in the analysis of the correlation with immune checkpoint genes, LMRS exhibited a similar landscape to the risk gene DGKZ. Moreover, PDCD1 (PD1), as an immune inhibitory receptor, is expressed in various types of immune system cells, particularly cytotoxic T cells (58). The expression of PDCD1 and CTLA4 is commonly associated with improved immunotherapeutic efficacy (59, 60). We observed a significant positive correlation between the two classical immune checkpoints (PDCD1 and CTLA4) with LMRS ( Figures 6B, C ). These findings suggest that higher LMRS might be associated with improved responsiveness to immune checkpoint blockade therapy.

During the course of patient response to immune checkpoint blockade therapy, the support of products of certain essential genes is required. These essential genes play a crucial role in antigen presentation and the cytotoxic capability of T cell CD8+, directly impacting the efficacy of immune checkpoint blockade therapy (61). We compiled 18 essential genes for immune checkpoint blockade therapy (62) and found a significant positive correlation between LMRS and the expression of multiple essential genes ( Figure 6D ). Subsequently, we incorporated a cohort of immune therapy (GSE67501), which studied the transcriptional landscape of 11 ccRCC patients before receiving nivolumab (anti-PDCD1) treatment and carefully documented the response of each patient to immunotherapy after nivolumab treatment (7 non-responsive samples and 4 responsive samples) ( Supplementary Table S12 ). We found that the expression of PDCD1 in the responsive group was higher than in the non-responsive group, although not statistically significant ( Figure 6E ). Meanwhile, the LMRS of the responsive group was significantly higher than the non-responsive group (T test, p value = 0.030). In conclusion, LMRS can be used to predict the response of ccRCC patients to immune checkpoint blockade therapy. Patients with higher LMRS are more likely to benefit from immunotherapy.

Current treatments for ccRCC often involve cytokine therapy and targeted drug therapy, yet they fail to generate a universal and enduring complete remission response (63). To explore the guidance significance of LMRS in anticancer drug therapy, we employed oncoPredict to predict the drug sensitivity (IC50 values) of 449 drugs included in GDSC for ccRCC samples in the TCGA-KIRC cohort ( Supplementary Tables S13 , S14 ). Lower IC50 values indicate higher sensitivity of the samples to anticancer drugs. Our study identified a significant correlation between LMRS and IC50 values of 348 drugs, with 90.23% showing a significant negative correlation ( Supplementary Table S15 ; Figure 7A ). Considering the target information of the drugs, we identified 9 drugs targeting metabolic pathways. Among them, LMRS exhibited a significant negative correlation with the IC50 values of 6 drugs targeting metabolic pathways (BX-912, AICA ribonucleotide, DMOG, OSU-03012, daporinad, CAP-232), and only a significant positive correlation with AGI-6780 ( Figure 7B ).

Sunitinib, an orally active small molecule tyrosine kinase inhibitor, with potent anti-angiogenic and antitumor activities, is widely used as a first-line treatment for renal cell carcinoma (64). Additionally, other drugs such as axitinib are also commonly used in the first-line treatment of ccRCC patients (6, 7). We found a significant negative correlation between LMRS and the IC50 values of sunitinib, axitinib, temsirolimus, rapamycin, and pazopanib ( Figure 7C ). While LMRS exhibited a negative correlation with the IC50 value of sorafenib, it was not statistically significant.

The structural heterogeneity of the tumor microenvironment in ccRCC patients leads to different clinical outcomes and postoperative recurrence risks, thus emphasizing the urgent need for personalized treatment strategies (65, 66). We further compared the drug sensitivity between the High-LMRS and Low-LMRS groups, discovering that the High-LMRS group showed greater sensitivity to 178 drugs such as gemcitabine (pyrimidine antimetabolite targeting DNA replication) and irinotecan (topoisomerase I inhibitor targeting DNA topology) ( Figure 7D ). In contrast, the Low-LMRS group demonstrated increased sensitivity to 86 drugs, including vincristine (microtubule assembly inhibitor) and WYE-125132 (mTOR inhibitor). Statistical analysis of the targeted pathways for each drug revealed that the sensitive drugs in the High-LMRS group primarily targeted kinases and DNA replication ( Figure 7E ). Conversely, the sensitive drugs in the Low-LMRS group mainly targeted RTK signaling, EGFR signaling, and PI3K/MTOR signaling pathways. In summary, these findings indicate that the LMRS we constructed could offer guidance for targeted therapy, facilitating the development of personalized treatment regimens for ccRCC.