Fig. 1 illustrates the workflow of the bioinformatics analysis. Seven microarray datasets related to LUSC, including GSE19188, GSE2088, GSE30219, GSE33479, GSE33532, GSE21933, and GSE8569, were sourced from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), with detailed dataset information provided in Supplementary Table 1. The probe IDs in the series matrix files were annotated to gene symbols according to the corresponding platform files. For genes with multiple probes, the expression values were averaged and normalized using the “limma” package (Version 3.56.2, Bioconductor, USA, https://www.bioconductor.org/) in R software (Version 4.3.1, R Foundation for Statistical Computing, Vienna, Austria) to minimize discrepancies across datasets. Among these datasets, GSE19188, GSE2088, GSE30219, and GSE33479, comprising 150 patients with LUSC and 122 normal controls, were designated as the training dataset. External validation was performed using three additional GEO datasets: GSE33532, GSE21933, and GSE8569.

Fig. 1.

Identified differentially expressed genes (DEGs) between normal and lung squamous cell carcinoma (LUSC) individuals from four Gene Expression Omnibus (GEO) datasets. (A) Principal component analysis (PCA) of four GEO datasets before batch correction. (B) PCA of four GEO datasets after batch correction. (C) Volcano plot illustrating DEGs between normal and LUSC patients. (D) Heatmap depicting DEGs between normal and LUSC patients.

To integrate the datasets (GSE19188, GSE2088, GSE30219, and GSE33479) into a unified training set, batch effects were corrected using the Combat function from the “sva” package (Version 3.48.0, Bioconductor, USA). DEGs between patients with LUSC and normal controls were identified using the “limma” package. Genes with a $|$log2FC$|$ greater than 0.585 and an adjusted p-value below 0.05 were considered significantly differentially expressed. Volcano plots and heatmaps were generated to visualize the DEGs between patients with LUSC and normal controls.

The LM22 signature matrix, which includes 22 types of immune cells, was employed to estimate the relative abundances of infiltrating immune cells for each sample in the merged training dataset through the CIBERSORT algorithm (https://cibersort.stanford.edu). The total proportion of the 22 immune cell types in each sample was normalized to 1. The type and distribution of infiltrating immune cells were analyzed using the “ggplot2” package (Version 3.4.4, R Foundation for Statistical Computing, Vienna, Austria) in R software and presented as barplots and boxplots. Significant differences in the relative abundances of immune cell types between patients with LUSC and normal controls were assessed using the Wilcoxon test, with a p-value 0.05. Spearman correlation analysis of the 22 immune cell types in 150 patients with LUSC was conducted using the “corrplot” package (Version 0.92, R Foundation for Statistical Computing, Vienna, Austria), with significance determined at a p-value 0.05.

To identify the significant co-expression module associated with macrophages, the top 50% of genes exhibiting the highest variance were selected for analysis using the WGCNA package (Version 1.72.1, Bioconductor, USA), ensuring robust results [21, 22]. The optimal soft-thresholding power was determined using a weighted adjacency matrix, which was then converted into a topological overlap matrix (TOM). A minimum module size of 50 genes was set to ensure biologically relevant modules. Modules were assigned distinct colors based on their TOM dissimilarity measure, calculated as 1 minus the TOM value, using hierarchical clustering. The eigengenes of each module were used to represent the overall gene expression profiles within the module. The significance of each module was evaluated for its association with corresponding immune cells. Additionally, gene significance was assessed by evaluating the strength of the relationship between individual genes and immune cells. Higher gene significance values indicated stronger correlations between specific genes and immune cells [22]. Ultimately, 567 key genes in the blue module associated with M1 macrophages were identified.

The R package “VennDiagram” was used to identify 142 overlapping genes between DEGs and blue module-associated genes. To investigate the biological functions of these genes, GO functional enrichment and KEGG pathway analyses were conducted using the R packages “clusterProfiler” (Version 4.8.3, Bioconductor, USA), “org.Hs.eg.db” (Version 3.17.0, Bioconductor, USA), and “enrichplot” (Version 1.20.3, Bioconductor, USA). Statistical significance was determined with an adjusted p-value threshold of less than 0.05. Results were visualized using bar plots and circle plots, generated with the “ggplot2” (Version 3.4.4, R Foundation for Statistical Computing, Vienna, Austria) and “circlize” (Version 0.4.15, R Foundation for Statistical Computing, Vienna, Austria) packages in R, respectively.

The 142 overlapping genes were then used as input features for four machine learning models: SVM, RF, GLM, and XGB, aimed at identifying significant predictive genes associated with LUSC. The 150 LUSC samples were randomly split into two datasets: a training dataset (60%, N = 90) and an internal validation dataset (40%, N = 60). These models were executed with default parameters, and the “caret” R package (Version 6.0.94, R Foundation for Statistical Computing, Vienna, Austria) was employed to optimize model parameters via grid search. The “DALEX” package (Version 2.4.3, R Foundation for Statistical Computing, Vienna, Austria) was subsequently used to analyze residual distributions and assess feature importance for the four models [23]. To evaluate the predictive accuracy of the models, the internal validation dataset was used for area under the curve (AUC) analysis via 5-fold cross-validation. Higher AUC values indicate superior model performance [24]. The most effective machine learning model was identified, and the top five variables from this model were determined as key predictive genes linked to LUSC. Finally, receiver operating characteristic (ROC) curves were generated using three external datasets (GSE33532, GSE21933, and GSE8569) to validate the diagnostic value of the signature derived from the best-performing machine learning model.

A nomogram model was constructed using the “rms” package (Version 6.6.0, R Foundation for Statistical Computing, Vienna, Austria) to predict the likelihood of LUSC. Scores were assigned to five key genes, and their sum formed the “total score”, which was used as an indicator of LUSC risk severity. Calibration curve analysis and DCA were also performed to assess the predictive accuracy and clinical utility of the nomogram model.

The discriminative ability of the predictive model to distinguish between patients with LUSC and healthy controls was validated using three independent external datasets (GSE33532, GSE21933, and GSE8569). ROC analysis was conducted, and the “pROC” R package (Version 1.18.5, R Foundation for Statistical Computing, Vienna, Austria) was employed to generate and visualize the ROC curves.

Forty-eight pairs of LUSC tissue samples were collected from Tianjin Medical University Cancer Institute and Hospital, with approval from the Ethics Committee (No. bc2023152). Informed consent was obtained from all subjects involved in the study. The study was carried out in accordance with the ethical guidelines of the Declaration of Helsinki. Immunohistochemical (IHC) analysis was performed on a cohort of 48 individuals diagnosed with LUSC following the manufacturer’s protocol (#PV-9000, ZSGB-BIO, Beijing, China). Tissue sections were dewaxed and hydrated using a gradient ethanol approach. After antigen retrieval in citric acid buffer (pH 6.0) with microwave heating, endogenous peroxidase was blocked for 10 minutes. Sections were then incubated overnight at 4 °C with primary antibodies: ZIC2 rabbit polyclonal antibody (1:100, HA500195, Huabio, Hangzhou, Zhejiang, China), RTN1 rabbit polyclonal antibody (1:50, 15048-1-AP, Proteintech, Wuhan, Hubei, China), MMP7 rabbit polyclonal antibody (1:50, 10374-2-AP, Proteintech, Wuhan, Hubei, China), KLRB1 mouse monoclonal antibody (1:500, 67537-1-lg, Proteintech, Wuhan, Hubei, China), and CXCL13 rabbit monoclonal antibody (1:50, HA722117, Huabio, Hangzhou, Zhejiang, China). After washing, sections were treated with a reaction enhancer and enzyme-linked secondary antibody (sheep anti-rabbit/mouse IgG) (#PV-9000, ZSGB-BIO, Beijing, China). Staining was visualized using DAB chromogen, and counterstaining was performed with hematoxylin. Digital images at $\times$20 magnification were obtained using a BX51 microscope (OLYMPUS, Tokyo, Japan). The evaluation of ZIC2, RTN1, MMP7, KLRB1, and CXCL13 staining followed specific scoring criteria: a positive stained area 5% was scored 0; 5%–25% received 1 point; 26%–50% equaled 2 points; 51%–75% was given 3 points; and areas exceeding 75% were assigned 4 points. For color intensity, no staining was scored 0 points; yellow staining received 1 point; brown staining was assigned 2 points; and tan staining received 3 points. The immunostaining score was calculated by multiplying the color intensity and positive stained area scores.

Statistical analyses were performed using R software (Version 4.3.1). For comparisons between two groups, unpaired Student’s t-tests were applied for normally distributed variables, and Wilcoxon rank-sum tests were used for non-normally distributed variables. Spearman correlation analysis was conducted to estimate correlation coefficients. Statistical significance was defined as p 0.05 in two-tailed tests, unless otherwise specified.

To identify differential expression between lung squamous cell carcinoma (LUSC) and normal samples, four GEO datasets (GSE19188, GSE2088, GSE30219, and GSE33479) underwent batch correction and were subsequently merged into a large sample pool to form the training cohort. Principal component analysis (PCA) and boxplot analyses confirmed that batch effects were effectively minimized across the four original datasets (Fig. 1A,B; Supplementary Fig. 1). Differential expression analysis identified a total of 1032 DEGs, including 466 upregulated and 566 downregulated genes, between 122 control samples and 150 LUSC samples (Fig. 1C; Supplementary Table 2). Heatmaps depicted the top 20 genes with the most significant expression changes, either upregulated or downregulated (Fig. 1D).

Next, we explored the immune landscape in LUSC by assessing the abundance of 22 immune cell types using the CIBERSORT algorithm (Fig. 2A). Significant differences between control and LUSC samples were analyzed using the Wilcoxon rank-sum test. Boxplots and statistical results revealed that T cells CD4 memory resting, monocytes, mast cells resting, T cells CD4 memory activated, and M1 macrophages were the top five immune cell types showing significant differences between control and LUSC groups (Fig. 2B; Supplementary Table 3). Correlation analyses of the 22 types of infiltrating immune cells in patients with LUSC are shown in Fig. 2C.

Fig. 2.

Analysis of immune characteristics in LUSC. (A) Barplot showing the distribution of 22 types of infiltrating immune cells in control and LUSC samples. (B) Differential expression analysis of 22 types of infiltrating immune cells between control and LUSC samples. (C) Correlation analysis among 22 types of infiltrating immune cells in LUSC patients. *p 0.05 , **p 0.01 and ***p 0.001.

To gain insight into the genes associated with M1 macrophages, the WGCNA algorithm was applied to identify gene modules strongly associated with immune cell types. The optimal soft-thresholding parameters were determined as $\beta$ = 3 and R² = 0.9 for constructing a scale-free network (Fig. 3A). Using the dynamic tree-cutting algorithm, seven distinct co-expression modules, each represented by a unique color, were generated (Fig. 3B). Genes from these modules were then evaluated for their co-expression and adjacency relative to the 22 types of immune cells. By integrating the results of WGCNA and immune cell differential analysis, the blue module, consisting of 567 genes, was found to have a significantly stronger association with M1 macrophages (R = 0.37, p = 7.8 $\times$ 10^-20), and was selected for further investigation. These genes within the blue module were positively correlated with M1 macrophages (Fig. 3C,D).

Fig. 3.

Identification of M1 macrophage-related genes through WGCNA. (A) The process for selecting the soft threshold power. (B) Cluster dendrogram of genes grouped into co-expression modules, with various modules represented by distinct colors. (C) Correlation analysis between module eigengenes and 22 types of infiltrating immune cells. (D) Scatter plot showing the relationship between module membership in the blue module and gene significance for M1 macrophages.

To identify M1 macrophage-related DEGs, the DEGs were intersected with the genes from the WGCNA blue module. A Venn diagram revealed that 143 M1 macrophage-related DEGs were obtained (Fig. 4A; Supplementary Table 4). GO analysis indicated that these genes are involved in leukocyte- or lymphocyte-mediated immunity, antigen processing and presentation of peptide antigens, assembly of major histocompatibility complex (MHC) class I or II protein complexes, regulation of immune response signaling pathways, and positive regulation of immune cell-cell adhesion (Fig. 4B; Supplementary Table 5). KEGG pathway analysis revealed that these genes are associated with signaling pathways in rheumatoid arthritis, asthma, Fc epsilon RI signaling, chemokine signaling, and Tumor Necrosis Factor Signaling (TNF) signaling (Fig. 4C; Supplementary Table 6).

Fig. 4.

Biological functional analysis of differentially expressed M1 macrophage-related genes. (A) Venn diagram illustrating the overlapping and unique differentially expressed M1 macrophage-related genes. (B) Circular visualization of Gene Ontology (GO) analysis highlighting the biological functions associated with differentially expressed M1 macrophage-related genes. (C) Barplot depicting the signaling pathways enriched in the differentially expressed M1 macrophage-related genes.

To explore the potential of these M1 macrophage-related DEGs for LUSC diagnosis, four machine learning models—SVM, XGB, RF, and GLM—were constructed based on the expression profiles of 143 M1 macrophage-related DEGs in the training dataset to identify diagnostic genes for LUSC. Residual distribution analysis revealed that SVM and XGB exhibited smaller residuals compared to RF and GLM (Fig. 5A,B). The top 10 feature variables for each model were ranked by their root mean square error (RMSE), and their corresponding feature importance was visualized in Fig. 5C. Additionally, the performance of each machine learning model was evaluated using ROC curves generated through five-fold cross-validation on the testing datasets. The results showed that the SVM model achieved the highest AUC value (AUC_SVM = 0.995, AUC_RF = 0.988, AUC_XGB = 0.983, and AUC_GLM = 0.622) (Fig. 5D). These findings demonstrate that the SVM model was the most effective for distinguishing patients with LUSC. Consequently, MMP7, RTN1, ZIC2, KLRB1, and CXCL13 were identified as novel predictive genes (5-gene signature) for LUSC due to their high significance in the SVM model.

Fig. 5.

Development and evaluation of the machine learning models based on 143 M1 macrophage-related genes. (A) The cumulative residual distribution of Support Vector Machine (SVM), Extreme Gradient Boosting (XGB), Random Forest (RF), and Generalized Linear Model (GLM) models. (B) The residual boxplots of the four machine learning models. (C) The important features of the four machine learning models. (D) Receiver operating characteristic (ROC) analysis of these models using 5-fold cross-validation in the internal testing dataset.

To assess the predictive power of the 5 hub genes, a nomogram was developed to estimate the risk associated with LUSC (Fig. 6A). The calibration curve showed excellent agreement between actual and predicted risks (Fig. 6B). DCA further validated the high accuracy and clinical utility of the nomogram (Fig. 6C). The predictive performance of the 5 hub genes was also evaluated in three independent external testing datasets. The AUC values for GSE33532, GSE21933, and GSE8569 were all 1.000, as demonstrated by the ROC curves (Fig. 6D–F). These results confirm that the diagnostic signature developed in this study exhibits robust performance in distinguishing LUSC from normal individuals across both internal and external validation cohorts.

Fig. 6.

Validation of the 5-gene diagnostic signature. (A) Nomogram of 5 hub genes for predicting the risk of LUSC. (B,C) Calibration curve (B) and decision curve analysis (DCA) (C) for validating the predictive accuracy of the nomogram. (D–F) ROC curve of the 5 hub genes in external GSE33532 (D), GSE21933 (E), and GSE8569 (F) datasets.

To validate the predictive performance of M1 macrophage-associated biomarkers, IHC experiments were performed to evaluate the protein expression levels of candidate genes in tissue samples. IHC results show ZIC2 is predominantly expressed in the nuclei of lung cancer cells, with elevated levels in cancer tissues. RTN1 is mainly found in the membrane of lung cancer cells, but its expression is downregulated in these tissues. MMP7, a secreted protein, is present in the nucleus, cytoplasm, and extracellular matrix, with elevated expression in lung cancer tissues. KLRB1, a transmembrane protein, is downregulated in lung cancer tissues. CXCL13 is expressed in the cytoplasm and extracellular matrix, with upregulation in lung cancer tissues. The results showed that the protein expression of those 5 hub genes was largely consistent with their corresponding mRNA expression levels in both the training dataset and the TCGA dataset (Fig. 7A,B; Supplementary Fig. 2; Supplementary Table 2). Furthermore, the predictive efficacy of the five hub genes was assessed based on the IHC results, yielding an AUC value of 0.992 in the ROC curve analysis (Fig. 7C). These results confirm the reliability of the data mining outcomes and highlight their potential for LUSC diagnosis.

Fig. 7.

Immunohistochemistry validation of the predictive performance of M1 macrophage-associated biomarkers. (A) Representative images of immunohistochemical staining for Zinc finger protein ZIC 2 (ZIC2), Reticulon-1 (RTN1), Matrix metalloproteinase-7 (MMP7), Killer cell lectin-like receptor subfamily B member 1 (KLRB1) and C-X-C motif chemokine 13 (CXCL13) in normal lung tissues and LUSC tissues (scale bar, 50 µm). (B) Statistical analysis of immunohistochemistry results in normal lung tissues and LUSC tissues. Data are mean $\pm$ SD. n = 48. *p 0.05, and ***p 0.001. (C) ROC curve of the 5-gene signature in the immunohistochemical (IHC) data.