Patients with no survival time available and follow-up time of less than 1 month or more than 120 months were excluded, then mRNA data and clinical information of 257 CC patients were downloaded from the TCGA database. The clinical statistical information of the TCGA all set is shown in Additional file 1: Supplementary Table S1. Another dataset GSE44001 consisting of 300 CC patient with associated prognostic information was obtained from the Gene Expression Omnibus (GEO) database.

First, we extracted 13 m5C regulatory factors from the TCGA expression matrix. Based on consistent clustering of these 13 genes, 257 CC samples were clustered by the “NMF” method, which was used to select the standard “brunet” option for 50 iterations. The number of clusters k was set at 2 to 10 and the average contour width of the common member matrix was determined by the “NMF” package. The minimum member of each subclass was set to 10. According to cophenetic, rss and silhouette, the optimal number of clusters was determined. KM analysis was used to analyze the difference in prognosis between different subtypes of the patients using the “survival” package and heatmaps were drawn using the “pheatmap” package.

In order to identify the immune infiltration differences between different m5C modification subtypes, we used “MCPcounter” package to evaluate the score of 10 immune cells and the score of 28 immune cells were evaluated by the “single sample Gene set enrichment analysis (ssGSEA)” method in the “gene set variation analysis (GSVA)" package (Charoentong et al., 2017). Besides, we analyzed the mRNA level differences of 13 m5c-related genes between different subtypes.

The “limma” package was used to calculate the DEGs between different m5C modification subtypes, and the filter was applied according to the thresholds FDR <0.05 and |log₂FC| > log₂ (1.5). Furthermore, 601 upregulated DEGs and 113 downregulated DEGs were analyzed by the “WebGestalt” package for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation.

Under the premise there is no preference in the distribution of clinical characteristics of the grouped samples, a total of 257 patients in TCGA all set were randomly divided into training set (n = 128) and test set (n = 129). The TCGA training set and TCGA test set were evaluated by Chi-square test, the sample information is shown in Additional file 1: Supplementary Table S2. In addition, we also used the TCGA all set (n = 257) and GSE44001 set (n = 300) as validation set for subsequent verification. In the training set, a univariate Cox regression analysis was conducted by the “survival coxph function” package using the 714 DEGs and survival data, p < 0.01 was selected as the threshold for filtering. Next, we used the “glmnet” package for the least absolute shrinkage and selection operator (LASSO) regression to further compress the screened genes to reduce the number of genes, and finally a novel gene signature was established. LASSO retains the advantages of subset shrinkage, and is a biased estimate for processing data with multicollinearity. Based on the LASSO regression results, we developed a prognostic risk score formula, which was calculated as follows: Risk score  ( patient ) =  Σ i C o e f f i c i e n t ( mRNAi )   ×  Expression  ( mRNAi )

We calculated the risk score of each sample depending on the signature gene and drew the risk score distribution of the sample in the TCGA training set. Furthermore, we used the “timeROC” package to perform receiver operating characteristic (ROC) analysis to explore the prediction accuracy of 1 year, 3 years and 5 years survival rates. Finally, we calculated the risk score and divided the samples with risk score greater than zero into the high-risk group and samples with risk score less than zero into the low-risk group to draw KM curves. To determine the robustness of the signature, we used the same coefficient to perform the same analysis used for the TCGA test set, TCGA all set and external validation data set GSE44001.

We selected the corresponding gene expression profiles of these samples for ssGSEA via the “GSVA” R package to observe the relationship between the risk score and the KEGG pathway. We calculated the score of each sample in different KEGG pathway and obtained the ssGSEA score of each sample. Next, we calculated the correlation between these pathways and the risk score and selected pathways with a correlation greater than 0.4. Finally, the top 18 KEGG pathways were selected and clustered according to their enrichment score.

To identify the independence of the gene signature in clinical parameters, we used univariate and multivariate Cox regression to analyze the hazard ratio (HR), 95% confidence interval (CI) of HR and p value in the TCGA all set. We systematically analyzed the clinical information of TCGA patient records, including age, T stage, N stage, FIGO stage, grade, chemotherapy and risk score. According to the results of univariate and multivariate Cox analyses, we constructed a nomogram with the T stage and risk score for predicting survival outcomes (1 year, 3 years and 5 years). Then we performed the calibration curve by the “rms” package to determine the consistency between the actual survival rates and the nomogram-predicted rates. In order to evaluate the reliability of the nomogram, we performed DCA (decision curve analysis) using the “rmda” package. DCA analysis is a method that can assess whether the nomogram improves clinical decision-making. This method can tell us whether it is beneficial to use the model to make clinical decisions, or which of the two models will lead to better decisions.

Fresh adjacent normal tissues and CC tissues were obtained from the Chinese Academy of Medical Sciences and the CAMS & PUMC Medical College. All patients were not treated preoperatively and signed informed consent forms provided by the Cancer Hospital, CAMS & PUMC. The normal surgical margin tissue and the morphology of the primary tumor area were immediately excised from each patient by an experienced pathologist and stored in liquid nitrogen. The study was approved by the Ethics Committee of the Cancer Institute (Hospital), CAMS & PUMC (20/207–2,403).

The human CC cell line SiHa was provided by the Cell Resource Center, IBMS, CAMS/PUMC. The cell lines were cultured in DMEM medium supplemented with 10% fetal bovine serum (Invitrogen, San Diego, CA) at 37 °C and 5% CO2 in a humidified incubator. Human specific siRNA sequences are shown in Additional file 1: Supplementary Table S3. The transfection method was described in our previous article (Cao et al., 2018).

Total RNA was extracted using RNApure Tissue & Cell Kit (Cwbiotech, Beijing, China). Isolated RNA was used as a template for reverse transcription reactions using HiFiScript cDNA Synthesis Kit (Cwbiotech, Beijing, China). Real-time quantitative PCR analysis was performed using SYBR^® Fast qPCR Mix (TaKaRa, Shiga, Japan) and a CFX96 Real-Time System (Bio-Rad, California, United States of America). The primer sequences are shown in Additional file 1: Supplementary Table S4. GAPDH served as the internal control.

Western blotting and IHC analysis were performed as described previously (Cao et al., 2018). The antibodies used were as follows: anti-FNDC3A antibody (Abcam, Cambridge, United Kingdom), anti-VEGFA antibody (Proteintech, Wuhan, China), anti-OPN3 antibody (Affinity Biosciences, Cincinnati, United States), and anti-CPE antibody (Proteintech, Wuhan, China). The IHC quantization analysis was calculated by ImageJ software and statistically analyzed in three random fields. Data are shown as mean ± SEM.

Cells were inoculated into 96-well plates at a concentration of 2000 cells per well. According to the manufacturer’s directions, cell viability was determined by the Cell Counting Kit-8 (CCK-8, Dojindo, Japan). The absorbance was measured at 450 nm by an automatic microplate reader (BioTek, Winooski, United States). Measurements were taken every 24 h for seven consecutive days.

SiHa cells treated with siRNA were plated in 6-well plates at a density of 500 cells per well. After overnight incubation, the cells were cultured for 14 days to form colonies, fixed with methanol and stained with crystal violet. The data represent the mean ± SD of three independent experiments.

700 μL DMEM medium supplemented with 20% serum was added to the lower chamber of the Transwell plates, and 1×10⁵ suspended cells were added to the upper compartment. For the invasion experiment, 50 μL Matrigel was added to the membrane of the upper chamber. The Transwell plates were incubated in a carbon dioxide incubator for 16 h in the migration experiment and for 24 h in the invasion experiment. Then, we removed the chamber, washed the cells once with PBS, fixed the cells with solution (methanol: acetone = 1:1) for 30 min, and then stained them with 0.5% crystal violet for 30 min. The chamber was washed with PBS, and then the upper cells were carefully removed, sealed with neutral gum and photographed for counting.

R software 3.5.3 and SPSS 22.0 software (SPSS Inc. Chicago, United States) were used for all statistical analyses. p < 0.05 was taken as the probability value to establish statistical significance. Chi-square test was used for statistics of multiple categories, Student’s t-test was used to determine the significance of differences between two groups, and ANOVA was used for comparisons among more than two groups.

To clearly illustrate the process of our research, a flow chart is shown in Figure 1.

We extracted mRNA levels of 13 m5C regulatory factors from the expression matrix of the TCGA. Then, 257 CC samples were clustered by “NMF” package. The optimal number of clusters was determined according to cophenetic, rss and silhouette analyses, the optimal number of clusters was 2 (Figures 2A,B). The expression levels of m5C methylation-related genes in the C1 and C2 subtypes were significantly different (Figure 2C). The KM curve revealed that overall survival (OS) rates of the C1 and C2 subtypes were significantly different (p < 0.05), and the prognosis of the C1 group was worse than that of the C2 group (Figure 2D).

Due to the significant difference in the prognosis of CC patients with two m5C modification subtypes, we next explored the difference in immune cell infiltration between C1 and C2 subtypes. The ssGSEA score suggested that levels of activated CD8 T cells, central memory CD4 T cells, CD56 bright natural killer cells, macrophages, MDSCs and neutrophils were markedly different between two subtypes, and the MCPcounter score indicated that the infiltration levels of CD8 T cells, NK cells, neutrophils, endothelial cells and fibroblasts were significantly different (Figures 3A,B). Furthermore, we also analyzed the expression of 13 genes between two subtypes. The expression of eight genes in C1 and C2 subtypes was substantially different; but no difference was found in NOP2, NSUN4, NSUN7, TRDMT1 and DNMT3A expression (Figure 3C). The above results indicated that there are significant differences in the immune infiltration of C1 and C2 subtypes, while the expression of most m5C regulatory factors in the two is also different.

DEGs between C1 and C2 subtypes consisted of 601 upregulated and 113 downregulated genes. The volcanic map of representative DEGs is represented in Figure 4A. Detailed information about these DEGs is represented in Additional file 2: Supplementary Table S5. We selected 50 genes with the most prominent changes in expression (upregulated and downregulated), as shown in the heat map (Figure 4B).

Next, 601 upregulated DEGs were analyzed by GO function and KEGG pathway enrichment annotation. For the GO function analysis of the upregulated DEGs, 493 gene ontologies were annotated to biological process, 88 to cellular component and 81 to molecular function with significant differences (p < 0.05). The top 10 annotations are shown in Figures 4C–E. For the KEGG pathway enrichment analysis of upregulated DEGs, the top 10 KEGG pathways annotated are shown in Figure 4F, including adherens junction, ECM-receptor interaction, amoebiasis, focal adhesion, human papillomavirus infection, PI3K-Akt signaling pathway, and pathways in cancer. More detailed information can be found in Additional file 3: Supplementary Table S6. For CC downregulated DEGs, the results of GO function and KEGG pathway enrichment analysis are shown in Additional file 4: Supplementary Table S7.

First, 257 samples in the TCGA-CC dataset were divided into a training set and a test set. The training set consisted of 128 samples, and the test set consisted of 129 samples. The statistical results showed that our groups had no preference, and there was no significant difference between the training set and the test set (Additional file 1: Supplementary Table S2).

Using the training set data, a univariate Cox proportional hazards regression analysis was conducted by the “survival coxph function” package for DEGs (714 genes) and the survival data, and p < 0.01 was selected as the threshold for filtering. Finally, 27 genes were selected, and the univariate Cox analysis results are shown in Figure 5A. Next, we used LASSO regression to further compress these 27 genes to reduce the number of genes in the risk signature. We performed 10-fold cross validation to construct the model and analyzed the confidence intervals under each lambda, as shown in Figure 5B. We selected the following four final genes with lambda = 0.09,009,939: FNDC3A, VEGFA, OPN3 and CPE.

The KM curves all showed that a higher mRNA level of those four genes indicated worse prognosis in the TCGA training set (Additional file 5: Supplementary Figure S1, p < 0.05). The final 4-gene signature formula was as follows:

Risk score = 0.3,250,335*FNDC3A (mRNA level)＋0.2,821,988*VEGFA (mRNA level)＋0.3,133,706*OPN3(mRNA level)＋0.1,857,458*CPE (mRNA level).

We calculated the risk score of each sample according to the mRNA level of the signature gene in the training set, and the proportion of deaths in the high-risk group was significantly higher than in the low-risk group. This demonstrated that the risk score is a critical prognostic factor. Consequently, with the increase in risk score, the mRNA levels of FNDC3A, VEGFA, OPN3 and CPE were upregulated (Figure 5C). To investigate the diagnostic accuracy of the risk signature, the AUC of the time-dependent receiver operating characteristic (ROC) curves was computed. The AUC values of the signature for predicting 1 year, 3 years and 5 years survival rates were 0.74, 0.76 and 0.80 (Figure 5D). The KM curve suggested that patients with higher risk score had worse prognosis than those with lower risk score (Figure 5E, p < 0.001).

To determine the robustness of the model, we used the TCGA test set (n = 129) with the same model and coefficient as the training set to calculate the risk score of each sample and drew the risk score distribution of the patients. Then, we performed ROC analysis on the prognostic classification of the risk score and analyzed the classification efficiency of 1 year, 3 years and 5 years survival rates. Finally, we divided the samples with risk score into high-risk group (n = 65) and low-risk group (n = 64) to perform the KM survival analysis. Importantly, the results of the above analysis were consistent with the performance of the TCGA training set (Additional file 6: Supplementary Figure S2 ).

In addition, we determined the robustness of the signature in the TCGA all set (Figures 6A–C) and the external validation dataset GSE44001 (Figures 6D,E ). The proportion of deaths in the high-risk group was significantly higher than in the low-risk group, which was consistent with the performance of the TCGA training set. The AUC values of the signature showed that the risk score was a good prognostic factor. Finally, the KM curve also showed that there were consistent differences between the high-risk and low-risk groups.

Furthermore, we performed KM survival analysis according to age, grade, stage, recurrence and chemotherapy treatment in the TCGA all set. For the patients in age < = 60, T1+T2 Stage, T3+T4 Stage, N0 Stage, N1 Stage, M0 Stage,Ⅰ+Ⅱ Stage, Ⅲ+Ⅳ Stage, G1+G2, G3+G4, Recurrence _ No and Chemotherapy _ Yes subgroup, the OS interval of patients in the high-risk group was significantly shorter than that of patients in the low-risk group (p < 0.05). Only in the age >60, M1, Recurrence _ Yes or Chemotherapy _ No subgroup was the OS interval of patients not different between the high- and low-risk groups (p > 0.05). The above findings further showed that the risk signature still has good predictive ability among different clinical subgroups (Figure 7 ).

To observe the relationship between risk score and the KEGG pathway, we selected the gene expression profiles corresponding to these samples for ssGSEA using the “GSVA” R package. The score of each sample in different KEGG pathway were calculated, and the ssGSEA score of each KEGG pathway corresponding to each sample were obtained. The correlation between these pathways and risk score was further calculated, and the function with correlation greater than 0.4 was selected (Figure 8A). Fifteen pathways were positively correlated with the risk score, and three pathways were negatively correlated with the risk score. The top 18 KEGG pathways were selected and clustered according to their enrichment score, as shown in Figure 8B. We can find that KEGG_MTOR_SIGNALING_PATHWAY, KEGG_ECM_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, KEGG _TGF_BETA_SIGNALING_PATHWAY, KEGG_ADHERENS_JUNCTION, KEGG_ PATHWAYS_IN_CANCER and other tumor-related pathways were activated with increasing risk score.

To identify the independence of the 4-gene signature in clinical parameters, we used univariate and multivariate Cox regression to analyze the related HR, 95% CI of HR and p value in the clinical information of the TCGA all set. The clinical data of patients were analyzed systematically, including age, T stage, N stage, FIGO stage, grade, chemotherapy and risk score (Figures 9A,B). In the TCGA all set, the risk score was an independent prognostic factor. Therefore, the 4-gene signature has good predictive performance and clinical application value.

According to the results of univariate and multivariate Cox regression analyses, we constructed a nomogram with clinical features, T stage and risk score (Figure 9C). We found that the risk score had the greatest impact on survival prediction, indicating that the risk score is indispensable in the nomogram. Furthermore, we used the calibration curve to evaluate the prediction accuracy of the signature, as shown in Figure 9D. The prediction calibration curves of the three calibration points for 1 year, 3 years and 5 years survival rates were close to the standard curves, indicating that the signature had good prediction performance. In addition, DCA showed that the benefits of the risk score and nomogram were significantly higher than those of the extreme curves. The nomogram curve was higher than that of the risk score, which indicated that the nomogram had good reliability (Figure 9E).

In order to explore the difference in protein expression of signature genes, we used IHC to detect 6 normal and 21 tumor tissues. The IHC quantization analysis was calculated by ImageJ software and statistically analyzed in three random fields. The results showed that the protein expression levels of the four signature genes in tumor were higher than those in normal (Figure 10A). At the same time, we conducted qRT-PCR experiments in 10 normal and 12 tumor tissues to explore the differences in the transcription level of those four genes. It showed that the mRNA levels of the four model genes were significantly increased in tumor tissues, indicating that the expression of model genes in tumor tissues may be abnormally activated (Figure 10B).

To clarify the functional role of FNDC3A, VEGFA, OPN3 and CPE in CC cells, we applied human-specific siRNA to decrease their protein expression (Figure 11A). The CCK-8 assay was applied to detect cell proliferation. The downregulation of FNDC3A, VEGFA or CPE expression significantly suppressed the proliferation (Figure 11B) and colony formation capacity of SiHa cells (Figure 11C). Transwell assays were applied to detect the invasion and migration ability of SiHa cells in vitro, and the number of cells that passed through the polycarbonate membrane was smaller in the FNDC3A, VEGFA or CPE siRNA group than in the negative control group, indicating that FNDC3A, VEGFA or CPE could significantly promote the invasion and migration of SiHa cells (Figure 11D ). Downregulation of OPN3 expression had no significant effect on the proliferation, colony formation ability, invasion and migration of SiHa cells.