All the raw data of series matrix files and clinical information from the GSE41271 (19), GSE42127 (20) and GSE32863 (21) datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/; accessed on 3 June 2022). The three matrix files originated from the GPL6884 Illumina HumanWG-6 v3.0 Expression Bead Chip platform. After a series of normalizations using the InSilicoMerging R package, the batch effect among the three microarray datasets was removed using empirical Bayes methods and the ‘combat’ algorithm. Finally, the three datasets were merged into one big matrix file. Next, samples from patients with squamous carcinoma, missing age of patients, OS of <1 month and information with ambiguity were eliminated. Finally, a total of 155 LUAD samples and 58 normal samples were used for the following analysis of differentially expressed genes (DEGs). A total of 97 LUAD cases had survival information that met the study criterion (GSE41271, 61 cases; GSE42127, 36 cases), which were included for the following WGCNA and survival-related prognostic analysis (Table I). Flow charts of the study's designed strategy for selecting patients are presented in Figs. 1 and 2, respectively.

DEGs in LUAD and normal samples were screened using the limma package (version 3.40.6) in R software. Volcano plots and heat maps were created to visualize DEGs. Genes with an adjusted P<0.05 and |log₂fold change (FC)|>1 were set as the filtration criterion for DEGs.

For survival analysis, the ‘survival’ package of R software (version 4.1) was used to integrate the survival time, survival status and gene expression data in order to further assess the significant prognosis-related genes using univariate regression analysis, and P<0.05 was considered to indicate a statistically significant difference.

The ‘WGCNA’ package was used to construct a WGCNA (18). First, according to the gene expression profiles, the standard deviation (SD) of each gene was calculated, and genes in the top 25% with the smallest SD were removed. The goodSamplesGenes function in the R software package was used to eliminate outlier samples and the WGCNA package (version 3.9) was used to analyze the original data.

To further analyze the module, the dissimilarity of module eigengenes was calculated, a cut line for the module dendrogram was selected and certain modules were merged. In addition, modules with distances of <0.25 were merged and 13 co-expressed modules were obtained. Module-membership (MM) was evaluated as the connectivity between gene expression values, while the gene-significance (GS) represented the correlation between each gene and OS. Genes with larger GS values had a greater influence on traits and genes with larger MM values were more correlated with modules. Genes in the key modules with |MM|>0.8 and |GS|>0.2 were selected as significant genes in the module.

A Venn diagram was used to display the overlapping genes among DEGs, prognosis-related genes and WGCNA. The overlapping genes were considered the candidate genes for further analysis.

For GO functional enrichment analysis, the GO annotations of genes in the R package org.Hs.eg.db (version 3.1.0) was used as the background to map the genes to the background set using the R package Cluster Profiler (version 3.14.3). The minimum gene set was set to 5 and the maximum to 5,000. P<0.05 was considered to indicate a statistically significant difference.

For gene set functional enrichment analysis, KEGG API was used (https://www.kegg.jp/kegg/rest/keggapi.html) to obtain the latest gene annotations. The KEGG pathway enrichment analysis was performed using the R software package Cluster Profiler (version 3.14.3).

The STRING online tool (https://string-db.org/cgi/) is a system that searches for known and predicted protein-protein interactions (PPIs). Such interactions include both direct physical interactions and indirect functional correlations between proteins. The 17 genes overlapping in the Venn diagram were selected to build a PPI network using the STRING database. Next, Cytoscape (version 3.9.1) was used to visualize the network, followed by the Cytohubba plug-in to explore hub genes in the network. The top 5 genes with the most connected traits were identified using the maximum clique centrality algorithm. Subsequently, the five candidate genes and clinical parameters, such as age, sex and stage, were included in the univariate and multivariate Cox regression analysis, and P≤0.05 was considered to indicate a statistically significant difference.

For the purpose of assessing the prognostic value of FOXM1 and CENPF expression in LUAD, the Kaplan-Meier plotter was used to examine the relationship between FOXM1, CENPF and OS in LUAD. Furthermore, a Pearson correlation analysis between FOXM1 and CENPF was conducted in clinical samples from the GEO database. Furthermore, to validate the correlation between FOXM1 and CENPF, an analysis with the GEPIA (http://gepia.cancer-pku.cn/) public web database was performed based on The Cancer Genome Atlas database (22). Based on its larger sample size and credible analysis results, this website was used to validate the correlation between FOXM1 and CENPF.

The GSEA software (version 3.0) was obtained from the GSEA (23). According to gene expression levels, genes were divided into high (≥50%) and low expression groups (<50%). Background GSEA signature gene set expression was obtained from the Molecular Signatures Database (23,24), and the c2.cp.kegg.v7.4.symbols.gmt subset was downloaded to evaluate related pathways and molecular mechanisms. Classification into groups was performed based on gene expression profiles and phenotypes, setting the minimum gene set to 5, the maximum gene set to 5,000 and 1,000 resampling. P<0.05 was considered to indicate a statistically significant difference.

A total of 22 normal adjacent tissues and 22 LUAD tissues (12 males and 10 females; mean age, 62.6 years; age range, 41–82 years), collected from The Second Hospital of Hebei Medical University from 2018 to 2023 (approval no. 2022-R676) were fixed with 4% paraformaldehyde for 24 h and then embedded in paraffin. Samples were cut into 5-µm sections following deparaffinization and dehydration. H&E staining was performed according to standard protocols. IHC was conducted to investigate FOXM1 and CENPF protein expression in LUAD and normal tissues. Tissues were incubated with FOXM1 (cat. no. ab207298; dilution, 1:200; Abcam) and CENPF (cat. no. 28568-1-AP; dilution, 1:200; Proteintech) antibodies overnight at 4°C. The next day, tissues were incubated with a secondary antibody (cat. no. PV-6001; ready-to-use; OriGene Technologies, Inc.) for 1 h at 37°C. The nuclei were then stained with hematoxylin following DAB kit (cat. no. ZLI-9018; OriGene Technologies, Inc.) staining. The results were observed under a bright-field microscope (Nikon Eclipse CI; Nikon Corporation).

For IF, samples were prepared using CENPF (cat. no. 28568-1-AP; dilution, 1:3,000; Proteintech) and FOXM1 (cat. no. ab207298; dilution, 1:200; Abcam) antibodies following the manufacturer's instructions. In the fluorescent images, CENPF appeared green and FOXM1 red. The clinicopathological data of the subjects are listed in Table SI. The results were observed under a bright-field microscope (Nikon Eclipse CI).

In the previous bioinformatics analysis, FOXM1 and CENPF were found to be highly expressed in LUAD as compared with normal tissue, and to be significantly associated with poor prognosis. Therefore, aiming to further investigate the role of deregulated expression, gene expression was knocked down by silencing FOXM1 and CENPF in A549 cells. The A549 human LUAD cell line was purchased from Haixing Biosciences and cultured in DMEM/F-12 (Gibco; Thermo Fisher Scientific, Inc.) medium containing 10% fetal bovine serum (cat. no. F7524; MilliporeSigma) in an incubator with a humidified atmosphere with 5% CO₂ at 37°C. Small inhibitory (si)RNA sequences targeting FOXM1 and CENPF were constructed by and purchased from Guangzhou RiboBio Co., Ltd. A total of 10 µl siRNA (100 nM) was transfected in each well of six-well plates using riboFECT CP Transfection Reagent (Guangzhou RiboBio Co., Ltd.), following the manufacturer's instructions. The sequences of si-NC, si-FOXM1 and si-CENPF (Guangzhou RiboBio Co., Ltd.) are listed in Table II.

Total RNA was isolated from the A549 cells using RNAiso Plus (Takara Bio, Inc.). Subsequently, 1 µg RNA was used for cDNA conversion. cDNA was obtained by RT at 42°C for 15 min and 95°C for 3 min. PCR was conducted with SYBR Green (cat. no. FP205; Tiangen Biotech Co., Ltd.) according to the manufacturer's instructions. The thermocycling conditions for PCR using a CFX96-Real-Time System (Bio-Rad Laboratories, Inc.), were as follows: Initial denaturation at 95°C for 15 min, followed by denaturation at 95°C for 10 sec, and annealing and extension at 60°C for 32 sec for a total of 40 cycles. Furthermore, the expression relative to GAPDH was calculated using the 2^−∆∆Cq method (25). The specific primer sequences are listed in Table III.

A549 cells were planted overnight in 96-well plates at a density of 3,000 cells/well and transfected with si-NC, si-FOXM1 orsi-CENPF the next day. Cells were further incubated for 24, 48 and 72 h in a cell incubator. Subsequently, 20 µl CCK-8 reagent (Wuhan Servicebio Technology Co., Ltd.) was added to each well and incubated at 37°C for 2 h in a cell incubator. The absorbance of each well was measured at a wavelength of 450 nm using a microplate reader.

For the colony-formation assay, 500 cells/well were incubated in a six-well plate overnight and transfected with si-NC, si-FOXM1 orsi-CENPF the next day. After 48 h of incubation, the medium was replaced and cells were incubated for 10 days. When colonies were observed (>50 cells), cells were fixed with 4% paraformaldehyde for 20 min and stained with 0.1% crystal violet for 20 min at room temperature. Finally, colony numbers were evaluated using ImageJ software (version 1.54d; National Institutes of Health).

Transwell assay inserts (pore size, 8 µM; Corning, Inc.) were used to evaluate the migration capability of A549 cells. Cells were harvested following transfection for 24 h. A total of 2×10⁴ cells were transfected with si-NC, si-FOXM1 or si-CENPF and seeded into the upper chamber in 100 µl medium with 2% fetal bovine serum. Furthermore, 600 µl medium containing 20% fetal bovine serum was added to the lower chamber and the plates were incubated for an additional 24 h. Migrated cells were fixed with 4% polyformaldehyde for 30 min and stained with 0.1% crystal violet for 20 min at room temperature. The results were evaluated using ImageJ software (version 1.54d; National Institutes of Health).

A549 cells (3.5×10⁵ cells/well) were seeded in six-well plates and transfected with si-NC, si-FOXM1 or si-CENPF. When the cell confluence had reached 90%, the cell monolayer was scratched with a 10-µl pipette tip to generate a line-shaped wound, then debris cells were washed away with phosphate-buffered saline. The scraped monolayer was incubated in medium containing 1% fetal bovine serum for an additional 48 h. Scratched fields were selected randomly and cell migration distances were further calculated.

The χ² test and non-parametric test were used for count data. Unpaired Student's t-test or ANOVA were used for measurement data and the Student-Newman-Keuls method was used as the post-hoc test. Kaplan-Meier curves were drawn for survival analysis by log-rank test. Univariate and multivariate Cox regression analysis was used to explore the independent risk factors for clinicopathological data and protein expression levels. Statistical analysis was performed using SPSS 26.0 (IBM Corp.). P<0.05 was considered to indicate a statistically significant difference.

Genes were screened using an adjusted P<0.05 and |log2FC|>1, and 1,018 DEGs were identified, including 416 upregulated and 602 downregulated genes (Table SII). The visualization results of the normalization of original data and DEGs are presented in Figs. 3 and 4 (FOXM1, logFC=1.187657, P=7.80×10⁻¹²; CENPF, logFC=1.453464, P=4.70×10⁻¹⁵).

Genes that were significantly associated with prognosis were screened using univariate regression analysis with log-rank P<0.05, likelihood P<0.05 and 95% confidence interval <1 as the criteria. A total of 2,923 genes met the screening criteria. Specific information of prognosis-related genes selected by univariate regression analysis is presented in Table SIII [FOXM1, hazard ratio (HR)=1.461737, 95% CI (1.178551–1.812968), P=0.000528; CENPF, HR=1.260067, 95% CI (1.070011–1.48388), P=0.005448].

Clustering analysis was performed according to the expression matrix and clinical characteristics of GSE41271 and GSE42127. The clinical characteristics of the patients with LUAD are presented in Table SIV. A total of 97 LUAD samples with complete OS data and corresponding series matrix files were used to determine the modules with highly correlated genes by WGCNA. The clinical variables of sex, age, stage and OS were analyzed using WGCNA. A sample dendrogram and trait heatmap are presented in Fig. 5A. A power value of β=30.89 was selected as the soft threshold parameter to build a scale-free network (Fig. 5B), with the mean connectivity at 32,744.67 (Fig. 5C). Similar expression patterns were gathered into the same module, and the modules with a cutting height difference of <0.25 were merged. A total of 13 co-expression modules were explored after preprocessing hierarchical clustering (Fig. 5D). The correlations among 13 modules are presented in the heatmap of Fig. 5E. Next, the hierarchical clustering and adjacency relationships between modules and clinical traits were analyzed. Among the 13 modules, modules with a significance of P<0.05 were selected for further analysis, and the royal blue module exhibited the strongest negative correlation with OS; thus, it was considered the key module for further exploration. The royal blue module included a total of 103 genes and is presented in Table SV. Based on the cut-off criteria (|MM|>0.8 and |GS|>0.2), a total of 16 genes with high connectivity in the clinically significant module were explored (Fig. 5F; Table SVI).

DEGs, prognosis-related and module-trait genes (royal blue module from WGCNA) were overlapped and a total of 17 genes (CDC20, PRC1, CCNF, KIF20A, C19orf48, CENPF, KIF2C, TTK, NCAPG, C16orf59, KIF14, AUNIP, FOXM1, POLQ, CDT1, DSP and CCNB1) were screened out, as indicated in the Venn diagram (Fig. 6A).

GO and KEGG enrichment analysis of 17 genes was performed and the top five sorted by P-value (P<0.05) were as follows. GO: ‘Mitotic sister chromatid segregation’ (P=2.13×10⁻¹³), ‘sister chromatid segregation’ (P=5.26×10⁻¹³), ‘nuclear chromosome segregation’ (P=7.66×10⁻¹²), ‘mitotic nuclear division’ (P=1.13×10⁻¹¹), ‘cell cycle process’ (P=1.99×10⁻¹¹), ‘microtubule cytoskeleton’ (P=4.93×10⁻⁸), ‘cytoskeletal part’ (P=4.93×10⁻⁸) and ‘protein kinase binding’ (P=5.51×10⁻⁴). KEGG: ‘Cell cycle’ (P=4.04×10⁻⁴), ‘oocyte meiosis’ (P=1.36×10⁻²) and ‘cellular senescence’ (P=1.43×10⁻²). Detailed terms of GO terms in the categories Biological Process, Cellular Component and Molecular Function, as well as KEGG pathways, are presented in Table SVII.

A total of 17 genes were used to establish a PPI network and the top five hub genes were selected (Fig. 6B). The five hub genes were FOXM1, CENPF, KIF2C, KIF20A and NCAPG. The expression of FOXM1 and CENPF was not only found to be statistically significant in the univariate regression analysis but also in the multivariate regression analysis for overall survival, as presented in Table IV. Thus, as independent predictors of prognosis for patients with LUAD, FOXM1 and CENPF were selected as the candidate genes for further experiments.

FOXM1 and CENPF exhibited a higher expression in LUAD vs. normal tissue (Fig. 4C and D). Kaplan-Meier survival curves indicated that higher expression of FOXM1 and CENPF resulted in poorer OS of patients with LUAD (Fig. 6C). Furthermore, correlation analysis revealed that FOXM1 expression was significantly positively associated with that of CENPF expression, both in the GEO datasets and GEPIA (Fig. 6D).

According to the expression levels of FOXM1 and CENPF, the filtered data matrix based on Fig. 2 was divided into high (≥50%) and low expression (<50%) groups and ssGSEA analysis was performed. The significant pathways and molecular functions are presented in Fig. 6E. FOXM1: ‘CELL-CYCLE’ [Enrichment Score (ES)=0.6773, Normal P-value (NP)=0.0000], ‘DNA_REPLICATION’ (ES=0.7756, NP=0.0000), ‘PROGESTERONE_MEDIATED_OOCYTE_MATURATION’ (ES=0.4814, NP=0.0000), ‘P53_SIGNALING_PATHWAY’ (ES=0.6226, NP=0.0000), ‘CELL_ADHESION_MOLECULES_CAMS’ (ES=−0.5472, NP=0.0309). CENPF: ‘OOCYTE_MEIOSIS’ (ES=0.5350, NP=0.0000), ‘DNA_REPLICATION’ (ES=0.7120, NP=0.0020), ‘PYRIMIDINE_METABOLISM’ (ES=0.5223, NP=0.0138), ‘CELL_ADHESION_MOLECULES_CAMS’ (ES=−0.5047, NP=0.0462), ‘MAPK_SIGALING_PATHWAY’ (ES=−0.3332, NP=0.0423).

IHC was performed to verify the expression of FOXM1 and CENPF in LUAD and normal lung tissue. These results indicated that both FOXM1 and CENPF were more highly expressed in LUAD compared with normal lung tissue (Fig. 7).

The RT-qPCR results indicated that si-FOXM1#2 and si-CENPF#2 exhibited significant knockdown efficiency (Fig. 8A). To further investigate the potential effect of FOXM1 and CENPF on cell proliferation, CCK-8 and cell colony formation assays were performed to evaluate cell proliferation of A549 cells transfected with si-FOXM1 and si-CENPF. The CCK-8 assay results suggested that cell proliferation in the si-FOXM1 and si-CENPF groups was decreased compared with the si-NC group (Fig. 8B). The colony-formation assay indicated that the colony-formation rate was significantly decreased in the si-FOXM1 and si-CENPF groups compared with the si-NC group (Fig. 8C).

Wound-healing and Transwell assays were performed to determine the relative migration ability of A549 cells. It was found that wound closure was significantly delayed in the si-FOXM1 and si-CENPF groups as compared with the si-NC group (Fig. 8E), which was also observed in the Transwell assay (Fig. 8D).