Comprehensive analysis of LD-related genes signature for predicting prognosis and immunotherapy response in clear cell renal cell carcinoma

Abstract

Background

Lipid droplets (LD) in renal clear cell carcinoma (ccRCC)play a crucial role in lipid metabolism and immune response modulation. The purpose of this study was to create a LD-related signature to predict prognosis and guide the immunotherapy and targeted therapy in ccRCC patients.

Methods

We conducted a comprehensive analysis using transcriptional profiles and clinical data obtained from The Cancer Genome Atlas (TCGA). LD-related genes were identified from existing literature and the GeneCards database, and differentially expressed genes were determined. Sequentially, we conducted Cox regression analysis and Lasso regression analysis, to establish a prognostic risk model. The performance of the risk model was evaluated using Kaplan–Meier (KM) analysis and time-dependent receiver operating characteristic (ROC) analysis. Additionally, gene set enrichment analysis (GSEA), ESTIMATE, CIBERSORT, and immunophenoscore (IPS) algorithm were used to assess the tumor microenvironment (TME) and treatment response.

Results

We constructed a risk signature with four LD-related genes in the TCGA dataset, which could be an independent prognostic factor in ccRCC patients. Then, patients were classified into two risk groups and exhibited notable differences in overall survival (OS), progression-free survival (PFS), and TME characteristics. Furthermore, we developed a comprehensive nomogram based on clinical features, which demonstrated good prognostic predictive value. According to the results of GSEA analysis, immune-related pathways were found to be significantly enriched in the high-risk group. Additionally, the high-risk group displayed high levels of immune cell infiltration, TMB and IPS scores, indicating better efficacy of immune checkpoint inhibitors (ICIs). Finally, high-risk demonstrated reduced IC50 values compared to the low-risk counterpart for specific targeted and chemotherapeutic drugs, suggesting that the patients receiving these targeted drugs in high-risk group had better treatment outcomes.

Conclusions

Our findings suggested that the LD-related gene signature could potentially predict the prognosis of ccRCC patients. Additionally, it showed promise for predicting responses to immunotherapy and targeted therapy in ccRCC patients. These insights might potentially have guided the clinical management of these patients, but further validation and broader data analysis are needed to confirm these preliminary observations.

Introduction

Renal cell carcinoma (RCC), a common malignant tumor of renal cells and renal tubular epithelial cells, constitutes 3% of all malignant diseases and is one of the most lethal malignancies in the urinary system [1]. Clear cell renal cell carcinoma (ccRCC), as the predominant subtype, accounts for 75–80% of all RCC cases and displays high invasiveness and recurrence rates [1,2,3].The incidence and mortality rates of ccRCC have seen a rapid rise in recent decades. Currently, surgical resection combined with adjuvant systemic therapy is the primary treatment method for ccRCC patients, but many patients experience tumor recurrence or metastasis after surgery [4]. In recent years, immune checkpoint inhibitors (ICIs) such as PD-L1, PD-1, and CTLA-4 have made progress in the treatment of ccRCC [5]. Nevertheless, for metastatic RCC, about 75% of patients develop resistance to immune checkpoint blockade (ICB) therapy, and even those who initially respond well eventually face disease progression [6, 7]. Therefore, finding new reliable prognostic biomarkers is crucial for developing immunotherapy.

Lipid droplets (LDs) are enveloped by a phospholipid monolayer and accompanied by LD surface proteins and function as organelles accountable for storing neutral lipids, mainly comprising cholesterol esters (CEs) and triglycerides (TGs) [8]. Previous studies have consistently demonstrated the notable accumulation of LDs in specific cancer cells, suggesting their involvement in tumor proliferation, invasion, metastasis, and chemoresistance across various malignancies [9,10,11]. ccRCC is characterized by an anomalous buildup of lipid droplets within the cytoplasm. The lipid droplets in ccRCC tissues contain significantly higher levels of free cholesterol and esterified cholesterol, approximately 8 times and 35 times higher than those in normal kidney tissues, respectively. Lipid accumulation is closely associated with disease progression [12,13,14]. Additionally, ccRCC, as one of the solid tumors characterized by extensive immune infiltration, manifests substantial leukocyte infiltration and displays a close correlation with immune responses [15]. LDs play an active role in immune metabolism and immune signal transduction as primary respondents of innate immunity. Furthermore, they serve as central hubs for metabolic-immune system integration, consequently contributing to anti-tumor immunity [16,17,18]. Therefore, LDs may be implicated in immune cell infiltration within KIRC tissues and possess therapeutic potential for tumor management. However, the gene set associated with LDs in KIRC has not been systematically investigated.

In the study, we established novel characteristics of LD-related genes based on the TCGA cohort. Subsequently, and then developed a prognostic model to explore the clinical relevance and potential prognostic implications of these genes in ccRCC. Our primary objective was to assess the clinical prognosis of KIRC and explore possible implications for immune-related and targeted therapies.

Materials and methods

Exploration of the lipid droplets associated genes

A systematic search was conducted on the PubMed database using the search term "lipid droplet-associated protein," with the final search executed on November 27, 2023. Publications containing the keyword "LD" were gathered, and factors linked explicitly to LDs or associated with organelles interacting with LDs, including mitochondria, endoplasmic reticulum, and the Golgi apparatus, were selected. No restrictions were applied to publication dates. From this search, we obtained 58 genes. Additionally, we retrieved 104 LD-associated genes from GeneCards (https://www.genecards.org/) using the keywords "LD-associated genes" and "LD-associated protein," with a relevance score ≥ 7 as the filtering criterion. Finally, we compiled a set of 162 genes associated with LDs from both databases.

Data collection

Clinical and transcriptome data from 539 ccRCC and 72 healthy kidney specimens were obtained from The Cancer Genome Atlas (TCGA) data portal. After removing duplicate samples, a total of 530 KIRC samples were collected. Among these, 513 samples with overall survival time (OS) greater than 0 days were randomly divided into a training set of 333 cases and a testing set of 180 cases, ensuring similar distributions of clinical features (Table S4). Additionally, the predictive performance of the model was validated using the E-MTAB-1980 dataset, which includes 101 samples, and the ICGC dataset, which includes 90 samples.

Identification and enrichment analysis of differentially expressed genes

LD-associated differentially expressed genes (DEGs) between the tumor and normal samples in the TCGA-KIRC cohort were identified using the "limma" R package, with thresholds of |log2 fold change|> 1 and P < 0.05. These DEGs were then subjected to Gene Ontology (GO) [19] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [19] enrichment analysis using the R package "clusterProfiler" [20].

Development and validation of a risk model based on LD-associated Genes in KIRC

We performed univariate Cox regression analysis to assess the association between LD-associated genes and the prognosis of KIRC. LD-related genes with negligible contribution were removed through least absolute shrinkage and selection operator (LASSO) analysis [21]. Next, we developed a risk model using multivariate Cox regression analysis. The risk score was calculated using the formula, ${\sum }_{i=1}^{n}expi*coefi$, where expi and coefi represented gene expression levels and regression coefficients, respectively. Subsequently, four genes were selected to construct an LD-related gene signature, and their differential expression was validated using the TCGA database, UALCAN and HPA databases (https://www.proteinatlas.org/; https://ualcan.path.uab.edu/index.html). Then, KIRC patients were stratified into high-risk and low-risk groups according to the median risk score. We confirmed the effectiveness of the risk model in separating groups using principal component analysis (PCA). The predictive performance of the signature was evaluated by time-dependent receiver operating characteristic (ROC) curves and Kaplan–Meier (K-M) curves.

Relationship between risk model and clinical characteristics

To assess the association between the risk model and clinical/pathological factors (age, sex, stage, grade, and T), we performed univariate and multivariate Cox regression analyses. Additionally, we analyzed differences in risk scores across various clinical features. Stratified analyses were conducted to evaluate the prognostic significance of the risk score in different clinicopathological subgroups.

Development and evaluation of the nomogram

We developed a nomogram to predict OS in KIRC based on independent prognostic factors. To assess its accuracy, calibration curves at 1, 3, and 5 years and ROC curves were used. Decision curve analysis (DCA) was performed to compare the net benefit of the comprehensive nomogram against a model with only clinical variables.

Gene set enrichment analysis (GSEA)

GSEA was conducted using the "clusterProfiler" R package to delineate the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between different risk groups in the TCGA-KIRC dataset [22]. An adjusted p-value < 0.05 was considered statistically significant.

Immune infiltration analysis

The ESTIMATE algorithm was used to calculate stromal, immune, and estimate scores for all TCGA-KIRC cases [23]. ssGSEA and Microenvironment Cell Populations (MCP) methods were applied to analyze immune-related functions and immune cell infiltration between high- and low-risk groups [24, 25]. We also compared immune checkpoint expression between the two risk groups to assess potential utility in immunotherapy.

Somatic variant analysis and mutation landscape

Somatic mutation data for KIRC patients were obtained from the TCGA database to calculate individual TMB scores [26]. We then compared TMB differences between high- and low-risk groups and conducted survival analysis. The mutation landscape of KIRC patients was visualized using a waterfall plot with the "maftools" R package [27].

IPS analysis

The Immunophenoscore (IPS) is a recognized predictor of immune checkpoint inhibitor (ICI) response [28, 29]. Using IPS data from The Cancer Immunome Atlas (TCIA) (https://tcia.at/), we assessed potential differences in immunotherapy response between high- and low-risk groups.

Comparisons of drug sensitivity

We used the "pRRophetic" R package to evaluate the risk model's predictive value for sensitivity to targeted and chemotherapeutic drugs [30, 31].

Statistical analysis

Statistical analyses were conducted using R software (version 4.2.2). Wilcoxon's signed-rank test was utilized to assess distinctions between two groups, while Kaplan–Meier curves and log-rank tests were employed to scrutinize survival disparities among various risk groups. A two-sided p-value < 0.05 was considered statistically significant.

Result

Identification of genes associated with lipid droplets based on the TCGA dataset

A total of 162 genes associated with lipid droplets (LDs) were identified through literature search and database analysis. Table S1 displays the genes associated with lipid droplets retrieved from the literature. We utilized the "limma" R package to conduct differential expression analysis on these 162 genes in KIRC. P-value < 0.05 was considered indicative of significant differential expression and set a fold change threshold of |log2 FC|> 1. The analysis identified 134 genes with significant expression changes, comprising 59 upregulated genes and 75 downregulated genes. (Table S2). The data was visualized using a heatmap (Fig. 1A). GO enrichment analysis demonstrated a close correlation between genes associated with lipid droplets and several functions, including lipid localization, lipid storage, lipase activity, carboxylic ester hydrolase activity, O-acyltransferase, lipid droplets and endocytic vesicles (Fig. 1B-D). KEGG enrichment evaluation verified that these genes were related to the PPAR signaling pathway, endocytosis and fat digestion absorption (Fig. 1E). The findings indicated that the gene set we identified was associated with lipid metabolism, further validating the role of lipid droplets and lipid metabolism in the progression of ccRCC.

Fig. 1

Functional enrichment analysis based on DEGS: A Differential expression heatmap of 134 DEGs. GO analysis (B) BP: biological processes, C MF: molecular function. D CC: cellular components. E KEGG pathway analysis

Construction of a risk model

First, we conducted univariate Cox regression analysis on the 134 DEGs in the training dataset, identifying 64 genes significantly associated with OS (P < 0.05). (Table S3) LASSO Cox regression analysis was applied to select these survival-related genes (Fig. 2A-B), followed by univariate and multivariate Cox regression analysis for these genes (Fig. 2C-D). Finally, a risk model based on four genes associated with lipid droplets was constructed in the train datase. The risk score was computed by considering both the coefficients and the expression levels of the four genes: (0.820 * expression of AUP1) + (0.264 * expression of SERPINF1) + (0.717 * expression of DGAT2) + (-0.286 * expression of METTL7A). Subsequently, these patients were stratified into distinct risk groups based on a median risk score of -0.1087. The heatmap and risk curve showed a significant increase in the mortality rate of patients in the high-risk group in the TCGA-KIRC and E-MTAB-1980 cohorts (Fig. 2E-F). The results suggested that the risk score might have potential as a prognostic biomarker for indicating the survival outcomes of ccRCC patients.

Fig. 2

Construction of LD-related genes risk model: A LASSO coefficient profiles of the 71 genes in the train dataset. B Selection of the optimal parameter (lambda) in the LASSO model. C Univariate and (D) multivariate Cox regression analysis of the four LD-associated genes in the risk mode. Risk scores and survival status of KIRC patients among diverse risk groups of (E) TCGA cohort and (F) E-MTAB-1980 cohort

Evaluation of the risk signature of LD-related genes

The process of selecting lipid droplet-associated genes for bioinformatics analysis was conducted. Primarily, PCA analysis of train dataset showed the signature had good performance in clustering. (Fig. 3A) Subsequently, the prognostic accuracy of the risk model was validated using entire datasets, train datasets, test datasets, ICGC datasets and MTAB-1980 datasets. The K-M plots demonstrated that ccRCC patients with a high-risk score exhibited reduced Progression-Free-Survival (PFS) in the TCGA-KIRC cohort (Fig. 3B-D), as well as lower OS in the TCGA-KIRC (Fig. 3E-G). Moreover, excellent results were observed in the area under the ROC curve (AUC) at 1-, 3-, and 5-year intervals. The values in the entire cohort were 0.764, 0.704, and 0.713 (Fig. 3H); in the train cohort, they were 0.761, 0.715, and 0.702 (Fig. 3I); in the test cohort, they were 0.751, 0.721, and 0.741 (Fig. 3J). The model's predictive ability was evaluated using ROC and Kaplan–Meier curves in the two datasets. Within the datasets of E-MTAB-1980 and ICGC, Kaplan–Meier curves elucidate inferior survival outcomes within the high-risk cohort, and the results of the time-dependent ROC curve analysis also indicate that the model possesses a predictive capability. (Figure S1) The findings indicate that our model exhibited high sensitivity and specificity in diverse datasets, enabling accurate prognosis prediction for patients with KIRC. Therefore, we further visualized the expression differences and survival probabilities of these genes in TCGA-KIRC samples. We found significant disparities in both aspects for these four genes (Figure S2). Additionally, we obtained validation from the UALCAN and HPA websites showing the differential expression and prognostic significance of these four genes in KIRC tissues (Figure S3).

Fig. 3

Validating the performance of the risk model: A PCA analysis of risk model analyses in the TCGA dataset. B-D PFS curves, E–G Kaplan– Meier curves and H-J the AUCs of the time-dependent ROC curves for the entire dataset, the train dataset and the test dataset

Relationships between the risk model and clinical features

Univariate and multivariate Cox regression analyses were conducted on the entire dataset's risk score and clinical features (Figs. 4A-B). To address potential confounding factors, variables such as age, gender, tumor stage, T, M, and N were included as comprehensively as possible in the multivariate Cox regression model. The results showed that our risk model could independently predict KIRC prognosis. Furthermore, the AUC result of 0.704 demonstrates the good predictive ability of the risk score (Fig. 4C). Box plots revealed that Stage, Grade, T, M and N had high risk scores (Figure S4). To assess the clinical advantage of the LD-related genes signature, four independent prognostic indicators (risk score, stage, age, T) derived from the multivariate Cox analysis were further explored in the nomogram analysis. The results showed that as the total risk score increased, the prognosis worsen (Fig. 4D). Calibration curves indicated that the nomogram performed similarly to the ideal model in predicting prognosis probability (Fig. 4E). Moreover, the time-varying ROC curve underscored the precision of the nomogram in forecasting prognosis (Fig. 4F). The application of decision curve analysis (DCA) was employed to evaluate the clinical applicability of the nomogram and revealed that compared to single clinical characteristics, the nomogram can generate more net benefits (Fig. 4G). Additionally, we further conducted stratified analyses to validate the prognostic value of the risk scores in subgroups characterized by different clinical features, demonstrating the robustness of our findings. The survival analysis of distinct subgroups according to the clinical features showed significant differences between two risk groups except for G1 + 2, and the patients in high-risk groups had poorer OS (Figs. 5A-M).

Fig. 4

Development and assessment of the nomogram. A Univariate and (B) multivariate Cox regression analyses of the risk score and clinical characteristics. C Comparison of ROC curves among clinical feature and risk score. D The nomogram and the (E) calibration curve analyses were performed to predict 1-, 3-, and 5-year OS according to risk score. F The time-dependent ROC curves of the nomogram. G DCA of the clinical usefulness of the constructed nomogram

Fig. 5

Kaplan–Meier curves of clinical subgroups. A age < 65, B age ≥ 65, C male, D female, E Grade 1 + 2, F Grade 3 + 4, G Stage I + II, H Stage III + IV, I T I + II, J T III + IV, K M1, L M0, M N0

Exploration the potential mechanisms between the two risk groups by GSEA

The GO assessment results showed enrichment of terms such as " humoral immune response mediated by circulating immune”,"mmunoglobulin production”,"T-cell receptor complex", "antigen binding " and " humoral immune response” (Fig. 6A). KEGG results revealed a high pathway enrichment including "cytokine-cytokine receptor interaction", "hematopoietic cell lineage", "allograft rejection", "graft versus host disease " and "type I diabetes mellitus” (Fig. 6C). These results suggested that the high-risk group was notable enriched in immune-related processes.

Fig. 6

Gene set enrichment analyses between the high- and low- risk groups. GO enrichment analyses in the high- risk group (A) and the low- risk group (B). KEGG pathway analyses in the high- risk group (C) and the low- risk group (D)

Immune landscape variations between the two-risk groups

The GSEA outcomes indicated that the risk model had a strong correlation with the immune response. To confirm this conclusion, we conducted immune analysis employing to investigate the tumor microenvironment in the TCGA-KIRC cohort by the ESTIMATE and CIBERSORT algorithms., ESTIMATE analysis revealed that the high-risk group had higher StromalScore, ImmuneScore and ESTIMATEScore (Fig. 7A). In addition, we systematically evaluated 22 types of immune cells and immune functions using ssGSEA. The deconvolution algorithm CIBERSORT showed that the high-risk group had significantly higher proportions of Plasma cells, T cells CD8, T cells follicular helper, T cells regulatory, NK cells activated and Macrophages M0 (Fig. 7B). Furthermore, the majority of immune Function aspects ADCS, APC co-stimulation, CCR chemotaxis, Cytolytic activity, T_cell_co-inhibition, and co-stimulation—exhibited higher scores within the high-risk group, whereas the Type II IFN Response displayed elevated scores in the low-risk group (Fig. 7C). We quantified cytotoxic lymphocytes, CD8 + T cells, B cells, endotheliocytes, mononuclear cells, neutrophils, myeloid dendritic cells and fibroblasts using MCF counter. B lineage, CD8 T cells, and Cytotoxic lymphocytes and Fibroblasts increased in the high-risk group, while Endothelial cells, Neutrophils, and NK cells increased in the low-risk group (Fig. 7E-K). Moreover, the risk score was positively correlated with B lineage, CD8 T cells, Cytotoxic lymphocytes, and fibroblasts, but negatively correlated with myeloid dendritic cells, NK cell endotheliocytes, and neutrophils (Fig. 7D). These findings suggested significant disparities in the immune landscape between the high- and low-risk groups.

Fig. 7

The immune landscape between the two risk groups. A ImmuneScore, StromalScore and ESTIMATE score, between two risk groups. Immune cells (B) and immune functions (C) infiltration between the two risk groups. D Correlation between the risk scores and the immune cells. The different expression of B lineage cell (E), CD8 + T cell (F), cytotoxic lymphocytes (G), endothelial cells (H), fibroblasts (I), neutrophils (J), and NK cells (K) between the two risk groups by MCP counter

Correlations between somatic mutation and risk model

TMB stands as an emerging biomarker that is progressively employed for predicting patient prognosis. We divided patients with TMB information into high-risk and low-risk groups, further classifying samples as H-TMB and L-TMB based on the median TMB score. (Table S5) The analysis showed that patients in the high-risk group had higher TMB (Fig. 8A) and poorer OS (Figs. 8B). The combination of TMB and risk score seemed to have good risk stratification (P < 0.001, Fig. 8C) The mutation information and frequencies of the 15 genes in two risk groups were shown in Figs. 8D-E. Among these genes, the mutation rates of VHL (Von Hippel-Lindau Tumor Suppressor), PBRM1 (Polybromo 1), and TTN (Titin) exceeded 10% in both groups, which was consistent with previous research in KIRC. In addition, four prevalent checkpoint genes (LAG3, CTLA-4, TIGIT, and PDCD1) exhibited a significant elevation in the high-risk group and were positively correlated with the risk score (Figs. 8F-G). The findings indicated that individuals in the high-risk group might experience potential benefits from immunotherapy.

Fig. 8

Somatic mutational analyses between high- and low- risk groups. A TMB difference between the two risk groups. B Survival analysis of OS between H-TMB and L- TMB groups. C Survival analysis of OS stratified by TMB and risk groups. The 15 frequently mutated genes in the high-risk group (D) and low-risk group (E). F The correlationp between (F) risk scores and the expression levels of the immune checkpoints. G The different expression of tbetween the immune checkpoints between the two risk groups

Comparisons of drug sensitivity

To ascertain whether increased TMB within the high-risk group corresponds to improved outcomes with ICIs treatment, we evaluated the immunotherapy response in the two risk groups using the IPS algorithm.The patients in the high-risk group with checkpoint inhibitor treatment had notable therapeutic advantages (Figs. 9A-D). We also analyzed the therapeutic efficacy of targeted drugs and common chemotherapy drugs in the risk model to investigate the clinical utility of LDs in precise ccRCC treatment. By comparing the IC50 of the drug in the ccRCC sample, individuals in the high-risk group showed higher sensitivity to chemotherapy drugs (such as 5-fluorouracil and cisplatin) and targeted drugs such as alpelisib, dabrafenib and entinostat (Figs. 9E-L). Therefore, the risk model possibly helped predict sensitivity to ICIs and targeted therapy in KIRC.

Fig. 9

Drug sensitivity analysis between the low- and high-risk groups. (A-D) Comparisons of the IPS in the two risk groups. The sensitivity analysis of 5-luorouracil (E), cisplatin (F), axitinib (G), alpelisib (H), cediranib (I), dabrafenib (J), entinostat (K) ibrutinib between the low- and high-risk groups

Discussion

The accumulation of LD in tissues other than adipose tissue has been recognized as a novel cancer hallmark [11]. LDs can interact with various organelles through membrane contact sites, regulating oxidative stress, preventing toxic substance accumulation in the endoplasmic reticulum, participating in metabolic regulation, and so on [32,33,34]. Previous reports have indicated that LD content is higher in cancer cells and tissues, such as prostate cancer, breast cancer and colorectal cancer, compared with normal cells and tissues [11]. Research has shown that the interaction between oncogenic and lipid metabolism pathways regulates cancer cells' LD balance [35]. Highly proliferative cancer cells increase the number of enzymes involved in lipid and cholesterol biosynthesis, and the abundance of lipid droplets storing excess lipids and cholesterol is associated with tumor aggressiveness [36]. A growing body of evidence highlighted the close correlation between tumor immunity and lipid metabolism in ccRCC [37]. For example, studies have indicated that fatty acid metabolism-related mRNA signatures are crucial in predicting survival prognosis and immunotherapy outcomes in ccRCC patients [38]. However, the exact role of lipid metabolism in ccRCC remains unclear. KIRC is an immunogenic tumor and its tumor immunity is closely associated with lipid metabolism [36, 39, 40]. Lipid droplets, important organelles in lipid metabolism and immune response, may contribute to ccRCC pathogenesis by regulating kidney tumor immune responses. However, there is still insufficient evidence regarding the functional role and clinical relevance of LD-related genes as a module signature in the context of precision medicine for KIRC. This study aimed to establish a LD-related signature as a promising biomarker and indicator for risk stratification, prognostic prediction, and evaluating responses to immunotherapy and targeted therapy in ccRCC patients.

Firstly, we identified 162 LD-related genes from medical literature and online databases, and then conducted GO and KEGG analyses based on 134 differentially expressed LD-related genes in TCGA database. The GO enrichment analysis revealed that these DEGs primarily participated in lipid storage, lipid metabolism, and the PPAR signaling pathway, which confirmed their association with lipids.

Risk score is commonly employed as a method to develop meaningful signatures. We constructed a risk signature with 4 LD-related genes based on the train data set. Furthermore, the risk model was verified in test dataset, entire dataset, ICGC cohort and E-MTAB-1980 cohort. The results showed that the low-risk group exhibited superior OS and PFS compared to the high-risk group, and the AUC values showed good prognostic performance. Both univariate and multivariate Cox regression analyses indicated that the T stage, clinical stages, age and risk score were predictive factors for individual survival outcomes in KIRC patients. The ROC curves also showed the accuracy of these predictive factors. Therefore, we developed a comprehensive nomogram by incorporating the signature with two clinical indicators (age and pathological stage) for accurate predictions. The results showed that the nomogram had good clinical applicability in estimating the survival rate of KIRC patients. Time-dependent ROC analysis demonstrated the favorable predictive performance of the risk features for the 1-year, 3-year, and 5-year survival of KIRC patients in the TCGA database. The DCA analysis suggested the obtained nomogram may better guide clinical decision-making.

In constructing the 4-gene risk signature, we observed high expression of AUP1 and SERPINF1 in KIRC samples, while DGAT2 and METTL7A were lowly expressed. After conducting a thorough literature review and analysis, we looked into the mechanisms of four genes present in tumor tissues. Our findings revealed that AUP1 plays a crucial role in regulating the synthesis of cholesterol esters and fatty acids (FAs) in KIRC cells. It does so by targeting and regulating sterol O-acyltransferase 1, which ultimately leads to the progression of KIRC [37]. The DGAT2 involved in lipid synthesis is modulated by RNASET2, consequently leading to an elevation in both triglyceride synthesis and the formation of lipid droplets within ccRCC cells [41]. Significant associations were found between the low expression of METTL7A and poor prognosis in KIRC [42]. SERPINF1 (PEDF) exhibits a negative correlation with anti-angiogenesis and a positive correlation with low tumor grade and pT stage [43]. These results indicate that these genes play important roles in lipid metabolism and tumor progression. However, more research is necessary to establish the precise mechanisms by which these genes influence clear cell renal cell carcinoma. Therefore, we further analyzed the molecular mechanisms between the two risk groups to speculate on potential mechanisms in KIRC. Based on KEGG and GO annotation, our findings imply that LD-related genes have the potential to impact the progression and immune response of KIRC via immune pathways.

Several studies have shown that lipids can act as intracellular signaling molecules that affect the functions of different immune cells, which has implications for the connection between lipid metabolism and tumor immunity. For instance, the reprogramming of lipid metabolism can indirectly influence immune cell functions, enhancing tumor immunotherapy [44,45,46,47]. ccRCC was recognized as an immunogenic tumor, with some patients potentially benefiting from ICIs, and the combination treatment of antiangiogenics and targeted immunotherapy has been recognized as a first-line treatment option [48, 49]. Based on these considerations and GSEA analysis results, we conducted an immunological analysis, uncovering patients in the high-risk group experienced decreased survival rates with increasing immune scores, consistent with prior research [50, 51]. CD8 + T cells, plasma cells, regulatory T cells and follicular helper T cells demonstrated elevated scores in the high-risk group, all of which are commonly associated with poor prognosis in ccRCC patients [52,53,54]. These findings further support the worse prognosis associated with the high-risk group, which also showed heightened scores for diverse immune functions, including T cell co-inhibition and co-stimulation, CCR and immunocheckpoints, with plasma cells and T cells serving essential roles in tumor immunity [55]. These conclusions suggested that the tumor microenvironment (TME) played a crucial role in ccRCC treatment and tumor progression [56]. Recent studies on Necroptosis-Pyroptosis genes and necrosis-associated miRNA profiles in ccRCC further demonstrated this point [57, 58].

Previous studies have indicated that TMB serves as a predictive biomarker for immune checkpoint inhibitor therapy [59, 60]. The elevated TMB found in the high-risk group in our study suggested a higher probability of recognition by the immune system, which was consistent with the elevated CD8 + T cell and plasma cell infiltration.. Additionally, the IPS-PD-1 + CTLA4 and IPS-CTLA4 scores were notably elevated in the high-risk group. The expression of CTLA4 was positively correlated with risk scores and had elevated in the high-risk group. As such, these findings suggested that high-risk patients might have a better potential to benefit from immunotherapy. Additionally, patients in the high-risk group were more sensitive to alpelisib, dabrafenib and entinostat and were more likely benefit from these target therapy. In summary, our risk model could help predict KIRC patient sensitivity to targeted therapies and immunotherapies.

However, this study has several limitations. Firstly, our experiments involved a relatively small sample size and lacked cross-referencing to diverse multi-omics datasets. Therefore, in the future, we need to use larger datasets for further data analysis and mining. Secondly, the exact mechanisms by which these four genes affect the prognosis of ccRCC patients remain unclear. Further research is needed to explore the correlations of these genes to validate the accuracy of ccRCC occurrence and prognosis prediction models, including functional analyses using animal models with knockout of candidate genes to determine whether the loss/inhibition of lipid droplet-associated factors alters the metabolic profile, proliferation, and metastatic capability of cells, as well as testing changes in sensitivity to commonly used ccRCC therapeutic drugs in cell lines with gene knockout or overexpression. Additionally, it is essential to recognize that this study is based on a retrospective design using public datasets, which introduces several potential biases. For instance, samples in public datasets may not fully represent the target population, leading to selection bias. Different analytical tools and algorithms can yield varying results from the same data, and updates to public datasets may include new data or error corrections, resulting in data processing bias. Lastly, our study lacks validation experiments. Given these limitations, we look forward to delving deeper and referencing different multi-omics databases in our future work, as well as conducting further validation experiments to analyze the specific mechanisms of the related genes in ccRCC to verify the reliability of our current findings.

Conclusion

In summary, we constructed a prognostic risk model using four LD-related genes. This model assisted in predicting the survival of ccRCC patients and demonstrated potential application value in forecasting responses to targeted therapy and immunotherapy. However, comprehensive validation experiments and additional multi-omics data are needed to confirm these results.