We downloaded the processed mRNA expression data of ovarian serous cystadenocarcinoma (OV) cohort from the TCGA database (https://portal.gdc.cancer.gov/). A total of 377 specimens with complete survival information were collected. The mutation data was obtained by downloading the SNP data of OV cohort while the GISTIC2 copy number data of this cohort was acquired from GDAC FireBrowse (http://firebrowse.org/). Predicted neoantigens for the OV cohort were retrieved from a previous analysis of the TCGA dataset (24). The GSE63885 (n=75 samples) and GSE17260 (n=110 samples) datasets were downloaded from the GEO database for gene expression profile and survival analyses of the corresponding patients in external validation. The combat function in the R SVA library was used for batch effect correction of merged datasets.

The flowchart of the study design is as indicated in Additional Figure S1 . The previously published metabolism-related genes (25) were retrieved and used for subsequent NMF clustering. Then, the R library “survival” was used for performing Cox regression analysis to evaluate the correlation between all candidate genes and overall survival (OS) of OV patients. Subsequently, the NMF package was used to perform an unsupervised NMF clustering method. The same candidate genes were used to apply the NMF clustering method to the two GEO external validation datasets. The k value at which the correlation coefficient begins to decrease was selected as the optimal number of clusters. SubMap analysis (Gene Pattern) is a method for evaluating the similarity of classification between patient cohorts based on expression profiles; this method was used to determine whether the subcategories identified in the above datasets were relevant. Then, based on the t-SNE method, the mRNA expression data of the above-mentioned metabolic genes were used to verify the subtype assignment.

The MCP algorithm was used to evaluate eight immune (immune cell populations: T cells, CD8+T cells, natural killer cells, cytotoxic lymphocytes, B cell lines, monocyte cell lines, and marrow cells, dendritic cells and neutrophils) and two non-immune stromal cell populations (endothelial cells and fibroblasts). In addition, the ssGSEA analysis, which calculates an enrichment score indicating the degree to which genes in a specific gene set are coordinated in a single sample, was also used to estimate immune infiltration. The GSVA-R software package was used to estimate six other immune cell populations, including regulatory T cells (Treg), helper T cells 1 (Th1), helper T cells 2 (Th2), and helper T cells 17 (Th17), central memory T cells and effective memory T cells (Tem), helper T cell (TFH), gamma delta T cells (Tgd). In addition, the ESTIMATE algorithm was used to calculate the immune score and matrix score, which can reflect the richness of the genetic characteristics of the matrix and immune cells.

The LIMMA library in R was used to analyze differentially expressed genes (DEGs) among OV subtypes; the selection was based on the condition of corrected P<0.05 and log₂FC greater or equal to 2. Only genes whose expression were significantly different in all three possible comparisons were considered subtype-specific genes. The top 30 genes with the largest log₂FC value in each subcategory were further selected to build a prediction model, and a 90-gene classifier was generated. Then, NTP algorithm was used to predict the sample subtypes based on the expression of the 90 genes in validation GEO datasets, followed by comparison of the classification results based on the NMF algorithm.

By measuring the similarity of gene expression profiles between the subtypes obtained in the present study and previously published SubMap analysis (gene pattern) of gene profiles in patients with melanoma immunotherapy, we indirectly predicted the efficacy of immunotherapy on the predicted subtypes. In addition, based on the largest pharmacogenomics database (GDSC Cancer Drug Sensitivity Genomics Database, https://www.cancerrxgene.org/), we used the R library “pRRophetic” to predict the chemotherapy sensitivity of each tumor sample. The regression method was used to obtain the IC50 estimated value of each specific chemotherapy drug treatment, and the GDSC training set was used to perform 10 cross-validation tests to test the regression and prediction accuracy.

Tissue chips were purchased from Shanghai Outdo Biotech Co.,Ltd. The chip number is HOvaC070PT01 include 69 cases. Immunohistochemical assay used the CXCL11(Affinity, Cat#DF9917) and PXDN (Cloud-Clone Corp, Cat#PAM070Hu01).The dilution of CXCL11 antibody is 1:200 and that of PXDN antibody is 1:100. The total score is the product of “staining intensity score” and “staining positive rate score” which represents CXCL11 and PXDN expression. According to the total score > 5 or < 5, they were divided into high expression group or low expression group. Use of all human samples was approved by the committee for ethical review of research involving human subjects at Shantou University and Jinan University.

Use the tissue slice digital scanner or imaging system to collect the scanned documents or images on the immunohistochemical section, automatically read the tissue measurement area by using the Seville image analysis system, and analyze and calculate the number of weak, medium and strong positive cells in the measurement area (negative without coloring, 0 point; weak positive light yellow, 1 point; medium positive brownish yellow, 2 points; strong positive brownish brown, 3 points), the total number of cells and the positive cumulative optical density IOD value, Positive pixel area, tissue area mm². The following results were calculated to reflect the degree of positivity. The following indexes can be selected according to the section to evaluate the intensity of positive cells.

All statistical analyses were performed using the R language software (version 4.0). Survival analysis was performed by Kaplan-Meier method, and the log-rank test was used for comparison. Univariate Cox proportional hazards regression model was used to estimate the hazard ratio of univariate analysis. All statistical tests were two-sided, and p<0.05 was retained as statistical significance treshold.

A list of 2752 previously published metabolism-related genes (these genes encode all known human metabolic enzymes and transport proteins) were used as input in the NMF analysis. After merging of TCGA and 2752 metabolism-related genes and correction of the batch effects using the Combat PCA algorithm ( Figure 1A ), Cox regression analysis was applied and a total of 177 prognostic-related candidate genes were identified ( Additional File S1 ). Next, the NMF consensus clustering method was used for OV classification based on the expression profiles of the above 177 candidate genes. The clustering was performed on the TCGA dataset, and after a comprehensive consideration, k=3 was selected as the optimal number of clusters based on the cophenitic coefficient ( Figure 1B ). For k=3, we performed dimensionality reduction analysis by t-SNE and found that the distribution of the subtypes was largely consistent with the two-dimensional t-SNE distribution pattern ( Figure 1C ). We observed significant prognostic differences among the three subtypes in the TCGA data set ( Figure 1D ). Compared with C2 and C3, C1 had a shorter median survival time ( Figure 1D ). Subsequently, the OV samples from the GEO datasets were used in independent verification, using the aforementioned k=3 classification, which also revealed three different molecular subtypes. Obvious differences between the three subtypes were observed in the GSE63885 ( Figure 1E ) and GSE17260 ( Figure 1F ) datasets, and the overall survival (OS) of the C1 subtype was significantly shorter than that of C2 and C3 subtypes.

Considering that the subtyping of OV samples was based on metabolism-related genes, we used the ssGSEA algorithm to quantitatively analyze the metabolic process. We uncovered multiple metabolic pathways that were significantly different among the three subtypes. Such as, the pathways predominant in the C1 subtype were the inositol phosphate metabolism (some specific metabolic genes: ALDH6A1, CDIPT, IMPA1, IMPA2), reactome inositol phosphate metabolism (some specific metabolic genes: CALM1, IMPA1, IMPA2, INPP1), wp melatonin metabolism and their effects (some specific metabolic genes: AANAT, ACHE, ADRB1, APOE) were significantly higher than the other two subtypes. Also, there were some pathways, certain differences between the C2 and C3 subtypes, such as kegg glyoxylate and dicarboxylate metabolism (some specific metabolic genes: ACO1, ACO2, AFMID, CS), kegg pyrimidine metabolism (some specific metabolic genes: AK3, CAD, CANT1, CDA) and wp pyrimidine metabolism (some specific metabolic genes: CAD, CDA, CMPK1, CMPK2) ( Figure 2A ). We further explored the immune characteristics between different subtypes. First, the immune process was quantified using the ssGSEA algorithm and differential analysis was performed to discover subtype-specific immune characteristics. The results showed that the immune characteristics of the three subtypes were clearly classifiable, and differences in the enrichment scores of B cells, Th1 cells, Th2 cells, TFH, Tgd, NK CD56dim cells, macrophages, mast cells, normal mucosa, blood vessels, lymph vessels were significant between C1 and other subtypes ( Figure 2B ). In addition, C1 had higher cell cycle, RTK RAS, WNT, and angiogenesis pathway activation scores than C2 and C3 ( Figure 2C ). Moreover, the ESTIMATE algorithm was used to calculate the immune scores and significant differences in the immune scores were observed between the three groups with the immune score of C3 being significantly lower than that of C1 and C2 ( Figure 2D ). The immune score of C1 was the highest ( Figure 2D ). Meanwhile, the stromal scores showed a tendency similar to that of the immune score ( Figure 2E ).

Due to the significant differences in the immune scores between the subtypes, the MCP and ssGSEA algorithms were further used to conduct in-depth research on immune infiltration. The abundance of 16 immune-related cell types was calculated using MCP counter and ssGSEA algorithm and displayed in the heatmap ( Figure 3A ). The results showed that there were significant differences in the level of immune infiltration between the three subtypes ( Figure 3A ), suggesting that the immune characteristics of the three subtypes are reproducible. In addition, the enrichment score of T cells, B lineage, monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells, fibroblasts, Tcm, Tem and TReg in C1 were higher than in C2 and C3 subtypes.( Figure 3B ). We further studied the association between the subtypes and the expression of 12 potentially targetable immune checkpoint genes, which were selected based on current inhibitor drugs in clinical trials or drugs approved for specific cancer types, and the results indicated remarkable differences in the expression of these genes between the three subtypes, except for FGL1, CTAG1B and MAGEA4 ( Figure 3C ). The immune checkpoint genes TNFRSF9 showed higher expression in C1 compared to the other two subtypes ( Figure 3C ).

We further explored the relationship between different subtypes and clinical symptoms in order to reveal their clinicopathological characteristics. Based on the TCGA dataset, the results of the chi-square test showed that the correlation between clinicopathol ogical characteristics and OV subtypes in the TCGA cohort was not obvious ( Figure 4A ). However, in validation set 1 ( Figure 4B ), the proportion of chemotherapy and residuals in the three groups was consistent with the results of KM analysis in Figures 1E–G whereas in validation set 2, the grade was higher in C1 ( Figure 4C ). In addition, we used GSEA to analyze the differences in signaling pathways between the subtypes. GO analysis results showed that the C1 subtype was mainly enriched in the extracellular matrix structural constituent and collagen binding signaling pathways while the C2 subtype was mainly enriched in the structural constituent of ribosome and ribosomal subunit signaling pathways ( Figure 4D ). The C3 subtype was mainly enriched in microtubule bundle formation, plasma membrane and cell projection cytoplasm signaling pathway ( Figure 4D ). KEGG analysis results showed that the C1 subtype was mainly enriched in the ECM-receptor interaction and focal adhesion signaling pathways ( Figure 4E ). The C2 subtype was mainly enriched in the ribosome and oxidative phosphorylation signaling pathways while the C3 subtype was mainly enriched in basal cell carcinoma and hedgehog signaling pathway signal pathway ( Figure 4E ).

The instability of the genome promotes the accumulation of mutations in cancer cells and leads to the rapid evolution of the cancer genome in response to the tumor microenvironment and treatment-induced related stresses during evolution (26). Detecting and characterizing these tumor somatic mutations has become an important means to analyze the occurrence and development of tumors. In order to study the differences in the frequency of somatic mutations among OV subtypes, and to observe different mutation patterns among OV subtypes, the somatic mutation data from the TCGA database was analyzed. The mutation characteristics among the three subtypes were shown in Figure 5A . Most of the mutations in the three subtypes were nonsense mutation, but C2 has fewer mutated genes than C1 and C3 ( Figure 5A ). On the other hand, in a multi-group comparison, the results showed that the C2 subtype had a higher number of neoantigens than the C1 and C3 subtypes but not significant ( Figure 5B , P>0.05). Moreover, the tumor mutation burden in the C2 subtype was slightly higher than that of the C1 and C3 subtypes but not significant ( Figure 5C , P>0.05).

After comprehensively considering accuracy and clinical application potential, the top 30 genes with the largest log₂FC values in each subtype were selected for the development of subtype classification. Therefore, a 90-gene classifier was generated and visualized ( Figure 6A ). Subsequently, the 90-gene classifier was used to repeat the subtype prediction in the GEO dataset. The results showed a consitency between the NMF and the NTP predicted clusters, confirming that the subtypes were predictable based on the 90-gene classifier ( Figure 6B ). In addition, based on the pRRophetic oftware package, we predicted the sensitivity of each OV subtypes to anticancer drugs. The results showed that the C2 subtype had a lower IC50 for BMS.708163 than C1 and C3, while the C1 subtype had a lower IC50 for AZD6482 than other two subtypes ( Figure 6C ). Moreover, we further based on the data set of melanoma immunotherapy to predict the sensitivity of the three subtypes to anti-tumor immunotherapy. The results showed that the C2 subtype was more sensitive to immunotherapy ( Figure 6D ).

Through the previous subtype analysis, we divided the samples into three subgroups: C1, C2 and C3. We performed difference analysis among subgroups (C1_vs_C2( Figure 7A ), C3_vs_C1( Figure 7B ) and C3_vs_C2( Figure 7C ) respectively) (|log₂FC| > 0.585 and. p value < 0.01). According to the Wayne analysis ( Figure 7D ), 1111 genes were differentially expressed between C2 and C1 groups; 1064 genes were differentially expressed between C3 and C1 groups; 1127 genes were differentially expressed between C3 and C2 groups. In the three differential analysis, 67 genes were differentially expressed at the same time. Finally, we analyzed the survival of these 67 genes and found that the expression of two genes was very significant in the survival analysis (p value < 0.01). The sequence is CXCL11( Figure 7E ) and PXDN ( Figure 7F ).

We used immunohistochemistry to stain the tissue chip, and then scored according to the staining intensity and staining area. Figure 8A shows the staining corresponding to different scores. After comparing the score with the patient information, we found that the expression of PXDN was inversely proportional to the tumor size ( Figure 8B ) and H-score increases with the increase of pathological grade positively correlated with the pathological grade ( Figure 8D ). The expression of CXCL11 was inversely proportional to the tumor size ( Figure 8C ), while the H-score was no correlation with the pathological grade ( Figure 8E ). Furthermore, it was found that there was an inverse correlation between tumor size and pathological grade ( Figure 8F ).