$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
This study used only publicly available, de-identified datasets and did not involve direct human or animal experimentation; therefore, additional ethics committee approval and informed consent were not required.
Data download and processing
RNA-sequencing data and corresponding clinical information for lung squamous cell carcinoma (LUSC) were obtained from The Cancer Genome Atlas (TCGA) database through the Genomic Data Commons data portal under the TCGA-LUSC project. The TCGA-LUSC expression matrix used in this study was based on FPKM values. Gene expression values were transformed and normalized before downstream analyses. Clinical variables included age, sex, tumor stage, pathologic TNM stage, grade, survival time, and survival status when available. A total of 489 TCGA-LUSC cases were initially retrieved, and 381 patients with complete expression profiles and overall survival information were included in prognostic model construction and internal evaluation.
Independent validation datasets were downloaded from the Gene Expression Omnibus (GEO) database on January 3, 2026. GSE30219 was based on the GPL570 platform and included 307 LUSC patients with available expression and overall survival information. GSE37745 was also based on the GPL570 platform and included 196 LUSC patients with available expression and overall survival information. GSE57148 was based on the GPL11154 platform and included 91 normal lung tissues and 98 lung tissues from patients with chronic obstructive pulmonary disease (COPD), totaling 189 samples. GSE57148 was used to identify COPD-associated differentially expressed genes, whereas GSE30219 and GSE37745 were used as independent external validation cohorts.
For GEO datasets, probe annotation was performed using the corresponding R annotation package and platform annotation files. Probe identifiers were converted to official gene symbols. When multiple probes mapped to the same gene, the probe with the highest average expression value was retained to represent that gene. No missing model genes were detected in the validation datasets after gene-symbol matching. TCGA-LUSC, GSE30219, and GSE37745 were analyzed as general LUSC cohorts because patient-level COPD comorbidity status was not confirmed in the annotations used for this analysis.
Because TCGA and GEO datasets were generated using different expression platforms, cross-platform normalization and batch-effect correction were performed using standard R-based preprocessing methods before model application. The risk model was trained in the TCGA-LUSC cohort and then independently evaluated in each external GEO cohort rather than by directly merging all cohorts. Differential expression analysis was performed using the limma R package. COPD-associated differentially expressed genes in GSE57148 were screened using |log2FC| > 0.263 and P < 0.05. This log2FC threshold corresponds to approximately 1.2-fold change and was used as an exploratory screening criterion to retain potentially relevant PANoptosis-associated genes. A total of 277 PANoptosis-associated genes were curated from previously published apoptosis-, pyroptosis-, necroptosis-, and PANoptosis-related studies and are provided in Supplementary Table 1.
Functional enrichment analysis of genes
To clarify the functional implications of the selected COPD-associated PANoptosis genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler R package and the org.Hs.eg.db annotation package. GO biological process, cellular component, and molecular function categories were evaluated. KEGG pathway enrichment analysis was conducted to identify signaling pathways associated with the selected genes. P values were adjusted for multiple testing using the Benjamini-Hochberg false discovery rate method where applicable. Enrichment terms with P < 0.05 were considered statistically significant in this exploratory analysis. Enrichment plots were generated using ggplot2.
Unsupervised clustering analysis of PANoptosis-associated gene expression patterns
To explore molecular heterogeneity associated with PANoptosis-associated gene expression in LUSC, consensus clustering was performed using the ConsensusClusterPlus R package. LUSC samples were clustered according to the expression profiles of survival-associated PANoptosis genes. Hierarchical clustering was applied with Pearson correlation distance. The maximum number of clusters was set to six, and 1,000 resampling iterations were performed to evaluate clustering robustness. The optimal cluster number was determined by evaluating the consensus matrix, cumulative distribution function curve, delta area plot, and biological interpretability of the resulting groups. Based on these criteria, k = 2 was selected for downstream analysis. Survival differences between the two molecular groups were assessed using Kaplan-Meier analysis and the log-rank test.
Analysis of immune microenvironment differences between subtypes
To compare immune microenvironment features between molecular subtypes, immune cell infiltration was estimated using the CIBERSORT deconvolution algorithm with the LM22 leukocyte signature matrix. The analysis was performed in R using the e1071 and preprocessCore packages. CIBERSORT permutation P values were recorded to assess the reliability of deconvolution estimates. Because this study was exploratory and based on retrospective transcriptomic data, immune-cell differences were interpreted as computationally inferred immune-infiltration patterns rather than direct cellular measurements.
The ESTIMATE algorithm was used to calculate stromal score, immune score, ESTIMATE score, and tumor purity for each tumor sample. GSVA was applied to estimate pathway-level enrichment scores based on selected gene sets. Group-wise differences in immune-cell fractions, immune checkpoint genes, HLA family genes, and ESTIMATE-derived scores were evaluated using nonparametric tests. For multiple immune-related comparisons, Benjamini-Hochberg correction was applied where appropriate; analyses reported using nominal P values were interpreted as exploratory. Spearman rank correlation analysis was used to evaluate associations between gene expression and immune-related markers, with both correlation coefficients and P values reported where applicable. HLA transcript differences were interpreted as antigen-presentation-related transcriptional changes rather than direct functional evidence of enhanced antigen-presentation capacity.
Establishing a prognostic signature related to PANoptosis-associated genes
The TCGA-LUSC cohort with complete expression profiles and overall survival information was used for prognostic model construction. Among the 489 initially retrieved TCGA-LUSC cases, 381 patients with complete overall survival data were included in the prognostic analysis. These patients were randomly divided into a training cohort and an internal testing cohort at a 7:3 ratio. Stratified randomization was performed according to survival status to maintain a comparable distribution of survival events between the training and testing cohorts.
In the training cohort, univariate Cox proportional hazards regression was first used to evaluate the association between each candidate PANoptosis-associated gene and overall survival. Genes with P < 0.05 were considered candidate prognostic genes and were subsequently entered into LASSO Cox regression using the glmnet R package. Ten-fold cross-validation was used to select the optimal penalty parameter and reduce overfitting. Based on the LASSO Cox regression coefficients and corresponding normalized gene expression values, an individualized risk score was calculated for each patient using the formula:
Risk score = Σ(coefi × Xi)
where coefi represents the regression coefficient of each selected gene and Xi represents the normalized expression value of the corresponding gene. The final prognostic model contained 12 genes: GSDMD, IL18, NFKBIA, PIK3CA, IL1B, BIRC3, MCL1, PSMB10, LMNA, CFLAR, IL1R1, and AKT3. The complete coefficient-based risk-score equation is provided in Supplementary Table 2.
The median risk score in the training cohort was used as the cutoff to classify patients into high-risk and low-risk groups. The same risk-score formula was applied to the internal testing cohort and to the external validation cohorts GSE30219 and GSE37745. Kaplan-Meier survival analysis, log-rank testing, and time-dependent receiver operating characteristic curve analysis were used to evaluate model performance. Because the validation datasets were generated using microarray platforms and did not contain confirmed COPD comorbidity annotation, external validation was interpreted as retrospective evaluation in independent LUSC cohorts rather than validation in clinically confirmed COPD-comorbid LUSC patients.
Drug sensitivity prediction analysis
Drug sensitivity was estimated using the pRRophetic R package, which predicts drug response from tumor gene expression profiles based on pharmacogenomic reference data from the Genomics of Drug Sensitivity in Cancer database. Predicted half-maximal inhibitory concentration (IC50) values were calculated for each patient sample. Expression matrices were processed according to pRRophetic input requirements, and batch-effect correction was performed using the standard pRRophetic-compatible workflow. Predicted IC50 values were reported on the pRRophetic output scale.
Eight candidate agents, including sorafenib, gefitinib, bleomycin, bosutinib, etoposide, lenalidomide, camptothecin, and methotrexate, were evaluated as an exploratory drug-sensitivity panel. Differences in predicted IC50 values between high- and low-risk groups were compared using the Wilcoxon rank-sum test. These results were interpreted as computational drug-sensitivity estimates rather than measured clinical chemotherapy response or experimentally confirmed drug resistance.
Statistical analysis
All statistical analyses were conducted using R software. Continuous variables between two groups were compared using the Wilcoxon rank-sum test, whereas comparisons among more than two groups were performed using the Kruskal-Wallis test when appropriate. Overall survival was defined as the primary survival endpoint. Kaplan-Meier survival curves were generated to compare survival differences between groups, and statistical significance was assessed using the log-rank test. Univariate and multivariate
Cox proportional hazards regression analyses were used to evaluate prognostic associations between clinical variables, risk group, and overall survival. Variables with clinical relevance or statistical significance in univariate Cox analysis were considered for multivariate Cox regression. The proportional hazards assumption was assessed using Schoenfeld residuals. Missing clinical variables were handled using complete-case analysis for Cox regression and nomogram construction. Collinearity among clinical variables was evaluated before multivariate modeling.
Time-dependent ROC curves were used to assess the predictive performance of the risk model for 1-, 3-, and 5-year overall survival. A nomogram was constructed using variables retained in the multivariate model or variables with sufficient clinical availability. Calibration plots were used to compare predicted and observed overall survival probabilities. Decision curve analysis was performed as an exploratory assessment of potential net benefit across selected threshold probabilities.
Spearman rank correlation analysis was used to evaluate associations between gene expression and immune-related features. Correlation coefficients and P values were reported where applicable. For multiple comparisons, Benjamini-Hochberg false discovery rate correction was applied where appropriate. Analyses reported using nominal P values were considered exploratory. A two-sided P value < 0.05 was considered statistically significant.