
==== Front
Discov Oncol
Discov Oncol
Discover Oncology
2730-6011
Springer US New York

39230769
1265
10.1007/s12672-024-01265-w
Analysis
Molecular subtypes and nomogram for predicting the prognosis of cervical cancer based on a matrix-immune signature
Liao Yuanyuan 1
Huang Qidan 1
Shen Guqun 2
Muhanmode Yalikun 2
Luo Xiaolin 1
Li Fen 2
Wen Mengke 2
Liu Jihong liujh@sysucc.org.cn

1
Huang He huangh@sysucc.org.cn

1
1 https://ror.org/0400g8r85 grid.488530.2 0000 0004 1803 6191 Department of Gynecological Oncology, Sun Yat-sen University Cancer Center, Guangzhou, 510000 China
2 https://ror.org/015tqbb95 grid.459346.9 0000 0004 1758 0312 The Second Department of Gynecology, Affiliated Tumor Hospital of Xinjiang Medical University, Urumqi, 830011 China
4 9 2024
4 9 2024
12 2024
15 40523 11 2023
22 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Cervical cancer is a kind of tumor related to chronic HPV infection. Currently, the treatment of cervical cancer is guided mainly by clinicopathological factors. The role of tumor microenvironment in the prognosis and treatment of cervical cancer has been ignored. We aimed to use bioinformatics to identify the molecular subtypes in cervical cancer and construct a predictive nomogram combining a matrix-immune signature (MIS) and clinicopathological factors to support treatment decisions. Two cervical cancer subtypes with different prognoses were identified based on matrix- and immune-genes in TCGA-CESC. The MIS was developed using Cox regression and Lasso algorithm and verified in the Cancer Genome Characterization Initiative (CGCI) using time-dependent receiver operating characteristic (ROC) curve analysis. Multivariable analysis identified lymph node metastases, lymphovascular space invasion, and the MIS as independent prognostic factors, which were used to construct the predictive nomogram. The areas under the ROC curve of the model were 0.872, 0.879, and 0.803 for the 1-, 3-, and 5-year periods, respectively. The C-index was 0.845. Calibration curves confirmed the excellent prognosis prediction of the nomogram. The nomogram indicted a 3-year survival rate of > 90% in patients with a total score > 110.1. The constructed predictive nomogram has significant implications for prognostic assessment and treatment selection in cervical cancer.

Keywords

Cervical cancer
Tumor microenvironment
Molecular subtypes
Prognosis
Nomogram
Bioinformatics
http://dx.doi.org/10.13039/501100001809 National Natural Science Foundation of China 8216110322 8216110322 8216110322 8216110322 8216110322 8216110322 8216110322 8216110322 8216110322 Liao Yuanyuan Huang Qidan Shen Guqun Muhanmode Yalikun Luo Xiaolin Li Fen Wen Mengke Liu Jihong Huang He issue-copyright-statement© Springer Science+Business Media, LLC 2024
==== Body
pmcIntroduction

Cervical cancer is a kind of tumor related to chronic inflammation. Persistent infection with high risk human papilloma virus (HPV) is the major risk factor for cervical cancer and precancerous lesions [1]. Currently, the treatment of cervical cancer is guided mainly by clinicopathological factors. However, a proportion of patients undergo recurrence. For patients with recurrent or metastatic disease, the prognosis is poor [2]. The tumor microenvironment (TME) plays a very important role in the occurrence and development of cervical cancer. The TME contains various components, such as the extracellular matrix (ECM) and immune cells, which can interact with each other and affect the immune response [3]. The immune response within the established TME is an important factor determining the prognosis of cancer. Specific types of immune cells or molecules determine whether tumor promoting or anti-tumor immune responses dominate in the TME [4]. In addition, it was reported that the ECM remodeled by tumor cells forms an unfavorable microenvironment that prevents immune cell penetration or serves as a barrier to anti‐cancer agents, which is known as immunological destruction of advanced cold tumors [5]. Thus, it is critical to explore the status of the TME in cervical cancer to provide suggestions for immunotherapy and to predict patient prognosis.

For patients with programmed cell death 1 ligand 1 (PD-L1)-positive recurrent or metastatic cervical cancer, immune checkpoints inhibitors (ICIs) combined with chemotherapy showed fine effects in the KEYNOTE-826 clinical trial [6]. Recent studies found that a treatment response of ICIs was also observed in PD-L1-negative tumors [7, 8]. Therefore, it is not suitable that the expression of PD-L1 serves as the main indication for the use of ICIs.

In this study, we employed the Estimation of Stromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) algorithm and weighted gene correlation network analysis (WGCNA) to screen hub matrix- and immune-related prognostic genes in patients with cervical squamous cell carcinoma (CESC) from The Cancer Genome Atlas (TCGA) database. Two clusters were identified with different prognoses and TMEs. of A risk assessment model based on the differentially expressed genes between the clusters was generated and named as the matrix-immune signature (MIS). Its prognostic efficacy was assessed and validated in a Cancer Genome Characterization Initiative (CGCI) cohort. Afterwards, we integrated clinical characteristics and the MIS to develop a nomogram model, aiming to provide a reliable and efficient quantitative method to guide further treatment.

Materials and methods

CESC data collection and preprocessing

Gene expression sequencing data (counts and fragments per kilobase of transcript per million mapped reads (FPKM) values), clinical data, and mutation data of 307 patients with CESC from the TCGA database were downloaded from the UCSC XENA database (https://xena.ucsc.edu). After excluding patients who lacked clinical information or whose survival time was 0, 282 CESC samples were retained for further analysis. Validation data included clinical information and survival data of patients with CESC, obtained from the HIV + Tumor Molecular Characterization Project—Cervical Cancer (HTMCP—CC) database (https://portal.gdc.cancer.gov/). The immune-related gene set was retrieved from the IMMPORT database (https://www.immport.org/), while the ECM-related genes were collected from previous study [9]. The microarray data of the CESC tissues and normal cervical tissues were extracted from Gene Expression Omnibus (GSE63514) (https://pubmed.ncbi.nlm.nih.gov/).

WGCNA identification of ECM- and immune-related genes associated with CESC prognosis

We used the expression matrix of patients with CESC to calculate Immune scores and stromal scores based on ESTIMATE algorithm via the R package ‘‘estimate’’ [10]. WCGNA was performed to investigate co-expressed gene modules that correlated with immune/stromal scores. Modules with the maximal correlation coefficient were considered as key ECM- and immune-related gene modules [11]. Univariate Cox regression analysis was applied to screen the relationships between ECM- and immune-related genes associated with patient prognosis using the survival package of R.

Protein interaction and regulatory network analysis

The STRING protein–protein interaction database was used to construct a PPI network of ECM-/immune-related genes. An interaction score of 0.9 was regarded as the cut-off criterion and the PPI was visualized. The top 10 core hub genes were further explored through the CytoHubba app in the Cytoscape software [12].

Consensus clustering analysis

Based on the 19 identified prognostic ECM- and immune-related genes, consensus clustering analysis was carried out to identify different CESC microenvironment patterns. The R package “ConsensuClusterPlus” was utilized to measure the similarity between and within each group via the Euclidean distance with 50 repetitions [13]. The optimum cluster number (k) and level of consensus stability was determined according to the CDF plots and the CDF curve, respectively. OS difference between the clusters was obtained using the R package “survival”. The DEGs between the clusters were estimated using the “DESeq2” package in R according to the following criteria: log fold change (logFC) > 1, and an adjusted P-value (adj.P-value) < 0.05 [14].The DEGs were then subjected to gene ontology GO analysis. The Chi-squared test or Fisher’s exact test were used to compare and analyze the statistical significance between two clusters in terms of their clinicopathological features.

Pathway enrichment analysis and the TME landscape

The CIBERSORT algorithm was adopted to estimate the relative fraction of 22 immune cell types [15]. Single-sample gene-set enrichment analysis (ssGSEA) was performed to quantify the enrichment levels of immune signatures using R package “GSVA” [16].

Single nucleotide polymorphism (SNP) analysis and prediction of ICI therapy response

To analyze SNPs and calculate the TMB, we used the maftools package to analyze frequently mutated genes in the patients in the two clusters [17]. The potential ICI response was predicted using the TIDE algorithm [18]. In addition, we compared the differential expression of immune checkpoints in the two clusters. This was also used to infer the treatment response to ICIs.

Development and validation of a CESC prognostic signature

Univariate Cox regression analysis was carried out to estimate the prognostic significance of the DEGs identified between the clusters, with a threshold of P < 0.05. Using the R package “glmnet”, these DEGs were subjected to LASSO-penalized Cox regression analysis to construct the prognostic model by narrowing down the candidate genes [19]. Then, we used a multivariate Cox proportional hazards regression model to identify further independent prognostic factors and establish the prognostic MIS. The formula to calculate the risk score of the risk model was as follows: MIS = ∑ n coef(i) ∗ mRNA(i)expression. ROC curves for 1, 3, and 5 years were plotted as time-dependent ROC curves to determine their accuracy [20]. In addition, the model was tested using the external test set (CGCI-HTMCP-CC) according to the regression coefficients of the genes in the model, and a further time-dependent ROC curve was drawn for validation.

Construction of a clinical prediction model based on the MIS

To demonstrate the independent predictive value, we used univariate and multivariate Cox analyses to estimate the predictive power of the MIS combined with clinicopathological features of patients with regard to OS. According to the clinical significance and statistical value, the clinical prediction nomogram was constructed. The discrimination of the nomogram was assessed using the C-index. Time-dependent ROC curves for 1, 3, and 5 years were plotted. A calibration curve was generated to assess the performance of the nomogram by comparing the predicted values of the nomogram with the actual observed survival data. We used the R package “pRRophetic” to estimate the chemotherapeutic response determined by the half maximal inhibitory concentration (IC50) of each patient with CESC on the Genomics of Drug Sensitivity in Cancer (GDSC) website [21].

Statistics

All statistical analysis was performed in RStudio (v4.1.3) and GraphPad Prism (v9). Values were presented as mean ± SEM. Compare the differences between the groups using T-test and Wilcoxon test. *p < 0.05; **p < 0.01, ***p < 0.001 were defined to be statistically significant and p > 0.05 was ns. (non-significant).

Results

Identification of hub matrix- and immune-related genes in cervical cancer

The expression profiles of 282 patients with CESC were obtained from TCGA, and 1027 matrix-related genes and 2483 immune-related genes were acquired from public databases. We used ESTIMATE techniques to calculate the ImmuneScore and StromalScore for each patient. To examine clusters of co-expressed matrix-related genes that were significantly related to the StromalScore, we used WGCNA, which identified a total of six modules, each of which was labeled with a distinct color. Among them, the turquoise module was most associated with the StromalScore and contained 446 matrix-related genes (Fig. 1A, B).Univariate Cox regression analysis was performed using these genes, which identified 69 prognostic genes that were filtered into a protein–protein interaction (PPI) network. Utilizing the Molecular Complex Detection Algorithm (MCODE), six significant modules from the PPI network complex were discovered. The top 10 hub gene cluster was obtained using Cytohubba in Cytoscape (Fig. 1C). We identified 10 hub immune-related genes in the same way. Among them, TNF (tumor necrosis factor) was not only a matrix-related gene, but also an immune related gene. Finally, we obtained 19 hub genes.Fig. 1 A Cluster trees of co-expressed matrix-related genes that were significantly related to the StromalScore. B Division of gene modules, heat map of the relationship between building blocks and characteristics. C Hub genes identified using Cytoscape and CytoHubba

Consensus clustering reveals two distinct tumor subgroups in cervical cancer

To delineate the immune and matrix heterogeneity within cervical cancer, the 19 hub genes were subjected to unsupervised consensus clustering analysis. The cumulative distribution function (CDF) curves of consensus matrix indicated that when k = 2, the interference between subgroups was minimal and the distinction was significant (Fig. 2A–C). The heatmap showed the different expressions of matrix- and immune-related genes and clinicopathological characteristics between Cluster 1 and Cluster 2 (Fig. 2D). We found that matrix-related genes were highly expressed in Cluster 1 and lowly expressed in Cluster 2. The expression of immune-related genes showed the opposite result. There was no statistically significant difference in distribution of clinicopathological features between the two clusters (Table 1). Most notably, the Cluster 1 subtype had a worse overall survival (OS) compared with the Cluster 2 subtype (P < 0.05) (Fig. 3A) and was significantly associated with a reduced disease control rate (DCR) (Fig. 3B).Fig. 2 A Consensus clustering matrix displaying the two CESC sample clusters with k = 2. B Cumulative distributive function for k = 2–9. C Delta graph displaying the change in the CDF curve’s area from k = 2–9. D Heat-map showing the 19 hub genes expression profiles in Cluster 1 and Cluster 2, and the associations between clinicopathological characteristics. CESC cervical squamous cell carcinoma and endocervical adenocarcinoma, CDF cumulative distribution function, N node

Table 1 Clinicopathological characteristics for two CESC sample clusters

Characteristic	Levels	Overall	Group			
			C1	C2	P-value	
n		282	156	126		
Age, n (%)	 ≤ 60	230 (81.6%)	133 (57.8%)	97 (42.2%)	0.0749	
	 > 60	52 (18.4%)	23 (44.2%)	29 (45.8%)		
Stage, n (%)	NA	6 (2.1%)	4 (66.7%)	2 (33.3%)	0.7047	
	Stage I	154 (54.6%)	82 (53.2%)	72 (46.8%)		
	Stage II	63 (22.3%)	33 (52.4%)	30 (47.6%)		
	Stage III	39 (13.8%)	24 (61.5%)	15 (38.5%)		
	Stage IV	20 (7.1%)	13 (0.65%)	7 (0.35%)		
Grade, n (%)	NA	27 (9.6%)	17 (63%)	10 (37%)	0.583	
	G1	16 (5.7%)	8 (50%)	8 (50%)		
	G2	124 (44%)	72 (58.1%)	52 (41.9%)		
	G3	115 (40.8%)	59 (51.3%)	56 (48.7%)		
Node, n (%)	NA	103 (36.5%)	73 (70.9%)	30 (29.1%)	0.7537	
	N0	125 (44.3%)	57 (45.6%)	68 (54.4%)		
	N1	54 (19.1%)	26 (48.1%)	28 (51.9%)		
LVSI, n (%)	NA	137 (48.7%)	82 (59.9%)	55 (40.1%)	0.666	
	Absent	68 (24.1%)	36 (52.9%)	32 (47.1%)		
	Present	77 (27.3%)	38 (49.4%)	39 (50.6%)		
CESC cervical squamous cell carcinoma and endocervical adenocarcinoma, LVSI Lymphovascular space invasion, NA Not Available

Fig. 3 A Kaplan–Meier survival curves between the two CESC sample clusters based on overall survival. B Stacked bar chart showed the number of patients with progressive disease (PD) and the total number of patients with complete response (CR), partial response (PR) and stable disease (SD) in the two CESC sample clusters. P values in the stacked bar chart were calculated using the chi-squared test. C Gene set variation analysis (GSVA) between Cluster 1 and Cluster 2. KEGG, Kyoto Encyclopedia of Genes and Genome

Gene set variation analysis (GSVA) of biological pathways between C1 and C2 clusters

To explore the biological behaviors among C1 and C2 clusters, we performed GSVA enrichment analysis, and the results were shown in a heatmap. Cluster 1 was markedly enriched in stromal and carcinogenic activation pathways, such as ECM receptor interaction, the TGF beta signaling pathway, and gap junctions. Cluster 2 presented enrichment pathways associated with full immune activation including antigen processing and presentation, natural killer cell mediated cytotoxicity, the B cell receptor signaling pathway, and the T cell receptor signaling pathway (Fig. 3C).

Differences in the tumor immune microenvironment between C1 and C2 clusters

We analyzed the abundance of immune cells between Cluster 1 and Cluster 2 using the CIBERSORT algorithm. There was a lower abundance of immune cells in Cluster 1, including dendritic cells, B cells, CD4 + T cells, and macrophages, compared with those in Cluster 2 (Fig. 4A).Fig. 4 A Box plots showing variation in immune cell infiltration among clusters. B The expression of immune checkpoint inhibitors (ICIs) in Cluster 1 and Cluster 2. C Violin plots showing variation in Tumor Immune Dysfunction and Exclusion (TIDE) score among clusters. D Violin plots showing the tumor mutation burden (TMB) of the patients among the clusters

Prediction of the immune therapy response between C1 and C2 clusters

The expression of immune checkpoints (ICs) in the two clusters were analyzed, and the results showed that compared with Cluster 1, most immune checkpoints were highly expressed in Cluster 2, such as programmed cell death 1 (PDCD1), cytotoxic T-lymphocyte associated protein 4 (CTLA4), lymphocyte activating 3 (LAG3), Hepatitis A virus cellular receptor 2 (HAVCR2), and inducible T cell costimulator (ICOS) (Fig. 4B). We also used tumor mutation burden (TMB) and the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to predict immunotherapeutic response between the two clusters. The results demonstrated that patients in Cluster 1 had a higher TIDE score and lower TMB than those in Cluster 2 (Fig. 4C, D). These findings suggested that patients in the Cluster 2 subgroup might be more sensitive to immunotherapy and obtain better therapeutic effects.

Gene mutations in the two clusters

We summarized the incidence of copy number variations and somatic mutations between the two clusters. Among 148 samples, 125 experienced mutations in Cluster 1, with a frequency of 84.46%. The gene mutation frequency in Cluster 2 was 91.45%. It was found that TNN (encoding tenascin N) exhibited the highest mutation frequency, followed by PIK3CA (encoding phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha) in both Cluster 1 and Cluster 2 (Fig. 5).Fig. 5 A, B Waterfall plots of tumor somatic mutations established by Cluster 1 (A) and Cluster 2 (B)

Differential expression analysis and gene ontology (GO) enrichment analysis

We performed differential expression analysis among two clusters and identified 1230 significant differentially expressed genes (DEGs) (Fig. 6A). GO analysis showed enrichment of biological processes, such as regulation of cell adhesion, cytokine production, and T cell activation. The cellular component was enriched in collagen-containing extracellular matrix, external side of plasma membrane, and endoplasmic reticulum lumen. The molecular functions were enriched in extracellular matrix structural constituent, immune receptor activity, and signaling receptor activator activity (Fig. 6D).Fig. 6 A Differential gene expression in Cluster 1 and Cluster 2. B, C Least absolute shrinkage and selection operator (LASSO) regression analysis, and the number of variables corresponding to the optimal λ value was 3. D Gene ontology (GO) analysis of differentially expressed genes (DEGs) in the two clusters. E Multivariate Cox stepwise regression analysis showing that four genes are independent prognostic factors

MIS based on DEGs and testing in the CGCI data

First, least absolute shrinkage and selection operator (LASSO) and multivariate Cox analyses of 32 prognostic DEGs were used to seek the optimum prognostic signature (Fig. 6B, C). We obtained four genes to construct a risk scoring system in the training sets using multivariate Cox regression analysis: MIS = (− 0.145406843672574 * expression of ITGA5) + (− 0.100834783610163 * expression of CXCL8) + (0.108030183332151* expression of CD6) + (0.0800602829915149* expression of SULT1E1) (Fig. 6E; ITGA5 (encoding integrin subunit alpha 4), CXCL8 (encoding C-X-C motif chemokine ligand 8), CD6 (encoding CD6 molecule), and SULT1E1 (encoding sulfotransferase family 1E member 1). Next, the levels of the four risk genes were measured in CESC tissues and normal cervical tissues in the Gene Expression Omnibus (GSE63514) database. Compared with those in the corresponding normal tissues, the expression levels of ITGA5 and CXCL8 were upregulated in tumor tissues (p < 0.05), and the expression levels of SULT1E1 and CD6 was downregulated (p < 0.05) (Fig. 7A). Additionally, the 1-, 2-, and 3-year survival rates of the MIS in the training sets were represented by AUC values of 0.79, 0.814, and 0.8 (Fig. 7B). In testing sets, there are only 1- and 2-year survival rates, which were represented by AUC values of 0.665 and 0.802 (Fig. 7C). Overall, the MIS showed moderate sensitivity and specificity.Fig. 7 A The expression of ITGA5, SULT1E1, CXCL8, and CD6 in the normal group (n = 24) and tumor group (n = 28). B Time‐dependent receiver operating characteristic (ROC) curve of the training set (The Cancer Genome Atlas (TCGA)-CESC. C Time‐dependent ROC curve of the testing set from the Cancer Genome Characterization Initiative (CGCI)

Construction and validation of the nomogram

To better predict the survival of patients with CESC, we constructed a nomogram that integrated the prognostic MSI model and clinical determinants. Considering the clinicopathological features and MIS, we developed a multivariate Cox regression risk forest plot for the overall cohort, which revealed that lymph node metastasis, lymphovascular space invasion (LVSI), and MIS were independent factors that affected the prognosis in CESC (Table 2). As shown in Fig. 8A, with each item assigned a score based on the actual condition, patients could get a total score to predict their survival rate at 1 y, 3 y, and 5 y. As the total score increased, the survival probability increased. The C-index was 0.845 (95% confidence interval (CI) 0.706–0.984). The applicable AUC values of the receiver operating characteristic (ROC) curve to predict the 1-, 3-, 5- year survival rates were 0.872, 0.879, and 0.803, respectively (Fig. 8B). The calibration curve showed that there was a slight difference between the actual risk and the anticipated risk (Fig. 8C). These results indicated that the nomogram model has good accuracy to predict the survival of patients with CESC. Table 2 Univariate and multivariate Cox analysis

Characters	Total (N)	Univariate analysis		Multivariate analysis		
		Hazard ratio (95% CI)	P Value	Hazard ratio (95% CI)	P Value	
Age	282	1.8 (1–3)	0.035	1.4 (0.45–4.1)	0.59	
 ≤ 60	188					
 > 60	94					
Grade	256	0.87 (0.51–1.5)	0.61			
1 ~ 2	140					
3	116					
Stage	276	2.4 (1.5–4)	 < 0.001	0.58 (0.16–2.1)	0.41	
I ~ II	217					
III ~ IV	59					
N	179	2.8 (1.4–5.6)	0.0048	2.4 (1–5.9)	0.047	
0	125					
1	54					
LVSI	145	9.1 (2.1–39)	0.003	10 (1.9–54)	0.0067	
Absent	68					
Present	77					
Histological	282	1 (0.66–1.5)	0.96			
Adenosquamous	29					
Squamous cell carcinoma	233					
Others	20					
MIS	282	1.3 (1.2–1.3)	 < 0.001	1.5 (1.2–1.9)	 < 0.001	
CI, confidence interval, LVSI Lymphovascular space invasion, MIS matrix-immune signature

Fig. 8 A Nomogram integrating the matrix-immune signature (MIS) and clinical characteristics. B Time‐dependent ROC curve of the nomogram. (C) Calibration of the nomogram for 1‐y, 3‐y, and 5‐y survival. LVSI, lymphovascular space invasion; DFS, disease free survival, OS, overall survival

We regarded the total score with a 3-year survival rate of 90% as the cutoff point to divide patients into high- and low-risk groups. Kaplan–Meier (K-M) analysis showed a significant difference in survival between the two groups (P < 0.001) (Fig. 9A). For patients receiving surgical treatment, we also found that patients with lower scores (Sgroup1) had poor prognosis (Fig. 9B).Fig. 9 A Kaplan–Meier survival curves between the low- and high-risk groups based on overall survival. B Kaplan–Meier survival curves between Sgroup1 and Sgroup2 based on overall survival. C Differential chemotherapeutic responses in low- and high-risk groups of patients. IC50, half maximal inhibitory concentration

We also investigated the response to chemotherapy in high- and low-risk patients with CESC, and found that the estimated IC50 values of Doxorubicin, Cytarabine, AUY922, BMS-509744, BMS-754807, Obatoclax Mesylate, OSU-03012, and PAC-1 were significant different between the high- and low-risk patients with CESC: high-risk patients with CESC showed decreased sensitivity to the eight chemotherapies (Fig. 9C).

Discussion

Cervical cancer is a tumor associated with chronic HPV infection [1]. The TME is a complicated system and is increasingly recognized to regulate tumor cell proliferation, migration, and immune escape ability, thereby promoting the occurrence and development of tumors in multiple cancer types [22]. The aberrant TME, especially immature and leaky vessels, prevents the penetration and accumulation of chemotherapeutics and results in the failure of chemotherapy to treat gynecologic cancer [23]. Additionally, cancer associated fibroblasts (CAFs) are the most common cells in the TME. Excessive secretion of ECM proteins by activated CAFs can cause tumor-associated fibrosis, forming a stromal barrier against ICI and immune cell penetration [24].

Immune cells and the surrounding ECM, as two major components in the TME, play an important role in tumor progression and influence the treatment effect of immune therapy by interacting with each other or tumor cells [25]. To explore the TME characteristics of cervical cancers, we screened 19 hub matrix- and immune-related genes and used them to identify two subgroups with diverse prognosis, immune, and molecular features. We evaluated the treatment response of immune therapy in patients with cervical cancer by calculating the TIDE score of each subgroup and supposed that the matrix-enrich cluster, with a higher TIDE score, was not sensitive to ICIs. This suggested that integrate matrix- and immune-related genes could be used to conduct molecular classification of cervical cancer to provide a predictive advantage in precision immunotherapy and prognosis for patients.

We then performed differential expression analysis on the two clusters and identified a four gene signature (ITGA5, CXCL8, SULT1E1, and CD6) to establish a risk assessment signature, which exhibited favorable predictive accuracy at different time points in time-dependent AUC. Meanwhile, we used an external database to verify that patients with higher scores are associated with poorer prognosis. Among the genes, ITGA5 and CXCL8 are high risk genes. ITGA5, belonging to the integrin alpha chain family, is significantly overexpressed in various tumors and participates in tumor progression by regulating tumor cell growth, migration, and invasion [26–28]. ITGA5 contributes to drug resistance by promoting vasculogenic mimicry formation in hepatocellular carcinoma [29]. Kyung et al. performed RNA sequencing in patients with locally advanced cervical cancer treated with platinum-based Concurrent chemoradiotherapy (CCRT) and found that ITGA5 was upregulated in the poor responders [30]. CXCL8 has been demonstrated as a key inflammatory factor in the TME. CXCL8 secreted from cancer‐associated stromal cells promotes metastasis and chemoresistance in multiple types of cancer. Recent studies found high systemic CXCL8 correlates with reduced clinical benefit of PD-L1 blockade. Blocking CXCL8 can reverse the immunosuppressive microenvironment of tumors and enhance the therapeutic effect of ICIs [31–33]. More effort should be made to explore the drug resistance mechanisms of ITGA5 and CXCL8 in cervical cancer.

Furthermore, lymph node metastasis, LVSI, and MIS were identified as independent prognostic factors by multivariate Cox analysis and these factors were chosen to establish a nomogram with a 3-year AUC of 0.879, which can display the impact of various factors on CESC prognosis intuitively. The calibration curves also suggested that this model had performed well for prognosis prediction. The nomogram model incorporating clinicopathological factors displayed a significantly higher AUC value, which improved the prediction of CESC prognosis compared with using the MIS alone. Among the TCGA database, patients with a score of over 110 had a 3-year survival rate of over 90%. Similar results were obtained from the subgroup analyses of patients who had accepted surgical treatment with early stage disease. These results indicated that the prognosis of patients with early stage disease is not only decided by clinicopathological factors, but also is dependent on the TME. Controversy remains as to whether to give additional adjuvant chemotherapy following radiotherapy for patients with CESC with intermediate-risk factors after surgery. Components such as chemokines, immune cells, and stromal cells in the TME are associated with efficacy of ICI therapies [34]. It might be better to combine clinicopathological features and TME characteristics to provide guidance as to whether to use ICIs for postoperative and advanced cervical cancer.

Limitation

Our study faces certain limitations that should be acknowledged. Primarily, our analysis is retrospective, based on pre-existing data, and its outcomes have not been validated through further clinical trials. To enhance the credibility of our findings, there is a pressing need for further efforts, including the collection of fresh clinical samples and the execution of detailed molecular studies to solidify our findings through cohort research. Moreover, a detailed investigation into specific genes such as ITGA5 and CXCL8 is warranted, which may unveil potential mechanisms of occurrence and drug resistance in cervical cancer.

Conclusion

Through bioinformatics, we identified two CESC subtypes with distinct tumor microenvironments and prognosis. We developed and initially validated a MIS associated with the prognosis of CESC patients. Moreover, by integrating the MIS with clinicopathological factors, we constructed a nomogram model that might offer insights for prognostic assessment and postoperative adjuvant treatment. Further studies on gene and protein expression are required to validate the model before its clinical application.

Acknowledgements

Not applicable.

Author contributions

YL and QH are co-first authors and contributed equally to this work. HH and JL are corresponding authors and contributed equally to this work. YL and QH contributed to the conception and design of the study. YL and QH performed the data analysis and wrote the manuscript. GS and YM contributed to the data acquisition and assessment. XL contributed to formal analysis. FL and MW contributed to prepare the figures and tables. HH and JL designed the study and revised the manuscript. All authors reviewed and approved the final manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, 8216110322.

Data availability

The data analyzed in current study can be download from related websites.

Code availability

Not applicable.

Ethics approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Yuanyuan Liao and Qidan Huang have attribute equally to this work.
==== Refs
References

1. Kjær SK Frederiksen K Munk C Iftner T Long-term absolute risk of cervical intraepithelial neoplasia grade 3 or worse following human papillomavirus infection: role of persistence J Natl Cancer Inst 2010 102 19 1478 1488 10.1093/jnci/djq356 20841605
Kjær SK, Frederiksen K, Munk C, Iftner T. Long-term absolute risk of cervical intraepithelial neoplasia grade 3 or worse following human papillomavirus infection: role of persistence. J Natl Cancer Inst. 2010;102(19):1478–88.20841605 10.1093/jnci/djq356
2. Siegel RL Miller KD Wagle NS Jemal A Cancer statistics CA Cancer J Clin 2023 73 1 17 48 10.3322/caac.21763 36633525
Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics. CA Cancer J Clin. 2023;73(1):17–48.36633525 10.3322/caac.21763
3. Bejarano L Jordāo M Joyce JA Therapeutic targeting of the tumor microenvironment Cancer Discov 2021 11 4 933 959 10.1158/2159-8290.CD-20-1808 33811125
Bejarano L, Jordāo M, Joyce JA. Therapeutic targeting of the tumor microenvironment. Cancer Discov. 2021;11(4):933–59.33811125 10.1158/2159-8290.CD-20-1808
4. Ozga AJ Chow MT Luster AD Chemokines and the immune response to cancer Immunity 2021 54 5 859 874 10.1016/j.immuni.2021.01.012 33838745
Ozga AJ, Chow MT, Luster AD. Chemokines and the immune response to cancer. Immunity. 2021;54(5):859–74.33838745 10.1016/j.immuni.2021.01.012
5. Yuan Z Li Y Zhang S Wang X Dou H Yu X Zhang Z Yang S Xiao M Extracellular matrix remodeling in tumor progression and immune escape: from mechanisms to treatments Mol Cancer 2023 22 1 48 10.1186/s12943-023-01744-8 36906534
Yuan Z, Li Y, Zhang S, Wang X, Dou H, Yu X, Zhang Z, Yang S, Xiao M. Extracellular matrix remodeling in tumor progression and immune escape: from mechanisms to treatments. Mol Cancer. 2023;22(1):48.36906534 10.1186/s12943-023-01744-8
6. Colombo N Dubot C Lorusso D Caceres MV Hasegawa K Shapira-Frommer R Tewari KS Salman P Hoyos UE Yañez E Pembrolizumab for persistent, recurrent, or metastatic cervical cancer N Engl J Med 2021 385 20 1856 1867 10.1056/NEJMoa2112435 34534429
Colombo N, Dubot C, Lorusso D, Caceres MV, Hasegawa K, Shapira-Frommer R, Tewari KS, Salman P, Hoyos UE, Yañez E, et al. Pembrolizumab for persistent, recurrent, or metastatic cervical cancer. N Engl J Med. 2021;385(20):1856–67.34534429 10.1056/NEJMoa2112435
7. Feun LG Li YY Wu C Wangpaichitr M Jones PD Richman SP Madrazo B Kwon D Garcia-Buitrago M Martin P Phase 2 study of pembrolizumab and circulating biomarkers to predict anticancer response in advanced, unresectable hepatocellular carcinoma Cancer 2019 125 20 3603 3614 10.1002/cncr.32339 31251403
Feun LG, Li YY, Wu C, Wangpaichitr M, Jones PD, Richman SP, Madrazo B, Kwon D, Garcia-Buitrago M, Martin P, et al. Phase 2 study of pembrolizumab and circulating biomarkers to predict anticancer response in advanced, unresectable hepatocellular carcinoma. Cancer. 2019;125(20):3603–14.31251403 10.1002/cncr.32339
8. Fujita T Amano H Nakamura M Hirano S Nakamura S Remarkable response to immune checkpoint inhibitor monotherapy in an EGFR-mutant pulmonary adenocarcinoma patient with 0% expression of PD-L1 J Thorac Oncol 2023 18 9 e93 e94 10.1016/j.jtho.2023.05.025 37599053
Fujita T, Amano H, Nakamura M, Hirano S, Nakamura S. Remarkable response to immune checkpoint inhibitor monotherapy in an EGFR-mutant pulmonary adenocarcinoma patient with 0% expression of PD-L1. J Thorac Oncol. 2023;18(9):e93–4.37599053 10.1016/j.jtho.2023.05.025
9. Naba A Clauser KR Hoersch S Liu H Carr SA Hynes RO The matrisome: in silico definition and in vivo characterization by proteomics of normal and tumor extracellular matrices Mol Cell Proteomics 2012 11 4 M111014647 10.1074/mcp.M111.014647
Naba A, Clauser KR, Hoersch S, Liu H, Carr SA, Hynes RO. The matrisome: in silico definition and in vivo characterization by proteomics of normal and tumor extracellular matrices. Mol Cell Proteomics. 2012;11(4):M111014647.10.1074/mcp.M111.014647
10. Yoshihara K Shahmoradgoli M Martínez E Vegesna R Kim H Torres-Garcia W Treviño V Shen H Laird PW Levine DA Inferring tumour purity and stromal and immune cell admixture from expression data Nat Commun 2013 10.1038/ncomms3612 24113773
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013. 10.1038/ncomms3612.24113773 10.1038/ncomms3612
11. Langfelder P Horvath S WGCNA: an R package for weighted correlation network analysis BMC Bioinform 2008 10.1186/1471-2105-9-559
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008. 10.1186/1471-2105-9-559.10.1186/1471-2105-9-559
12. Chin CH Chen SH Wu HH Ho CW Ko MT Lin CY cytoHubba: identifying hub objects and sub-networks from complex interactome BMC Syst Biol 2014 8 4 S11 10.1186/1752-0509-8-S4-S11 25521941
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(4):S11.25521941 10.1186/1752-0509-8-S4-S11
13. Wilkerson MD Hayes DN ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking Bioinformatics 2010 26 12 1572 1573 10.1093/bioinformatics/btq170 20427518
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.20427518 10.1093/bioinformatics/btq170
14. Love MI Huber W Anders S Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 Genome Biol 2014 15 12 550 10.1186/s13059-014-0550-8 25516281
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.25516281 10.1186/s13059-014-0550-8
15. Li T Fu J Zeng Z Cohen D Li J Chen Q Li B Liu XS TIMER2.0 for analysis of tumor-infiltrating immune cells Nucleic Acids Res 2020 48 W1 W509 W514 10.1093/nar/gkaa407 32442275
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–14.32442275 10.1093/nar/gkaa407
16. Hänzelmann S Castelo R Guinney J GSVA: gene set variation analysis for microarray and RNA-seq data BMC Bioinform 2013 10.1186/1471-2105-14-7
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013. 10.1186/1471-2105-14-7.10.1186/1471-2105-14-7
17. Mayakonda A Lin DC Assenov Y Plass C Koeffler HP Maftools: efficient and comprehensive analysis of somatic variants in cancer Genome Res 2018 28 11 1747 1756 10.1101/gr.239244.118 30341162
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.30341162 10.1101/gr.239244.118
18. Jiang P Gu S Pan D Fu J Sahu A Hu X Li Z Traugh N Bu X Li B Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response Nat Med 2018 24 10 1550 1558 10.1038/s41591-018-0136-1 30127393
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8.30127393 10.1038/s41591-018-0136-1
19. Friedman J Hastie T Tibshirani R Regularization paths for generalized linear models via coordinate descent J Stat Softw 2010 33 1 1 22 10.18637/jss.v033.i01 20808728
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22.20808728 10.18637/jss.v033.i01
20. Heagerty PJ Lumley T Pepe MS Time-dependent ROC curves for censored survival data and a diagnostic marker Biometrics 2000 56 2 337 344 10.1111/j.0006-341X.2000.00337.x 10877287
Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.10877287 10.1111/j.0006-341X.2000.00337.x
21. Geeleher P Cox NJ Huang RS Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines Genome Biol 2014 15 3 R47 10.1186/gb-2014-15-3-r47 24580837
Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. 2014;15(3):R47.24580837 10.1186/gb-2014-15-3-r47
22. de Visser KE Joyce JA The evolving tumor microenvironment: from cancer initiation to metastatic outgrowth Cancer Cell 2023 41 3 374 403 10.1016/j.ccell.2023.02.016 36917948
de Visser KE, Joyce JA. The evolving tumor microenvironment: from cancer initiation to metastatic outgrowth. Cancer Cell. 2023;41(3):374–403.36917948 10.1016/j.ccell.2023.02.016
23. Li M Wang Y Zhang L Liu Q Jiang F Hou W Wang Y Fang H Zhang Y Cancer cell membrane-enveloped dexamethasone normalizes the tumor microenvironment and enhances gynecologic cancer chemotherapy ACS Nano 2023 17 17 16703 16714 10.1021/acsnano.3c03013 37603464
Li M, Wang Y, Zhang L, Liu Q, Jiang F, Hou W, Wang Y, Fang H, Zhang Y. Cancer cell membrane-enveloped dexamethasone normalizes the tumor microenvironment and enhances gynecologic cancer chemotherapy. ACS Nano. 2023;17(17):16703–14.37603464 10.1021/acsnano.3c03013
24. Luo H Xia X Huang LB An H Cao M Kim GD Chen HN Zhang WH Shu Y Kong X Pan-cancer single-cell analysis reveals the heterogeneity and plasticity of cancer-associated fibroblasts in the tumor microenvironment Nat Commun 2022 13 1 6619 10.1038/s41467-022-34395-2 36333338
Luo H, Xia X, Huang LB, An H, Cao M, Kim GD, Chen HN, Zhang WH, Shu Y, Kong X, et al. Pan-cancer single-cell analysis reveals the heterogeneity and plasticity of cancer-associated fibroblasts in the tumor microenvironment. Nat Commun. 2022;13(1):6619.36333338 10.1038/s41467-022-34395-2
25. Cox TR The matrix in cancer Nat Rev Cancer 2021 21 4 217 238 10.1038/s41568-020-00329-7 33589810
Cox TR. The matrix in cancer. Nat Rev Cancer. 2021;21(4):217–38.33589810 10.1038/s41568-020-00329-7
26. Wang T Yang J Mao J Zhu L Luo X Cheng C Zhang L ITGA5 inhibition in pancreatic stellate cells re-educates the in vitro tumor-stromal crosstalk Med Oncol 2022 40 1 39 10.1007/s12032-022-01902-w 36469173
Wang T, Yang J, Mao J, Zhu L, Luo X, Cheng C, Zhang L. ITGA5 inhibition in pancreatic stellate cells re-educates the in vitro tumor-stromal crosstalk. Med Oncol. 2022;40(1):39.36469173 10.1007/s12032-022-01902-w
27. Wang JF Chen YY Zhang SW Zhao K Qiu Y Wang Y Wang JC Yu Z Li BP Wang Z ITGA5 promotes tumor progression through the activation of the FAK/AKT signaling pathway in human gastric cancer Oxid Med Cell Longev 2022 10.1155/2022/8611306 36620085
Wang JF, Chen YY, Zhang SW, Zhao K, Qiu Y, Wang Y, Wang JC, Yu Z, Li BP, Wang Z, et al. ITGA5 promotes tumor progression through the activation of the FAK/AKT signaling pathway in human gastric cancer. Oxid Med Cell Longev. 2022. 10.1155/2022/8611306.36620085 10.1155/2022/8611306
28. Li XQ Zhang R Lu H Yue XM Huang YF Extracellular vesicle-packaged CDH11 and ITGA5 induce the premetastatic niche for bone colonization of breast cancer cells Cancer Res 2022 82 8 1560 1574 10.1158/0008-5472.CAN-21-1331 35149589
Li XQ, Zhang R, Lu H, Yue XM, Huang YF. Extracellular vesicle-packaged CDH11 and ITGA5 induce the premetastatic niche for bone colonization of breast cancer cells. Cancer Res. 2022;82(8):1560–74.35149589 10.1158/0008-5472.CAN-21-1331
29. Shi Y Shang J Li Y Zhong D Zhang Z Yang Q Lai C Feng T Yao Y Huang X ITGA5 and ITGB1 contribute to Sorafenib resistance by promoting vasculogenic mimicry formation in hepatocellular carcinoma Cancer Med 2023 12 3 3786 3796 10.1002/cam4.5110 35946175
Shi Y, Shang J, Li Y, Zhong D, Zhang Z, Yang Q, Lai C, Feng T, Yao Y, Huang X. ITGA5 and ITGB1 contribute to Sorafenib resistance by promoting vasculogenic mimicry formation in hepatocellular carcinoma. Cancer Med. 2023;12(3):3786–96.35946175 10.1002/cam4.5110
30. Kim KH Chang JS Byun HK Kim YB A novel gene signature associated with poor response to chemoradiotherapy in patients with locally advanced cervical cancer J Gynecol Oncol 2022 33 1 e7 10.3802/jgo.2022.33.e7 34783210
Kim KH, Chang JS, Byun HK, Kim YB. A novel gene signature associated with poor response to chemoradiotherapy in patients with locally advanced cervical cancer. J Gynecol Oncol. 2022;33(1):e7.34783210 10.3802/jgo.2022.33.e7
31. Wang X Xu F Kou H Zheng Y Yang J Xu Z Fang Y Sun W Zhu S Jiang Q Stromal cell-derived small extracellular vesicles enhance radioresistance of prostate cancer cells via interleukin-8-induced autophagy J Extracell Vesicles 2023 12 7 e12342 10.1002/jev2.12342 37387557
Wang X, Xu F, Kou H, Zheng Y, Yang J, Xu Z, Fang Y, Sun W, Zhu S, Jiang Q, et al. Stromal cell-derived small extracellular vesicles enhance radioresistance of prostate cancer cells via interleukin-8-induced autophagy. J Extracell Vesicles. 2023;12(7):e12342.37387557 10.1002/jev2.12342
32. Walle T Kraske JA Liao B Lenoir B Timke C von Bohlen UHE Tran F Griebel P Albrecht D Ahmed A Radiotherapy orchestrates natural killer cell dependent antitumor immune responses through CXCL8 Sci Adv 2022 8 12 4050 10.1126/sciadv.abh4050
Walle T, Kraske JA, Liao B, Lenoir B, Timke C, von Bohlen UHE, Tran F, Griebel P, Albrecht D, Ahmed A, et al. Radiotherapy orchestrates natural killer cell dependent antitumor immune responses through CXCL8. Sci Adv. 2022;8(12):4050.10.1126/sciadv.abh4050
33. Liu H Zhao Q Tan L Wu X Huang R Zuo Y Chen L Yang J Zhang ZX Ruan W Neutralizing IL-8 potentiates immune checkpoint blockade efficacy for glioma Cancer Cell 2023 41 4 693 710.e8 10.1016/j.ccell.2023.03.004 36963400
Liu H, Zhao Q, Tan L, Wu X, Huang R, Zuo Y, Chen L, Yang J, Zhang ZX, Ruan W, et al. Neutralizing IL-8 potentiates immune checkpoint blockade efficacy for glioma. Cancer Cell. 2023;41(4):693-710.e8.36963400 10.1016/j.ccell.2023.03.004
34. Nagarsheth N Wicha MS Zou W Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy Nat Rev Immunol 2017 17 9 559 572 10.1038/nri.2017.49 28555670
Nagarsheth N, Wicha MS, Zou W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol. 2017;17(9):559–72.28555670 10.1038/nri.2017.49
