2.1 Study design and data sets

Four matrix datasets (GSE9750 [9], GSE63514 [10], GSE122697 [11] and GSE7803 [12]) containing data on the tissues from patients with normal cervical epithelium and CESC were selected for analysis. The GSE9750, GSE63514 and GSE122697 datasets were used as the training group, while the GSE7803 dataset was used as the validation group. The details of the training and validation groups are presented in Table 1. Fig. 1 illustrates the flowchart of the study design. The raw data for the above datasets were sourced from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

Fig. 1.The flowchart of the study design. GEO: Gene Expression Omnibus; CESC: cervical squamous cell carcinoma; DEGs: differentially expressed genes; WB: Western blot; RT-qPCR: real-time quantitative polymerase chain reaction.

The top 30 differentially expressed genes (DEGs) significantly contributing to cervical squamous cell carcinoma (CESC) diagnosis prediction were identified using the random forest algorithm. Gene expression scores were then computed based on the expression data of these DEGs in all samples in the training group. A neural network model was created to obtain weight values of the genes with an importance score >2 and the greatest predictive contribution to CESC diagnosis. Based on the obtained gene expression scores and weight values, a molecular prognostic scoring system was constructed. The validity of the gene model was confirmed using the publicly available GSE7803 dataset, and ethical review was not required as the dataset is publicly available in the GEO database.

2.2 Identification of DEGs

Both the training and validation group data were processed using R software (version 4.1.1, Statistics Department of the University of Auckland, Auckland, New Zealand), which facilitated data standardization and normalization. The steps of merging training set data are as follows: 1. Read and merge the gene information of each training set; 2. Take log2 for the data with large values; 3. Merge the data; 4. Output the corrected data result. Duplicate gene probes were eliminated, and the consolidated data was used for subsequent analysis (Supplementary Fig. 1). The R software “Limma” [13] package was used to eliminate the batch effect in normal and tumor samples of the training group when identifying DEGs based on logFC >2 (log2 (Fold Change) >2) and adj. p.value < 0.05. A volcano plot was then created to visualize the identified DEGs. Finally, “pheatmap” and “ggplot2” packages of R software are used to draw heat map and volcano map of DEGs [14].

2.3 Enrichment analysis of DEGs

To obtain comprehensive information on the biological functions and signaling pathways associated with the significant differentially expressed genes (DEGs) in the training group dataset, pathway enrichment analysis and gene ontology annotation were performed using the Metascape online database (http://metascape.org) [15]. The results were then grouped into clusters based on the similarity of the significant terms, and the most prominent term was selected to present each cluster. Gene Ontology (GO) circles showing biological functions are drawn through the R software “clusterProfiler”, “org.Hs.eg.db”, “enrichplot” and “GOplot” package [16, 17, 18, 19].

2.4 Identification of DEGs for predictive diagnosis of CESC based on random forest algorithm and artificial neural network

The R software “randomForest” package was used to identify the top 30 genes that significantly impacted the diagnosis of CESC [20], which were then used to generate a gene score table reflecting their expression levels [21]. The expression values of the DEGs were converted into binary values of 1 or 0 based on specific conversion criteria. If the expression value of an upregulated gene in a sample was higher than the median expression value of that gene across all samples, it was assigned a value of 1; otherwise, it was assigned a value of 0. Conversely, if the downregulated gene was higher, it was assigned a value of 0; otherwise, it was assigned a value of 1. The R software “neuralnet” and “NeuralNetTools” [22] packages were used to develop a predictive model featuring a single input layer, a single hidden layer, and a single output layer. The hidden layer was configured with five hidden nodes, and the output layer was designed with two nodes using a softmax activation function. The cross-entropy error function was set, and the optimization process involved selecting the maximum weight value of the DEGs within the hidden layer with values optimized accordingly [23].

2.5 Development and validation of new gene signatures

The molecular prognostic scoring system is a novel scoring system that has proven successful in predicting patients with breast cancer [23] and was implemented in this study to predict a new gene signature for the diagnosis of cervical squamous cell carcinoma (CESC). The scoring formula for evaluating each differentially expressed gene (DEG) was computed as follows: system score = (gene score × gene weight). The total system score of the eight DEGs with the highest predictive CESC diagnostic ability was then calculated. The publicly available GSE7803 dataset was utilized to validate the model, and the accuracy was assessed through the area under the curve (AUC) value of the receiver operating characteristic (ROC) curve generated using the “pROC” package in R. An AUC value greater than 0.8 indicated favorable accuracy, while an AUC value >0.9 indicated outstanding accuracy of the model [24].

2.6 Cell culture

The human cervical epithelial cell line HCerEpC (JZ-004978) and human cervical cancer cell line Hela (CL-0101) were cultured in MEM medium containing 10% fetal bovine serum and maintained in a 5% carbon dioxide (CO₂) incubator at 37 °C.

2.7 Western blot (WB) and real-time quantitative polymerase chain reaction (RT-qPCR)

Cell lysis buffer was used to lyse the collected HCerEpC and Hela cells, the proteins were extracted, and the protein concentration was determined. The protein levels of MELK (A10794, abclonal), CRISP3 (14847-1-AP, proteintech), GINS2 (A9172, abclonal), DTL (A12150, abclonal), C1orf112 (PA5-55082, ThermoFisher), KIF14 (A10275, abclonal), SPINK5 (A20916, abclonal) and CELSR3 (sc-293381) were detected by western blot (Primary Antibody Dilution Buffer for Western Blot & Secondary Antibody Dilution Buffer for Western Blot: Beyotime Biotechnology). Development was performed with the aid of a supersensitive luminescent substrate solution (Biosharp) using a protein gel imaging analyzer (Bio-Rad). Western blot was analyzed by Image J software (ImageJ 1.48v, National Institutes of Health, Bethesda, MD, USA). Total RNA was extracted from each cell and reverse-transcribed into cDNA using primer sequences (Sangon Biotech) listed in Supplementary Table 1 and following the instructions of the detection kit (QIAGEN kit). The reaction procedures were as follows: PCR initial heat activation at 95 °C for 2 min; denaturation at 95 °C for 5 s, Combined annealing/extension at 60 °C for 30 s for a total of 40 cycles. β-Actin was used as an internal control. The mRNA expression of MELK, CRISP3, GINS2, DTL, C1orf112, KIF14, SPINK5 and CELSR3 was quantified using the 2^−ΔΔCt method [25].

2.8 Immuno infiltration analysis

The CIBERSORT R script v1.03 was used to calculate the infiltration score of each immune cell in the two groups of samples. The R software “corrplot” [26] and “vioplot” packages were used to analyze the correlation and difference of immune cells, following a correlation heat map and infiltration histogram were subsequently drawn.

2.9 Statistical analysis

All statistical analyses were performed using R language software (version 4.1.1). Statistical significance was set at p < 0.05.

3.1 Identification of DEGs

A total of 241 DEGs were identified in the training dataset, and their expression profiles are displayed through heat maps and volcano plots (Fig. 2A,B). The heat map revealed that approximately half of the DEGs were highly expressed in CESC, while the rest exhibited low expression, thereby suggesting that both oncogenes and tumor suppressor genes were present among the DEGs. logFC of DEGs are listed in Supplementary Table 2.

Fig. 2.Expression of all differential genes in the training group. (A) Heat map (Red: High expression; Blue: Low expression). (B) Volcano map (Green: Down-regulated; Red: Up-regulated). logFC: log fold change; Sig: Significance Level.

3.2 Enrichment analysis of DEGs in the training dataset

Metascape (v3.5, The team of Yingyao Zhou, San Diego, CA, USA) was used to perform an enrichment analysis to better understand the functional and metabolic pathways associated with the DEGs. The results are depicted in a clustering network diagram and show the top 20 clusters in which the DEGs were significantly enriched (Fig. 3A).

The most significantly enriched biological processes were the “mitotic cell cycle process” and “keratinized envelope formation”. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that the DEGs were primarily involved in the “vitamin D receptor pathway” and “Polo-like Kinase 1 (PID PLK1) pathway”. Furthermore, GO circle diagrams were drawn to identify the regulatory role of the DEGs in the enrichment analysis. The majority of DEGs were upregulated in “mitotic sister chromatid separation, sister chromatid separation and organelle fission”, while they were downregulated in “epidermal development and skin development” (Fig. 3B).

Fig. 3.Enrichment analysis of all differential expressed genes in the training group. (A) Cluster network plot (Different colored squares represent different signaling pathways). (B) GO circle plot (Different colored squares represent different biological functions).

3.3 Random forest algorithm identification of disease signature genes

The expression data of 241 DEGs were integrated into the random forest algorithm classifier (Fig. 4A,B). A screening process was conducted to identify the eight disease signature genes with importance scores above 2, which were SPINK5, CRISP3, MELK, CELSR3, GINS2, DTL, C1orf112 and KIF14. Except for CRISP3 and SPINK5, the remaining six genes demonstrated high expression levels in CESC and low expression levels in the normal cervical epithelium (Fig. 4C). The neural network showed that MELK, CRISP3, GINS2, DTL, C1orf112, KIF14, SPINK5 and CELSR3 had good CESC diagnostic predicting (Fig. 4D).

Fig. 4.Random forest algorithm for identification of disease signature genes. (A) Random forest. (B) Random forest algorithm plot of the top 30 DEGs. (C) Heat map of DEGs with importance scores greater than 2. (D) Neural network model.

3.4 Artificial neural network-based molecular prognostic scoring system

After converting the expression data of the eight disease-specific genes with importance scores greater than 2 into “gene scores”, the weight values of each gene were optimized using an artificial neural network algorithm. The molecular prognostic score was calculated by summing the systematic scores (“gene score × gene weight”) of these eight signature genes (Table 2). The molecular prognostic scores of 134 samples in the training group were then used as predictive values, with the presence or absence of CESC in patients as the true value. The AUC of the new gene model was 0.995, indicating the model’s excellent predictive ability (Fig. 5).

Fig. 5.ROC plot of the training group evaluated of model accuracy. AUC: area under the curve; CI: confidence interval.

3.5 Validation of the prediction model

To validate the prediction model developed in the training set, an independent dataset (GSE7803) was used to evaluate its ability to diagnose CESC in other individuals. The random forest algorithm was utilized to select the top 30 DEGs and disease signature genes with importance scores greater than 2 from the validation set. The results showed that these genes were identical to those in the training set, demonstrating the scalability and stability of the random forest algorithm. Subsequently, the “gene score” and “molecular prognostic score” of GSE7803 were calculated using the same methodology in the testing set. The validation model generated a ROC curve with an AUC of 1.000, indicating that the prediction model is highly valid and stable (Fig. 6).

Fig. 6.ROC plot of the testing group evaluated of model accuracy. AUC: area under the curve; CI: confidence interval.

Next, we investigated the expression levels of the eight CESC-characteristic genes, namely MELK, CRISP3, GINS2, DTL, C1orf112, KIF14, SPINK5 and CELSR3, in human cervical epithelial cells (HCerEpC) and human cervical carcinoma cells (Hela). The results revealed that the expression of the CRISP3 and SPINK5 protein and mRNA was high in HCerEpC and low in Hela cells, while the expression levels of MELK, GINS2, DTL, C1orf112, KIF14 and CELSR3 protein and mRNA were significantly lower in HCerEpC but significantly higher in Hela cells (Fig. 7A–I,8A–H, Supplementary Table 3, Supplementary Fig. 2), indicating that these genes may have a role in the development and progression of CESC and could serve as potential targets for further research or therapeutic intervention.

Fig. 7.The expression of eight disease-specific proteins in human cervical epithelial cell line HCerEpC and human cervical cancer cell line Hela was detected by WB. (A) The total. (B) MELK. (C) CRISP3. (D) GINS2. (E) DTL. (F) C1orf112. (G) KIF14. (H) SPINK5. (I) CELSR3. **: p < 0.01; ***: p < 0.001; ****: p < 0.0001.

Fig. 8.The expression of 8 disease characteristic genes in human cervical epithelial cells HCerEpC and human cervical cancer cells Hela was detected by RT-qPCR. (A) MELK. (B) CRISP3. (C) GINS2. (D) DTL. (E) C1orf112. (F) KIF14. (G) SPINK5. (H) CELSR3. *: p < 0.05; **: p < 0.01.

3.6 Immunoinfiltration analysis

The analysis of immune infiltration in two groups based on the predicted model revealed a higher degree of macrophage infiltration in squamous cell carcinoma of the cervix than in the normal cervical epithelium (Fig. 9A). The quantitative comparison of immune infiltration in the two groups showed significant differences, based on the prediction model, in immune cells such as T cells CD8, T cells CD4 naive, T cells CD4 memory activated, Macrophages M0, Macrophages M1 and Dendritic cells resting (p < 0.05, Fig. 9B). Notably, Macrophages M0 were negatively correlated with Dendritic cells resting (−0.41), T cells CD8 (−0.33), and activated T cells CD4 memory (−0.14), while they were positively correlated with Macrophages M2 (0.37) (Fig. 9C). Collectively, these results suggest that the predicted model could be useful in analyzing immune infiltration in cervical cancer and identifying potential targets for immunotherapy.

Fig. 9.Immunoinfiltration analysis. (A) Two groups of immune cell infiltration distribution. (B) Difference analysis of immune cells between two groups. (C) The correlation between immune cells.