A bioinformatics analysis of KIF transcriptomic subtyping reveals prognostic heterogeneity and differential therapeutic responses in cervical cancer

Introduction

Cervical cancer (CC) remains one of the most prevalent malignancies among women worldwide, with particularly high incidence and mortality in low- and middle-income countries. Although the implementation of human papillomavirus (HPV) vaccination and screening programs has led to declining or stabilized incidence and mortality in numerous developed regions, CC continues to impose a significant global health burden, with ~660,000 new cases and 350,000 deaths annually (1–3). Identifying patient populations that can benefit from chemotherapy, targeted therapy, or immunotherapy is therefore critical for improving survival outcomes.

Precision medicine guided by molecular subtyping has become a cornerstone of oncology. For instance, patients with advanced triple-negative breast cancer have been stratified into molecular subgroups, with subtype-guided therapy resulting in improved clinical outcomes (4). Recent advances in multi-omics technologies have enabled comprehensive molecular characterization of CC. Integrated multi-omics and single-cell analyses have revealed that productive HPV integration correlates with elevated E6/E7 protein expression, enhanced tumor aggressiveness, immune evasion, and increased disease progression from carcinoma in situ to advanced stages (5). Nevertheless, despite these technological advances, effective molecular biomarkers and subtyping strategies specifically for CC remain limited, particularly those based on kinesin family genes. Kinesin family members (KIFs) are microtubule-associated motor proteins essential for intracellular transport, mitosis, and cellular signaling (6). In cancer, dysregulation of KIFs frequently promotes cell proliferation, chromosomal instability, invasion, metastasis, and therapy resistance. For example, overexpression of KIF2C drives tumor progression and is associated with poor prognosis across multiple cancers, including CC (7–9). Although individual KIFs have been investigated in CC, a comprehensive and systematic analysis of the entire KIF family in this context remains lacking.

Given their critical roles in cell division, intracellular transport and signaling pathways, it was hypothesized that the expression patterns of KIFs may reflect underlying tumor heterogeneity and serve as a basis for molecular classification with prognostic and therapeutic relevance in CC. To test this hypothesis, KIF expression profiles in CC were systematically analyzed using multi-omics data from The Cancer Genome Atlas (TCGA). The objectives of the present study were to identify KIF-based molecular subtypes, characterize their multi-omics phenotypes, tumor microenvironment (TME) features and therapeutic responses, and validate these findings in independent cohorts. The present study provides a novel perspective on CC heterogeneity and offers potential biomarkers to guide personalized therapeutic strategies.

Materials and methods

Patient cohorts and multi-omics data acquisition

The study population consisted of three primary cohorts. TCGA Cohort: Clinical information, follow-up data, RNA-seq, DNA copy number, and somatic mutation data for patients with CC were obtained from the Genomic Data Commons http://portal.gdc.cancer.gov/). RNA-seq data were normalized using the Fragments Per Kilobase of Exon Model (FPKM) method. FPKM values were used exclusively for survival analyses and non-negative matrix factorization (NMF) subtyping. Prior to subtyping, FPKM values were transformed using log2(x+1) to approximate a normal distribution. Subsequently, expression values of each KIF gene were centered and scaled using the preProcess function from the caret package with the option method=c(‘center’, ‘scale’). South Korea and Mexico Cohorts: Gene expression and clinical follow-up data were retrieved from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE44001 (10) and GSE52903 (11), respectively. These cohorts included information such as age, Federation of Gynecology and Obstetrics (FIGO) stage, and survival outcomes, disease-free survival (DFS) for the South Korea cohort and overall survival (OS) for the Mexico cohort. Gene expression profiles were normalized according to the original methods reported for each dataset.

Immunotherapy datasets comprised three cohorts. The GIDE cohort included 41 patients with melanoma treated with anti-PD-1 monotherapy and 32 patients receiving combined ipilimumab and anti-PD-1 therapy, with RNA sequencing performed on tissue samples before or during early treatment (12). The Van cohort consisted of 39 patients with melanoma treated with ipilimumab, with RNA sequencing conducted prior to therapy (13). The IMvigor cohort included 298 patients with metastatic urothelial carcinoma treated with intravenous atezolizumab, with RNA sequencing of tumor samples performed before treatment as part of the phase II IMvigor210 clinical trial (14). Due to the limited availability of publicly accessible immunotherapy-treated CC cohorts with transcriptomic data, these well-characterized pan-cancer immunotherapy cohorts were used as surrogates to evaluate the predictive performance of the KIF subtype-based immunotherapy response models.

Non-negative matrix factorization

NMF decomposes a large non-negative matrix into two smaller non-negative matrices (15). In the present study, NMF was employed to reduce the dimensionality of high-dimensional KIF gene expression data. Prior to NMF, the input matrix was preprocessed as follows: (i) genes expressed in fewer than 50% of samples were excluded; (ii) genes not significantly associated with OS in patients with CC (P>0.05) were removed; and (iii) log2(expression+1) transformation was applied, followed by centering and scaling using the preProcess function from the caret package. This standardization reduces technical variability and improves comparability across samples. The selected KIF genes were then subjected to NMF using the nmf package in R, employing the Brunet method (16) with 30 iterations. The rank range was set from 2 to 7, and the optimal rank was determined based on the Brunet criterion, which considers the maximum change in cophenetic correlation, consensus matrix clustering performance, and the sample sizes within each category. KIF genes were pre-selected based on their association with OS to focus clustering on clinically relevant genes, thereby enhancing the biological interpretability of the resulting subtypes. To avoid potential circularity, using survival information for both gene selection and outcome evaluation, the prognostic value of the identified subtypes was subsequently validated in two independent external cohorts (GSE44001 and GSE52903), which were not used during gene selection or model development. The consistent survival differences observed in these external datasets support the robustness and generalizability of the KIF-based subtyping approach.

Survival analysis, concordance index and decision curve analysis

Survival analyses were performed using the Kaplan-Meier (KM) method and the Cox proportional hazards regression model. Differences between KM curves were evaluated using the log-rank test. Variables with P<0.10 in univariate Cox regression analysis, including age, FIGO stage, body mass index (BMI), and smoking history, were included in the multivariable Cox proportional hazards model. Stepwise selection was applied to determine the final model. Optimal cut-off values were identified using maximally selected rank statistics from the R package survminer (https://CRAN.R-project.org/package=survminer). Hazard ratios (HRs) and 95% confidence intervals (95% CIs) were calculated. The clinical utility of KIF subtypes was assessed by comparing the concordance index (C-index) and decision curves, with 95% CIs for C-index values estimated via bootstrap resampling using the R package CsChange (https://CRAN.R-project.org/package=CsChange).

Neural network model based on RNA expression of KIFs (KRNNM)

The TCGA cohort was randomly split into a training set (70%, n=213) and a validation set (30%, n=91). Expression values were log-transformed and normalized prior to analysis. KRNNM was constructed using the resilient backpropagation (RPROP) algorithm. The input layer contained neurons corresponding to the number of KIFs, while the output layer contained neurons corresponding to the number of KIF subtypes. Predicted KIF subtypes were compared with those obtained from NMF. Model performance was evaluated using metrics including accuracy, precision, sensitivity, specificity, and F1 score, derived from the confusion matrix. The area under the receiver operating characteristic curve (AUROC) was calculated using the R package pROC (https://CRAN.R-project.org/package=pROC).

Gene enrichment analysis and proteomics data analysis

Differential expression analysis between KIF subtypes was performed using DESeq2 on count matrices. P-values were adjusted using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR). Genes with FDR <0.05 and |log2(fold change)| >1 were considered significantly differentially expressed. Functional enrichment analysis was conducted using the clusterProfiler package in R (https://www.bioconductor.org/packages/clusterProfiler), focusing on Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Pathways with FDR <0.05 were considered significantly enriched. Level 4 proteomics data for the TCGA CC cohort were obtained from The Cancer Proteome Atlas (TCPA, http://www.tcpaportal.org/) and analyzed using reverse phase protein array technology, covering 218 proteins.

TME components and immune characteristics analysis

The xCell algorithm (17) was used to estimate the cellular composition of CC samples. Cell types with scores of zero in more than 50% of samples or with a standard deviation <0.1 were excluded, resulting in 48 retained cell types. Immune and stromal scores were calculated using the ESTIMATE algorithm. Immune characteristics included: (i) Expression of 84 immune-related genes, categorized into 12 types; and (ii) five immune feature scores, wound healing, IFN-γ response, overall lymphocyte infiltration, TGF-β response and monocyte/macrophage regulation, derived from Thorsson et al (18).

Somatic mutation and DNA copy number variation analysis

Somatic mutation data for patients with CC were obtained from the TCGA database. The R package maftools (https://www.bioconductor.org/packages/maftools) was used to construct a binary mutation matrix. Only non-synonymous mutations were considered, while synonymous mutations were treated as wild type. Fisher's exact test was applied to identify differentially mutated genes among KIF subtypes. Tumor mutation burden (TMB) was calculated based on all non-synonymous mutations. Oncogenic pathways analyzed included TP53, Wnt, RTK-RAS, PI3K, NOTCH, Cell Cycle, Hippo, MYC, NRF2 and TGF-β. DNA copy number variation (CNV) data, excluding germline mutations, were also retrieved from the TCGA database. GISTIC 2.0 analysis was employed to identify significant amplification and deletion sites in autosomes, using hg38 as the reference genome and a q-value threshold of <0.05. Visualizations were generated using the maftools package.

Prediction of antitumor treatment efficacy

Chemotherapy drug sensitivity was evaluated using the oncoPredict algorithm (https://CRAN.R-project.org/package=oncoPredict). Drug response data from the Genomics of Drug Sensitivity in Cancer (GDSC, http://www.cancerrxgene.org/) database, together with transcriptomic profiles of CC cell lines, including HeLa, SISO, ME-180, MS751, SiHa, Ca-Ski, C-33-A, dot2-4510, CAL-39, SKG-IIIa, HT-3, SW756 and C-4-I, were used to train a regression model. This model predicted the sensitivity of different KIF subtypes in patients with CC to chemotherapy agents. The potential response of patients with CC to immunotherapy was inferred using six computational methods: MHC-I association immunoscore (MIAS) (19), Immuno-predictive score (IMPRES) (20), T cell-inflamed gene expression profile (GEP) score (21), Tumor immune dysfunction and exclusion (TIDE) score (22), Immunophenoscore (IPS) (23) and subclass mapping (submap) (24). Predictions from these six algorithms were integrated into a composite scoring system. For each algorithm, subtypes with higher immune scores (indicating improved predicted response) received 1 point, and subtypes with lower scores received 0 points. A total score >3 (at least four out of six algorithms predicting response) was used as an empirical threshold to define a high probability of benefit. This threshold is exploratory and has not been optimized on training data, and its predictive performance was evaluated in three independent immunotherapy cohorts. The prediction model was subsequently applied to immunotherapy datasets to assign KIF subtypes and compare objective response rates (ORR). Treatment efficacy was assessed using RECIST 1.1 criteria (http://www.recist.com/), with complete response or partial response defined as an objective response, and progressive disease or stable disease classified as non-objective response. It should be noted that chemotherapy sensitivity predictions were derived from cell line models (GDSC database) and may not fully recapitulate in vivo drug responses in patients with CC due to differences in TME, pharmacokinetics, and inter-patient heterogeneity. Therefore, these predictions are exploratory and require further validation in preclinical models and clinical cohorts.

Immunohistochemistry

To independently validate the protein expression of the identified marker KIFs (KIF1A and KIF4A) in an external cohort, commercially available tissue microarrays (TMAs) containing tumor tissues and paired adjacent non-tumor tissues from 80 Chinese patients with CC were obtained (cat. no. CXC1601; Shanghai Zhuohao Pharmaceutical Company). Written informed consent had been obtained from all patients, and ethical approval for the collection and use of these samples was granted by the company's institutional review board. Immunohistochemistry (IHC) staining for KIF1A and KIF4A was performed in our laboratory on these TMAs. After staining, slides were evaluated as follows: five fields of view were randomly selected per slide under a microscope (magnification, ×200). Staining intensity was scored as 0 (very weak), 1 (weak), 2 (moderate), or 3 (strong), and the proportion of protein-positive cells was scored as 1 (0–25%), 2 (25–50%), 3 (50–75%), or 4 (>75%). The staining score for each field was calculated by multiplying the intensity score by the proportion score, and the mean of these scores was used as the final staining score for the slide. Cases were classified into high- and low-expression groups based on the optimal cutoff value. The resulting IHC scoring data were then used to develop and validate the KKIHC classification system as an independent confirmation of the bioinformatics-derived KIF subtypes.

Statistical analysis

Subtype-specific differentially expressed genes were identified as those consistently upregulated or downregulated in one subtype compared with the other two, with no significant expression differences observed between the other two subtypes, as determined by the Wilcoxon test. Multinomial logistic regression analysis was subsequently performed on the selected KIFs, and genes with P<0.001 were designated as marker KIFs. Missing data were present in several clinical variables: Serum CA125 levels (12% missing), serum SCC levels (15% missing), BMI (8% missing), and tumor grade (5% missing). Variables with less than 20% missingness were imputed using the Monte Carlo method with multiple imputation by chained equations (5 imputations, 10 iterations). The imputation model included all variables used in subsequent analyses to preserve multivariate relationships. Sensitivity analyses comparing results with and without imputation yielded similar conclusions, indicating that the imputation did not introduce substantial bias. Group differences in categorical variables were assessed using the Chi-square test or Fisher's exact test. For continuous variables, two-group comparisons were conducted using the Wilcoxon test, and comparisons among three or more groups were performed using the Kruskal-Wallis test. Correlation analyses were conducted using Spearman's rank correlation. Unless otherwise specified, all statistical tests were two-sided, and a P<0.05 was considered to indicate a statistically significant difference. All statistical analyses were performed using R software (version 4.2.2; R Foundation for Statistical Computing).

Results

Expression patterns and prognostic values of KIFs

The overall study design is illustrated in Fig. 1A. Expression levels varied substantially across KIFs, with certain members (for example, KIF1A) exhibiting pronounced inter-patient heterogeneity (Fig. 1B). The expression of most KIFs was positively correlated, although a few (for example, KIFC3 and KIF9) demonstrated negative correlations with other KIFs (Fig. 1C). A total of 23 KIFs were significantly associated with OS in patients with CC (Fig. 1D).

Figure 1. Expression patterns of KIFs and survival analysis. (A) Research workflow. (B) Expression levels of KIFs in patients with CC from the TCGA cohort. KIFs are arranged from left to right in ascending order of expression and color-coded according to their subfamilies. (C) Correlation analysis and log-rank test of KIF expression in the TCGA CC cohort. Correlations with FDR <0.0001 are displayed. Red lines represent positive correlations, indigo lines represent negative correlations, and line thickness indicates the strength of the correlation. Genes are color-coded according to their KIF subfamilies. Black outlines denote log-rank test P<0.05. (D) Forest plot of univariate Cox regression for KIF expression. KIFs are arranged from top to bottom in ascending order of hazard ratio values, with dot size representing expression levels. KIF, kinesin family genes; CC, cervical cancer; TCGA, The Cancer Genome Atlas; FDR, false discovery rate; CI, confidence interval.

Clinical characteristics and prognostic differences among KIF subtypes

A total of 304 patients with CC from the TCGA cohort were analyzed, with additional validation cohorts from South Korea (n=300) and Mexico (n=55) (Tables SI and SII). NMF clustering of the 23 prognostically significant KIFs identified three subtypes: C1 (n=130), C2 (n=140) and C3 (n=34) (Fig. 2A and B). Based on univariate analysis, variables with P<0.10, including age, FIGO stage, BMI and smoking history, were entered into the multivariable Cox proportional hazards model (Table SIII).

Figure 2. Identification of KIF subtypes and construction of the KRNNM model. (A) Cluster number analysis from rank=2 to rank=7; rank=2 exhibited superior clustering effectiveness. (B) Cophenetic coefficient analysis showing the largest change between rank=3 and rank=4, indicating rank=3 as optimal according to Brunet's criterion. (C) Kaplan-Meier curves and univariate Cox regression analysis of KIF subtypes in The Cancer Genome Atlas cohort. (D) Construction and validation of the KRNNM neural network model. The workflow includes: Standardization of training set data, conversion into an input matrix (n samples × 23 KIF gene expressions), neural network training to generate the prediction model, prediction of KIF subtypes for both training and test sets, and generation of a confusion matrix. (E) Kaplan-Meier curves for overall survival of KIF subtypes in the Mexico cohort. (F) Kaplan-Meier curves for disease-free survival of KIF subtypes in the South Korea cohort. KIF, kinesin family genes; KRNNM, neural network model based on RNA expression of KIFs.

C2 patients exhibited worse OS compared with C1 patients (log-rank test P=0.001), with multivariate analysis showing an increased risk of poor OS (HR=1.887, 95% CI: 1.052–3.384, P=0.033, Table I). C2 patients also had poorer disease-specific survival relative to C1 (log-rank test P=0.016; HR=1.903, 95% CI: 1.018–3.556, P=0.044, Table I).

Table I. Multivariate Cox model analysis of the prognostic value of KIF subtypes.

Table I.

Multivariate Cox model analysis of the prognostic value of KIF subtypes.

Model	HR	95LL	95UL	P-value
Multivariable analysis model 1a (overall survival)
C1	Ref.
C2	1.887	1.052	3.384	0.033
C3	1.757	0.764	4.041	0.185
Multivariable analysis model 2a (disease-specific survival)
C1	Ref.
C2	1.903	1.018	3.556	0.044
C3	2.078	0.875	4.934	0.097
Multivariable analysis model 3a (disease-free survival)
C1	Ref.
C2	1.318	0.523	3.318	0.558
C3	3.859	1.390	10.717	0.010
Multivariable analysis model 4a (progression-free survival)
C1	Ref.
C2	1.514	0.842	2.723	0.166
C3	2.892	1.447	5.777	0.003

a Additional factors adjusted for in multivariate model 1 include age, FIGO stage, BMI and smoking history; in multivariate model 2, the adjustments include age and FIGO stage; in multivariate model 3, the adjustment includes race; and in multivariate model 4, the adjustments include FIGO stage, BMI, and race.

Compared with C1, C3 patients had shorter DFS (log-rank P=0.008) and progression-free survival (PFS, log-rank P<0.001) (Fig. 2C). Multivariate analysis indicated that the C3 subtype was associated with shorter DFS (HR=3.859, 95% CI: 1.390–10.717, P=0.010) and shorter PFS (HR=2.892, 95% CI: 1.447–5.777, P=0.003, Table I).

Significant differences in BMI (P=0.025) and tumor grade (P=0.033) were observed among the three subtypes, with C3 patients exhibiting higher tumor grade, a feature generally associated with more aggressive tumor behavior (Table SIV).

Development and initial validation of KRNNM

A neural network prediction model (KRNNM) was constructed based on the expression of the 23 KIFs (Fig. 2D). The model achieved near-perfect performance on the training set (AUROC = 0.972, accuracy = 97.7%) and high accuracy on the validation set (AUROC = 0.930, accuracy = 89.0%) (Table SV). This model was applied to two independent CC cohorts. In the Mexico cohort, C2 patients exhibited significantly worse OS than C1 patients (HR=2.47, 95% CI: 1.09–6.87, P=0.044, Fig. 2E). In the South Korea cohort, C3 patients had significantly worse DFS than C1 patients (HR=2.91, 95% CI: 1.28–6.60, P=0.011, Fig. 2F).

The prognostic value of KIF subtypes was further assessed using C-index and decision curve analysis (DCA). In the TCGA cohort, the C-index for OS improved from 0.616 (95% CI: 0.536–0.696) with FIGO staging alone to 0.700 (95% CI: 0.631–0.769) when combined with KIF subtypes (P=0.001). In the Mexico cohort, the C-index increased from 0.518 (95% CI: 0.401–0.633) to 0.642 (95% CI: 0.539–0.747) (P=0.040). DCA indicated that incorporation of KIF subtypes increased the net clinical benefit for patients (Fig. S1A-F).

Genomic alterations in the C3 subtype are enriched in Wnt and TP53 pathways

The most frequently mutated genes across all KIF subtypes were TTN and PIK3CA (Fig. S2). TP53 and FSIP2 exhibited higher mutation rates in C3 patients, while PIK3R4 mutations were almost exclusively observed in C3 (P<0.05, Fig. 3A). Significantly higher mutation frequencies were observed in the Wnt and TP53 signaling pathways in C3 patients compared with other subtypes (P<0.05, Fig. 3B). No significant differences in TMB were detected among the three subtypes (Fig. 3C).

Figure 3. Somatic mutation profiles and DNA copy number variation of patients with cervical cancer across different KIF subtypes. (A) Genes with differential mutation frequencies among C1, C2 and C3 groups, showing only genes with an overall mutation rate >10%. (B) Mutation status of common oncogenic pathway gene sets across C1, C2, and C3 groups. (C) Tumor mutation burden across C1, C2, and C3 groups (ns, not significant). (D) From top to bottom: C1, C2 and C3 subtypes. The left panel displays the distribution of GISTIC scores across chromosomal loci, with red indicating significantly amplified regions and blue indicating significantly deleted regions. The right panel presents a bubble chart depicting sample alteration frequency and the number of genes involved in the ten chromosomal loci with the lowest false discovery rate values. ns, not significant.

Analysis of CNV revealed that C1 and C2 patients exhibited relatively high-level CNVs, whereas C3 patients displayed fewer significantly amplified or deleted loci (Fig. 3D). FNDC3B was the only gene commonly amplified across all three subtypes (Fig. S3A and B). Enrichment analysis of subtype-specific CNV genes indicated that C1-specific amplified genes were enriched in B-cell differentiation and TNF receptor binding, whereas C2-specific amplified genes were enriched in T-cell proliferation pathways (Fig. S3C).

Transcriptomics and proteomics analyses

Transcriptome analysis showed that YAP1 was upregulated in C2, while KIF1A was upregulated and IL1A downregulated in C3 (Fig. 4A). Functional enrichment analysis revealed that immune-related pathways, such as leukocyte chemotaxis, were upregulated in C2 but downregulated in C3. Hippo-YAP and VEGFR signaling pathways were enriched in C2, whereas Wnt signaling was enriched in C3 (Fig. 4B). Proteomics analysis revealed consistent alterations: in C2 patients, YAP and YAP_pS127, as well as VEGFR2, were upregulated, whereas LCK protein was downregulated in C3 patients (Fig. 4C).

Figure 4. Transcriptome and proteomics analyses. (A) Volcano plots illustrating differential gene expression across the transcriptome for patients classified as C1, C2 and C3 subtypes. Significantly upregulated genes are shown in red, significantly downregulated genes in blue, and non-significant genes in gray. The x-axis represents log2 FC in gene expression between groups, and the y-axis depicts -log10 FDR. Genes with the highest log2FC and FDR values are labeled to highlight prominent differential expression. (B) Enrichment analysis of differentially expressed genes across the transcriptome for C1, C2 and C3 patients. Pathways and biological processes are grouped into categories including oncogenic pathways, growth factor activity, cellular response to stress or stimuli, ECM interaction, immune activity, and inorganic ion transport. The heatmap illustrates relative enrichment within each subtype, with red indicating upregulation and blue indicating downregulation. (C) Proteomics results for patients with cervical cancer. The heatmap displays proteins specifically upregulated or downregulated in each subtype. Clinical variables, including grade, smoking history, FIGO stage, race, body mass index and age, are annotated at the top. The upper panel shows upregulated proteins, and the lower panel shows downregulated proteins. FC, fold-change; FDR, false discovery rate; ECM, extracellular matrix.

TME and immune characteristics

Within the TME, C1 patients exhibited higher infiltration of lymphoid cells, including CD8+ T cells, B cells and Th1 cells, whereas C2 patients displayed increased mesenchymal and endothelial cell populations. C3 patients demonstrated reduced lymphoid and myeloid cell infiltration (for example, CD4+ Tem cells, dendritic cells) and increased fibroblast content (Fig. S4).

C1 patients had lower stromal scores, whereas C3 exhibited lower immune scores; overall microenvironment scores were highest in C2 (Fig. 5A). HLA molecule expression was reduced in C3 compared with C1 patients (Fig. 5B). Cytotoxic genes, including PRF1, GZMA and GNLY, showed similar expression levels in C1 and C2, whereas inflammation-related genes, such as IL1A and TGFB1, were highly expressed in C2. TCF1 (HNF1A) expression was lower in C2 patients (Fig. 5B), and TGF-β response scores were elevated in both C2 and C3 patients (Fig. 5C).

Figure 5. TME and immune characteristics across different KIF subtypes. (A) Comparison of TME scores among KIF subtypes, including Stromal Score, Immune Score, and ESTIMATE Score. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. ns, not significant. (B) Differential expression of 49 immune-related genes among KIF subtypes. The heatmap highlights genes with significant expression differences (P<0.05) between subtypes. (C) Heatmap illustrating the relationship between KIF subtypes and immune feature scores. TME, tumor microenvironment; KIF, kinesin family genes; IFN, interferon.

Predicted therapeutic responsiveness difference among KIF subtypes

Predicted sensitivity to nine chemotherapy agents was assessed across KIF subtypes. C2 patients exhibited significantly lower predicted IC50 values for bleomycin, cetuximab and mitomycin C, whereas they were less sensitive to cisplatin and gemcitabine. C1 patients were less sensitive to 5-fluorouracil and docetaxel but more sensitive to cisplatin and gemcitabine, while C3 patients were less sensitive to cetuximab and mitomycin C but more sensitive to gemcitabine (Fig. S5).

Predicted immunotherapy responses differed among subtypes. C1 and C2 patients had significantly higher predicted MIAS and GEP scores than C3, and C1 patients exhibited lower predicted TIDE scores (Fig. 6A). For IMPRES and IPS scores, C1 and C2 included significantly more high-scoring patients than C3 (P<0.05, Fig. 6B). Submap analysis predicted C1 patients as likely responders to immunotherapy (Fig. 6C). Overall, C1 patients displayed the highest predicted probability of benefiting from immunotherapy, followed by C2, with C3 patients showing the lowest predicted response (Fig. 6D).

Figure 6. Immunotherapy sensitivity scores and predicted objective response rates. (A) MIAS, GEP and TIDE scores for different KIF subtypes, with statistical significance annotated. **P<0.01. ns, not significant. (B) Distribution of high-scoring patients according to IMPRES and IPS assessments. (C) Assessment of potential immunotherapy responsiveness using the Submap algorithm. (D) Summary of overall immune benefit scores for each subtype across all six prediction methods. (E-G) Comparison of objective response rates (CR/PR) and overall survival probabilities among KIF subtypes in three independent immunotherapy cohorts (IMvigor, GIDE and Van). Left panels display the proportion of patients achieving CR/PR vs. PD/SD, and right panels present Kaplan-Meier survival curves with pairwise comparisons of overall survival among subtypes. KIF, kinesin family genes; CR, complete response; PR, partial response; PD, progressive disease; SD, stable disease; CI, confidence interval.

Using the KRNNM model, KIF subtypes were predicted in three independent immunotherapy cohorts (IMvigor, GIDE, Van), and ORR were calculated as exploratory endpoints. In the IMvigor cohort (metastatic urothelial carcinoma), ORR was 33.9% for C1, 17.5% for C2, and 12.2% for C3, with significant differences among subtypes (P<0.05). C2 patients had worse OS than C1 (HR=1.483, 95% CI: 1.074–2.048, Fig. 6E). In the GIDE cohort (melanoma), ORR was 66.7% for C1, 55.9% for C2, and 38.9% for C3 (Fig. 6F). In the Van cohort (melanoma), ORR was 18.8% for C1, 23.5% for C2, and 0% for C3; C2 and C3 patients exhibited worse OS than C1 (HR=2.534 and 3.903, respectively, Fig. 6G). These results indicate that C1 patients had the highest predicted probability of immunotherapy response, followed by C2, with C3 patients predicted to have the lowest probability.

Association of KIF subtypes with clinical characteristics and prognosis in Chinese patients with CC

KIF4A was identified as a marker of C2 (OR for C1 and C3 vs. C2: 3.274 and 2.923, respectively, Table II), and KIF1A was a marker of C3 (OR for C1 and C2 vs. C3: 0.009 and 0.003, respectively, Table III). KIF1A exhibited specific high expression in C3, whereas KIF4A expression was significantly lower in C2 (Fig. 7A).

Figure 7. KIF subtypes in Chinese patients with CC. (A) Expression differences of KIF1A and KIF4A among different subtypes. (B) IHC scoring workflow for KIF1A and KIF4A on tissue microarrays. Samples were first classified by KIF1A expression, with high KIF1A indicating the C3 subtype. For samples with low KIF1A, classification was further refined by KIF4A expression: low KIF4A indicated C2, and high KIF4A indicated C1. Optimal cutoff values were applied to distinguish high vs. low expression. In the OS dataset, cutoff values were 10 for KIF1A and 9.6 for KIF4A. In the PFS dataset, cutoff values were 10 for KIF1A and 8.4 for KIF4A. (C) Kaplan-Meier curves for OS of KIF subtypes in the Chinese CC cohort. (D) Kaplan-Meier curves for PFS of KIF subtypes in the Chinese CC cohort. (E and F) Comparison of C-index between basic Cox models and Cox models incorporating KIF subtypes. The C-index with 95% CI is displayed for each model across various clinical variables. Statistical significance between models is indicated by asterisks. *P<0.05, **P<0.01 and ***P<0.001. Blue represents the base model, and yellow represents the subtype-extended model. Error bars indicate variability of the C-index, emphasizing comparative predictive accuracy. KIF, kinesin family genes; CC, cervical cancer; OS, overall survival; PFS, progression-free survival; IHC, immunohistochemical; ns, not significant.

Table II. Screening of KIF characteristic genes in C2 patients.

Table II.

Screening of KIF characteristic genes in C2 patients.

	C1			C3
Variables	β	Sth. error	OR	β	Sth. error	OR
KIF4A	5.317	0.223	3.274 (2.114–5.069)a	4.186	0.256	2.923 (1.769–4.831)a
KIF9	2.498	0.032	1.078 (1.011–1.148)b	0.043	0.172	1.113 (1.023–1.211)b
KIF22	3.923	0.016	1.064 (1.031–1.097)a	1.575	0.022	1.035 (0.992–1.079)

a P<0.001 and

b P<0.05.

Table III. Screening of KIF characteristic genes in C3 patients.

Table III.

Screening of KIF characteristic genes in C3 patients.

	C1			C3
Variables	β	Sth. error	OR	β	Sth. error	OR
KIF1A	−4.007	1.180	0.009 (0.001–0.089)a	−4.347	1.337	0.003 (0.000- 0.041)a

a P<0.001.

IHC staining of KIF1A and KIF4A was performed on 80 tissue samples, and the KKIHC classification system was developed: patients were first stratified by KIF1A expression (high → C3); among those with low KIF1A, high KIF4A indicated C1, and low KIF4A indicated C2 (Fig. 7B). Among the three subtypes, only serum SCC levels differed significantly (P=0.030, Table SVI).

Survival analysis demonstrated that C2 and C3 patients had poorer OS compared with C1 (log-rank P<0.05). Univariate Cox regression yielded HRs of 3.090 for C2 vs. C1 and 5.269 for C3 vs. C1 (Fig. 7C). Multivariate analysis confirmed that C2 and C3 subtypes were independently associated with poorer OS (HR=2.136 and 5.284, respectively, Table IV). Similarly, C2 and C3 patients had shorter PFS than C1 (log-rank P<0.05), with univariate HRs of 4.058 and 4.695, and multivariate HRs of 4.146 and 6.125, respectively (Fig. 7D, Table IV).

Table IV. Multivariate Cox model analysis of the prognostic value of KIF subtypes.

Table IV.

Multivariate Cox model analysis of the prognostic value of KIF subtypes.

Model	HR	95L	95H	P-value
Multivariable analysis model 1a (overall survival)
C1	Ref.
C2	3.371	1.043	10.895	0.186
C3	6.669	2.331	19.076	0.001
Multivariable analysis model 2a (progression-free survival)
C1	Ref.
C2	3.612	1.306	9.992	0.008
C3	3.764	1.606	8.822	0.000

a Additional factors adjusted for in multivariate model 1 and 2 include age, TNM, tumor size, CA125, SCC and HPV.

Integration of KIF subtypes into Cox models substantially improved the C-index for both OS and PFS. For example, the C-index for TNM stage increased from 0.758 to 0.835 for OS (P<0.001) and from 0.752 to 0.826 for PFS (P=0.005). Similar improvements were observed for models including age, tumor size, CA125, SCC and HPV status (Fig. 7E and F).

Discussion

In the present study, three distinct molecular subtypes of CC (C1, C2 and C3) were identified based on the expression profiles of 23 kinesin family genes. These subtypes exhibited marked differences in prognosis: C1 patients demonstrated the most favorable survival outcomes, C2 patients showed intermediate prognosis characterized by active but partially therapy-resistant immune features, and C3 patients had the poorest prognosis, with Wnt pathway activation and an immune-silent phenotype. Multi-omics analyses further revealed subtype-specific molecular alterations, including enrichment of Hippo-YAP and VEGFR signaling in C2 and frequent Wnt and TP53 mutations in C3, as well as predicted differential responses to chemotherapy and immunotherapy. These findings establish a framework for patient stratification and suggest potential therapeutic targets, thereby supporting the development of personalized treatment strategies for CC.

The biological and immunological distinctions among KIF subtypes help explain variations in prognosis and predicted immunotherapy efficacy, while also suggesting potential intervention points. C1 was characterized by downregulated expression of several KIFs, including KIF1A, and relatively low activity of oncogenic pathways. In C1 patients, genes located at significantly amplified CNV loci were enriched in pathways associated with negative regulation of immune cell apoptosis and differentiation of immature B cells. Functional enrichment analysis further indicated that most oncogenic pathways, including Wnt, VEGFR and PI3K-Akt, were comparatively downregulated, consistent with lower proliferative activity. Additionally, the TME in C1 patients was enriched with infiltrating immune cells, including CD8+ T cells and helper T cells, which may enhance tumor recognition and elimination, thereby strengthening antitumor immune responses (25,26). Most immune-related genes were significantly upregulated, suggesting active immune function. The combination of high lymphocyte infiltration and absence of prominent immunosuppressive signals indicates that C1 patients possess a relatively immune-active TME (27,28). These features likely contribute to the superior clinical outcomes and predicted immunotherapy efficacy observed in C1 patients (29,30). However, direct experimental evidence in CC cohorts remains limited, and the precise mechanisms underlying this favorable prognosis warrant further investigation.

Compared with C1, C2 exhibited distinct KIF expression patterns, including downregulation of KIF4A and upregulation of specific mitotic KIFs, such as KIFC2, which may contribute to aberrant activation of downstream signaling pathways. The Hippo pathway, for example, was highly enriched in C2 patients, accompanied by increased YAP activity and elevated YAP_pS127 protein levels. Excessive YAP activation is known to promote tumor cell proliferation in other cancers (31–33), suggesting a potential oncogenic role for YAP in C2 progression, although functional validation in CC models is required. Additionally, VEGFR signaling was enriched in C2, indicating potential suitability for anti-angiogenic therapy (34,35). Upregulation of VEGFR2 in C2 patients suggests activation of VEGF signaling, which may contribute to the aggressive phenotype observed in this subtype. Despite both C1 and C2 patients being predicted as high-probability responders to immunotherapy, C2 patients may derive comparatively less benefit. Mild increases in PD-L1 expression and immune cell infiltration were observed in C2, which could partially account for their response to immunotherapy (36–39); however, direct evidence in CC is lacking. The proportion of stromal cells, including fibroblasts, was higher in C2 patients. Cancer-associated fibroblasts (CAFs), potentially induced by pro-tumor signals such as IL-1, may contribute to immunotherapy resistance in C2, as suggested by previous studies (40–43). Immune-active pathway enrichment and elevated TGF-β response were also prominent in C2 patients. High levels of TGF-β, possibly secreted by CAFs, could suppress effector T-cell function while promoting regulatory T-cell activity, thereby attenuating antitumor immunity (44–48). Furthermore, reduced HLA gene expression in both C2 and C3 patients may further compromise immunotherapy efficacy (49–52). These complex molecular and microenvironmental alterations help explain why, despite exhibiting high immune activity, C2 patients display partial resistance to immunotherapy and have poorer clinical outcomes compared with C1 patients. Lower TCF1 expression in C2 suggests a shift of T cells toward an exhausted phenotype, which could account for diminished immunotherapy responsiveness despite the presence of cytotoxic effectors. Moreover, the combination of elevated TGF-β signaling and increased mesenchymal cell content is associated with transcriptomic enrichment of TGF-β pathways, potentially contributing to the observed partial immunotherapy resistance in C2.

Patients with the C3 subtype exhibited the poorest clinical prognosis, characterized by pronounced immunosuppression and highly aggressive tumor behavior. This subtype was defined by marked upregulation of KIF1A, which may act as a driver of downstream effects, including excessive activation of the Wnt signaling pathway. In C3, Wnt signaling was hyperactivated, potentially reflecting the higher mutation rate observed in this subtype, although this hypothesis requires experimental validation. Aberrant Wnt signaling is known to promote tumor growth, metastasis and maintenance of stemness in multiple cancers (53–56) and may similarly contribute to the aggressive phenotype of C3. Furthermore, Wnt activation can impair immune cell function and facilitate immune evasion (57,58). In the C3 subtype, protein levels of KU80, FASN and IGFBP2 were significantly upregulated. KU80 regulates FOXF2-mediated Wnt/β-catenin signaling and promotes colorectal cancer progression (59), FASN facilitates epithelial-mesenchymal transition (EMT) through the PRRX1/Wnt/β-catenin axis in salivary adenoid cystic carcinoma (60), and IGFBP2 activates Wnt signaling to promote EMT in solid tumors (61). Consistent with these molecular alterations, C3 patients exhibited lower predicted immunotherapy benefit, lower microenvironment scores, and reduced enrichment of immune-related pathways, reflecting an immune-silent phenotype. Downregulation of LCK protein in C3 further supports the lack of immune cell infiltration, while decreased HLA expression suggests impaired antigen presentation, potentially facilitating immune evasion. Additionally, C3 patients showed increased fibroblast infiltration, which can differentiate into CAFs under stimulation from TGF-β, IL-1, and Wnt signaling (40). Collectively, the extremely poor prognosis of C3 patients is likely associated with the combined effects of Wnt pathway activation, reduced lymphoid and myeloid cell populations, increased fibroblast content, and elevated TGF-β response. The immune-cold phenotype of C3, characterized by low lymphoid/myeloid infiltration and diminished HLA expression, is consistent with LCK downregulation and the enrichment of Wnt pathway mutations observed in genomic and proteomic analyses. Elevated TGF-β response scores in both C2 and C3 patients further indicate an immunosuppressive microenvironment, concordant with increased fibroblast infiltration in these subtypes.

The strength of the present study lies in the integration of multi-omics data, which allowed the linkage of genomic alterations (for example, Wnt and TP53 mutations in C3) to transcriptomic and proteomic changes (for example, KIF1A upregulation, LCK downregulation) and features of the tumor immune microenvironment. Importantly, KIF family genes themselves emerged as central hubs that may orchestrate these multi-layer alterations, rather than merely serving as downstream markers of tumor heterogeneity. While the current integration is primarily descriptive, the concordance across data layers supports the biological validity of the KIF subtypes and generates testable hypotheses for future mechanistic studies.

These findings indicate that KIF expression heterogeneity contributes directly to the distinct molecular and clinical characteristics of each subtype, supporting the use of KIF genes as both subtyping biomarkers and potential therapeutic targets. The present study provides a foundation for hypothesis-generating strategies in personalized treatment, although prospective validation is required before clinical application. If validated in future studies, KIF subtypes could inform clinical decision-making. For instance, patients with the C1 subtype, who demonstrate favorable prognoses and predicted benefit from immunotherapy, might be considered for immune-based therapies in prospective trials. Patients with the C2 subtype may require additional interventions targeting Hippo-YAP and VEGFR signaling, although this remains speculative and warrants preclinical and functional validation. Furthermore, the lower predicted IC50 values for certain chemotherapeutic agents in C2 patients suggest potential sensitivity to alternative regimens, such as bleomycin or cetuximab. C3 patients, characterized by aggressive Wnt-driven tumors and an immune-silent phenotype, may benefit from combination strategies aimed at inhibiting Wnt signaling and restoring immune function; however, these hypotheses await empirical testing. No direct clinical recommendations can be drawn from the present retrospective analysis. The use of KIF1A and KIF4A as biomarkers for prognosis and predicted therapeutic response is particularly promising, pending validation in CC-specific immunotherapy cohorts. High KIF1A expression, associated with the C3 subtype, and low KIF4A expression, characteristic of the C2 subtype, provide actionable insights for patient stratification, enabling the optimization of treatment plans according to individual molecular profiles. This approach not only enhances prognostic accuracy but also facilitates the development of targeted therapies, ultimately promoting more personalized and effective clinical management of CC. Two KIF-based subtyping strategies were developed in the present study: KRNNM and KKIHC. KRNNM, based on RNA-seq data, may be applied in research contexts, whereas KKIHC classifies patients using IHC-determined KIF1A and KIF4A expression. Both approaches remain exploratory and require prospective validation in independent clinical cohorts before clinical implementation. Their potential to improve prognostic precision or guide therapeutic decisions remains to be established.

Several limitations of the present study should be acknowledged. First, the KRNNM neural network model showed near-perfect performance on the training set, which may indicate a risk of overfitting given the limited sample size (n=304). Although the model architecture was simplified, regularization was introduced, and consistent prognostic separation was observed in two independent external cohorts, the model's performance metrics should be interpreted with caution, and prospective validation is needed. Second, the 23 KIF genes were pre-selected based on their association with OS before performing NMF clustering, which could introduce a risk of circular analysis. To mitigate this concern, the prognostic value of the identified subtypes in two external cohorts (GSE44001 and GSE52903) that were not involved in any step of gene selection or model development. The consistent survival differences observed in these independent datasets support the robustness and generalizability of the KIF-based subtyping approach. Third, due to the scarcity of publicly available CC cohorts treated with immune checkpoint inhibitors, melanoma and urothelial carcinoma cohorts were used as surrogates to infer immunotherapy response. These tumor types have distinct immune microenvironments and mutation burdens, which may limit the direct generalizability of the authors' predictions. Fourth, chemotherapy sensitivity predictions were based on cell line-derived models (GDSC database) and may not fully reflect in vivo drug responses due to differences in TME, pharmacokinetics, and inter-patient heterogeneity. These predictions should be interpreted as hypothesis-generating and require experimental validation. Fifth, the IHC validation cohort was relatively small (n=80), and the optimal cut-off values for KIF1A and KIF4A staining require prospective validation in larger independent cohorts. Sixth, all analyses are retrospective, and the proposed clinical utility of KIF subtyping remains hypothesis-generating. Future studies with larger, independent, prospective cohorts are urgently needed to validate the subtypes, the prediction models, and the therapeutic implications.

In conclusion, the present study highlights the critical role of three KIF-defined molecular subtypes in determining prognosis and treatment response among patients with CC. C1-C3 subtypes displayed distinct prognostic profiles: C2 patients exhibited a similarly immune-active state as C1 patients, whereas C3 patients demonstrated an immunosuppressive phenotype, suggesting potential resistance to immunotherapy, although this requires confirmation in CC-specific cohorts. C2 was characterized by activation of Hippo-YAP and VEGFR signaling, representing potential therapeutic targets. In C3, Wnt pathway-related genes, TP53 mutations, and KIF1A overexpression may serve as actionable targets. If validated in future prospective studies, integration of KIF subtypes into clinical practice could enhance personalized treatment strategies, offering new avenues to improve outcomes for patients with CC. However, the current findings are preliminary and require further validation. Notably, predicted immunotherapy responses derived from non-CC cohorts necessitate confirmation in dedicated CC immunotherapy trials before any clinical application can be considered.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by the National Natural Science Foundation of China (grant nos. 82403850 and 82203502), the Cross-Innovation Talent Project of Renmin Hospital of Wuhan University (grant no. JCRCGW-2022-002), the Wu Jieping Medical Foundation (grant no. 320.6750.2023–05-90) and the Beijing Bethune Charitable Foundation (grant nos. ZLZX027 and 2024-YJ-226-J-021).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

XZ summarized the findings and drafted the manuscript. XZ and JW analyzed and interpreted the data. JW and LL provided critical revisions and feedback to enhance the manuscript. JW, LL and XZ designed the study. HZha and JL analyzed and interpreted the data. HZho and YD collected the data. HC, BX and QS conceived and designed the study.. XZ and JW confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

TMA samples for the CC cohort were purchased from Guilin Fanpu Biotech, and the ethical approval for use of these samples was obtained through the Guilin Fanpu Biotech Ethics Committee [approval no. Fanpu-2018 (23)]. All research involving human data was conducted in accordance with the Declaration of Helsinki.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References