Construction and validation of a palmitoylation-related prognostic model for lung adenocarcinoma based on integrated bioinformatics and machine learning

Introduction

Lung adenocarcinoma (LUAD) is the predominant subtype of non-small cell lung cancer (NSCLC), originating from bronchial glands or alveolar epithelial cells. The subclassifications include minimally invasive adenocarcinomas (MIA) and invasive mucinous adenocarcinomas (IMA), among others (1). Lung cancer is the leading cause of cancer-related deaths worldwide, with LUAD accounting for over 40% of cases (2). It is the most common type of lung cancer among both male and female smokers and non-smokers (3). Molecular heterogeneity includes mutations in various driver genes, such as EGFR, ALK, ROS1, and KRAS (4). Although the application of targeted therapies and immune checkpoint inhibitors has significantly improved the prognosis for certain patients (5), the survival rate remains unsatisfactory, with an overall 5-year survival rate of approximately 15% (6). Tumor heterogeneity, drug resistance, and metastatic recurrence continue to pose significant challenges in clinical treatment.

Over the past decade, the prognostic assessment of LUAD has primarily relied on traditional clinical predictors such as tumor stage, histological grade, age, and gender. With the advancement of precision medicine, personalized treatment strategies have gradually become a research focus. In this context, genomic and transcriptomic analyses have provided critical insights into molecular subtyping and the identification of therapeutic targets for LUAD (6,7). However, treatment response and prognosis in LUAD are influenced by multiple factors, including tumor biological characteristics, patient immune status, and the tumor microenvironment (TME) (8). Recently, bioinformatics has highlighted the critical roles of metabolic reprogramming (7), immune evasion (8) and apoptosis (9) in LUAD by integrating data from the genome, epigenome, and proteome. Emerging biological processes such as cuproptosis (10), ribosome biogenesis (11), and ferroptosis (12), which are closely associated with the TME, have also been incorporated into prognostic signatures, offering diverse prognostic guidance for lung cancer treatment. Although the number of predictive model studies has increased in recent years, most models suffer from methodological limitations. These include the use of only a single algorithm for feature selection—leading to the loss of key features and model overfitting—as well as a lack of comparative analysis of algorithm combinations and independent external validation, resulting in insufficient generalizability and inconsistent model quality. Therefore, there is an urgent need to integrate novel prognostic biomarkers with clinical phenotypes to more accurately assess patient outcomes and guide personalized precision therapy (13).

Palmitoylation is a dynamic post-translational modification process that covalently attaches palmitic acid to specific amino acid residues of proteins through thioester bonds. This process can be classified into three categories: S-palmitoylation, N-palmitoylation, and O-palmitoylation (10). The enzymes responsible for catalyzing this reaction, known as palmitoyl transferases (PATs) belong to a family of zinc finger-like proteins that contain DHHC motifs (14). Palmitoylation plays a crucial role in various biological functions, including cell signaling, metabolic reprogramming, and immune regulation, primarily by influencing the membrane localization, stability, and interaction networks of proteins (15). A number of pro-oncogenic proteins and cancer suppressors are shown to undergo palmitoylation and are tightly linked to tumor development (16). On one hand, palmitoylation enhances the sustained activation of pro-survival signaling pathways, such as ERK and PI3K/AKT, by regulating the distribution and localization of pro-oncogenic proteins (e.g., RAS, EGFR) on the lipid membrane (17-19). On the other hand, the oncogenic SCRIB proteins promote the phosphorylation of YAP and TAZ, which consequently inhibit the activation of signaling transduction pathways, such as MAPK, thereby suppressing tumor growth (20). Other signaling transduction pathways, such as Wnt and PD-L1, are also regulated by protein palmitoylation, which affects tumor progression (21). Research has shown that members of the ZDHHC family regulate LUAD progression by targeting key proteins for palmitoylation. For example, ZDHHC21 mediates multisite palmitoylation of EGFR, promoting survival and migration of NSCLC cells; its inhibition enhances EGFR activity (19). ZDHHC9 promotes LUAD proliferation and immune evasion by stabilizing PD-L1 protein levels, while its deletion triggers PD-L1 degradation and activates anti-tumor immunity (22). Knocking down DHHC20, in turn, inhibits the growth of KRAS-mutant LUAD by blocking EGFR palmitoylation and reducing the activity of the PI3K signaling pathway, as well as the stability of the proto-oncoprotein MYC (23). These findings suggest that the palmitoylation network is not only involved in the pathogenesis of LUAD but may also serve as a crucial target for prognostic markers and drug development.

Research on the prognostic significance of LUAD palmitoylation-related genes (LPRGs) in the TME is still relatively limited. In this study, we systematically explored genes significantly associated with palmitoylation using bulk transcriptomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. A robust predictive model was constructed by integrating multiple machine learning algorithms, and a clinically applicable nomogram incorporating key clinical features was developed, demonstrating strong predictive performance. During model construction, core candidate genes were identified through diverse computational approaches to further elucidate potential molecular targets related to palmitoylation. Finally, comprehensive stratified analyses—including immune characteristics, mutational landscape, and drug sensitivity—were conducted by categorizing patients into high- and low-risk subgroups. As the first comprehensive study to reveal the multidimensional role of palmitoylation in LUAD, our work establishes a clinically relevant prognostic framework and offers novel therapeutic perspectives for targeted LUAD treatment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1389/rc).

Methods

Data sources

We obtained RNA sequencing data and clinical information for LUAD from TCGA database (https://portal.gdc.cancer.gov/). This dataset included raw counts and transcripts per million (TPM) values from 541 LUAD patients and 59 normal tissue samples, which were used for model construction. Corresponding clinical survival data were obtained from UCSC Xena (http://xena.ucsc.edu/). Patients with incomplete survival information or duplicate samples were excluded from further analysis. The stability and accuracy of the model were evaluated in validation cohorts comprising three datasets from GEO database (http://www.ncbi.nlm.nih.gov/geo/): GSE30219, GSE72094, and GSE31210. The detailed characteristics of these datasets are summarized in Table 1. Palmitoylation-related genes (PRGs) were retrieved from the GeneCards database (https://www.genecards.org/). Briefly, a gene set was initially acquired by querying the keyword ‘palmitoylation’ in GeneCards. The Relevance Score provided by GeneCards for each gene—a metric that quantifies the strength of evidence linking the gene to the biological process of interest—was used to refine the gene list. To ensure high confidence, we applied a data-driven filtering strategy based on the distribution of relevance scores, ultimately identifying 1,424 PRGs for subsequent analysis. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Identification of differentially expressed genes (DEGs) and the key modular genes by weighted gene co-expression network analysis (WGCNA)

Differential gene expressions between tumor and adjacent normal samples in the TCGA dataset were analyzed using the R package DESeq2 (version 1.46.0). DEGs were identified with a threshold of adjusted P value (padj) <0.05 and |log₂ fold change (FC)| >0.585. The ‘ggplot2’ and ‘pheatmap’ packages were employed to generate a volcano plot and a heatmap, respectively, for visualizing the DEGs. To identify key phenotype-related gene modules, WGCNA was performed using the core WGCNA R package (version 1.73). Initially, low-expression genes and outlier samples were filtered out. An appropriate soft-thresholding power (β) was selected to ensure a scale-free topology, with the criterion of a scale-free fit index >0.90 while maintaining a high mean connectivity. The adjacency matrix and topological overlap matrix (TOM) were subsequently calculated. Hierarchical clustering was then applied to group genes into co-expression modules, and a dendrogram with module assignments was constructed. Finally, module-trait associations were evaluated by correlating module eigengenes with clinical traits. Hub genes from modules significantly associated with key phenotypes were extracted as potential key regulators for further analysis.

Selection of feature hub genes and functional enrichment analysis

The intersection of DEGs, key module genes identified by WGCNA, and PRGs was determined. The resulting overlapping genes were defined as differentially expressed LPRGs and visualized using a Venn diagram. To explore the signaling pathways and biological functions associated with LPRGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the ‘clusterProfiler’ R package.

Machine learning for development and validation of prognostic LPRGs

After identifying the core LPRGs, we conducted a rigorous data cleaning process on the TCGA database once again. First, we excluded non-tumor samples. Then, we removed duplicate patient IDs from the data and deleted patients with unknown survival status or follow-up duration of less than one month. Lastly, we conducted a complete case analysis without employing any imputation procedures, resulting in a final dataset comprising patients with complete data for both predictor and outcome variables. Prognostic model construction was then performed using the ‘MIME’ package (version 0.0.0.9). The TCGA-LUAD cohort served as the training set, while the GSE30219, GSE72094, and GSE31210 datasets were utilized as external validation cohorts. An initial pool of 152 candidate features was subjected to univariate Cox regression analysis to identify LPRGs with significant prognostic value. Subsequently, the training set was used to conduct ten-fold cross-validation, evaluating a total of 117 algorithm combinations derived from 10 distinct machine learning methods. These included Supervised Principal Components (SuperPC), Random Survival Forest (RSF), Elastic Net (Enet), Stepwise Cox (StepCox), CoxBoost, Partial Least Squares Regression for Cox (plsRcox), Generalized Boosted Regression Modeling (GBM), Survival Support Vector Machine (Survival-SVM), Ridge regression, and least absolute shrinkage and selection operator (Lasso), alongside four potential variable selection filters (Lasso, StepCox, CoxBoost, and RSF). All predictive models were evaluated using Harrell’s C-index across both the training and external validation sets. Model performance was ranked based on the average C-index, leading to the identification of the optimal algorithm combination. The optimal cutoff value for the LPRGs signature was determined using the minimum P value method, stratifying patients into distinct risk subgroups. K-M survival analysis was then performed using the ’survminer’ and ’survival’ R packages, with overall survival (OS) differences between subgroups compared via the log-rank test. Additionally, time-dependent receiver operating characteristic (ROC) curves were generated using the ‘timeROC’ R package (version 0.4) to assess the predictive accuracy of the prognostic model at different time points. Simultaneously, we collected signatures of LUAD prediction models published over the past 2 years for comparison with our model, thereby clarifying the clinical benefit assessment of our model.

Construction and validation of predictive nomogram

Clinical data from the TCGA training cohort were retrieved and summarized in Table 2. A heatmap was used to visualize the distribution of risk scores and clinical characteristics—including T stage, N stage, sex, and age—across the high- and low-risk subgroups. To further evaluate the prognostic value of the model in a clinical context, univariate and multivariate Cox regression analyses were performed to verify whether the LPRG-based signature served as an independent prognostic factor, distinct from conventional clinical variables such as age, sex, and tumor stage. Subsequently, a clinical prediction nomogram was constructed based on the regression coefficients, integrating survival time, survival status, clinical stage, and risk score into a cumulative point system. This nomogram enables the estimation of individual probabilities of survival outcomes. Developed using the ‘regplot’ package, the nomogram offers an intuitive visualization of prognosis and provides quantifiable survival probability predictions for LUAD patients. To systematically evaluate the predictive performance of this tool, multiple validation approaches were employed: (I) calibration curves were plotted to assess the agreement between predicted probabilities and observed outcomes; (II) time-dependent area under the curve (AUC) analysis was conducted to evaluate the predictive accuracy of the model and individual variables at 1 to 5 years follow-up time points; (III) decision curve analysis (DCA) was performed to quantify the net clinical benefit across a range of threshold probabilities.

Hub prognostic biomarkers screening

To identify the core variables within the prognostic signature, the ‘MIME’ package was employed to re-analyze the previously identified significant prognostic features using an integrated approach involving eight machine learning algorithms. These included Lasso, Enet, Boruta, CoxBoost, RSF, eXtreme Gradient Boosting (Xgboost), StepCox, and Support Vector Machine-Recursive Feature Elimination (SVM-RFE). Candidate genes were ranked based on their importance, and those with the highest selection frequency were identified as the core PRGs most critically associated with LUAD prognosis. Subsequently, we examined the expression patterns of the five prognostic genes across risk subgroups and evaluated survival outcomes according to their high or low expression levels. Finally, gene set enrichment analysis (GSEA) was performed to investigate potential underlying pathways and mechanisms.

Analysis of tumor mutational burden (TMB) between LPRGs risk groups

We obtained TMB data from TCGA and analyzed the mutational landscape of LUAD patients across different risk subgroups using the maftools package. The results are presented in a waterfall plot displaying the top 20 most frequently mutated genes, while a violin plot is used to compare TMB scores between the two subgroups, serving as a metric for potential response to immunotherapy.

Immune characteristics and TME analysis

The ‘Immunedeconv’ R package (version 2.1.0) was employed to quantify the association between the palmitoylation modification signature and immune cell infiltration using six advanced algorithms—including QUANTISEQ, TIMER, CIBERSORT, XCELL, MCPCOUNTER, and EPIC—thereby enhancing the robustness of the findings. Furthermore, single-sample gene set enrichment analysis (ssGSEA) was applied to evaluate variations in immune cell characteristics, immune functional activity, and cancer hallmark pathways between LPRG-defined subgroups. The Tumor Immune Dysfunction and Exclusion (TIDE) framework (http://tide.dfci.harvard.edu) was utilized to assess potential responsiveness to immunotherapy and the likelihood of tumor immune escape. Additionally, the ESTIMATE algorithm was used to compute stromal, immune, ESTIMATE, and tumor purity scores for LUAD patients, enabling the validation of TME heterogeneity across distinct risk subgroups.

Prediction of drug sensitivity

To identify potential drugs for high-risk patients, we accessed a dataset of drug responses that included 809 tumor cell lines and 198 compounds from the Genomics Database for Cancer Drug Sensitivity (GDSC2, https://www.cancerrxgene.org/). Subsequently, we employed the ‘oncoPredict’ R software package (v1.2) to integrate tumor gene expression profiles from TCGA and systematically assess the differences in sensitivity to each type of drug between high- and low-risk groups, based on the drug’s half-maximal inhibitory concentration (IC₅₀) index.

Statistical analysis

All statistical analyses were conducted using R software (version 4.4.3). A P value of <0.05 was considered statistically significant for all tests.

Results

Identification of core modules and genes associated with LUAD

Figure 1 illustrates the workflow of this study. As detailed in the Methods section, we identified DEGs between LUAD tissues and normal lung tissues from the TCGA dataset. A total of 8,532 DEGs were identified, comprising 5,505 upregulated and 3,027 downregulated genes, which were visualized using volcano and heatmaps (Figure 2A,2B). Furthermore, to identify co-expression modules associated with LUAD, we employed WGCNA to construct a gene co-expression network. This network exhibited scale-free topology characteristics (R²=0.90), leading to the determination of an optimal soft threshold power of 7 (Figure 2C). Fifteen distinct modules were identified, and Pearson correlation coefficients were calculated between module eigengenes and sample traits. Module-trait relationship analysis revealed that the MEbrown module was significantly and negatively correlated with the LUAD phenotype (r=−0.76) (Figure 2D-2F). A total of 1,762 key genes were selected for subsequent analyses, providing crucial insights for elucidating the core molecular mechanisms of LUAD. Subsequently, we intersected the disease-related module eigengenes identified by WGCNA in LUAD, the DEGs from the differential expression analysis, and the palmitoylation-associated genes obtained via correlation scoring (Figure S1), resulting in 152 common DEGs (Figure 2G).

Figure 1 Study flowchart. DEGs, differentially expressed genes; FC, fold change; GEO, Gene Expression Omnibus; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden; WGCNA, weighted gene co-expression network analysis.

Figure 2 Identification of DEGs and key module genes in LUAD patients. (A) A gradient volcano plot was generated to illustrate the differential gene expression between tumor and adjacent normal tissues from the TCGA database. (B) Heatmap displaying the top 20 upregulated and downregulated DEGs and their clustering relationships with samples. Each row represents a differentially expressed gene, and each column corresponds to a sample from either the LUAD group or the normal group. (C) Analysis of the scale-free topology fit index and mean connectivity under different soft threshold powers. (D) A total of 15 co-expression modules were identified using the dynamic tree cut method, followed by the merging of similar modules. (E) Heatmap depicting the correlations between the identified lung cancer modules and phenotypic traits. (F) Scatter plot of gene significance against module membership reveals key module genes associated with LUAD. (G) Venn diagram shows the intersection among DEGs, palmitoylation-related genes, and genes obtained from WGCNA. DEGs, differentially expressed genes; LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas; WGCNA, weighted gene co-expression network analysis.

Functional enrichment of LPRGs

To further elucidate the potential biological mechanisms of LPRGs, we performed functional enrichment analyses, including GO, KEGG, and GSEA. GO analysis revealed that the DEGs were significantly enriched in various biological processes, including the regulation of phospholipid metabolic process, regulation of membrane potential, and glycerolipid metabolic process. For cellular components, significant terms included membrane microdomain, endocytic vesicle, and caveola. Molecular functions were primarily enriched in carboxylic ester hydrolase activity, lipase activity, and ion channel regulator activity (Figure 3A-3C). As shown in Figure 3D, KEGG pathway analysis indicated that these DEGs were predominantly involved in several key pathways, including glycerophospholipid metabolism, the AGE-RAGE signaling pathway in diabetic complications, and ECM-receptor interaction. These alterations in biological functions are associated with tumor invasion, metastasis, and remodeling of the TME, which are recognized pathological features of LUAD progression.

Figure 3 Pathway enrichment analysis. (A-C) Bar plot displaying the GO analysis results based on BP, CC, and MF. (D) Sankey bubble plot of KEGG pathway enrichment analysis. BP, biological process; CC, cellular component; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Development and validation of a prognostic model based on LPRGs

This study ultimately included 493 patients with LUAD, of whom 178 experienced an outcome event (death). A total of 152 differentially expressed LPRGs were included as candidate features in an integrated machine learning workflow to construct a robust and accurate prognostic model. The TCGA-LUAD cohort served as the training set, while the GSE30219, GSE72094, and GSE31210 datasets were used as external validation sets. This approach enhanced sample heterogeneity and reduced the risk of overfitting. Utilizing the modeling workflow from the mime1 package, univariate Cox regression initially identified 51 genes with prognostic significance as the final predicted variable. Subsequently, 117 algorithm combinations were employed to build candidate models, thereby mitigating potential bias inherent to any single algorithm. Based on the highest average C-index across the training and all validation sets, the ‘RSF + Ridge’ combination emerged as the optimal model (Figure 4A). Next, LUAD patients were stratified into high-risk (N=246) and low-risk groups (N=247) based on the median risk score derived from the model. Survival analysis across all cohorts indicated that the risk score had a significant impact on patient prognosis (P<0.05). The low-risk group consistently exhibited a significant survival advantage across all datasets, suggesting that a lower risk score is associated with a higher probability of survival and underscoring the favorable prognostic value of this signature (Figure 4B-4E, Figure S2). To evaluate the model’s robustness, time-dependent ROC analysis was performed on all four datasets. Using AUC >0.5 as the threshold for statistical significance, the AUC for 1-year survival prediction in the external validation sets was all over 0.7, except GES31210. The ROC curves demonstrated consistent predictive accuracy, aligning with the results from the training set (Figure 4F-4H). To further validate the relative advantage of our model, the ‘RSF + Ridge’ model was compared against 33 previously published LUAD prognostic models from the past 2 years. Hazard ratio (HR) analysis revealed that our model achieved a higher C-index across all datasets (Figure 5A, Table S1). Concurrently, its C-index and AUC also ranked among the highest compared to other models (Figure S3A,S3B). Integrating these findings, the ‘RSF + Ridge’ model demonstrated superior robustness and predictive performance and was therefore selected as the definitive prognostic signature model related to palmitoylation, comprising 51 feature genes.

Figure 4 Development of an LPRGs-based prognostic signature using machine learning. (A) Performance evaluation of 117 models through internal-external cross-validation, ranked by their mean C-index across validation sets. (B-E) Kaplan-Meier curves for overall survival of the LPRGs-based subgroups across the four datasets. (F-H) Time-dependent ROC curve analysis of the optimal model demonstrates consistent predictive accuracy at 1-, 3-, and 5-year overall survival. AUC, area under the curve; LPRG, LUAD palmitoylation-related gene; LUAD, lung adenocarcinoma; ROC, receiver operating characteristic; RSF, Random Survival Forest; TCGA, The Cancer Genome Atlas.

Figure 5 Development and validation of the prognostic nomogram. (A) Our model outperforms 33 previously published prognostic models for LUAD over the past two years. (B) Association analysis between the LPRGs-based risk subgroups and clinical characteristics. (C) A predictive nomogram constructed based on LPRGs and the clinical stage. The “Total points” axis corresponds to the probability of 1-, 3-, and 5-year overall survival, reflecting the prognostic outcome of LUAD patients under different total scores at these time points. (D) Calibration curves and (E) time-dependent ROC curves of the multi-parameter model demonstrate its predictive accuracy. (F) Decision curve analysis using net benefit assessment illustrates the clinical value of the model. ***, P<0.001. AUC, area under the curve; HR, hazard ratio; LPRG, LUAD palmitoylation-related gene; LUAD, lung adenocarcinoma; OS, overall survival; ROC, receiver operating characteristic.

Incorporation of clinical features to construct an LPRGs-based nomogram

The distribution of the risk score alongside various clinical features (gender, stage, and T, N classification) is displayed in a heatmap (Figure 5B). To further validate the independence of our model, we incorporated the risk score and patient clinical indicators (including gender, age, T stage, N stage, M stage, and overall stage) into both univariate and multivariate Cox regression analyses (Figure S4). The results demonstrated that the stage and risk score served as the most reliable and independent predictors of prognosis in LUAD patients [HR = 2.214, 95% confidence interval (CI): 1.781–2.752, P<0.001], confirming LPRGs as an independent prognostic factor across all datasets (Figure S3A,S3B). A nomogram was subsequently constructed based on the overall stage and the LPRGs signature (Figure 5C), enhancing the clinical applicability of the model. The calibration curves indicated a high degree of consistency between the nomogram-predicted and the observed 1-, 3-, and 5-year OS rates. And the calibration slope for the first year reached 0.9722 (P=0.0065), with an intercept of 0.1252 (P=0.834), confirming the model’s excellent calibration and clinical utility in LUAD prognosis, which validated its clinical utility for LUAD prognosis (Figure 5D and Table S2). Time-dependent AUC analysis showed that the nomogram, which incorporated clinical factors, exhibited a superior predictive ability with AUC values consistently above 0.7. Furthermore, the AUC values demonstrated increasing stability over time (Figure 5E). Additionally, DCA confirmed that the predictive performance of this nomogram was better than that of any single clinical factor (Figure 5F).

Screening of core prognostic genes

Subsequently, we applied eight individual machine learning algorithms to identify the top 19 most frequently selected genes, as well as the intersecting genes derived from all eight algorithms, to pinpoint key genes most critically associated with LUAD prognosis (Figure 6A). These genes were distributed across various chromosomes, excluding chromosomes 2, 10, 14, 18–21, and the sex chromosomes (Figure 6B). Among the top five candidate genes, TXN is recognized as a potential therapeutic target in lung cancer. It can enhance proliferation and antioxidant capacity, activate the ERK pathway, and promote immune escape, thereby being closely linked to tumor progression, poor prognosis, and resistance to immunotherapy (24,25). GPD1L and ATP8A2 have been widely implicated as playing tumor-suppressive roles in various malignancies, including LUAD (13,26,27). DNAJB4, a member of the heat shock protein family, may specifically influence lung cancer cell proliferation and tumorigenesis when downregulated (28). In contrast, the role of SCN2B in LUAD remains insufficiently studied. Kaplan-Meier survival analysis revealed that high expression of TXN and DNAJB4 was significantly associated with shorter OS in LUAD patients, while elevated expression of the other three genes also indicated a favorable prognosis (Figure S5). To investigate the functional roles of these core genes in LUAD, we examined their differential expression across various subgroups and performed GSEA. The analysis suggested that DNAJB4 might activate pathways such as endocytosis, focal adhesion, IgSF CAM signaling, and platelet activation. Conversely, ATP8A2 potentially suppresses pathways including the cell cycle and oxidative phosphorylation, influencing the TME and contributing to a hypoxic state (Figure 6C-6E and Figures S6,S7).

Figure 6 Exploration of differences in core genes. (A) Top 19 genes with the highest importance identified through 8 machine learning algorithms. (B) Chromosomal distribution of LPRGs visualized in a circular plot. (C) Relative expression levels of the five key genes across LPRG-based subgroups. Top 5 pathways enriched in GSEA for (D) ATP8A2 and (E) DNAJB4. GSEA, gene set enrichment analysis; LPRG, LUAD palmitoylation-related gene; LUAD, lung adenocarcinoma; RSF, Random Survival Forest; SVM-RFE, Support Vector Machine-Recursive Feature Elimination.

TMB mutation in LPRGs subgroups

To further investigate the genetic heterogeneity between the high- and low-risk subgroups, we independently analyzed the genomic mutation profiles of both groups using the TCGA-LUAD dataset (Figure 7A,7B). Significant differences in mutation spectra were observed. Among the 244 high-risk samples, 236 (96.72%) carried gene mutations, with a significantly higher mutation frequency compared to the low-risk group. TP53—a key tumor suppressor gene that coordinately regulates cell cycle, DNA repair, and apoptosis (29)—and TTN—a critical protein involved in modulating myofibrillar elasticity, were the most frequently mutated genes common to both subgroups. However, their distribution showed a marked disparity: in the high-risk group, TP53 mutations accounted for 61% and TTN for 57%, whereas in the low-risk group, their mutation rates were only 41% (30).

Figure 7 Mutational landscape of the LPRG subgroups. (A,B) The waterfall plot demonstrates that patients in the high-risk subgroup (B) exhibit a higher mutation frequency, whereas the low-risk subgroup (A) shows the opposite pattern. (C) Analysis of mutual exclusivity and co-occurrence among the top 20 mutated genes in the low-risk subgroup. (D) Differences in TMB scores between the high-risk and low-risk subgroups. ***, P<0.001. LPRG, LUAD palmitoylation-related gene; LUAD, lung adenocarcinoma; TMB, tumor mutation burden.

Furthermore, analysis of co-occurrence and mutual exclusivity among the top 20 mutated genes revealed that low-risk patients exhibited a greater tendency towards mutually exclusive mutations, while high-risk patients showed a significantly higher frequency of co-occurring mutations (Figure 7C). TMB analysis indicated a significantly higher mutation load in the high-risk subgroup compared to the low-risk subgroup (Figure 7D), underscoring the distinct genomic characteristics between the LPRGs-defined subgroups.

Characteristics of the immune microenvironment and predictive potential for immunotherapy response

TME is composed of various cell types, including tumor cells, stromal cells, and infiltrating immune cells. Given the close association between LUAD pathophysiology and the immune microenvironment, we further investigated the distribution of immune cells in LUAD. Using six distinct algorithms to quantify immune infiltration levels across subgroups, the results indicated that the infiltration of most immune cell types was negatively correlated with the risk score (Figure 8A, Figure S8). Notably, ssGSEA analysis revealed higher infiltration levels of activated CD4⁺ and CD8⁺ T cells in the high-risk subgroup, whereas the low-risk subgroup exhibited significantly higher infiltration of B cells, dendritic cells, and several other immune cells (Figure 8B). Furthermore, we assessed 13 immune functional signatures across the subgroups. The high-risk subgroup demonstrated significantly higher scores for gene sets related to APC co-inhibition, cytolytic activity, inflammation-promoting, and T cell co-stimulation (Figure 8C), suggesting that high-risk tumors are more likely to promote tumor progression and immune escape by enhancing immunosuppressive and inflammatory responses. Concurrently, LUAD patients with high LPRG scores showed elevated TIDE scores, indicating potentially reduced responsiveness to conventional immunotherapy and a greater propensity for immune escape (Figure 8D). Additionally, the expression levels of certain immune-activating or co-stimulatory molecules (such as TNFSF18 and HLA family genes) were significantly higher in the low-risk subgroup, suggesting that these patients might derive greater benefit from immune checkpoint blockade (ICB) therapy. These collective differences indicate that the high-risk subgroup is characterized by an immunosuppressive TME, whereas the low-risk subgroup exhibits a more active anti-tumor immune state (Figure 8E).

Figure 8 Characteristics of the tumor microenvironment in LPRGs-based subgroups. (A) Correlation analysis between immune cell proportions, assessed by six algorithms, and the risk score. (B-E) Analysis of immune cell levels (B), T-cell exclusion and dysfunction (C), TIDE scores (D), and immune checkpoint gene expression differences (E) across different subgroups. *, P<0.05; **, P<0.01; ***, P<0.001. LPRG, LUAD palmitoylation-related gene; LUAD, lung adenocarcinoma; TIDE, Tumor Immune Dysfunction and Exclusion.

Drug sensitivity analysis

We further investigated drug sensitivity across different risk subgroups. Among the 198 chemotherapeutic agents evaluated, 70.2% (139 compounds) exhibited significant differences in IC₅₀ values between the high- and low-risk subgroups (P<0.05), as indicated by drug sensitivity differential analysis. Figure 9 highlights the top 10 agents with the lowest P values. Compared to patients in the low-risk subgroup, those in the high-risk subgroup showed significantly lower predicted IC₅₀ values for therapeutic agents such as luminespib, MK-1775, SCH772984, and selumetinib. This suggests that patients in the high-risk subgroup may exhibit stronger treatment responses to these drugs.

Figure 9 Differences in IC₅₀ values for commonly used drugs between high and low-risk groups. ***, P<0.001. IC₅₀, half-maximal inhibitory concentration.

Discussion

Despite advancements in cancer treatment, LUAD presents significant challenges due to its clinically nonspecific early symptoms, late-stage diagnosis, and highly invasive and metastatic nature. Recently, prognostic models for patients with various cancers, including LUAD, have been developed using bioinformatics techniques (31,32). It has been reported that palmitoylation plays a crucial role in tumor development and progression, and genes associated with palmitoylation can serve as prognostic markers (33-35). However, to our knowledge, no study has yet integrated LUAD transcriptomic data with palmitoylation levels to construct a survival prediction model. Therefore, investigating palmitoylation to identify novel therapeutic targets is crucial for improving the diagnosis and treatment of LUAD.

In our study, we established an LPRGs-based scoring system for prognostic assessment in LUAD patients by focusing on DEGs that were significantly altered compared to normal lung tissues. Using WGCNA, we identified relevant module genes and intersected them with a PRG set, thereby defining 152 prognosis-associated DEGs. Subsequently, univariate Cox regression was applied to obtain 51 prognostic markers, which were then used to construct candidate models through 117 combinations derived from 10 machine learning algorithms. Based on the optimal model, patients were ultimately stratified into high-risk and low-risk subgroups. Analyses revealed that the prognostic model based on PRG expression consistently indicated poorer outcomes in the high-risk subgroup. The model was subsequently validated using GEO datasets as external cohorts. Following this validation, a risk prediction nomogram was developed by integrating the risk score and clinical stage to provide a more accurate assessment of mortality risk and to inform better treatment selection. Finally, we explored multiple dimensions of differences between the two subgroups, including genomic mutations, immune characteristics (such as immune infiltration levels and immune checkpoint expression), and drug sensitivity. These results integrate potential prognostic features of LUAD with the TME, offering a personalized prognostic assessment and therapeutic guidance for LUAD patients.

Exploration of the top five contributing genes (core genes) within LPRGs and their mechanisms involving palmitoylation identified TXN and DNAJB4 as risk factors, while ATP8A2, GPD1L, and SCN2B were recognized as protective factors for LUAD prognosis. TXN, a member of the thioredoxin family, is a key antioxidant and ferroptosis regulator (36). It maintains the cellular reductive environment through its antioxidant activity, which can indirectly regulate the palmitoylation mediated by DHHC enzymes and their substrates. This subsequently alters the membrane localization, stability, and activity of critical oncogenic signaling proteins such as RAS, EGFR, Wnt, and PD-L1, ultimately influencing tumor cell proliferation, metastasis, and immune evasion. Liu et al. demonstrated that TXN regulates c-Myc expression via the ERK1/2 and ERK5 signaling pathways, thereby promoting NSCLC progression. Both in vitro and in vivo knockdown of TXN significantly suppressed cell proliferation and invasion while enhancing apoptosis, whereas TXN overexpression reversed these malignant phenotypes (24). DNAJB4, a DNAJ/Hsp40 chaperone, has emerged as a novel prognostic biomarker in lung cancer. Chen et al. experimentally showed that loss of DNAJB4 can inhibit epithelial-mesenchymal transition (EMT) and reduce lung cancer metastasis. It also downregulates the formation of oncogenic complexes associated with EGFR, FAK, and STAT3 signaling pathways (37). Our model suggests ATP8A2 as a protective factor. Belonging to the P4-ATPase family, ATP8A2 actively translocates phosphatidylserine and phosphatidylethanolamine from the outer to the inner leaflet of the plasma membrane, generating and maintaining phospholipid asymmetry. Aberrant methylation of ATP8A2 has been identified in various cancers, including LUAD (38), a finding corroborated by Zhang’s research (39). GPD1L, a protein involved in metabolic regulation, has been implicated in the progression of several cancers such as lung, colorectal, and laryngeal cancer. We observed that decreased GPD1L expression in the high-risk subgroup correlated with poor prognosis, consistent with previous studies reporting its specific downregulation in tumor samples and its tumor-suppressive function (40). Fan et al. confirmed that, in addition to inhibiting proliferation, migration, and invasion, GPD1L also alleviates apoptosis and mitochondrial damage in vitro (13). Studies in renal cancer models found that knocking down GPD1L may remodel the cellular membrane phospholipid composition, thereby affecting the localization and activity of palmitoylating enzymes and their substrates on the membrane, consequently impeding the palmitoylation process and membrane targeting of its substrates (41). SCN2B encodes the sodium channel β2 subunit. Current evidence regarding its role in lung cancer progression is limited, marking it as a potential yet understudied target. It may promote migration and perineural invasion in certain cancers (e.g., prostate cancer) through mechanisms involving cell adhesion and non-conducting sodium channel signaling (42).

It is well known that tumor immune infiltration is markedly correlated with cancer prognosis. Consequently, the TME has progressively become a research focus. In our study, we observed activated functions related to APC co-inhibition and inflammation-promoting activity, along with higher TIDE scores in the high-risk subgroup. This suggests that high-risk individuals have a greater likelihood of developing immune escape mechanisms, potentially reducing the benefit from immunotherapy. Additionally, our pathway enrichment analysis revealed that PRGs were significantly enriched in cell cycle regulation pathways. Recent studies have shown that palmitoylation dynamically modifies cell cycle regulatory proteins and their associated pathways, influencing tumor cell proliferation, cell cycle checkpoints, and genomic stability. For instance, ZDHHC1-mediated palmitoylation of the tumor suppressor p53 enhances the binding stability of the p53-DNA complex and effectively suppresses its ubiquitination-dependent degradation pathway. This reinforces the activation of p53-dependent cell cycle checkpoints. This regulatory network ultimately induces tumor cell apoptosis and significantly inhibits malignant proliferation (43). Furthermore, ZDHHC13-catalyzed palmitoylation of the melanocortin 1 receptor (MC1R) has been identified as a core molecular switch for activating MC1R signaling cascades. This receptor mediates UVB radiation-induced G1 phase cell cycle arrest through a dual regulatory mechanism. This palmitoylation-dependent signaling activation system plays a significant role in regulating cellular senescence and inhibiting melanoma progression (44).

While the integration of multi-dimensional data and machine learning strengthens our findings, this study has limitations. Its retrospective nature necessitates validation in prospective clinical cohorts, and the functional roles of the identified biomarkers require further experimental investigation. Nonetheless, this work offers insightful clues and a framework for developing palmitoylation-based prognostic strategies in LUAD.

Conclusions

In summary, our study successfully employed bioinformatics and machine learning approaches to identify 51 prognosis-associated genes and develop a novel prognostic model based on PRGs. This model enables accurate prediction of patient outcomes and evaluation of immune responses in LUAD. The observed correlations between PRGs and both tumor mutation burden and immune cell infiltration suggest the crucial role of palmitoylation in modulating the TME. These findings provide new perspectives for prognostic assessment and open new avenues for developing personalized treatment strategies in LUAD.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1389/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1389/prf

Funding: The study was supported by the Sanming Project of Medicine in Shenzhen (grant No. SZZYSM202311001).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1389/coif). All authors report that the study was supported by the Sanming Project of Medicine in Shenzhen (grant No. SZZYSM202311001). The authors have no other conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Nicholson AG, Tsao MS, Beasley MB, et al. The 2021 WHO Classification of Lung Tumors: Impact of Advances Since 2015. J Thorac Oncol 2022;17:362-87. [Crossref] [PubMed]
2. Denisenko TV, Budkevich IN, Zhivotovsky B. Cell death-based treatment of lung adenocarcinoma. Cell Death Dis 2018;9:117. [Crossref] [PubMed]
3. Zappa C, Mousa SA. Non-small cell lung cancer: current treatment and future advances. Transl Lung Cancer Res 2016;5:288-300. [Crossref] [PubMed]
4. Jaiyesimi IA, Leighl NB, Ismaila N, et al. Therapy for Stage IV Non-Small Cell Lung Cancer With Driver Alterations: ASCO Living Guideline, Version 2023.3. J Clin Oncol 2024;42:e1-e22. [Crossref] [PubMed]
5. Thai AA, Solomon BJ, Sequist LV, et al. Lung cancer. Lancet 2021;398:535-54. [Crossref] [PubMed]
6. Spella M, Stathopoulos GT. Immune Resistance in Lung Adenocarcinoma. Cancers (Basel) 2021;13:384. [Crossref] [PubMed]
7. Zhu HE, Li T, Shi S, et al. ESCO2 promotes lung adenocarcinoma progression by regulating hnRNPA1 acetylation. J Exp Clin Cancer Res 2021;40:64. [Crossref] [PubMed]
8. Li Y, Zeng M, Qin Y, et al. The role of KRT18 in lung adenocarcinoma development: integrative bioinformatics and experimental validation. Discov Oncol 2024;15:841. [Crossref] [PubMed]
9. Liu Y, Ouyang Y, Feng Z, et al. RASGRP2 is a potential immune-related biomarker and regulates mitochondrial-dependent apoptosis in lung adenocarcinoma. Front Immunol 2023;14:1100231. [Crossref] [PubMed]
10. Zhang W, Qu H, Ma X, et al. Identification of cuproptosis and immune-related gene prognostic signature in lung adenocarcinoma. Front Immunol 2023;14:1179742. [Crossref] [PubMed]
11. Song Z, Wang Y, Zhu M, et al. Exploring ribosome biogenesis in lung adenocarcinoma to advance prognostic methods and immunotherapy strategies. J Transl Med 2025;23:503. [Crossref] [PubMed]
12. Huang Z, Chen X, Wang Y, et al. SLC7A11 inhibits ferroptosis and downregulates PD-L1 levels in lung adenocarcinoma. Front Immunol 2024;15:1372215. [Crossref] [PubMed]
13. Fan Z, Wu S, Sang H, et al. Identification of GPD1L as a Potential Prognosis Biomarker and Associated with Immune Infiltrates in Lung Adenocarcinoma. Mediators Inflamm 2023;2023:9162249. [Crossref] [PubMed]
14. Lemonidis K, Werno MW, Greaves J, et al. The zDHHC family of S-acyltransferases. Biochem Soc Trans 2015;43:217-21. [Crossref] [PubMed]
15. Resh MD. Palmitoylation of proteins in cancer. Biochem Soc Trans 2017;45:409-16. [Crossref] [PubMed]
16. Ko PJ, Dixon SJ. Protein palmitoylation and cancer. EMBO Rep 2018;19:e46666. [Crossref] [PubMed]
17. Bollu LR, Katreddy RR, Blessing AM, et al. Intracellular activation of EGFR by fatty acid synthase dependent palmitoylation. Oncotarget 2015;6:34992-5003. [Crossref] [PubMed]
18. Busquets-Hernández C, Triola G. Palmitoylation as a Key Regulator of Ras Localization and Function. Front Mol Biosci 2021;8:659861. [Crossref] [PubMed]
19. Runkle KB, Kharbanda A, Stypulkowski E, et al. Inhibition of DHHC20-Mediated EGFR Palmitoylation Creates a Dependence on EGFR Signaling. Mol Cell 2016;62:385-96. [Crossref] [PubMed]
20. Chen B, Zheng B, DeRan M, et al. ZDHHC7-mediated S-palmitoylation of Scribble regulates cell polarity. Nat Chem Biol 2016;12:686-93. [Crossref] [PubMed]
21. Liu Z, Xiao M, Mo Y, et al. Emerging roles of protein palmitoylation and its modifying enzymes in cancer cell signal transduction and cancer therapy. Int J Biol Sci 2022;18:3447-57. [Crossref] [PubMed]
22. Li Z, Jiang D, Liu F, et al. Involvement of ZDHHC9 in lung adenocarcinoma: regulation of PD-L1 stability via palmitoylation. In Vitro Cell Dev Biol Anim 2023;59:193-203. [Crossref] [PubMed]
23. Kharbanda A, Walter DM, Gudiel AA, et al. Blocking EGFR palmitoylation suppresses PI3K signaling and mutant KRAS lung tumorigenesis. Sci Signal 2020;13:eaax2364. [Crossref] [PubMed]
24. Liu X, Dong X, Hu Y, et al. TXN promotes tumorigenesis by activating the ERK1/2 and ERK5 signaling pathways regulating c-Myc in non-small cell lung cancer. Cell Signal 2025;125:111517. [Crossref] [PubMed]
25. Csiki I, Yanagisawa K, Haruki N, et al. Thioredoxin-1 modulates transcription of cyclooxygenase-2 via hypoxia-inducible factor-1alpha in non-small cell lung cancer. Cancer Res 2006;66:143-50. [Crossref] [PubMed]
26. Yang M, Pang D, Gong C, et al. GPD1L as a potential biomarker associated with Treg cell infiltration and lipid metabolism in clear cell renal cell carcinoma. J Transl Med 2025;23:872. [Crossref] [PubMed]
27. Ding L, Zhang H. Circ-ATP8A2 promotes cell proliferation and invasion as a ceRNA to target EGFR by sponging miR-433 in cervical cancer. Gene 2019;705:103-8. [Crossref] [PubMed]
28. Kim HY, Hong S. Multi-Faceted Roles of DNAJB Protein in Cancer Metastasis and Clinical Implications. Int J Mol Sci 2022;23:14970. [Crossref] [PubMed]
29. Wang Z, Strasser A, Kelly GL. Should mutant TP53 be targeted for cancer therapy? Cell Death Differ 2022;29:911-20. [Crossref] [PubMed]
30. Fullenkamp DE. Regulation of TTN as a mechanism of and treatment for heart failure. J Clin Invest 2025;135:e189335. [Crossref] [PubMed]
31. Piper Vallillo AJ, Viray H, Feldman J, et al. Management of Treatment Resistance in Patients With Advanced Epidermal Growth Factor Receptor-Mutated Lung Cancer: Personalization, Parsimony, and Partnership. J Clin Oncol 2024;42:1215-21. [Crossref] [PubMed]
32. Chmielecki J, Mok T, Wu YL, et al. Analysis of acquired resistance mechanisms to osimertinib in patients with EGFR-mutated advanced non-small cell lung cancer from the AURA3 trial. Nat Commun 2023;14:1071. [Crossref] [PubMed]
33. Qin G, Pang G, Wu S, et al. Construction and validation of a novel prognostic model with palmitoylation-related genes for glioblastoma. Transl Cancer Res 2024;13:6117-35. [Crossref] [PubMed]
34. Zeng S, Zeng L, Xie X, et al. Palmitoylation-related gene expression and its prognostic value in ovarian cancer: insights into immune infiltration and therapeutic potential. Discov Oncol 2024;15:802. [Crossref] [PubMed]
35. Zhu H, Hu H, Hao B, et al. Insights into a Machine Learning-Based Palmitoylation-Related Gene Model for Predicting the Prognosis and Treatment Response of Breast Cancer Patients. Technol Cancer Res Treat 2024;23:15330338241263434. [Crossref] [PubMed]
36. Wang J, Zhou J, Wang C, et al. Thioredoxin-1: A Promising Target for the Treatment of Allergic Diseases. Front Immunol 2022;13:883116. [Crossref] [PubMed]
37. Chen CH, Chang WH, Su KY, et al. HLJ1 is an endogenous Src inhibitor suppressing cancer progression through dual mechanisms. Oncogene 2016;35:5674-85. [Crossref] [PubMed]
38. Lv L, Huang Y, Li Q, et al. A Comprehensive Prognostic Model for Colon Adenocarcinoma Depending on Nuclear-Mitochondrial-Related Genes. Technol Cancer Res Treat 2024;23:15330338241258570. [Crossref] [PubMed]
39. Zhang X, Dong W, Zhang J, et al. A Novel Mitochondrial-Related Nuclear Gene Signature Predicts Overall Survival of Lung Adenocarcinoma Patients. Front Cell Dev Biol 2021;9:740487. [Crossref] [PubMed]
40. Liu T, Zhu H, Ge M, et al. GPD1L inhibits renal cell carcinoma progression by regulating PINK1/Parkin-mediated mitophagy. J Cell Mol Med 2023;27:2328-39. [Crossref] [PubMed]
41. Yao CH, Park JS, Kurmi K, et al. Uncoupled glycerol-3-phosphate shuttle in kidney cancer reveals that cytosolic GPD is essential to support lipid synthesis. Mol Cell 2023;83:1340-1349.e7. [Crossref] [PubMed]
42. Jansson KH, Castillo DG, Morris JW, et al. Identification of beta-2 as a key cell adhesion molecule in PCa cell neurotropic behavior: a novel ex vivo and biophysical approach. PLoS One 2014;9:e98408. [Crossref] [PubMed]
43. Tang J, Peng W, Feng Y, et al. Cancer cells escape p53's tumor suppression through ablation of ZDHHC1-mediated p53 palmitoylation. Oncogene 2021;40:5416-26. [Crossref] [PubMed]
44. Chen S, Zhu B, Yin C, et al. Palmitoylation-dependent activation of MC1R prevents melanomagenesis. Nature 2017;549:399-403. [Crossref] [PubMed]