We retrieved the Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) database’s microarray data set and associated clinical information for cervical cancer, and 306 cervical cancer samples were utilized for further investigation. The Gene Expression Omnibus (GEO) database was searched for, and the microarray data set for the GSE30759 chip was acquired; it comprises data on 48 cervical cancer patients. Only cervical cancer patients with comprehensive data, such as status and overall survival time, were gathered. Transcripts per kilobase million (TPM value) of the microarray data set (FPKM value) were translated to preserve compatibility between the two databases. The original data of the GSE30759 chip obtained from the GEO database was unified, and then the gene probes were annotated.

Combined TCGA and GEO databases to assess cervical cancer’s tumor immune microenvironment level. Based on 547 gene expression values, the proportion of cervical cancer samples with mixed cell types in 22 different immune cell subpopulations was evaluated using the CIBERSORT algorithm to obtain the ICI score matrix for different immune cells. Used the R package “ConsensuClusterPlus” for analysis and iterated 1,000 times to obtain stable classification results. Analyzed the distinctive characteristics of the stromal and immune cell transcription patterns in cervical cancer and deduced tumor cellularity and tumor purity using the ESTIMATE algorithm.

Using the R package “limma,” differentially expressed genes (DEGs) were found, where the absolute folding change>1 and the critical value standard P<0.05 were considered to be statistically different. The data were split into different genomic clusters according to the DEGs value using unsupervised clustering. Gene features are positively correlated with gene cluster A, whereas gene characteristics are negatively correlated with gene cluster B. The first principal component (PC1) produced by principal component analysis (PCA), after dimensionality reduction analysis using the Boruta approach, was chosen as the feature score, two PC1 values were obtained for each sample, the optimal cutoff value was taken as the dividing line and divided into high ICI score group or low ICI score group. ICI scores were generated for each instance using the gene expression rank index approach, using the formula:

Gene ontology (GO) analysis was carried out on the two gene clusters to understand their biological roles better. Gene Set Enrichment Analysis (GSEA) based on two ICI score grouping, P<0.05, was determine statistically different.

The TCGA database gathered data on somatic mutations in cervical cancer. Data on somatic mutations were examined and categorized using the ICI score system. The “maftools” package was used to identify the 20 mutated genes in both mutation data sets.

The Genomics of Drug Sensitivity in Cancer (GDSC) website (https://www.cancerrxgene.org/) provided information on the drug’s predicted half-maximal inhibitory concentration (IC50), and variations in drug sensitivity to Cisplatin, 5-Fluorouracil, Vinblastine, Vincristine, Paclitaxel, Docetaxel, Gemcitabine, Cyclophosphamide, Topotecan, Epirubicin, Olaparib, Dactinomycin were evaluated among ICI score groups.

Cervical cancer tissue samples were collected from patients before receiving immunotherapy at Wuhan University Central South Hospital. Tissues were fixed with formalin and then embedded with paraffin. All samples were stored at room temperature 20°C-25°C. All specimens were diagnosed and agreed upon by at least two pathologists based on pathological features. Finally, this study included 16 cases of cervical cancer tissues that had received immunotherapy. The Zhongnan Hospital of Wuhan University ethical committee approved the research (ethics number: 2020029).

Cervical cancer tissues were fixed with formalin, followed by paraffin embedding to obtain 4um thick sections of tumor tissue. Paraffin-embedded portions were dewaxed with xylene and ethanol, hydrated and blocked paraffin sections. Sections were treated with primary antibodies, which were then incubated with them overnight in a moist box kept at 4°C. After that, PBS was used to wash the parts three times for a total of eight minutes each. The sections were then treated with secondary antibodies for 30 minutes at room temperature before being rinsed with PBS. Visualized with 3,5-diaminobenzidine (DAB) substrate kit and restrained with hematoxylin for 2 min. We divided the staining intensity of CD8A, CXCL10 and GZMB into 1, 2, 3 or 4. The stronger the staining, the higher the score. Two pathologists evaluated the expression of CD8A, CXCL10, and GZMB in the tissue microarrays in a blind manner.

Xylene was used to dewax paraffin-embedded (FFPE) tissue sections (4 m) before they were hydrated with a series of graded ethanol solutions in deionized water. Blocking was performed by adding 3% H₂O₂ and left for 25min at room temperature protected from light, then rinsed three times with PBS for 5min each. 3% BSA was added and incubated for 30min. Serial staining was performed with the following antibodies, CD3, CD4, CXCL10, CD8A, and GZMB, and incubation was maintained overnight at 4°C. Next, the sections were shaken in PBST (pH 7.4) for three washes, each lasting 8 min. Anti-rabbit polymeric horseradish peroxidase-labeled secondary antibody was subsequently added and incubated for 50 min, protected from light. Tyrosine-CY3 was then introduced and continued to work for 20 minutes. Three PBST washes were then performed, each lasting for five minutes. DAPI staining solution was added to stain nuclei for 10 min. Sections were sealed with anti-fluorescence quenching sealer. Under the excitation of UV light, blue light represents DAPI-stained nuclei, and red light, green light, and pink light represent the corresponding fluorescein-labeled nuclei, respectively.

Univariate Cox regression analysis was used to assess each DEG’s expression levels. Multivariate Cox regression analysis was used to identify the three key genes, CST7, IL1B, and ITGA5, with the greatest predictive power. Patients were divided into high and low risk groups based on the results of the median risk assessment. Kaplan-Meier curves were used to describe overall survival rates, and the precision of the model’s predictions was evaluated using the AUC of the ROC curve.

Data used in this study were publically available and downloaded from the Gene Expression Omnibus and Genome Sequence Archive (GSE168652) (31), included two samples: tumor and normal), and Seurat (32) was used to perform the next analysis. Low-quality cells were filtered according to the original articles and our requirement: cells with less than 200 expressed genes or > 15% of mitochondrion-derived counts or< 0.8 of the number of genes detected per UMI were removed. All cells that passed QC were normalized to identify the top 2000 high-variable genes. Mitochondrial score and ribosome score and cell cycle were regressed in the ScaleData function to remove batch effects. Then the top 20 principal components (PCs) were selected by the RunPCA function and were set in the RunUMAP function and RunTSNE function to reduce dimensions and in cell clustering. Subsequently, the main cell clusters were identified with the FindClusters function with the method “original Louvain algorithm”. Method “MAST” in FindAllMarkers function was used to identify marker genes of different clusters, DEGs were selected by P<0.05, and cell types of clusters are first identified by singleR and Cell Markers database, then they are manually corrected by the markers in the original article. The different immune cells and non-immune cell expression matrices and metadata were imported into monocle2 (33). The high variable genes identified by monocle2 were used for pseudo-trajectory analysis to identify cell differentiation trajectories and identify dynamical expression genes. The dynamical expression heatmaps were shown by the plot-pseudotime-heatmap function. Then we imported the expression matrix and metadata into CellChat (34), CellChat provides a database of ligand receptors in humans and mice to identify cell-cell actions in single-cell datasets. We first identified the different numbers and strengths of cell interactions between immune cells and non-immune cells to find the main sources and outputs of different types of cells, and we used Euclidean distance and information flow to infer the conservative and specific cell-cell communication signal pathways and ligand-receptor pairs that mediate cell communication of cell groups(thresh<0.05).

In DMEM conditions containing 10% FBS, cells from the Hela and Siha cell lines were cultured. We purchased si-IL1B and si-ITGA5 from GenePharma. The following is the corresponding assay sequence of IL1B and ITGA5: si-IL1B: 5’-GATGTCTGGTCCATATGAA-3’; si-ITGA5: 5’-ACGAACCTCTTCTGTGATGGA-3’. We obtained the pGL4.10-CST7 promoter plasmid from Obio Technologies, Inc. Lipofectamine 2000 was used for cell transfection.

Cells (3000/well) were inoculated on 96-well plates for the MTT assay, and absorbance was recorded at various time intervals to gauge cell viability.

Cells (2000/well) were inoculated in 6-well plates for the clonogenic assay, fixed after 14 days, and stained with crystal violet.

Thirty thousand cells were inoculated in the top chamber (Corning) of the Transwell migration test using serum-free media to promote migration to the bottom chamber. Then, cells were added in transwell chambers with or without matrigel covered in the chambers for invasion andmigration detection, respectively. The cells were fixed and stained after 24 hours. Using NIH ImageJ software, migratory cells were counted in three randomly chosen fields in each well.

Following the manufacturer’s instructions, we extracted total RNA using the RNeasy mini kit (Cat. #74,101, Qiagen). The amount of RNA was then measured, and cDNA was produced using the reverse transcription. Real-time reverse transcription-PCR (qRT-PCR) was used to evaluate the cDNA, and the iQ™ SYBR^® Green Super Mix (Bio-RAD) was used. The following were the primer sequences:

After placing the cells in RIPA buffer (MilliporeSigma), which contains phosphatase inhibitor and protease inhibitor, and lysing them on ice for 30 minutes, the supernatant was collected by high-speed centrifugation. Protein extract samples were run through SDS-PAGE gels to separate them. They were then put on PVDF membranes, covered with 5% skim milk, and subjected to primary and secondary antibody incubation.

The Wilcoxon test compared two sets of data, while the Kruskal-Wallis test was used when there are more than two data groups. The Kaplan-Meier curve was used to visualize the results of the survival analysis (log-rank test). Built heat maps using the “pheatmap” package. The link between the ICI score subgroups and the incidence of somatic mutations was examined using a chi-square test and Spearman correlation analysis. Further analysis of the results of the CIBERSORT method was found statistically meaningful with a P<0.05. Each analysis was performed using the R program (version 4.1.1).

We searched public databases and merged 354 cervical cancer date from the TCGA and GEO databases. Based on the ESTIMATE and CIBERSORT algorithms, evaluated the immune cell subsets of cervical cancer and the immune characteristics of cervical cancer. According to the results of the CIBERSORT algorithm, unsupervised clustering of 22 immune cell types led to the discovery of three separate ICI clusters ( Supplementary Figure 1 ; Figure 1A ). The three separate ICI clusters’ overall survival rates were significantly different from one another(P=0.019) ( Figure 1B ); among them, the prognosis for Cluster B is the best, while Cluster C is the worst. To investigate potential relationships between immune cell infiltration and immune cell clusters, we examined the immune cell composition in three distinct clusters ( Figure 1C ). We discovered significant differences in immune scores between the three clusters but no differences in stromal scores. The relationship between the different immune cell infiltrates was shown using a heat map ( Figure 1D ). T cells CD8 and T cells CD4 memory activated demonstrated greater levels of activation in cluster B. In Cluster C, there are more Dendritic cells resting, NK cells activated, Mast cells resting, and T cells CD4 memory resting. Additionally, in order to more thoroughly examine the innate immunotherapy potential of various immune cell clusters, we analyzed the level of prominent immune checkpoints in these three ICI clusters. It showed that CTLA4 ( Figure 1E ), PD-L1 ( Figure 1F ), and PD-1 ( Figure 1G ) had higher expression levels in cluster B.

To compare the variations in transcriptional expression across the three distinct ICI clusters, we utilized the “limma” package to obtain 119 DEGs from 3 ICI clusters. Unsupervised clustering was used to separate the samples from the TCGA and GEO datasets into distinct gene groups, which were given gene cluster A and B ( Figure 2A ; Supplementary Figure 2 ), genes in gene cluster A are connected with characteristics in a positive way, whereas genes in gene cluster B are correlated with traits in a negative one. The Kaplan-Meier curve was used to compare the overall survival between the two gene clusters. There was no significant difference in the overall survival between the two gene clusters (P=0.441) ( Supplementary Figure 3A ). However, when comparing the median survival time between the two, we found that the median survival time was higher for gene cluster A. We performed GO enrichment analysis on two gene clusters, and it showed that gene cluster A was mainly enriched in immune cell activation and immune signal transmission pathway, and it was an important part of the cell transmembrane signal transmission process ( Figure 2B ). In gene cluster B, it showed strong enrichment characteristics during the development and differentiation of the epidermis and matrix ( Figure 2C ). The infiltration level of 22 immune cells in two gene clusters was examined using the CIBERSORT algorithm to learn more about the potential roles of different gene clusters in the immunological milieu ( Figure 2D ). The majority of the immune cells in gene cluster B have been discovered to have a higher active degree of immune infiltration. This implied that gene cluster B was in a state of immune activation. In addition, in gene cluster B, immune-related genes CTLA4, PD-L1, and PD-1 all had higher levels of expression ( Figures 2E–G ), which indicates that immunotherapy can exert a better effect in gene cluster B.

We performed a PCA analysis to quantify the combined factors of immune cell infiltration in cervical cancer. Recorded the sum of the individual scores of each sample as the ICI scores of two gene clusters to construct an ICI scoring system based on ICI phenotypic characteristics. Sankey diagrams depict the pattern of various gene clusters, ICI scores, and patient survival rates for cervical cancer ( Figure 3A ). When the survival prognosis of the high and low ICI score groups was compared using the Kaplan-Meier curve, the overall survival of the high ICI score group was greater (P=0.002) ( Figure 3B ). Used GSEA to analyze the function of the KEGG cell pathway, immune cell activation, immune cell killing and antigen presentation, and other immune cell working pathways were the main enriched pathways in the high ICI score group; growth of stromal cells was the main enrichment pathway in the low ICI score group ( Figure 3C ). To more effectively assess the ICI score’s effectiveness in predicting the prognosis of cervical cancer, we combined different clinical traits and performed a correlation analysis of the clinical traits between different ICI score groups ( Supplementary Figure 3B ), it showed that for patients with different grades, the patients with the worse grade had higher ICI scores (P=0.04). However, when we examined the relationship between various age groups and ICI scores, we discovered no statistically significant difference (P=0.23) ( Supplementary Figure 3C ). We also analyzed the survival prognosis in these two ICI score groups with different clinical characteristics. In patients under 60 years of age, overall survival was higher in patients with high ICI scores (p=0.005) ( Figure 3D ), and the same findings were also seen in the G1-2 patient group with high ICI scores (P=0.024) ( Figure 3F ). However, in the group of people older than 60 years ( Figure 3E ) and the G3 stage group ( Figure 3G ), there were no significant differences between the two ICI score groups.

Somatic mutations often occur in nature; when the tumor mutation burden (TMB) increases, it can lead to the production of new antigens and induce anti-cancer responses in immune cells, which also provides opportunities for immunotherapy, high TMB suggests benefit for immunotherapy of tumors (35). Somatic mutations (SNVs) data for cervical cancer were obtained from the TCGA database to investigate the potential relationship between tumor mutation burden in the immune microenvironment of cervical cancer and ICI scores. Exploring the difference in TMB between the two ICI score groups revealed that the group with higher ICI score had a higher TMB index (p=0.0017) ( Supplementary Figure 4A ) and a positive correlation between the two variables (Spearman coefficient: R=0.19, P=0.0013) ( Supplementary Figure 4B ). In addition, overall survival was more satisfactory in the high mutation burden group compared to the low mutation burden group (P=0.042) ( Supplementary Figure 4C ). In order to do a stratified survival prognosis analysis, we integrated the ICI scores and TMB. The survival outcomes of these four different subgroups were significantly different (P=0.007) ( Supplementary Figure 4D ).

The “maftools” package was used to identify the top 20 genes most likely to be mutated in each ICI score group to determine the genes with the most significant mutation risk in each ICI score subgroup ( Supplementary Figures 4E, F ), including missense, nonsense. Additionally, we discussed a possible connection between TMB and ICI scores ( Table 1 ).

Immune checkpoint-related genes such as PDCD1(PD-1), GZMB, CD274(PD-L1), and CTLA4, as well as immune activation-related genes CD8A and CXCL10, were used as characteristic genes to assess the immune activity of the two ICI score groups. The expression levels of these genes were compared. The results revealed that the expression of all these genes was higher in the high ICI score group ( Figure 4A ). To investigate the potential correlation between ICI scores and immunotherapy, we included 16 cervical cancer patients who had received PD1, PDL1, and VEGF immunotherapy, and obtained cervical tissues from these patients prior to receiving immunotherapy. The Response Evaluation Criteriain Solid Tumours (RECIST) efficacy assessment method was used to evaluate the immunotherapy effect of these patients, among which 5 patients were sensitive to immunotherapy and 11 patients were resistant to immunotherapy. IHC results revealed that patients who responded well to immunotherapy had higher expressions of CD8A, CXCL10, and GZMB ( Figure 4B ). Next, we performed a further evaluation using mIHC, in which patient1 showed immunotherapy resistance and patient2 and patient3 showed immunotherapy sensitivity. It can be seen that in the immunotherapy-sensitive group, the expressions of CD45, CD3 and CD4 were all increased, in addition, the expressions of CD8A, CXCL10 and GZMB were also increased ( Figure 5 ).

We evaluated the IC50 of Cisplatin, 5-Fluorouracil, Vinblastine, Vincristine, Paclitaxel, Docetaxel, Gemcitabine, Cyclophosphamide, Topotecan, Epirubicin, Olaparib, Dactinomycin in cervical cancer, and found that the drug sensitivity was higher in the high ICI score group ( Supplementary Figure 5 ). This indicates that ICI score has great potential in predicting chemotherapy drug sensitivity.

Univariate COX analysis was used to filter important genes from 119 DEGs, and 21 genes associated with the prognosis of cervical cancer were eliminated. Then, three important genes were found using the multivariate COX analysis, namely CST7, IL1B, and ITGA5 ( Supplementary Figure 6A ). These three key genes were used to build a prognostic model ( Supplementary Figures 6B–D ), and the results of the 1, 3, 5year ROC curves created using the model demonstrated that they could make accurate predictions ( Supplementary Figures 6E, F ).

We then performed cellular experiments to investigate the role of CST7, IL1B and ITGA5 in cervical cancer cells. We treated CST7 with plasmid overexpression and verified the knockdown efficiency of IL1B and ITGA5 and the overexpression efficiency of CST7 by qRT-PCR and western blot. The findings demonstrated that si-IL1B and si-ITGA5 had substantial knockdown efficiencies, and that CST7 overexpression was clearly present ( Figures 6A–F ). According to Transwell data, the capacity of cervical cancer cells to migrate was dramatically reduced when IL1B and ITGA5 were knocked down and CST7 was overexpressed in comparison to the control group ( Figures 6G–L ), similarly, the invasive ability of cervical cancer cells was significantly down-regulated ( Figures 6M–R ). MTT results showed that knockdown of IL1B and ITGA5 and overexpression of CST7 significantly reduced the proliferation of cervical cancer cells ( Figures 7A–F). Additionally, the findings of the clonogenic assay demonstrated that overexpressing CST7 and knocking down IL1B and ITGA5 decreased the capacity of cervical cancer cells to proliferate ( Figures 7G–L ).

To explore the potential role of IL1b, CST7, and ITGA5 in the immune microenvironment, we downloaded single-cell data of cervical cancer in GEO to validate our results and explore a more precise transcriptional landscape. The data set GSE168652 includes two samples, one from a tumor and one from adjacent normal tissue. After filtering low-quality cells, a total of 24498 cells (Normal: 11394; Tumor: 13104) remained to participate in downstream analysis. We used t-distributed stochastic neighbor embedding (t-SNE) to visualize the global distribution of cells ( Figures 8A–C ). Tumor tissue and normal tissue showed heterogeneity, with a distinct distribution of cells of different origins ( Figure 8A ). On this basis, unsupervised cluster analysis identified 13 clusters, and then we divided the cells into 7 main cell types according to the differentially expressed genes of each cluster and the cell-specific marker genes provided by the original article, including epithelial cells, smooth muscle cells, fibroblasts, endothelial cells, endometrial stromal cells, t cells and macrophages ( Figures 8B, C ). Then we showed the expression levels of these three key genes, CTS7, IL1B, and ITGA5, in different types of cells on the map ( Figures 8D–G ). All three genes were indeed expressed in immune cells. In particular, CST7 is highly expressed in T cells and IL1B is highly expressed in macrophages. Interestingly, ITGA5 is more widely expressed in non-immune cells, such as endothelial and epithelial cells. However, the expression of ITGA5 is concentrated in the epithelial subgroup Epithelial_02, suggesting that there may be some underlying knowledge in this cluster. The analysis of cell communication also shows that there are frequent interactions among cells expressing these three genes ( Figures 8H, I ).

We use monocle2 to explore the temporal order of immune cells and epithelial cells. We further subdivided T cells into 5 T cell subtypes according to the differentially expressed genes and classical marker genes of T cell subsets, including cytotoxic CD8+ T cells, helper T cells, regulatory T cells, naive T cells, and central memory CD8+ T cells. Four nodes and nine branches appeared on the pseudo-trajectories of T cells. Consistent with the results of pseudo-trajectories analysis, naive T cells first appeared and mainly gathered in state1, differentiated into other T cell subtypes, besides, helper T cells and regulatory T cells interacted with cytotoxic T cells along pseudo-trajectories ( Figures 9A–G ). However, macrophages show three cell states and differentiate from one starting point to two directions ( Supplementary Figures 7A–E ), which is consistent with our known knowledge. Under the influence of various cytokines, macrophages are polarized into M1/M2 macrophages, M1 macrophages are mainly involved in pro-inflammatory responses, and M2 macrophages are mainly involved in anti-inflammatory responses (36–38). Epithelial cells in tumor tissues exhibit complex diversity. We extracted Epithelial_02, a subpopulation highly expressing ITGA5, to reconstruct differentiation trajectories ( Supplementary Figures 7H–J ). Epithelial_02 is a malignant tumor cell that highly expresses some genes that promote tumor migration, invasion, and metastasis, such as the S100 family and CEACAM family genes (39–41). The expression patterns of 3 key genes were different in the pseudo-trajectories of the 3 types of cells, but IL1B and ITGA5 were more similar ( Figures 9H–J ; Supplementary Figures 7F–G , Supplementary Figures 7K–L ), suggesting that CST7 and the other two work in a different mode of action in cancer.