Patients who met the following inclusion criteria were retrospectively analyzed: i) Patients diagnosed with LUSC according to World Health Organization criteria (15); ii) those who underwent next-generation sequencing (NGS) analysis at West China Hospital (Sichuan, China) between January 2018 and December 2020 (identification period). The observation period was measured from treatment initiation until radiographically confirmed disease progression or the administrative censoring date of January 2023; iii) aged >18 years; and iv) with tumor cell content ≥20% in formalin-fixed, paraffin-embedded (FFPE) tissue sections (3-µm thickness), as determined by an experienced pathologist on H&E-stained slides. No additional exclusion criteria were applied. Additionally, a cohort of patients with LUSC who received first-line chemotherapy were also retrospectively included to evaluate the performance of the genomic signature developed in the present study. All samples were analyzed using a 56-gene panel (LungCore; Burning Rock Biotech, Ltd). This panel was part of standard routine care and covered the entire exon regions of 56 lung cancer-related genes which was previously described in our published article (16). This was a retrospective study, which utilized existing clinical records, and all patient data and samples were rigorously de-identified prior to analysis. The full study protocol, which included a waiver of the requirement for informed consent, was reviewed and approved by the Ethics Committee of West China Hospital (Sichuan, China) and conducted according to the local ethical guidelines (approval no. 2022-1849).

Genomic DNA was extracted from FFPE tumor tissues using QIAsymphony^® DSP DNA Mini Kit (cat. no. 937236; Qiagen GmbH), according to the manufacturer's instructions. A total of 200 ng of DNA was used for the preparation of the NGS library. DNA integrity was verified by 1.5% agarose gel electrophoresis. The concentration of DNA was quantified using the Qubit dsDNA HS Assay Kit on a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Inc.). DNA was fragmented using the Covaris M220 system (Covaris, LLC). The fragmented DNA was end-repaired and 3′-adenylated in a combined reaction, followed by adapter ligation using T4 DNA ligase (Burning Rock Biotech, Ltd.). Both enzymatic reactions were carried out in PCR tubes with thermal cycling under a heated lid (85°C). DNA fragments ranging from 200 to 400 bp were purified using the Agencourt AMPure XP Kit (Beckman Coulter, Inc.), followed by hybridization with capture probes; specifically, 50 ng/µl of capture probe baits were co-incubated with DNA fragments in a PCR instrument at 65°C for 16–24 h. This was followed by hybrid selection with magnetic beads and PCR amplification with DNA polymerase (Burning Rock Biotech, Ltd.). PCR amplification was performed with an annealing temperature of 60°C for 12 cycles with Illumina-provided i5/i7 primers. Library quality was assessed using the Qubit dsDNA HS Assay Kit (cat. no. Q32854; Invitrogen; Thermo Fisher Scientific, Inc.) and the Agilent 2100 Bioanalyzer System (Agilent Technologies, Inc.). Indexed libraries were primarily sequenced on the Miseq (Illumina, Inc.) platform. A subset of 66 libraries was sequenced on the MiniSeq (cat. no. SY-410-1003; Illumina, Inc.) platform due to instrument availability. The loading concentration of the final library was 14 pM for Miseq and 1.5 pM for MiniSeq. All sequencing was performed using 150-bp paired-end reads.

Sequencing data were mapped to the human reference genome (hg19) using Burrows-Wheeler aligner 0.7.10 (https://sourceforge.net/projects/bio-bwa/files/) (17). Local realignment and variant calling were performed using the Genome Analysis Toolkit version 3.2 (18) and VarScan software 2.4.3 (Genome Institute at Washington University) (19). Variants were filtered using the VarScan filter pipeline, requiring a minimum depth of 100×. From FFPE tumor tissue, at least five supporting reads for short insertions and deletions (indels) and eight for single nucleotide variations (SNVs) were required. Variants with a population frequency of >0.1% were considered as common single-nucleotide polymorphisms and excluded according to the Exome Aggregation Consortium (https://gnomad.broadinstitute.org), 1000 Genomes Project (https://www.internationalgenome.org/), dbSNP (https://www.ncbi.nlm.nih.gov/snp/) and ESP6500SI–V2 (https://genome.ucsc.edu/cgi-bin/hgTables?db=hg19&hgta_group=varRep&hgta_track=evsEsp6500&hgta_table=evsEsp6500&hgta_doSchema=describe+table+schema) databases. The remaining variants were annotated using ANNOVAR (2016-02-01 release) (20) and SnpEff v3.6 (21). DNA translocation was detected using Factera v.1.4.3 (22). Indels, copy number variations (CNVs) and large genomic rearrangements (LGRs) were identified as previously described (23,24).

Data clustering was performed using NMF based on Euclidean distance (25,26). Non-synonymous SNVs, indels, CNVs and structural variations were included in the NMF clustering analysis. Patients harboring only synonymous SNVs or undetectable genomic alterations were excluded from the analysis. The R package 'NMF' (version 0.22.0) was implemented in R (version 3.4.0; R Foundation for Statistical Computing) to estimate the optimal factorization rank using Lee and Seung's algorithm with 2:8 ranks (27). A marked decrease in the cophenetic correlation coefficient was observed between ranks 4 and 5 (28), therefore rank 4 was selected as the optimal value, thus resulting in four subgroups.

Sequencing data from patients with LUSC were downloaded from The Cancer Genome Atlas (TCGA; http://portal.gdc.cancer.gov/) database using the R package ‘TCGAbiolinks’. The clinical characteristics of patients, including age, sex tumor stage, progression-free survival (PFS) and patient outcomes were also collected. The prognostic value of the NMF-based model was validated in the TCGA LUSC cohort.

All analyses were performed using R statistics packages (R v3.4.0; Posit Software). The differences between the two groups were compared using Fisher's exact test for categorical variables. For comparisons involving more than two groups, the Kruskal-Wallis test was first applied. When the overall test was significant, pairwise post-hoc comparisons were performed using the Wilcoxon rank-sum test with Holm correction to control for multiple testing. PFS was assessed using the Kaplan-Meier method, and differences between survival curves were compared using the log-rank test. P<0.05 was considered to indicate a statistically significant difference.

A total of 317 patients with LUSC were included in the present study. The majority of patients were men (89.9%; 285/317), aged ≥60 years (67.2%; 213/317) and were former or current smokers (57.4%; 182/317). Histologically, 55.2% (175/317) of tumors were classified as keratinizing squamous cell carcinoma (KSCC), 28.7% (91/317) as non-KSCC (NKSCC), while 16.1% (51/317) were of unknown subtype. The distribution of tumor stage varied, with 3.2% (10/317) of patients diagnosed with stage I, 6.6% (21/317) with stage II, 28.7% (91/317) with stage III, 27.8% (88/317) with stage IV, while 33.8% (107/317) of cases were of unknown stage. Among the 317 patients, 35 received first-line chemotherapy. Of these, 91.4% (32/35) were men and 88.6% (31/35) were former or current smokers. Histologically, 65.7% (23/35) were diagnosed with KSCC, 25.7% (9/35) with NKSCC and 8.6% (3/35) were of unknown subtype. In terms of disease stage, 5.7% (2/35) had stage I disease, 25.7% (9/35) stage II, 34.3% (12/35) stage III, 28.6% (10/35) stage IV and 5.7% (2/35) had unknown stage. The detailed clinical characteristics of patients are presented in Table I. The two cohorts showed no statistically significant differences in terms of age (P=0.101), sex (P=1.000), subtype (P=0.557), T stage (P=0.948), N stage (P=0.180) and M stage (P=0.341). However, marginally significant differences were observed in terms of smoking history (P=0.048) and tumor stage (P=0.049). In addition, the clinical characteristics of patients from TCGA database and those who received radiochemotherapy are shown in Tables SI and SII.

Genomic alterations, including non-synonymous SNVs, indels, CNVs, gene fusions and LGRs, were detected in 98% (311/317) of patients (Fig. 1A). TP53 was the most commonly altered gene (91%; 287/317) followed by PIK3 catalytic subunit (PIK3CA; 46%; 145/317), cyclin dependent kinase inhibitor 2A (CDKN2A; 24%; 77/317) and EGFR (17%; 53/317). CNVs were identified in 209 cases, with the most prevalent CNVs detected in PIK3CA (38%; 120/317), fibroblast growth factor (FGF) receptor 1 (FGFR1; 15%; 48/317) and genes located at chromosome 11q13 (11%; 35/317), including cyclin D1 (CCND1; 11%; 35/317), FGF3 (9%; 28/317), FGF4 (8%; 25/317) and FGF19 (10%; 33/317). SNVs and indels were detected in 305 patients, and the highest frequency was recorded in TP53 (90%; 285/317), followed by CDKN2A (22%; 71/317), PIK3CA (14%; 45/317), PTEN (14%; 43/317), EGFR (11%; 36/317) and Erb-B2 receptor tyrosine kinase 4 (ERBB4; 10%; 32/317). Additionally, nine fusion events were identified, including EML4-ALK fusions in two patients, and SEC61G-EGFR, CASC21-MYC, FGFR1-ADAM9, FGFR3-TACC3, C12orf66-KIT and FGFR3-TACC3, FGFR1-IGFBPL1 fusions in one patient each. Two patients carried LGRs, including RB transcriptional corepressor 1 (RB1) in one case and TP53 in the other.

NMF-based clustering was performed to classify patients according to their genomic profile. Excluded from the analysis were two patients harboring synonymous SNVs and six with undetectable genomic alterations. Therefore, NMF analysis was conducted for a total of 309 patients with LUSC. These patients were classified into the following four distinct clusters (Fig. 1A): i) Cluster 1 (C1) consisted of 129 patients characterized by TP53 alterations; ii) Cluster 2 (C2) included 116 patients with PIK3CA amp; iii) Cluster 3 (C3) composed of 33 patients characterized by amp of CCND1, FGF19, FGF3 and FGF4; and iv) Cluster 4 (C4) included 31 patients characterized by amp of FGFR1, kinase insert domain receptor (KDR), KIT proto-oncogene (KIT) and platelet derived growth factor receptor α (PDGFRA). Among the four clusters, C3 and C1 exhibited the highest and lowest CNV frequency, respectively (Fig. 1B).

The association between NMF-based molecular clusters and survival outcomes was subsequently investigated. Among the 317 patients with LUSC, 35 patients who received first-line chemotherapy and had available PFS data were included in a retrospective cohort. A single patient was excluded due to the undetectable genetic alterations, leaving 34 patients for NMF analysis. Of the aforementioned patients, 13, 10, 5 and 6 were allocated into the C1, C2, C3 and C4 clusters, respectively (Fig. 2A). No statistically significant difference in PFS was observed among the four clusters (P=0.12; Fig. 2B). Pairwise comparisons further confirmed that PFS did not differ significantly among C2, C3 and C4 (all P>0.5; Fig. 2B). Given their comparable chemotherapy outcomes and the shared feature of gene amplification, C2, C3 and C4 were therefore merged in a data-driven manner for subsequent analyses. Notably, patients in C1 exhibited significantly shorter PFS compared with those in C2/3/4 [3.5 vs. 12.5 months; P=0.018; hazard ratio (HR) =0.35; 95% confidence interval (CI): 0.14–0.87; Fig. 2C). Additionally, tumor stage distribution was similar among the clusters, thus indicating a comparable background in tumor stage. (P=0.415; Fig. S1A; and P=0.293; Fig. S1B).

The LOF alterations in TP53 are known contributors to chemotherapy resistance in several types of cancer (29). In the present study, the prognostic effect of TP53 LOF in 35 patients with LUSC treated with first-line chemotherapy was explored. In the present chemotherapy cohort, all 29 patients (83%) harboring TP53 mutations were annotated as LOF according to the OncoKB database (https://www.oncokb.org/gene/TP53#tab=FDA). However, patients with TP53 mutations showed no significant difference in PFS compared with TP53-wild-type patients (P=0.835; Fig. S2). Considering that truncating alterations (frameshift, splice-site, nonsense mutations or copy number deletions) are expected to result in complete functional loss through nonsense-mediated mRNA decay or production of severely truncated proteins, whereas missense mutations may exert heterogeneous effects including partial LOF, dominant-negative or gain-of-function (30,31), TP53 LOF was defined in the present study exclusively as truncating alterations and did not include missense mutations in the LOF category. Using this definition, the results revealed that patients with TP53 LOF (n=15) exhibited significantly shorter median PFS (mPFS) compared with those without TP53 LOF (4.8 vs. 17 months; P=0.029; HR=0.37; 95% CI: 0.15–0.93; Fig. 3A). No significant difference was observed in the distributions of TP53 LOF across all clusters (P=0.389; Fig. S3A), nor specifically among C2, C3 and C4 clusters (P=0.288) and the distribution of TP53 LOF was also similar between C1 and C2/3/4 clusters (P=0.728; Fig. S3B). A comparison of TP53 background can be performed across these clusters. Furthermore, the combined predictive value of TP53 LOF and NMF-based molecular clustering in patients with LUSC was evaluated. The results demonstrated that patients in C1 with TP53 LOF exhibited worse prognosis (mPFS=2.4 months) compared with C1 without TP53 LOF (mPFS=7.1 months), C2/3/4 with TP53 LOF (mPFS=5.3 months) and C2/3/4 without TP53 LOF (mPFS not reached; Fig. 3B). Additionally, to refine the clustering approach, the status of TP53 alterations and gene amps were considered, thus resulting in the classification of 34 patients into four subtypes. Notably, patients without TP53 LOF, but with amp of at least one of the nine specified genes (PIK3CA, CCND1, FGF19, FGF3, FGF4, FGFR1, KDR, KIT and PDGFRA) displayed the longest mPFS (not reached) compared with those with TP53 LOF but without Amp(9G), who exhibited the shortest PFS (mPFS=2.6 months; Fig. 3C). Consistent results were observed when patients with Amp(9G) and TP53 LOF were merged with those without Amp(9G) and TP53 LOF into a single subgroup (Fig. 3D). Multivariable Cox analysis demonstrated that the Amp(9G) combined with TP53 LOF-based grouping was independently associated with PFS after adjustment for tumor stage (Fig. S4).

To validate the aforementioned findings, DNA sequencing data and clinicopathological characteristics of patients with LUSC were downloaded from TCGA database. The analysis included only patients who received chemotherapy, resulting in a validation cohort of 99 patients (TCGA-LUSC cohort). The results showed that TP53 LOF was associated with shorter PFS compared without TP53 LOF (19.0 vs. 61.6 months; P=0.014; HR=0.48; 95% CI: 0.26–0.87; Fig. 4A). NMF analysis further stratified the TCGA-LUSC cohort into the same four clusters (Fig. S5). Furthermore, evaluation of the combined TP53 LOF and Amp(9G) revealed that patients with Amp(9G) but without TP53 LOF had favorable mPFS (61.6 months), while those with TP53 LOF and no Amp(9G) had worse mPFS (18.7 months; Fig. 4B). The Amp(9G) and TP53 LOF signature remained an independent prognostic factor for PFS in multivariable analysis after controlling for tumor stage (Fig. S6).

To further investigate the feasibility of this signature in predicting the efficacy of first-line chemotherapy in LUSC, a total of 24 patients with LUSC treated with first-line radiochemotherapy were respectively analyzed. These patients were allocated into three groups based on the Amp(9G) and TP53 LOF status. Therefore, TP53 LOF combined with Amp(9G) was not associated with PFS (P=0.419; Fig. S7), indicating that its prognostic value could be specific to patients receiving first-line chemotherapy alone, and not to those treated with concurrent first-line radiochemotherapy.