Identification and experimental verification of necroptosis-related prognostic gene signature and characterization of tumor immune infiltration in lung squamous cell carcinoma

Datasets and data processing

This investigation employed a retrospective design, with experimental data derived from methodologies previously established (Sun et al., 2022). Transcriptomic profiles (RNA-Seq FPKM values), demographic variables, and survival outcomes of TCGA-LUSC cases were retrieved from the UCSC Xena repository (https://xena.ucsc.edu/; complete metadata) (Goldman et al., 2020). Clinical and pathological parameters associated with NRG expression in LUSC were systematically summarized in Table 1. To strengthen analytical reliability, additional LUSC cohorts were incorporated from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), including GSE19188 (24 LUSC samples), GSE41271 (80 LUSC samples), GSE42127 (43 LUSC samples), GSE81089 (67 LUSC samples), GSE157009 (249 LUSC samples), GSE30219 (61 LUSC samples), GSE73403 (69 LUSC samples), GSE8894 (75 LUSC samples), and GSE74777 (107 LUSC samples) (Barrett et al., 2013). Independent validation cohort with comparable clinicopathological characteristics. These datasets were selected because they provide comprehensive gene expression data and clinical information, making them suitable for external validation. A curated set of 159 NRGs was extracted from the KEGG database, with complete expression profiles provided in Table S1 (Kanehisa & Goto, 2000). Validation was conducted through cross-referencing with published literature and empirical evidence concerning necroptotic signaling pathways. All computational procedures were implemented in R statistical software (version 4.0.3).

NRGs differential expression, survival and mutation analysis

Comprehensive bioinformatics evaluation identified 35 differentially expressed NRGs (DENRGs) with significant differential expression in LUSC (Table S2). Expression profiles of tumor and adjacent normal lung tissues were compared using the limma package with stringent cutoffs (—log₂FC—>1, adjusted P-value < 0.05) (Liu et al., 2021). Visualization of results employed several R packages: ComplexHeatmap for hierarchical clustering, Enhanced Volcano for differential expression mapping, ggplot2 for quantitative profiling, and survminer for Kaplan–Meier curve construction (Ito & Murphy, 2013; Gu, 2022). Survival outcomes were further assessed using the log-rank test in conjunction with Cox proportional hazards regression. In addition, mutational features of NRGs in LUSC were systematically characterized with the maftools package (Mayakonda et al., 2018).

Functional annotation and pathway analysis

To clarify the biological relevance of differentially expressed NRGs, functional enrichment analyses were conducted. Gene Ontology (GO) annotation delineated major biological processes, molecular functions (MF), and cellular components (CC) linked to these genes (The Gene Ontology Consortium, 2019). KEGG pathway analysis further mapped the associated systemic functional networks. Both analyses were graphically represented and interpreted using ggplot2 to provide an integrated view of enrichment outcomes.

Development and validation of necroptosis-related prognostic signature

Candidate NRGs with prognostic significance were identified through univariate and multivariate Cox regression analyses, yielding independent biomarkers for model construction. Genes selected at λ.1se were entered into a multivariate Cox proportional hazards framework, and only statistically significant variables (P < 0.05) were preserved in the final model. A Least Absolute Shrinkage and Selection Operator (LASSO) regression approach incorporating three genes was applied to establish a prognostic signature. Risk scores were calculated as NRGs_score = Σ(Expi × coefi), and patients were classified into high- or low-risk groups based on the median cutoff. Kaplan–Meier survival analysis and ROC curve evaluation were employed to assess predictive accuracy (Kamarudin, Cox & Kolamunnage-Dona, 2017). A nomogram integrating clinical factors was further generated to predict 1-, 3-, and 5-year survival. All analyses were conducted using dedicated R packages, and independent GEO cohorts served as external validation datasets.

Cell lines and cell culture

The BEAS-2B human bronchial epithelial cell line and the NCI-H520 human LUSC cell line were obtained from Sangon Biotech (Shanghai, China). BEAS-2B cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM; Gibco) containing 10% fetal bovine serum (FBS; Gibco), whereas NCI-H520 cells were maintained in Roswell Park Memorial Institute 1640 medium (RPMI-1640; Gibco) supplemented with 10% FBS (Gibco). All cultures were incubated at 37 °C in a humidified environment with 5% CO₂.

RNA extraction and RT-qPCR assay

mRNA expression of CAMK2A, CHMP4C, and PYGB was assessed by RT-qPCR. Total RNA was isolated with TRIzol reagent (Ambion) following standard procedures, and complementary DNA (cDNA) was synthesized using 5 × HiScript qRT SuperMix II (Vazyme Biotech). Quantitative amplification was carried out on an Mx3000P QPCR system (Agilent Technologies, Santa Clara, CA) with ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech) under the manufacturer’s recommended conditions. The cycling protocol consisted of an initial denaturation at 95 °C for 10 s, followed by 40 cycles of 95 °C for 15 s and 60 °C for 60 s for annealing/extension. Relative expression levels were calculated using the 2^−ΔΔCt method, with GAPDH as the internal reference. Each sample was analyzed in triplicate to ensure reproducibility. Primer sequences for GAPDH and the target genes were as follows: GAPDH: Sequence (5′→3′).

Forward: TCAAGAAGGTGGTGAAGCAGG, Reverse: TCAAAGGTGGAGGAGTGGGT, 115 bp.

CAMK2A: Sequence (5′->3′).

Forward: CAGAAGTGCTGCGGAAGG, Reverse: ATGCGTTTGGATGGGTTA, 235 bp.

CHMP4C: Sequence (5′->3′).

Forward: GTTGGCTTTGGTGATGACTT, Reverse CTGGTTTTCTATTTGGCTGT, 148 bp.

PYGB: Sequence (5′->3′).

Forward: ACGGCTATGGAATCCGCTAT, Reverse: TGTTGACGGTGTTGTTCTTGT, 255 bp.

Clinical specimen analysis and immunohistochemical validation

A total of 21 paired LUSC and adjacent non-tumorous tissues were collected from surgical resections at Liuzhou People’s Hospital between 2019 and 2021. This study was approved by the Medical Ethics Committee of Liuzhou People’s Hospital (Reference No. KY2021-025-02) and was performed according to the Declaration of Helsinki, and written informed consent was obtained from all participants. After routine tissue preparation, immunohistochemical staining was carried out with antibodies against CAMK2A, CHMP4C, and PYGB. Quantification of staining intensity was achieved by calculating integrated optical density (IOD) values through Image-ProPlus6.0 software. The corresponding clinical parameters are presented in Table S3.

Immune cell infiltration, immune checkpoints, TMB and MSI analysis

The R packages “GSVA,” “immunedeconv,” “estimate,” “ggplot2,” “pheatmap,” and “ggstatsplot” were applied to evaluate the relationship between the five prognostic NRGs and immune cell infiltration, using three advanced algorithms: ssGSEA, ESTIMATE, and CIBERSORT (Chen et al., 2022; Chen et al., 2018). These approaches, together with immune profiling of malignant tumor tissues, provide essential insights into the tumor microenvironment (TME) and its relevance for prognosis and therapeutic response. ssGSEA, an adaptation of the Gene Set Enrichment Analysis (GSEA) algorithm, generates enrichment scores for each sample–gene set pair rather than across grouped samples and pathways (https://www.genepattern.org/modules/docs/ssGSEAProjection/4). ESTIMATE calculates tumor purity and quantifies stromal and immune cell infiltration from expression data, with results accessible across TCGA tumor types and platforms. CIBERSORT deconvolutes bulk RNA-sequencing data using a reference gene expression matrix to estimate the relative abundance of 22 immune cell types, thereby enabling detailed characterization of immune landscapes and their association with disease features or therapeutic outcomes. In addition, stromal, immune, and ESTIMATE scores were examined in parallel with sixty immune checkpoint molecules, 150 marker genes from five immune pathways, TMB, and MSI. All statistical analyses and data visualizations were conducted using R version 4.0.3 (Ito & Murphy, 2013; Chen et al., 2022; Chen et al., 2018; Sturm, Finotello & List, 2020; R Core Team, 2020).

Statistical analysis

Multiple analytical strategies were applied for statistical evaluation. Differential expression was determined by fold-change (FC) values, whereas survival outcomes were examined using hazard ratios (HR) and log-rank tests. Associations between variables were quantified with Spearman’s or Pearson’s correlation analyses, selected according to distributional properties of the data. Statistical significance was defined as P < 0.05, derived from conventional hypothesis testing or log-rank comparisons.

Screening of DENGRs in LUSC and normal lung tissues

The workflow of the study is presented in Fig. 1. A total of 159 NRGs were retrieved from the KEGG database for subsequent analysis of differential expression between LUSC and normal lung tissues. Transcriptomic profiling using the UCSC Xena database identified 5,357 DEGs, which were illustrated by a heat map (Fig. 2A) and a volcano plot (Fig. 2B). Functional enrichment analysis with KEGG and GO revealed that these DEGs were primarily associated with pathways related to the cell cycle, DNA replication, and immune response (Figs. 2C, 2D). Integration of the 159 NRGs with the 5,357 DEGs yielded 35 DENRGs through Venn diagram analysis, comprising 15 genes upregulated and 20 genes downregulated relative to normal controls, which were subsequently subjected to further investigation (Figs. 3A–3C).

Analysis of mutation and functional enrichment for 35 DENRGs

Mutation profiling and functional enrichment of 35 DENRGs in LUSC revealed distinct genetic and biological characteristics. Mutations were present in 31.71% of samples (156/492), with missense mutations representing the predominant type. Single nucleotide polymorphisms (SNPs) accounted for the majority of alterations, with C > T transitions as the most frequent single nucleotide variant (SNV) category. Among the NRGs, NLRP3, TLR4, STAT4, and PYGM exhibited the highest mutation frequencies (Figs. S1A, S1B). Functional annotation through GO and KEGG analyses delineated the biological relevance of these genes (Fig. 3D, Fig. S2). Enriched BP terms were primarily associated with apoptotic and necroptotic pathways, cellular responses to peptides, and I-kappaB kinase/NF-kappaB signaling (Fig. S2A). The encoded proteins were mainly localized within endosome membranes, membrane microdomains, membrane rafts, and the cytosolic compartment, as reflected in CC terms. MF enrichment highlighted tumor necrosis factor receptor superfamily binding, death receptor binding, ubiquitin-like protein ligase binding, and phospholipase A2 activity. Z-score clustering of GO terms indicated enrichment in necroptotic process, programmed necrotic cell death, necrotic cell death, and positive regulation of NF-kappaB transcription factor activity (Fig. S2B). KEGG pathway analysis further indicated associations with necroptosis, viral infections (Influenza A, Measles, Hepatitis, Epstein-Barr virus), NOD-like receptor signaling, TNF signaling, apoptosis, and Th1/Th2 cell differentiation (Fig. 3D). KEGG z-score clustering confirmed significant enrichment of the 35 DENRGs in pathways such as “Necroptosis” (hsa04217) and “Influenza A” (hsa05164), characterized by elevated z-scores and low P.adjust values (Fig. S2C).

Construction and validation of a prognostic nomogram incorporating NRGs in LUSC

A prognostic model for LUSC grounded in necrosis-related genes was constructed by combining LASSO regression with multivariate Cox regression analyses. This integrative strategy yielded an optimized three-gene signature, defined by CAMK2A, CHMP4C, and PYGB, with a penalty coefficient (λ) of 3 (Figs. 4A, 4B). The risk score was formulated as: Risk Score = (0.1269  × CAMK2A) + (0.2268  × PYGB) + (0.1256  × CHMP4C). Patients were subsequently categorized into high- and low-risk cohorts according to the median risk score. Figure 4C illustrates gene expression profiles, distribution of risk scores, and corresponding survival outcomes across the two subgroups. A clear positive association was observed between elevated risk scores and both higher mortality risk and shorter survival duration. Kaplan–Meier analysis revealed markedly reduced OS in the high-risk group compared with the low-risk group, with median OS times of 3.0 and 5.7 years, respectively (HR = 1.519, P = 0.0027). Time-dependent receiver operating characteristic (ROC) analysis further validated the prognostic performance of the signature, with area under the ROC curve (AUC) values of 0.588, 0.621, and 0.636 for 1-, 3-, and 5-year OS prediction, respectively (Figs. 4D–4E). To assess whether the signature provided independent prognostic information in LUSC, both univariate and multivariate Cox regression analyses were applied (Figs. 5 and 6). In the multivariate model, the risk score retained significance as an independent prognostic factor (HR = 2.76, P = 0.0012), while the pTNM stage showed only borderline significance (HR = 1.60, P = 0.01) (Figs. 5A, 5B). Based on these results, a nomogram incorporating the risk score and clinicopathological variables was established. The concordance index (C-index) of 0.627 (P < 0.05) reflected moderate predictive performance (Fig. 5C). Within the full cohort, the nomogram demonstrated reliable prediction of 1-, 3-, and 5-year OS, with estimates closely aligned to the ideal reference values (Fig. 5D).

Univariate Cox regression analysis demonstrated the prognostic relevance of three NRGs (Figs. 6A, 6B). Elevated CAMK2A expression correlated with reduced overall survival (OS) (HR = 1.35, P = 0.029). Similarly, increased CHMP4C expression was linked to shorter OS (HR = 1.47, P = 0.005) and disease-specific survival (DSS) (HR = 1.70, P = 0.016). PYGB overexpression was associated with inferior OS (HR = 1.62, P = 0.001), progression-free survival (PFS) (HR = 1.65, P = 0.003), and DSS (HR = 2.05, P = 0.001) (Figs. 6A, 6B). Validation using GEO datasets confirmed these associations (Fig. 6C). High CAMK2A expression consistently predicted poor OS in four cohorts (GSE19188: P = 0.034; GSE41271: P = 0.010; GSE42127: P = 0.007; GSE81089: P = 0.012). CHMP4C expression correlated with inferior OS in two cohorts (GSE81089: P = 0.0084; GSE157009: P = 0.0031). PYGB was a robust indicator of poor OS in four cohorts (GSE19188: P = 0.034; GSE30219: P= 3e−04; GSE73403: P = 0.0079; GSE81089: P = 0.00096), and was further predictive of disease-free survival (DFS) in GSE30219 (P = 0.014) and recurrence-free survival (RFS) in GSE8894 (P = 0.0084). Collectively, consistent overexpression of CAMK2A, CHMP4C, and PYGB across multiple datasets indicated an unfavorable prognostic impact in LUSC across diverse survival metrics.

Mutation analysis of three NRGs in LUSC

To clarify the mutational characteristics of the three NRGs in LUSC, mutation profiles were analyzed in 18 TCGA-LUSC samples. All cases exhibited genetic alterations in CAMK2A, CHMP4C, and PYGB. The mutation spectrum comprised missense, nonsense, and multi-hit variants, with missense mutations representing the predominant type (Fig. 7A). Among the three genes, CAMK2A displayed the highest alteration frequency (50%), followed by PYGB (39%) and CHMP4C (11%). Copy number variation (CNV) assessment demonstrated extensive heterozygous deletion in CAMK2A, while CHMP4C and PYGB showed pronounced heterozygous amplification, with gene-specific CNV frequencies quantified. mRNA levels of CAMK2A and PYGB exhibited strong positive associations with CNV status, whereas CHMP4C showed only a modest correlation (Fig. 7B). Single nucleotide polymorphisms (SNPs) were the most frequent nucleotide alterations, with C>A (n = 12) and C>T (n = 6) transitions dominating the mutational landscape (Fig. 7C). Methylation profiling revealed a significant inverse relationship between methylation status and transcript levels of CHMP4C and PYGB, in contrast to CAMK2A, where no such correlation was observed. These results point to gene-specific epigenetic regulation that selectively influences the transcriptional activity of CHMP4C and PYGB (Fig. 7D).

Analysis of mRNA and protein expression of three NRGs in LUSC

Differential expression of the three NRGs between LUSC and normal lung tissues in TCGA was verified previously (Fig. 3C). Subsequent validation using RT-qPCR and IHC was performed to evaluate their expression in clinical LUSC and normal lung samples. Quantitative RT-qPCR analysis assessed CAMK2A, CHMP4C, and PYGB mRNA levels in BEAS-2B bronchial epithelial cells and NCI-H520 LUSC cells. CAMK2A expression was significantly higher in BEAS-2B cells than in NCI-H520 cells (P < 0.01) (Fig. 8A). Conversely, CHMP4C expression was markedly increased in NCI-H520 cells, with minimal expression in BEAS-2B cells (P < 0.001) (Fig. 8B). PYGB also demonstrated elevated expression in NCI-H520 cells compared with BEAS-2B cells (P < 0.05) (Fig. 8C). Analysis of 21 paired LUSC tumor and adjacent normal tissues revealed that CHMP4C and PYGB were predominantly localized to the nucleus and cytoplasm of tumor cells, with positive immunoreactivity indicated by brown staining (Figs. 8E, 8F). In contrast, their expression was weak or undetectable in normal tissues (Figs. 8E, 8F). CAMK2A displayed the opposite pattern, with higher expression in normal tissues and reduced levels in tumor tissues (Fig. 8D). Immunohistochemical evaluation confirmed markedly increased CHMP4C and PYGB expression in LUSC compared with adjacent non-tumor tissues (P < 0.001), while CAMK2A expression was significantly enriched in normal tissues (P < 0.001) (Fig. 8G). The protein expression profiles aligned with mRNA expression data from the TCGA-LUSC cohort (Fig. 3C). Elevated CAMK2A, CHMP4C, and PYGB expression correlated with unfavorable prognosis, highlighting their value as prognostic indicators. Notably, CAMK2A was consistently higher in normal lung tissue, suggesting a potential protective role in maintaining pulmonary homeostasis that becomes attenuated during carcinogenesis. In contrast, CHMP4C and PYGB were overexpressed in LUSC and associated with poor outcomes, implying roles in tumor progression and metastatic potential. The distinct expression patterns and prognostic associations highlight the need for further mechanistic and translational investigations to clarify their contribution to LUSC pathogenesis and therapeutic sensitivity.

Verification of the prognostic model

The predictive model was externally validated using the GSE74777 dataset. Patients were stratified into high- and low-risk groups based on the median risk score (Figs. 9A–9C). Heatmap visualization indicated an association between risk scores and the three NRGs (Fig. 9A). ROC analysis demonstrated consistent predictive capacity for 1-, 3-, and 5-year survival across all four datasets, with AUC values exceeding 0.6 (Fig. 9D). Kaplan–Meier curves showed a trend toward reduced overall survival in the high-risk group, although statistical significance was not reached (Fig. 9E). Multivariate Cox regression further evaluated the influence of risk score alongside clinical parameters including T stage, gender, N stage, and age in LUSC patients (Fig. 9F). Collectively, the model maintained reliable performance in the GSE74777 cohort.

The correlation between the expression levels of three NRGs and the clinical features in patients with LUSC

Analysis of TCGA-LUSC data revealed distinct prognostic patterns for the three NRGs. CAMK2A expression was significantly higher in patients older than 65 compared with those at stages I–II (P < 0.05) (Table 2). CHMP4C expression was elevated in stage III–IV patients relative to stage I–II (P < 0.05) (Table 3) and was also increased in individuals with progressive disease compared with those classified as stable disease (SD), partial response (PR), or complete response (CR) (P < 0.05) (Table 3). By contrast, PYGB expression exhibited no significant association with the clinical features of LUSC (P < 0.05) (Table 4).

The three NRGs were associated with tumor immune infiltration in LUSC

In LUSC, the expression profiles of three NRGs display significant associations with clinical characteristics, while tumor-infiltrating lymphocytes function as independent indicators of tumor stage, grade, and lymph node involvement. Using TCGA data, the associations between the expression of these prognostic NRGs and immune infiltration were systematically examined. The tumor microenvironment, composed of malignant, stromal, and infiltrating immune cells, plays a decisive role in determining clinical outcomes, with immune infiltration recognized as an independent predictor of sentinel lymph node status and overall survival across multiple cancers. To evaluate these interactions, the ESTIMATE algorithm in R was applied to calculate immune, stromal, and ESTIMATE scores and to correlate them with the expression of CAMK2A, CHMP4C, and PYGB in LUSC. CAMK2A expression correlated positively with stromal score (R = 0.39, P = 4.7e−19), immune score (R = 0.15, P = 6.4e−4), and ESTIMATE score (R = 0.28, P = 3.2e−10), suggesting its involvement in promoting stromal and immune activity. In contrast, CHMP4C expression exhibited negative correlations with all three metrics, including immune score (R = −0.25, P = 1.5e−8), stromal score (R = −0.24, P = 8.0e−8), and ESTIMATE score (R = −0.24, P = 1.2e−5), reflecting an opposing relationship with stromal and immune components. PYGB displayed weaker associations, with a minor negative correlation with the immune score (R = −0.15, P = 1.2e−3) and no statistically significant associations with stromal or ESTIMATE scores (P = 0.82 and P = 0.09, respectively).

Given the strong association of the three NRGs with immune infiltration, their immunological relevance in LUSC was further examined using TCGA data and single-sample Gene Set Enrichment Analysis (ssGSEA) analysis (Fig. 10B). CAMK2A displayed an immunostimulatory pattern, showing positive correlations with NK cells (R = 0.470, P < 0.001), macrophages (R = 0.344, P < 0.001), Th1 cells (R = 0.327, P < 0.001), Tem cells (R = 0.295, P < 0.001), along with fifteen additional immune cell types. In contrast, CHMP4C exhibited an immunosuppressive profile, with negative correlations involving Tregs (R = −0.283, P < 0.001), B cells (R = −0.255, P < 0.001), and pDCs (R = −0.255, P < 0.001). PYGB showed a dual regulatory influence, positively correlated with NK cells (R = 0.155, P < 0.001), but negatively correlated with helper T cells (R = −0.249, P < 0.001), cytotoxic cells (R = −0.223, P < 0.001), and CD8⁺ T cells (R = −0.234, P < 0.001). Collectively, the analysis revealed distinct immunomodulatory patterns, with CAMK2A associated with pro-inflammatory responses, PYGB reflecting mixed regulatory activity, and CHMP4C predominantly linked to immunosuppression. The robustness of these results was further confirmed through CIBERSORT, which validated the immune landscape characterized by ssGSEA and demonstrated concordance across analytical methods (Fig. 10C). Continued investigation into the correlation between these NRGs and immune infiltration in LUSC is warranted to consolidate their functional associations with tumor-infiltrating immune cells.

Immune-related genes TMB and MSI analysis of the three NRGs

To delineate the relationship between NRGs and immune infiltration in LUSC, mRNA expression patterns of CAMK2A, CHMP4C, and PYGB were systematically correlated with immune-related gene sets, including chemokines, chemokine receptors, major histocompatibility complex (MHC) molecules, immunoinhibitors, and immunostimulators, across 32 TCGA cancer types (Fig. 11A). CAMK2A displayed broad positive correlations, with distinct clusters of strong associations in several cancers, whereas CHMP4C and PYGB exhibited heterogeneous profiles characterized by both positive and negative associations. CAMK2A was most closely linked with chemokine receptors such as CXCR4 and CCR5, while CHMP4C and PYGB were more strongly associated with MHC class I/II components and immunoinhibitory molecules, suggesting context-dependent functions in immune cell trafficking and antigen presentation.

Subsequent evaluation of the three NRGs in relation to immune checkpoint genes in LUSC (Fig. 11B) revealed that CAMK2A expression correlated positively with SIGLEC15, CTLA4, HAVCR2, IGSF8, PDCD1, PDCD1LG2, and TIGIT. CHMP4C was positively associated with IGSF8 but negatively correlated with CD274, CTLA4, HAVCR2, ITPRIP1, LAG3, PDCD1, PDCD1LG2, and TIGIT. PYGB demonstrated positive correlation with IGSF8 and negative correlation with CTLA4 and TIGIT. Collectively, these patterns indicate a regulatory involvement of CAMK2A, CHMP4C, and PYGB in immune checkpoint pathways in LUSC. Given the relevance of TMB and MSI as biomarkers for immunotherapy response, their associations with NRGs were further assessed. CAMK2A exhibited a significant correlation with TMB, while CHMP4C and PYGB showed no such association (Fig. 11C). None of the three NRGs demonstrated a correlation with MSI.