3.1 Identification of subtypes based on GSTTKs gene

The flow chart of this study was shown in Fig. 1. A total of 1310 GSTTKs were collected and enrolled in our study (Supplementary Table S3). The Univariate Cox regression analysis suggested 61 GSTTKs seemed to have substantial prognostic value (Supplementary Table S4). The tight association between 61 GSTTKs and immune cells was demonstrated as shown in Fig. 2A. Furthermore, GSTKKs score presented substantial positive correlations with immune score, stromal score, and estimation score, as well as a significant negative connection with tumor purity (Fig. 2B–E). This suggested that GSTTKs may play a crucial role in tumor microenvironment and tumor immunity.

Fig. 1

The flow chart of our research

Fig. 2

The relationship between prognostic GSTKKs and immune landscape. A The correlation between 61 prognostic GSTKKs and immune cells, Red indicated positive correlations, and blue indicated negative correlations. Asterisks denoted p-value (**p < 0.01, *p < 0.05). Blank cells represented no statistical significance of the correlation. (B–E). Correlation of GSTKKs score with immune score (B), stromal score (C), estimate score (D), and tumor purity (E) generated by estimate algorithm

Three molecular subtypes constructed based on the expression amounts of the 61 GSTKKs by the NMF algorithm were exhibited in Fig. 3A. The UMAP and PCA analysis also confirmed that the three types were robust (Fig. 3B, Supplementary Fig. S1A). Previous studies have reported that a larger value of the contour coefficient indicates a higher matching relationship between the target and the cluster in which he or she is located, illustrating when it is close to 0 that he should be on the boundary points of the two clusters. Thus, patients with a contour coefficient > 0 were included in the study to achieve more precise patient stratification [29]. By comparison, there was a strong correlation between GSTKKs-based subtypes and published classical molecular subtypes, especially C2 was similar to Immunoreactive (Fig. 3D).

Fig. 3

Identification of GSTKKs-based subtypes. A The consensus map after NMF clustering revealed three clusters with no overlap between clusters. B The UMAP algorithm displayed the two-dimensional principal component diagram of three subtypes, with each point representing a single sample. C The pie chart showed variations of clinical characteristics between the three subtypes by Fisher’s exact test. D The links between GSTKKs-based subtypes and immune subtypes (TCGA 2014). E Kaplan–Meier curve of OS according to GSTKKs-based subtypes in TCGA-OV cohort. F Kaplan–Meier curve of OS according to immune subtypes in TCGA-OV cohort

Consent with the previous result, the K-M analysis displayed significant prognoses differences among the three types of samples (Fig. 3E). The survival curve of classical molecular subtypes also displayed similar results (Fig. 3F). Among them, C1 got the worst prognosis, C2 possessed the most excellent prognosis, and C3 seemed to have an intermediate prognosis.

To further validate the accuracy and stability of NMF results, NTP approach was performed on three independent cohorts based on the gene signature (Supplementary Table S5), including GSE32062, GSE53963, and GSE140082 (Fig. 4A–C). Similarly, K-M analysis revealed that C1 had the worst prognosis and C2 had the best prognosis in each cohort. These results indicated that the subtypes based on GSTKKs were repeatable and practical.

Fig. 4

Validation of GSTKKs-based subtypes in three independent cohorts by the nearest template prediction (NTP) approach. A-C From left to right: verification of the three GSTKKs-based subtypes in GSE32062, GSE53963, and GSE140082 cohorts with NTP analysis, Kaplan–Meier curves of OS according to the GSTKKs-based subtypes in GSE32062, GSE53963, and GSE140082 datasets, and proportion of different stages in the subtypes

3.2 Analysis of clinical characteristics related to subtypes prognosis

The relevance between clinical features of the TCGA cohort and the GSTKKs-based subtypes was manifested in Fig. 3C. The independent prognostic elements of TCGA-OV and three GEO cohorts were displayed in Supplementary Fig. S1C-E, G. In the TCGA cohort, C1 subtype was an independent prognostic factor; In the GSE140082 cohort, age, stage, and C2 subtype were independent prognostic factors; There were none both in GSE53963 and GSE32062 cohorts. The proportion of each stage and grade in the three subtypes were also performed (Fig. 4, Supplementary Fig. S1F, H). Consistent with previous studies, advanced patients were the majority in each subtype. To sum up, patients with favorable prognoses were mainly distributed in the C2 subtype, while poor prognoses were mainly distributed in the C1 subtype.

3.3 Biological characteristics of GSTKKs-based subtypes

GSEA enrichment analysis was utilized to characterize metabolic pathways and specific biological processes among the three subtypes (Fig. 5A). Proliferation-associated pathways were enriched in C1, such as Hedgehog signaling pathway. Immune-related pathways were enriched in C2, such as autoimmune thyroid disease. Metabolism pathways were enriched in C3, such as Th1 and Th2 cell differentiation Thus, molecular features of C1 were defined as proliferation, C2 as immune, and C3 as immune and metabolism. Further gene set variation analysis (GSVA) also validated the above molecular features (Fig. 5B), C1 was significantly associated with proliferative activity, C2 was mainly related to immune pathways, and C3 was significantly enriched in metabolic processes.

Fig. 5

Biological characteristics of the three subtypes. A From top to bottom, GSEA enrichment analysis respectively revealed activated marker pathways of C1, C2, and C3. The FDR of the biological function was < 0.05. B Heatmap based on the score of each subtype in 50 Hallmark gene sets. The higher the score, the higher the pathway activity

3.4 Genomic alterations analysis of GSTKKs-based subtypes

As illustrated in Fig. 6A, the landscapes of gene mutations and copy number alterations in each subtype were explored. Top 20 frequently mutated genes were exhibited in Fig. 6B. Among these genes, TTN and FLG showed differences in the three subtypes. Their mutation frequency in C1 was highest, suggesting that TTN and FLG might play a critical role in tumor development. Notably, while some mutations were not statistically different, such as LRP2 and BRCA1, which separately had the highest mutation rate in C3 and C2, these genes still deserve attention.

Fig. 6

The mutational landscape of three subtypes. A Mutation landscape of top 20 frequently mutated genes (FMGs) in the three clusters. B The mutation frequency of top 20 FMGs among three subtypes. P values are shown as *P < 0.05. C The differences of fraction of genome altered (FGA), fraction of genome gained (FGG), and fraction of genome lost (FGL) in clinical characteristics. D The burden of copy number gain or loss in arm and focal levels. P values are shown as *P < 0.05; **P < 0.01; ***P < 0.001

The association between copy number variants and the clinical characteristics of patients was further identified (Fig. 6C). Fraction of genome altered (FGA) and fraction of genome lost (FGL) did not differ in clinical characteristics. Grade2 and Grade3 were associated with fraction of genome gained (FGG) only when graded diagnoses were made (P < 0.001). FGA was a measure of genomic instability and represented more specific mechanisms of chromosome alterations. After that, the genomic variants of the three isoforms were further analyzed (Fig. 6D). C1 had fewer deletions and amplifications at the focal level, whereas at the chromosome arm level, there was no significant difference among subtypes. Further evaluation of the mutational landscape found that there was no difference among the three subtypes of HRD, new antigen and TMB (Supplementary Fig. S2B).

Altogether, in the analysis of genomic variants, differences brought about by mutations and copy number alterations might have contributed to the diversified outcomes of the three heterogeneous molecular subtypes.

3.5 Immunocyte infiltration landscape of GSTKKs-based subtypes

Because the three subtypes were enriched on significantly different biologically relevant pathways, further exploration of the immune landscape of each isoform was of great value in determining the best treatment. As illustrated in Fig. 7A-D, the immune score, stromal score and estimated score were the highest in C2, while the tumor purity was the lowest. This might proclaim that C2 was more prone to immunotherapy. The immune infiltrate status of the three subtypes was assessed using eight algorithms (Supplementary Fig. S2A). Respectively, C1, C2, and C3 were seemed as the moderate immune infiltrating tumor, immune-hot tumor, and immune-cold tumor.

Fig. 7

Analysis of the proportion of immune cells in three subtypes. A–D The immune score (A), stromal score (B), estimate score (C), and tumor purity (D) in three subtypes. E Estimated proportion of immune cells among the three subtypes. P values are shown as *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, and ns was the abbreviation of no significance. F The heatmap illustrated the correlation between three subtypes and clinical characteristics, and the infiltration abundance of 28 immune cell subsets evaluated by ssGSEA algorithm

Assessing the proportion of immune cell components in each subtype using the CIBERSORT algorithm illustrated that all three subtypes were certified with a higher proportion of macrophages and T cells (Fig. 7E). Meanwhile, the proportion of immune cells was the highest in C2, such as M1 macrophages and activated CD8 memory T cells. This may be related to better immunotherapeutic efficacy. Figure 7F revealed the correlation between subtypes and immune cells, C2 displayed a significant positive correlation in most immune cells, suggesting that C2 might be associated with immune activation, harboring a more vital anti-tumor killing ability. C3 was deemed immune evasion because the infiltration was lowest in almost all immune cells and was significantly negatively correlated with most immune cells. C1 was in an intermediate state, which only significantly negatively correlated with activated B cells and natural killer T cells.

3.6 Immune checkpoints analysis of GSTKKs-based subtypes

In co-inhibited and co-stimulated ICPs, C2 was expressed at the highest levels in most cases (Supplementary Fig. S3A, B). Due to C2 had the highest expression level of HLA molecules, it was presumed that its patients had better antigen delivery ability (Supplementary Fig. S3C). In line with these results, shown in a heatmap of 27 ICP molecules expression (Supplementary Fig. S3D), C2 was expressed at the highest levels, followed by C1 and C3. Also, in both leukocyte fraction and TIS score, the C2 subtype was significantly different from the other two groups (P < 0.001), which means C2 is more responsive to immunotherapy (Supplementary Fig. S3E, F).

The use of immune checkpoint inhibitors (ICIs) has revolutionized the treatment paradigm for various cancers. Thus, enrichment scores for the seven anti-tumor immune cycle steps were calculated using the ssGSEA algorithm, with all steps showing higher scores (P < 0.001) in C2 but lower levels in C1 and C3 (Fig. 8A, B). The TME was most abundant in C2 according to the results of Fges enrichment, indicating that C2 would respond positively to immunotherapy. Notably, CAFs (cancer associated fibroblasts) was considerably enriched in C1 (Fig. 8C), which can secrete a range of growth factors, cytokines, extracellular matrix, etc. These factors are crucial for promoting tumorigenesis, proliferation, tumor angiogenesis, invasion, and metastasis. This might also imply a poorer prognosis in C1. Additionally, the results of Submap demonstrated that C2 subtype might achieve excellent clinical effectiveness from an immunotherapy response (Fig. 8D). Overall, the evidence presented above all indicated that C2 was an immunological subtype, and immunotherapy for C2 patients might be a more potent weapon.

Fig. 8

Characteristics of immune circulation GSTKKs-based subtypes and clinical treatment cohorts. A The butterfly diagram illustrated the distribution of metabolic pathway and cancer immune cycle among the three subtypes. B Enrichment scores for the seven anti-tumor immune cycle steps were calculated with ssGSEA algorithm. P values are shown as *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001. C The radar map displayed the proportion of the immune-related characteristics and immune molecules in the three subtypes. D Submap analysis of the three subtypes, including GSE100797 cohort with detailed CAR-T therapy information, GSE126044 cohort with detailed anti-PD1 therapy information and GSE115821 with detailed anti-PD1 and anti-CTLA4 therapy information. E Analysis of the response of three subtypes to clinical drugs in the treatment cohorts, including cisplatin, paclitaxel, and platinum

3.7 Pharmacotherapy prediction for GSTKKs-based subtypes

The outcomes of the clinical medication treatment cohort were displayed in Fig. 8E. Compared with C2, C1 was more sensitive to cisplatin, and C3 patients might benefit from platinum treatment. According to the CMap analysis results, 29 drugs harbored individualized therapeutic potential for subtypes (Supplementary Fig. S2D), and the molecular pathways and genes they target were exhibited (Supplementary Fig. S2E). In addition, based on the ridge regression model of pRRophetic package, candidate drugs for different subtypes were developed by calculating half-maximal inhibitory concentration (IC50) and quantifying drug sensitivity data (Supplementary Fig. S2C, Supplementary Fig. S4A-C). Altogether, all these candidates may bring better efficacy to specific OV patients.