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 log₂^(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.

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) log₂^{(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 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).

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).

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 |log₂(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.

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 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.

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.

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.

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).

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).

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).

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).

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).

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).

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).

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).

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).

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).

Predicted sensitivity to nine chemotherapy agents was assessed across KIF subtypes. C2 patients exhibited significantly lower predicted IC₅₀ 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).

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.

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).

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).

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).