Integration of single-cell and bulk transcriptomics to explore prognostic and immunotherapeutic characteristics of sodium overload-induced cell death in lung adenocarcinoma

Introduction

Background

Lung cancer remains the leading cause of cancer-related mortality, accounting for approximately 20% of all cancer deaths (1). Despite advances in treatment, the 5-year survival rate for lung cancer patients has shown only marginal improvement, primarily due to late-stage diagnoses that result in poor prognosis (2). Non-small cell lung cancer (NSCLC) constitutes approximately 85% of all lung cancer cases (3), with lung adenocarcinoma (LUAD) being the most prevalent histological subtype, representing approximately 40% of lung tumors (4). Although significant progress has been made in therapeutic strategies for LUAD, the 5-year survival rate remains low, ranging from 4% to 17% depending on disease stage and regional variations (5). Consequently, a comprehensive investigation into the molecular mechanisms underlying LUAD development and progression, particularly the role of sodium death-associated genes (NECSOs), is critical. Such research holds substantial promise for identifying novel therapeutic targets and prognostic biomarkers, offering new avenues for the diagnosis and treatment of LUAD.

Rationale and knowledge gap

Necrocide 1 (NC1) targets the transient receptor potential cation channel subfamily M member 4 (TRPM4), a non-selective monovalent cation channel that facilitates sodium ion influx, ultimately inducing cell necrosis. NC1 is closely linked to energy depletion and sodium death. Human TRPM4 gain-of-function mutations associated with arrhythmias exhibit sensitivity to NC1 or 2-deoxy-D-glucose-induced NECSO activation (6). Previous studies have demonstrated that pharmacological inhibition of TRPM4 mitigates cardiac ischemia-reperfusion injury and reduces arrhythmia incidence. TRPM4 is widely expressed in cardiac tissue and contributes to various physiological and pathological processes. Notably, TRPM4 exhibits voltage-dependent properties, with its activity markedly regulated by membrane potential, profoundly influencing cardiomyocyte electrical activity through membrane depolarization, likely mediated by sodium influx. It plays an essential role in maintaining the structural integrity and electrophysiological function of the heart (7-9). Under inflammatory conditions, TRPM4 promotes vascular endothelial cell migration and reactive oxygen species production (10). In oncology, TRPM4 represents a promising therapeutic target (11). For instance, TRPM4 inhibitors suppress colorectal cancer cell proliferation and disrupt cell cycle progression (12). In prostate cancer, TRPM4 enhances cell survival, migration, cell cycle transition, and adhesion (13). Additionally, TRPM ion channels, including TRPM4, are significantly implicated in bladder cancer development, positioning them as potential novel therapeutic targets (14).

Objective

This study integrates single-cell and bulk transcriptome data to investigate NECSOs and their interactions with the tumor microenvironment (TME), aiming to develop a prognostic model that identifies key targets linked to sodium death. This approach seeks to provide a novel theoretical framework for the prognosis and precision treatment of LUAD. Through single-cell TME analysis, we characterize immune features and examine dynamic changes in NECSOs using pseudotime analysis. We explore the mechanisms underlying differential prognoses in patients with varying NECSO risk profiles through multiple approaches, including enrichment scoring, consensus clustering, prognostic modeling, and immune infiltration analysis. Additionally, we construct a nomogram to accurately predict survival outcomes in LUAD patients. The findings of this study are expected to yield new prognostic biomarkers and therapeutic targets for early-stage LUAD, enhancing clinical prognosis prediction and treatment strategies. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-919/rc).

Methods

Data collection and processing

Single-cell RNA sequencing data were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSM5699782 and GSM5699783. Bulk transcriptome and copy number variation (CNV) data were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/), Encompassing 501 LUAD clinical samples and 59 normal clinical samples. For validation, two GEO datasets, GSE31210 (n=226) and GSE50081 (n=127), were used. Gene set enrichment analysis (GSEA) utilized the h.all.v2023.1.Hs.symbols.gmt and c2.cp.kegg.symbols.gmt datasets, sourced from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb). This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

TRPM4 co-expression

Using a correlation filter of |r| >0.5 and a P value threshold of P=0.001, we identified 40 genes co-expressed with TRPM4.

Single-cell analysis

Single-cell data processing was performed using the “Seurat” R package. Initially, quality control was conducted on the single-cell data, followed by normalization. We used “FindVariableFeature” to identify 2,000 highly variable genes and applied scaleData to ensure equal weighting of each gene. Dimensionality reduction was performed using “RunPCA”. The “DoubletFinder” R package was employed to remove doublets, while batch effects between samples were eliminated using the “Harmony” R package. We utilized “FindClusters” and “FindNeighbors” to construct cell clusters and identify marker genes for cluster annotation. The results were visualized using t-distributed stochastic neighbor embedding (t-SNE). Additionally, we calculated sodium death-associated scores for each cell using AddModuleScore and divided cells into high and low sodium death-associated expression groups based on the median score for differential analysis.

Weighted gene co-expression network analysis (WGCNA) analysis

WGCNA was employed to identify co-expression gene modules, examine relationships between modules and traits or phenotypes, and select network hub genes and highly correlated gene modules. The optimal soft threshold was determined using the pickSoftThreshold function and set to 3. Co-expression networks were constructed using a one-step method, with a minimum gene number of 20 per module. Module-trait relationships were analyzed by calculating Pearson correlation coefficients between each module and the tumor group. Subsequently, the module with the highest correlation to the tumor group was selected as the key module for further analysis.

Differential and prognostic analysis of NECSOs

Differential expression analysis was performed on genes within modules selected by WGCNA. Differentially expressed genes (DEGs) were identified using criteria of |log₂ fold change (FC)| >1 and false discovery rate (FDR) <0.05. Subsequently, univariate Cox analysis was conducted on these differential genes with a significance threshold of P<0.001.

Epithelial cell trajectory analysis

CytoTRACE (Cellular Trajectory Reconstruction Analysis using Gene Counts and Expression), a computational method for predicting cell differentiation states from single-cell RNA sequencing data, was used to determine the differentiation starting point of epithelial cells. Monocle, which employs principal component analysis and minimum spanning tree techniques, was utilized to accurately reconstruct the dynamic pathways of cells during development. We used the “Monocle” R package to simulate the dynamic development process and differentiation trajectory of epithelial cells, and to plot the changing trends of risk and protective genes along pseudotime.

Prognostic network and copy number analysis

Among the 14 prognostic genes we previously established, we constructed a prognostic network to investigate the relationship between NECSOs and tumors. Nine genes were identified as risk factors significantly associated with poor prognosis, while five favorable factors were significantly associated with good prognosis. CNVs of all 14 genes were calculated.

Consensus clustering analysis

Using the 14 previously identified prognostic genes, we conducted consensus clustering analysis with the “ConsensusCluster” R package to stratify LUAD patients into high- and low-sodium death-associated score groups. Kaplan-Meier survival analysis was performed to evaluate differences in overall survival (OS) between groups, and pathway enrichment analysis was conducted to identify pathways associated with each group.

Construction of prognostic model and validation

To mitigate overfitting, we applied least absolute shrinkage and selection operator (LASSO) regression for prognostic model construction. The risk score was calculated as: risk score = Σ (coefficient × gene expression). LUAD patients were stratified into high- and low-risk groups based on the median risk score. Kaplan-Meier survival curves were generated for these groups using the “survival” and “survminer” R packages. Receiver operating characteristic (ROC) curves at 1, 3, and 5 years were constructed using the “timeROC” package to evaluate the model’s predictive performance. The GSE31210 and GSE50081 datasets were used as external validation cohorts.

Nomogram construction

A prognostic nomogram was established by integrating the risk score with clinical parameters, including age, tumor (T)-stage, metastasis (M)-stage, node (N)-stage, tumor stage, and gender. Calibration plots and multi-ROC curves were used to evaluate the accuracy of the nomogram predictions.

Analysis of immune cell infiltration

Immune infiltration scores were calculated using seven algorithms: XCELL, Tumor IMmune Estimation Resource (TIMER), Quantitation of the Tumor Immune Contexture from human RNA-seq data (QUANTISEQ), Microenvironment Cell Populations-counter (MCPCOUNTER), Estimating the Proportions of Immune and Cancer cells (EPIC), Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts - Absolute (CIBERSORT-ABS), and Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT). The relationship between immune scores and risk scores was investigated. Furthermore, the single sample gene set enrichment analysis (ssGSEA) algorithm was employed to assess immune function in the high- and low-risk score groups. Additionally, the expression of immune checkpoint genes was analyzed in the high- and low-risk groups. Through evaluation of immune cell fractions, immune functions, immune scores, and immune checkpoint gene expression, the immune characteristics and differences between the high- and low-risk groups were explored.

Cell communication analysis

We performed cell communication analysis using the “cellchat” R package. CellChat infers biologically meaningful cell-cell communications by assigning a probability value to each interaction and conducting permutation tests. The built-in human cell-cell interaction database (CellChatDB.human) was employed to identify overexpressed ligand-receptor pairs, which were subsequently mapped to protein-protein interaction networks to validate interactions. Active signaling pathways were identified by calculating the combined expression intensity of all ligand-receptor pairs within specific pathways. Network centrality scores were computed to quantify the communication roles of cell populations.

Statistical analysis

All statistical analyses were performed using R software (version 4.2.1). The Wilcoxon test was employed to compare differences between groups. Log-rank tests were used to compare Kaplan-Meier survival curves. Univariate Cox analyses were conducted to identify independent prognostic factors. All P values were two-sided, and values less than 0.05 were considered statistically significant.

Results

TRPM4 co-expression

Co-expression analysis of TRPM4 in LUAD was conducted using a correlation coefficient threshold of |r| >0.5 and a P value cutoff of P=0.001, identifying 40 co-expressed genes.

Identification of NECSOs

Single-cell sequencing data were obtained from GSM5699782 and GSM5699783. The “Seurat” R package was utilized to convert 10× platform data into Seurat objects, followed by quality control filtering with the criteria: nFeature_RNA >500, nFeature_RNA <7,500, percent.mt <15, and nCount_RNA >1,000. The parameter nExp was set to 0.08 × ncell²/10,000, PN to 0.25, and PK to 0.09. Doublet cells were removed using the “DoubletFinder” R package, and batch effects were corrected with the “Harmony” R package. Cell clustering was performed using the “FindNeighbors” and “FindClusters” functions, resulting in 20 distinct cell clusters. Marker genes for each cluster were identified using the “FindAllMarkers” function, enabling the annotation of nine cell types, which were visualized via t-SNE (Figure 1A). A heatmap was generated to display the top five expressed genes for each cell type (Figure 1B). The “AddModuleScore” function was used to calculate the NESCO score for each cell based on the expression levels of 40 genes co-expressed with TRPM4 in LUAD, revealing elevated scores in epithelial cells compared to other cell types (Figure 1C). Violin plots were used to illustrate the distribution of sodium death-associated scores across cell types (Figure 1D). Cells were stratified into high- and low-score groups based on sodium death-associated scores, and differential expression analysis was conducted. The Wilcoxon rank-sum test was applied to determine the statistical significance of DEGs (P_adj<0.05), with other parameters set to default values, yielding 512 upregulated and 410 downregulated genes. GSEA indicated upregulation of the oxidative phosphorylation pathway and MYC targets V1 pathway in the high-score group (Figure 1E,1F). Activation of the oxidative phosphorylation pathway was closely associated with energy depletion, which may trigger TRPM4 activation, leading to sodium influx and subsequent cell swelling and death.

Figure 1 Identification of sodium death-associated genes via single-cell analysis. (A) t-SNE plot of cell subpopulations. (B) t-SNE plot annotated with marker genes. (C) t-SNE plot of sodium death-associated scores calculated using the “AddModuleScore” function. (D) Violin plot illustrating sodium death-associated scores across different cell subpopulations. (E,F) GSEA analysis showing pathway enrichment in the high-score group, including oxidative phosphorylation and MYC targets V1 pathways. GSEA, gene set enrichment analysis; NECSO, sodium death-associated gene; t-SNE, t-distributed stochastic neighbor embedding.

Identification of tumor-associated modules via WGCNA

WGCNA was performed on 932 DEGs to identify key genes associated with tumor development. Outlier samples were initially removed (Figure 2A). The “PickSoftThreshold” function was employed to determine an appropriate soft threshold, which was set to 3 (Figure 2B,2C). The minimum module size was established at 20 genes, and a topological overlap network was constructed. Pearson correlation analysis was subsequently conducted to evaluate the association between each module and tumor traits. The brown module exhibited a significant correlation with tumor characteristics (Figure 2D). A total of 126 genes were identified within the brown module. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed on these 126 genes, revealing significant enrichment in pathways related to phagosome, antigen processing and presentation, and pathogen clearance through antigen presentation. Enrichment was observed in pathways associated with viral myocarditis. Sodium-induced cell death in cardiovascular diseases may be linked to energy depletion and ionic imbalance in cardiomyocytes. The viral myocarditis pathway likely reflects metabolic and cell death mechanisms in stressed cardiomyocytes. Notably, Staphylococcus aureus infection was found to influence the bone marrow hematopoietic microenvironment through chronic inflammation (Figure 2E). Staphylococcus aureus infection pathway may be associated with hypoxic conditions in the TME, particularly as Staphylococcus aureus pulmonary infections can induce hypoxia and cyanosis, triggering immune responses that overlap with TME dynamics. These genes may regulate hypoxic stress and inflammatory responses, potentially impacting the progression of LUAD.

Figure 2 Screening prognostic genes in key modules. (A) Removal of outlier samples. (B) Soft threshold calculation via topological overlap network analysis. (C) Dendrogram of gene clustering. (D) Correlation between module genes and clinical features. (E) KEGG enrichment analysis of identified genes. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Differential expression analysis, univariate Cox analysis, and CNV analysis

Differential expression analysis was performed on 126 genes to compare normal and tumor samples, using the criteria |log₂FC| >1 and FDR <0.05, yielding 28 DEGs. Univariate Cox regression analysis was subsequently conducted with a significance threshold of P<0.05, identifying 14 prognostic genes. A prognostic network diagram was constructed to visualize these genes (Figure 3A). Of these, LST1, SEC11C, FBP1, ALDH2, and CD52 were classified as protective genes, with higher expression associated with improved prognosis. Additionally, CNVs for the 14 prognostic genes were calculated and visualized (Figure 3B).

Figure 3 Prognostic and copy number variation analysis. (A) Network prognostic diagram of 14 genes. (B) Copy number variations of sodium death-associated genes. CNV, copy number variation; NECSO, sodium death-associated gene.

Cell developmental trajectory analysis

Given the highest scores observed for epithelial cells, pseudotime analysis was conducted on this cell type. Epithelial cells were extracted from the GSM5699782 dataset, and their differentiation potential was assessed using CytoTRACE to determine the differentiation starting point (Figure 4A). Subsequently, pseudotime analysis was conducted using the ‘Monocle’ R package, incorporating Epithelial cells into the differentiation trajectory to reveal their evolutionary pattern (Figure 4B,4C). A heatmap illustrating the expression of 14 prognostic genes was generated (Figure 4D), alongside trend plots for six genes exhibiting significant changes (Figure 4E,4F). The risk genes RANBP1, GAPDH, PKM, LDHA, and TPI1 were found to be upregulated in the later stages of epithelial cell pseudotime, whereas the protective genes LST1, FBP1, and SEC11C were downregulated in these stages. LUAD is characterized by malignant proliferation and uncontrolled growth of glandular epithelial cells, with risk genes playing a critical role in the transformation of epithelial cells into LUAD cells during the later stages of differentiation. The expression of the 14 prognostic genes was also visualized at the single-cell level (Figure 4G).

Figure 4 Cell developmental trajectory analysis. (A) CytoTRACE analysis of differentiation potential in epithelial cells. (B,C) Pseudotime trajectory analysis of epithelial cells. (D) Heatmap displaying dynamic expression of 14 sodium death-associated genes, with low expression indicated in blue and high expression in red. (E,F) Trend plots of eight genes with significant changes. (G) Expression of 14 genes in single-cell data. CytoTRACE, Cellular Trajectory Reconstruction Analysis using Gene Counts and Expression; t-SNE, t-distributed stochastic neighbor embedding.

Cluster analysis

Consensus clustering analysis was conducted on 14 prognostic genes using the ConsensusClusterPlus R package. The optimal number of clusters was determined by evaluating the cumulative distribution function (CDF) and consensus heatmap, with k=2 selected, partitioning samples into two subgroups (Figure 5A,5B). Kaplan-Meier survival analysis revealed significant differences in clinical and prognostic characteristics between these subgroups (Figure 5C). KEGG pathway analysis indicated enrichment in pathways associated with the cell cycle, motor proteins, myocyte cytoskeleton, and hematopoietic cell lineage (Figure 5D). Sodium death-associated characteristics were closely linked to energy metabolism, myocardial infarction, and hypoxia, suggesting a potential mechanism involving myocyte cytoskeleton and hematopoietic function. GSEA was performed using the org.Hs.eg.db, clusterProfiler, and enrichplot R packages, with the c2.cp.kegg.symbols.gmt dataset as the reference. GSEA results showed upregulation of pathways, including cell cycle, DNA repair, oocyte meiosis, proteasome, and spliceosome, in the high-expression group of NECSOs, while pathways such as asthma, cell adhesion, and drug metabolism via cytochrome P450 were upregulated in the low-expression group (Figure 5E,5F).

Figure 5 Results of consensus clustering analysis. (A) Consensus heatmap for k=2. (B) Cluster CDF at k=2–9. (C) Kaplan-Meier survival curves for NECSO high and low-expression groups. (D) Bubble plot of KEGG enrichment analysis. (E) GSEA-enriched pathways upregulated in the NECSO high-expression group. (F) GSEA-enriched pathways upregulated in the NECSO low-expression group. CDF, cumulative distribution function; GSEA, gene set enrichment analysis; NECSO, sodium death-associated gene.

Construction of a prognostic model using LASSO regression

To mitigate overfitting, LASSO regression was applied to 14 prognostic genes, resulting in the selection of eight genes (FBP1, PKM, SEC11C, LST1, GAPDH, UBE2S, LDHA, SEC61G) for further analysis (Figure 6A,6B). A risk score was constructed as follows: risk score = FBP1 * (−0.0432666385825583) + PKM * 0.0476905297804326 + SEC11C * (−0.115629620330642) + LST1 * (−0.0913973369455636) + GAPDH * 0.032007089152869 + UBE2S * 0.0754427469488148 + LDHA * 0.314256579208699 + SEC61G * 0.13224704174282. The TCGA cohort was utilized as the training set, while GSE31210 and GSE50081 served as validation sets. LUAD patients were stratified into high- and low-risk groups based on the median risk score. Kaplan-Meier survival analysis was performed, demonstrating that the high-risk group exhibited significantly worse OS compared to the low-risk group in each cohort (P<0.05) (Figure 6C-6E). Five risk genes (PKM, GAPDH, UBE2S, LDHA, SEC61G) were highly expressed in the high-risk group, whereas three protective genes (FBP1, SEC11C, LST1) were lowly expressed in the high-risk group, as visualized in heatmaps (Figure 6F-6H). The area under the curve (AUC) for the TCGA cohort was 0.681 at 1 year, 0.675 at 3 years, and 0.620 at 5 years. For the GSE31210 cohort, the AUC was 0.711 at 1 year, 0.605 at 3 years, and 0.654 at 5 years. For the GSE50081 cohort, the AUC was 0.683 at 1 year, 0.700 at 3 years, and 0.683 at 5 years (Figure 6I-6K). These results indicate robust predictive performance of the model.

Figure 6 LASSO regression modeling of sodium death-associated genes. (A,B) LASSO regression coefficients for 14 NECSO genes. (C-E) Survival curves for the TCGA, GSE31210, and GSE50081 cohorts. (F-H) Risk score distribution, patient survival status, and NECSO gene expression heatmap for the prognostic model. (I-K) ROC curves for each cohort. AUC, area under the curve; LASSO, least absolute shrinkage and selection operator; NECSO, sodium overload; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Construction of nomogram

To enhance the accuracy of the diagnostic model, a prognostic model was constructed by integrating risk scores with clinical features, including age, tumor-node-metastasis (TNM) staging, tumor stage, and gender (Figure 7A). Calibration curves demonstrated a high concordance between the nomogram predictions and actual outcomes, closely aligning with the 45-degree line (Figure 7B). ROC curves were generated for 1-, 3-, and 5-year predictions, with area under the curve (AUC) values of 0.703, 0.714, and 0.651, respectively (Figure 7C-7E).

Figure 7 Construction of nomogram. (A) Nomogram for predicting 1-, 3-, and 5-year survival in the TCGA cohort. (B) Calibration plot demonstrating prediction accuracy. (C-E) Multi-index ROC curves for predicting 1-year (C), 3-year (D), and 5-year (E) survival rates in the TCGA cohort. AUC, area under the curve; M, metastasis; N, node; OS, overall survival; ROC, receiver operating characteristic; T, tumor; TCGA, The Cancer Genome Atlas.

Analysis of the immune cell infiltration

B cells memory, T cells memory resting, macrophages M2, dendritic cells resting, and Mast cells resting were found to be more abundant in the low-risk group, whereas T cells memory activated, NK cells resting, macrophages M0, macrophages M1, and mast cells activated were more prevalent in the high-risk group (Figure 8A). Seven algorithms, including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT, were employed to calculate immune scores, revealing a negative correlation between immune scores and risk scores (Figure 8B). Chemokine receptor (CCR), human leukocyte antigen (HLA), and T cell co-stimulation scores were observed to be higher in the low-risk group, while major histocompatibility complex I (MHC I) scores were higher in the high-risk group (Figure 8C). Expression levels of immune checkpoint genes were analyzed in the high- and low-risk groups, with most genes showing higher expression in the low-risk group (Figure 8D). Notably, a subset of genes, including CD274, IFNG, IL23A, LDHB, LDHC, PVR, TNFRSF9, TNFSF4, TNFSF9, and YTHDF1, exhibited increased expression in the high-risk group.

Figure 8 Immune status analysis of high- and low-risk groups. (A) Differences in immune infiltration between high- and low-risk groups. (B) Correlation between immune infiltration and risk scores calculated using seven algorithms. (C) Differences in immune function between high- and low-risk groups assessed via ssGSEA. (D) Expression differences of immune checkpoint-related genes between high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. APC, antigen-presenting cell; CCR, chemokine receptor; HLA, human leukocyte antigen; MHC, major histocompatibility complex; ssGSEA, single sample gene set enrichment analysis.

Cell communication analysis

Cell communication analysis was conducted using the “CellChat” R package. Compared to epithelial cells with low NECSO expression, epithelial cells with high NECSO expression exhibited upregulation in the COMPLEMENT, MK, and EGF signaling pathways (Figure 9A-9C). In the high NECSO expression group, epithelial cells demonstrated stronger signaling through the COMPLEMENT pathway to macrophages, monocytes, and dendritic cells, as well as enhanced signaling through the MK pathway to immune cells. Additionally, epithelial cells in the high NECSO expression group were found to receive signals from macrophages and dendritic cells via the EGF signaling pathway. In the COMPLEMENT pathway, the key ligand-receptor pair was identified as C3-(ITGAM + ITGB2). Key genes in this pathway, including C3 and C5, were expressed at higher levels in epithelial cells with high NECSO expression compared to those with low NECSO expression, indicating activation of the complement pathway (Figure 9D). It is hypothesized that epithelial cells with high NECSO expression may have an increased risk of malignant transformation. In the MK pathway, the key ligand-receptor pair was MDK-SDC4. Key genes, including MDK, SDC1, SDC4, NCL, and ITGB1, were expressed at higher levels in epithelial cells with high NECSO expression compared to those with low NECSO expression, whereas the key ligand MDK was nearly absent in the low NECSO expression group (Figure 9E). In the EGF pathway, the key ligand-receptor pair was AREG-(EGFR + ERBB2). While AREG expression levels were similar between the high and low NECSO expression groups, EGFR and ERBB2 were highly expressed in the high NECSO expression group but nearly absent in the low NECSO expression group (Figure 9F). It is speculated that epithelial cells with higher sodium death-associated scores may promote growth, proliferation, and differentiation through elevated expression of (EGFR + ERBB2) receptor genes, potentially leading to malignant transformation.

Figure 9 Cell communication analysis. (A-C) Interactions and roles of different cell types in the COMPLEMENT, MK, and EGF signaling pathways. (D-F) Roles of key ligand-receptor pairs in the COMPLEMENT, MK, and EGF signaling pathways and expression of key genes in these pathways across different cell types. EGF, epidermal growth factor; MK, Midkine; NECSO, sodium death-associated gene.

Analysis of NECSOs

In the prognostic model established, UBE2S, SEC61G, SEC11C, PKM, GAPDH, and LDHA exhibited higher expression in tumor samples compared to normal samples, whereas LST1 and FBP1 showed higher expression in normal samples relative to tumor samples (Figure 10A). ROC analysis was performed for the eight genes to evaluate their prognostic relevance, with LDHA, GAPDH, and PKM identified as the top three genes based on ROC performance (Figure 10B-10D). Subsequently, random forest analysis was conducted on these eight genes, revealing that LDHA, GAPDH, and PKM contributed most significantly to tumor development (Figure 10E,10F). These three genes play a critical role in NECSO-related mechanisms in LUAD.

Figure 10 Identification of sodium death-associated genes. (A) Expression differences of eight risk genes between normal and tumor samples. (B-D) Genes most closely associated with prognosis identified through screening. (E,F) Genes contributing most significantly to tumor development identified via random forest analysis. ***, P<0.001. AUC, area under the curve.

Discussion

Currently, multiple strategies are employed to combat LUAD, including chemotherapy, surgical resection, radiotherapy, and molecular targeted therapies (15). Despite these approaches, the prognosis for LUAD patients remains poor, underscoring the urgent need to develop prognostic models for LUAD prediction and to identify potential therapeutic targets. Immunotherapy, as an emerging treatment modality, has shown promise, with prior studies demonstrating that inhibitors of immune checkpoints cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), programmed cell death protein 1 (PD-1), and programmed cell death ligand 1 (PD-L1) induce durable remission in subsets of patients with metastatic disease (16). Additionally, delta-HE has been identified as a novel predictive and prognostic biomarker for PD-1/PD-L1 inhibitor therapy in non-small cell lung cancer (NSCLC) patients (17). Immunotherapy markedly enhances patients’ anti-tumor immune responses, necessitating further investigation into the immune characteristics of NECSO-derived models and NECSO-related signaling pathways within the TME.

In this study, single-cell analysis was utilized to elucidate the immune microenvironment associated with NECSO in LUAD. The TCGA cohort was employed as the training set, with two GEO datasets serving as validation sets. A prognostic model was constructed using LASSO regression, incorporating the genes FBP1, PKM, SEC11C, LST1, GAPDH, UBE2S, LDHA, and SEC61G. FBP1, a key gluconeogenesis enzyme, catalyzes the hydrolysis of fructose-1,6-bisphosphate to fructose-6-phosphate, promoting glucose synthesis. In hepatocellular carcinoma (HCC), the MAGE-TRIM28 complex accelerates the Warburg effect and HCC progression by degrading FBP1 (18). PKM, particularly the PKM1/PKM2 ratio, is altered in tumors compared to normal controls. PKM2 is highly expressed in ovarian cancer, inducing elevated glycolysis rates and contributing to malignant invasion and metastasis through epithelial-mesenchymal transition (EMT). PKM2 inhibitors suppress ovarian cancer cell migration and growth by disrupting the Warburg effect (19). SEC11C is closely associated with signal peptide cleavage (20). The LST1 protein functions as a transmembrane adapter protein that inhibits signal transduction and serves as a membrane scaffold to promote the formation of tunneling nanotubes (21). Combining immunotherapy with GAPDH targeting to induce ferroptosis in LUAD may offer a novel therapeutic strategy. Deleterious single-nucleotide polymorphisms (D39Y and S51Y) are located in functionally conserved domains of GAPDH. Moreover, post-translational modifications in the disordered region of GAPDH may enhance its glycolytic activity in LUAD (22,23). UBE2S promotes chemoresistance in HCC via PTEN-AKT signaling (24). LDHA enhances NLRP3 lactylation to induce pyroptosis in cardiomyocytes, thereby promoting myocardial ischemia/reperfusion injury (25). Similarly, elevated glycolysis remains a hallmark of cancer metabolism, largely driven by dysregulated metabolic enzymes. Lactate dehydrogenase A facilitates glycolysis by converting pyruvate to lactate (26). SEC61G promotes glycosylation, stabilization, and membrane presentation of immune checkpoint ligands (ICLs), facilitating immune evasion and tumorigenesis, and enhances cervical cancer proliferation via activation of the MAPK signaling pathway (27,28).

To investigate the relationship between the risk scores of our constructed model and immune characteristics, multiple algorithms, including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT, were employed to calculate immune scores, which exhibited a negative correlation with risk scores. The results revealed that CCR, HLA, and T cell co-stimulation were expressed at higher levels in the low-risk group, whereas MHC I was more prominently expressed in the high-risk group. Selective inhibition of CCRs that recruit immunosuppressive cells or enhance tumor migration and survival, while preserving receptors essential for normal immune cell trafficking, holds promise for significantly improving therapeutic outcomes. For instance, the CCR4-CCL17/CCL22 axis promotes an immunosuppressive microenvironment by recruiting regulatory T cells (Tregs), a key mechanism of tumor immune evasion. Inhibitors targeting this axis may enhance anti-tumor immune responses (29,30). The polymorphism and expression levels of HLA molecules are closely associated with tumor development, as they regulate tumor cell proliferation and suppress anti-tumor immune responses (31). Innate immune signaling mediated by TLRs and STING can modulate T cell receptor (TCR) signaling pathways, thereby influencing T cell functional states (32). In the context of MHC I, dendritic cells can elicit anti-tumor CD8 T cell responses through cross-presentation via the MHC (33).

To further explore the differential immune responses of epithelial cells with varying NECSO scores, cell communication analysis was conducted. The COMPLEMENT, MK, and EGF pathways were upregulated in epithelial cells with high NECSO scores. The complement system is traditionally recognized as an immune resistance mechanism against cancer, known to facilitate antibody-dependent cytotoxicity in immunotherapy (34). However, studies also indicate that complement activation within the TME may promote tumor growth and enhance metastasis (35). The MDK-SDC4 signaling pathway operates through the interaction of the heparin-binding growth factor MDK with its receptor SDC4, a transmembrane heparan sulfate proteoglycan (36). This pathway plays a central role in regulating cell proliferation, migration, and extracellular matrix remodeling, which are critical for tissue repair and regeneration but may also enhance tumor invasiveness and metastatic potential (37). In the EGF pathway, approximately half of NSCLC patients harbor epidermal growth factor receptor (EGFR) gene mutations. Compared to the EGFR high-expression group, the EGFR low-expression group exhibits lower rates of TP53 co-mutations and EGFR amplification but higher rates of EGFR subclonal mutations (38).

Although the role of NECSOs has been explored at the single-cell level, and the prognostic model derived from these genes demonstrates predictive accuracy for LUAD prognosis while characterizing immune features, further enhancements are required. Specifically, expanding the number of single-cell samples and cells, as well as increasing the sample size at the bulk transcriptome level, is critical to support more comprehensive and robust investigations.

Conclusions

In summary, a prognostic model based on NECSOs was constructed using the TCGA cohort as the training set and validated with the GSE31210 and GSE50081 cohorts. The model was further optimized by incorporating clinical data through nomogram construction. Single-cell analysis revealed that epithelial cells with high NECSOs expression exhibited upregulation in the COMPLEMENT, MK, and EGF signaling pathways.

Acknowledgments

We express our gratitude to all contributors who supported this study. The datasets utilized are publicly accessible from the following sources: Single-cell data were obtained from GSM5699782 and GSM5699783 (https://www.ncbi.nlm.nih.gov/geo/). Transcriptome data and CNV data were downloaded from TCGA (https://portal.gdc.cancer.gov/), comprising 501 LUAD clinical samples. Validation cohorts included GSE31210 (n=226) and GSE50081 (n=127) datasets from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The h.all.v2023.1.Hs.symbols.gmt and c2.cp.kegg.symbols.gmt datasets were acquired from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb). We sincerely thank these databases for providing critical data on LUAD for this study.

Footnote

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

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

Funding: The study was supported by the Discipline (Professional) Leader Cultivation Project of the Anhui Provincial University Young and Middle-aged Teacher Training Program (DTR2025029) and the Provincial Base for High-Level Cultivation of Basic Discipline Students.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-919/coif). The authors have no 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. Cao M, Chen W. Epidemiology of lung cancer in China. Thorac Cancer 2019;10:3-7. [Crossref] [PubMed]
2. Schabath MB, Cote ML. Cancer Progress and Priorities: Lung Cancer. Cancer Epidemiol Biomarkers Prev 2019;28:1563-79. [Crossref] [PubMed]
3. Liang H, Xu Y, Zhao J, et al. Hippo pathway in non-small cell lung cancer: mechanisms, potential targets, and biomarkers. Cancer Gene Ther 2024;31:652-66. [Crossref] [PubMed]
4. Xu JY, Zhang C, Wang X, et al. Integrative Proteomic Characterization of Human Lung Adenocarcinoma. Cell 2020;182:245-261.e17. [Crossref] [PubMed]
5. Hirsch FR, Scagliotti GV, Mulshine JL, et al. Lung cancer: current therapies and new targeted treatments. Lancet 2017;389:299-311. [Crossref] [PubMed]
6. Fu W, Wang J, Li T, et al. Persistent activation of TRPM4 triggers necrotic cell death characterized by sodium overload. Nat Chem Biol 2025;21:1238-49. [Crossref] [PubMed]
7. Guinamard R, Bouvagnet P, Hof T, et al. TRPM4 in cardiac electrical activity. Cardiovasc Res 2015;108:21-30. [Crossref] [PubMed]
8. Hu Y, Cang J, Hiraishi K, et al. The Role of TRPM4 in Cardiac Electrophysiology and Arrhythmogenesis. Int J Mol Sci 2023;24:11798. [Crossref] [PubMed]
9. Arullampalam P, Essers MC, Boukenna M, et al. Knockdown of the TRPM4 channel alters cardiac electrophysiology and hemodynamics in a sex- and age-dependent manner in mice. Physiol Rep 2023;11:e15783. [Crossref] [PubMed]
10. Sarmiento D, Montorfano I, Cerda O, et al. Increases in reactive oxygen species enhance vascular endothelial cell migration through a mechanism dependent on the transient receptor potential melastatin 4 ion channel. Microvasc Res 2015;98:187-96. [Crossref] [PubMed]
11. Borgström A, Peinelt C, Stokłosa P. TRPM4 in Cancer-A New Potential Drug Target. Biomolecules 2021;11:229. [Crossref] [PubMed]
12. Stokłosa P, Borgström A, Hauert B, et al. Investigation of Novel Small Molecular TRPM4 Inhibitors in Colorectal Cancer Cells. Cancers (Basel) 2021;13:5400. [Crossref] [PubMed]
13. Borgström A, Hauert B, Kappel S, et al. Small Molecular Inhibitors Block TRPM4 Currents in Prostate Cancer Cells, with Limited Impact on Cancer Hallmark Functions. J Mol Biol 2021;433:166665. [Crossref] [PubMed]
14. Ceylan GG, Önalan EE, Kuloğlu T, et al. Potential role of melastatin-related transient receptor potential cation channel subfamily M gene expression in the pathogenesis of urinary bladder cancer. Oncol Lett 2016;12:5235-9. [Crossref] [PubMed]
15. Miller KD, Goding Sauer A, Ortiz AP, et al. Cancer Statistics for Hispanics/Latinos, 2018. CA Cancer J Clin 2018;68:425-45. [Crossref] [PubMed]
16. Wang SJ, Dougan SK, Dougan M. Immune mechanisms of toxicity from checkpoint inhibitors. Trends Cancer 2023;9:543-53. [Crossref] [PubMed]
17. Hong G, Lee SI, Kang DH, et al. Delta-He as a Novel Predictive and Prognostic Biomarker in Patients With NSCLC Treated With PD-1/PD-L1 Inhibitors. Cancer Med 2025;14:e70826. [Crossref] [PubMed]
18. Jin X, Pan Y, Wang L, et al. MAGE-TRIM28 complex promotes the Warburg effect and hepatocellular carcinoma progression by targeting FBP1 for degradation. Oncogenesis 2017;6:e312. [Crossref] [PubMed]
19. Wang Y, Xu N, Ndzie Noah ML, et al. Pyruvate Kinase M1/2 Proteoformics for Accurate Insights into Energy Metabolism Abnormity to Promote the Overall Management of Ovarian Cancer Towards Predictive, Preventive, and Personalized Medicine Approaches. Metabolites 2025;15:203. [Crossref] [PubMed]
20. Liaci AM, Steigenberger B, Telles de Souza PC, et al. Structure of the human signal peptidase complex reveals the determinants for signal peptide cleavage. Mol Cell 2021;81:3934-3948.e11. [Crossref] [PubMed]
21. Weidle UH, Rohwedder I, Birzele F, et al. LST1: A multifunctional gene encoded in the MHC class III region. Immunobiology 2018;223:699-708. [Crossref] [PubMed]
22. Ouyang X, Zhu R, Lin L, et al. GAPDH Is a Novel Ferroptosis-Related Marker and Correlates with Immune Microenvironment in Lung Adenocarcinoma. Metabolites 2023;13:142. [Crossref] [PubMed]
23. John P, Sudandiradoss C. Structure, function and stability analysis on potential deleterious mutation ensemble in glyceraldehyde 3-phosphate dehydrogenase (GAPDH) for early detection of LUAD. Life Sci 2024;358:123127. [Crossref] [PubMed]
24. Gui L, Zhang S, Xu Y, et al. UBE2S promotes cell chemoresistance through PTEN-AKT signaling in hepatocellular carcinoma. Cell Death Discov 2021;7:357. [Crossref] [PubMed]
25. Fang L, Yu Z, Qian X, et al. LDHA exacerbates myocardial ischemia-reperfusion injury through inducing NLRP3 lactylation. BMC Cardiovasc Disord 2024;24:651. [Crossref] [PubMed]
26. Feng Y, Xiong Y, Qiao T, et al. Lactate dehydrogenase A: A key player in carcinogenesis and potential target in cancer therapy. Cancer Med 2018;7:6124-36. [Crossref] [PubMed]
27. Zeng K, Zeng Y, Zhan H, et al. SEC61G assists EGFR-amplified glioblastoma to evade immune elimination. Proc Natl Acad Sci U S A 2023;120:e2303400120. [Crossref] [PubMed]
28. Fan Y, Wang Y, Liu F, et al. SEC61G Promotes Cervical Cancer Proliferation by Activating MAPK Signaling Pathway. Dis Markers 2022;2022:7016079. [Crossref] [PubMed]
29. Kraus S, Kolman T, Yeung A, et al. Chemokine Receptor Antagonists: Role in Oncology. Curr Oncol Rep 2021;23:131. [Crossref] [PubMed]
30. Goenka A, Khan F, Verma B, et al. Tumor microenvironment signaling and therapeutics in cancer progression. Cancer Commun (Lond) 2023;43:525-61. [Crossref] [PubMed]
31. Liu DH, Mou FF, An M, et al. Human leukocyte antigen and tumor immunotherapy Int J Oncol 2023;62:68. (Review). [Crossref] [PubMed]
32. Imanishi T, Saito T. T Cell Co-stimulation and Functional Modulation by Innate Signals. Trends Immunol 2020;41:200-12. [Crossref] [PubMed]
33. MacNabb BW, Tumuluru S, Chen X, et al. Dendritic cells can prime anti-tumor CD8(+) T cell responses through major histocompatibility complex cross-dressing. Immunity 2022;55:982-997.e8. [Crossref] [PubMed]
34. Meri S, Magrini E, Mantovani A, et al. The Yin Yang of Complement and Cancer. Cancer Immunol Res 2023;11:1578-88. [Crossref] [PubMed]
35. Afshar-Kharghan V. The role of the complement system in cancer. J Clin Invest 2017;127:780-9. [Crossref] [PubMed]
36. Hashimoto M, Kojima Y, Sakamoto T, et al. Spatial and single-cell colocalisation analysis reveals MDK-mediated immunosuppressive environment with regulatory T cells in colorectal carcinogenesis. EBioMedicine 2024;103:105102. [Crossref] [PubMed]
37. Chang X, Wang C, Wang F, et al. Global research trends of tumor microenvironment in non-small cell lung cancer with epidermal growth factor receptor mutation: a bibliometric analysis from 2014 to 2023. Front Immunol 2025;16:1555216. [Crossref] [PubMed]
38. Torasawa M, Yoshida T, Shiraishi K, et al. Implications of EGFR expression on EGFR signaling dependency and adaptive immunity against EGFR-mutated lung adenocarcinoma. Lung Cancer 2025;202:108494. [Crossref] [PubMed]