The TCGA database (https://portal.gdc.cancer.gov/) included the gene transcripts and clinical details of 429 ovarian cancer tissue samples from patients with the disease. The patient’s clinical features encompassed survival status, time, tumor grade, age, etc. In the meantime, the GTEx database (https://www.gtexportal.org/home/) was accessed to download 88 normal ovarian samples. For validation, the gene expression profiling microarray datasets for ovarian cancer tissues were acquired from the GEO database (https://www.ncbi.nlm.nih.gov/geo). These datasets, GSE26712 (Zheng et al., 2019) and GSE102073 (Ye et al., 2021), each contained 153 and 84 ovarian cancer tissues, respectively. Additionally, ovarian cancer single-cell datasets EMTAB8107 (Ding et al., 2024), GSE118828 (Yu et al., 2022), GSE130000, and GSE154600 (Jiang et al., 2023) were downloaded to explore gene expression at the single-cell level.

Inflammation-related gene sets were obtained from the Molecular Signatures Database (MSigDB) (Castanza et al., 2023) (https://www.gsea-msigdb.org/gsea/msigdb/), including BIOCARTA_INFLAM_PATHWAY (v2023.2.), GOBP_CHRONIC_INFLAMMATORY_RESPONS (v2023.2), HALLMARK_INFLAMMATORY_RESPONSE (v2023.2.), and REACTOME_INFLAMMASOMES (v2023.2.).

With the help of the GSVA (Hänzelmann et al., 2013) R package (v2.0.4), we evaluated the gathered sets of gene sets associated with inflammation using single-sample gene set enrichment analysis (ssGSEA). ssGSEA is an extension of the GSEA method, primarily designed for individual samples where GSEA is not applicable. The algorithm uses the empirical cumulative distribution function to calculate enrichment scores (ES) and rank normalized gene expression values for a given sample. The prognostic correlation between gene sets associated with inflammation and patients with ovarian cancer was assessed simultaneously using the Cox proportional hazards model.

Consensus clustering was often used in cancer subtype classification studies using the Consensus ClusterPlus R package (v4.12.6) (Wilkerson and Hayes, 2010). This study’s ovarian cancer subtype classification was conducted based on the aforementioned inflammation-related gene set scores. The optimal clustering effect was determined by combining the consensus cumulative distribution function (CDF) with the Proportion of Ambiguous Clustering (PAC) score. In the CDF plot, the consensus matrix’s cumulative distribution function was displayed for different values of k (represented by colors), aiding in identifying the approximate maximum CDF value, where consensus and cluster confidence are maximized, resulting in the most reliable clustering analysis. In PAC analysis, lower PAC values indicate more ideal clustering effects.

Currently, the most widely used method for survival analysis is the Kaplan-Meier approach. The KM approach, as it is commonly called, was proposed by Kaplan and Meier. The Kaplan-Meier survival analysis compares the survival circumstances of two patient groups using a univariate analysis that integrates patients’ survival times and terminal states. The Kaplan-Meier survival curve, a commonly encountered representation, visually reflects survival differences under various conditions.

Based on the IOBR (Zeng et al., 2021) R package (v2.0), we employed built-in algorithms such as TIMER, CIBERSORT, MCPcounter, EPIC, and quanTIseq (Newman et al., 2015; Becht et al., 2016; Finotello et al., 2019; Racle et al., 2017; Li et al., 2016) to assess the abundance of immune cell infiltration in the tumor immune microenvironment of each ovarian cancer tissue.

The OncoPredict (Maeser et al., 2021) R package (v1.2) was created by Maeser et al. and was used to predict medication reactions in cancer patients. OncoPredict adapts tissue gene expression patterns to the semi-maximal inhibitory concentration (IC50) of drugs taken from cancer cell lines in the Genomics of Drug Sensitivity in Cancer (GDSC) database and the Cancer Cell Line Encyclopedia (CCLE) maintained by the Broad Institute. An unpaired t-test was used to assess the sensitivity of 198 drugs (between high-risk and low-risk groups). Set at p < 0.05 was the significance level.

The Cancer Immunome Atlas (TCIA) (Charoentong et al., 2017) database (https://tcia.at/) was used to download the Immunophenotype Scores (IPS) for CC. Subsequently, IPS were compared across different tumor groups to predict sensitivity to immunotherapy.

Using the limma (Ritchie et al., 2015) R package (3.60.4), differential gene expression analysis was performed on the TCGA data. This involved data preprocessing, normalization, and identifying significant differences in gene expression levels through linear modeling. Statistical thresholds were set (adj. P. Val. < 0.01 and |log2(FC)| > 1) to screen for significantly differentially expressed genes. Finally, to decipher the biological significance of these differential genes, enrichment analysis and functional annotation were performed.

The clusterProfiler (Wu T. et al., 2021) R package (v4.12.6) was used to perform enrichment analysis, which consisted of two steps: (1) Over-Representation Analysis (ORA) to investigate the functional enrichment of gene sets through Genomes (KEGG) analyses, Kyoto Encyclopedia of Genes, Gene Ontology (GO), and (2) Gene Set Enrichment Analysis (GSEA) to examine the enrichment of validated gene sets in KEGG pathways. These enrichment results revealed the gene sets’ biological functions and pathway associations.

In systems biology, co-expression gene modules and their correlation to phenotypes are identified using Weighted Gene Co-Expression Network Analysis (WGCNA) (Langfelder and Horvath, 2008). Gene co-expression networks are constructed by calculating gene-gene correlations and then converting the correlation matrix into a weighted matrix. Gene modules are then built based on the weighted matrix, and the eigengene for each module is computed. Subsequently, the correlation between module eigengenes and phenotypic data determines the module-phenotype relationship. Ultimately, gene modules associated with phenotypes of interest are identified, revealing underlying biological pathways and mechanisms. The WGCNA method aids in uncovering key modules and gene correlations within gene regulatory networks.

The TCGA dataset was initially used to construct a survival prognostic model using the multiCox and Least Absolute Shrinkage and Selection Operator (LASSO) techniques. Gene features were selected using LASSO, and multiCox was employed for multivariable Cox regression to establish the prognostic model. The model was then used for independent external validation datasets, and its ability to predict survival was evaluated using Kaplan-Meier survival analysis and time-dependent ROC curve analysis. These validation analyses verified the prognostic predictive efficacy of the model across different datasets, ensuring its reliability and generalizability.

Derived from the PCBC database’s mean-centered RNA-Seq data of PSCs (syn2701943) (Gul et al., 2023). A stem cell feature signature was found using the One-Class Logistic Regression (OCLR) machine learning approach, and it was confirmed using leave-one-out cross-validation. A Spearman correlation analysis was then used to compare the stem cell features and the normalized expression matrix of tumor samples. Finally, by scaling the Spearman correlation coefficient between 0 and 1, the mRNA Stemness Index (mRNAsi) was determined. A higher mRNAsi indicates a higher degree of tumor dedifferentiation and stronger stemness.

The following methods were used to process ovarian cancer single-cell sequencing data: We first converted the scRNA-seq data into a Seurat (Hao et al., 2024) object using the Seurat R package. We performed quality control (QC) by determining the percentage of ribosomal or mitochondrial genes and eliminating low-quality cells. FindVariableFeatures was used to determine the top 2000 genes exhibiting high variability. Furthermore, dimensionality reduction techniques were used to group approximately 2000 genes using Principal Component Analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP). We could identify marker genes for each cluster using the FindAllMarkers tool with |Log2FC| and min. 0.3 and 0.25 are the respective pct cutoff values. Different cell types were annotated using the SingleR (Aran et al., 2019) R package. Finally, we used the AddModuleScore function to compute the expression levels of prognostic model genes at the single-cell level.

Fenghuishengwu in China is where human ovarian surface epithelial cells (HOSE) are sourced (HOSE, CL0154). From the American Type Culture Collection, the SKOV3 cell line was acquired. Mycoplasma is routinely tested for in all cell lines. A complete medium, consisting of 1% double antibiotics and 10% fetal bovine serum (FBS), is used to cultivate both HOSE and SKOV3. When the cell confluence reaches 80%–90%, they are passaged.

Similarly, 5 × 10⁵ cells were cultured in each well of a 6-well plate. For transfection, HOSE and SKOV3 cells were treated with 15 nM of siRNA IL-6, siRNA TGF-β1 and siRNA NC (P4157, GenePharma, Shanghai, China), respectively, using Lipofectamine 3000 (L3000150, Thermo, New York, Waltham, MA, United States).

HOSE and SKOV3 cells were lysed to obtain total RNA using a TRIzol reagent (15,596,026, Invitrogen, New York, NY, United States). The qRT-PCR analysis was performed using the HiScript II One Step qRT-PCR SYBR Green kit (P131, Vazyme, Nanjing, China) and a Bio-Rad CFX96 PCR system (Bio-Rad, Hercules, CA, United States). RuiBiotech (Beijing, China) created and manufactured the primers used in this investigation. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal reference for assessing the target gene expression using the 2^−ΔΔCT method.

The total protein was extracted using a cell lysate solution, and the proteins were separated using 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). PVDF membranes were then used to hold the separated proteins (03010040001, Millipore, Billerica, MA, United States). After a 30-min blocking, the membranes were subjected to adding primary antibodies and incubated at 4°C. The primary antibodies were obtained from Abcam (Abcam, United States). Secondary antibodies (BA1054, 1:5000, Boster, Wuhan, China) were added to the membranes and incubated for 2 h at room temperature after washing. The membranes were then visualized using an ECL development kit (A38554, Invitrogen, New York, NY, United States) and photographed with a GE Las-4000 (GE Healthcare, Piscataway, NJ, United States). After conducting the experiment thrice, the gray values were obtained using Media Cybernetics’ ImageJ 1.8.0 program (Silver Spring, MD, United States). An internal reference was beta-actin.

The analysis was conducted using SPSS 26.0 and the R programming language. The measurement data was expressed using the standard deviation (x ± s). A one-way ANOVA was employed to compare the groups. Dunnett’s multiple comparisons were performed to determine whether the variance was uneven. Measurement data with a normal distribution were displayed as mean ± standard deviation, and t-tests were used to compare groups. The Mann-Whitney U test was performed using measurement data that was not normally distributed and displayed as the median and interquartile range to compare groups. Count data were expressed as rates, and group comparisons were conducted using the χ2 test.

First, we obtained four inflammatory-related gene sets from the MSigDB database, including BIOCARTA_INFLAM_PATHWAY, GOBP_CHRONIC_INFLAMMATORY_RESPONS, HALLMARK_INFLAMMATORY_RESPONSE, and REACTOME_INFLAMMASOMES. Most inflammatory-related genes had significantly different expression profiles in tumor and normal tissues, with most genes significantly elevated in tumor tissues, according to our analysis of the TCGA ovarian cancer dataset and the corresponding normal tissues from the GTEx database (Figure 1). This suggests a correlation between inflammatory phenotypes and tumor development.

Meanwhile, based on ssGSEA, we performed enrichment scoring of inflammatory-related gene sets for each ovarian cancer tissue (Figure 2A). Univariate COX regression analysis identified three inflammation-related gene sets as significantly and positively associated with better prognosis (HR < 1) in ovarian cancer patients. These gene sets are thought to be protective prognostic factors (Figure 2B). This suggests these inflammatory-related gene sets’ potential research value and clinical significance. We, therefore, performed consensus clustering analysis using the Consensus ClusterPlus R package based on the enrichment scores of these four inflammatory-related gene sets. The maximum number of clusters was set to 10, with 100 subsamples drawn, a sample proportion of 0.8, K-means as the clustering algorithm, and Euclidean distance as the metric. Ultimately, we performed nine clusters with k values ranging from 2 to 10. Through comprehensive evaluation using CDF curves and PAC analysis, we selected the ideal number of clusters as 2 (Figures 2C–E). At the same time, we found significant differences in inflammatory enrichment scores (Figure 2F) and overall survival rates (Figure 2G, log-rank p = 0.076) among patients with different inflammatory gene expression patterns. Chi-square tests for clinicopathological features revealed differences in the age distribution (with 65 years as the cutoff) and clinicopathological grading among patients in different groups (Figure 2H).

Previous studies have reported that different tissue types often have distinct immune infiltration microenvironments. Thus, we used five immune microenvironment analysis methods in this study—CIBERSORT, MCPcounter, quanTIseq, EPIC, and TIMER—for integrated evaluation and analysis of immune cell infiltration profiles to thoroughly examine the immune profiles among various subtypes. We found that the C2 subtype had much more infiltrating NK cells, B cells, macrophages, CD8⁺ T cells, and CD4⁺ T cells than the C1 subtype, with consistent results across various analysis methodologies. This suggests that the C2 subtype exhibits the biological characteristics of a so-called “hot” tumor immune microenvironment (Figure 3A). In the meanwhile, we found that the C2 subtype had far higher amounts of immunomodulators and cytokines expressed than the C1 subtype, based on a list of genes encoding immunomodulators and chemokines that we downloaded from the TISIDB database (Figure 3B).

With the deepening of research on the immune tumor microenvironment, substantial evidence suggests that tumors with different levels of immune cell infiltration have distinct response rates to chemotherapy and immunotherapy. Therefore, we examined the IC50 values of vinblastine, paclitaxel, docetaxel, and cisplatin. We found that, except cisplatin, the C2 subtype’s IC50 values were considerably higher than the C1 subtype’s, indicating a noticeably lower response rate of the C2 subtype to these three drugs (Figure 4A). Simultaneously, we utilized IPS, IPS-PD1/PD-L1/PD-L2, IPS-CTLA4, and IPS-PD1/PD-L1/PD-L2 + CTLA4 to assess differences in the response rates to immune checkpoint inhibitor therapy among different subtypes. We found that IPS, IPS-PD1/PD-L1/PD-L2 and IPS-CTLA4 were significantly higher in the C2 subtype compared to the C1 subtype (Figure 4B).

We hypothesize that the distinct immune infiltration microenvironments and responses to drug treatment among different subtypes are based on significantly different signaling pathways and biological differences. Therefore, we identified genes that were differently expressed among various subtypes using the limma R package (Figure 5A). We carried out an over-representation analysis (ORA), which included GO_BP/CC/MF (Figure 5B), and last, we performed a GSEA analysis. The results suggested multiple signaling pathways were significantly deregulated in C2 (Figures 5C, D). In summary, subtypes based on inflammatory gene enrichment scores exhibit distinct immune infiltration microenvironments, responses to drug treatment, and deregulated signaling pathways, warranting further investigation.

To conduct the WGCNA analysis, we first added the differentially expressed genes identified in the previous step. Four co-expression modules were obtained after setting the soft threshold β to three and the minimum number of genes in a module to 30 (Figures 6A, B). We utilized GO enrichment analysis and found that modules other than the gray module possessed distinct biological characteristics. Since previous analyses suggested that patients in the C2 subtype had a better prognosis, responsiveness to drug treatment, and deregulated signaling pathways, we hypothesized that genes significantly associated with C2 might be involved in tumor development, invasion, and resistance to drug treatment. Based on WGCNA co-expression network analysis, we found that the C2 subtype had the strongest positive correlation with the turquoise module (Figure 6C, Cor = 0.54), which contained 2323 genes. We could run functional enrichment analyses on these genes using thresholds (MM > 0.6 and GS > 0.3) to identify key genes inside the module (Figures 6D, E). The results showed that immune receptor activity and other biological processes were the primary roles of the module’s key genes (Figure 6F).

The key genes found in the previous step served as the prognosis model’s input genes. In this study, we used the TCGA dataset to train the model. Then, we assessed the model’s prognostic efficacy using independent external datasets based on time-dependent ROC curves and Kaplan-Meier survival analysis.

With the help of multivariate COX regression analysis, Least Absolute Shrinkage and Selector Operation (LASSO) (Figure 7A), and module key genes, we constructed a prognosis model in this study. The genes included in the model and their corresponding regression coefficients are shown (Figure 7B). Risk scores based on the prognostic model were simultaneously computed for all ovarian cancer samples in the training and validation sets. The samples were divided into low-risk and high-risk groups based on the median risk score. Our result shows the plotted Kaplan-Meier curves, which indicate substantial differences between the two groups (Figures 7C–E). Using a time-dependent ROC curve analysis, the prediction efficiency at 1, 3, and 5 years was also evaluated, and the results indicated that the prognostic model had good prediction performance.

Based on the CIBERSORT immune microenvironment analysis algorithm, we analyzed the TCGA dataset divided into high-risk and low-risk groups. We found that the low-risk group had higher M1 macrophages and CD8⁺ T cell infiltration abundances, while the high-risk group had significantly higher abundances of M2 macrophage infiltration (Figure 8A). This suggests an active anti-tumor immune response in the tissue immunological microenvironment of the high-risk group. Further analysis of T cell exhaustion markers and M2 macrophage markers showed that the high-risk group had significantly higher expression levels than the low-risk group (Figures 8B, C). Simultaneously, analysis of cellular stemness levels also indicated higher stemness scores in the high-risk group, suggesting more pronounced tumor stemness, i.e., dedifferentiation (Figure 8D). Since cellular stemness levels are negatively correlated with prognosis, these results imply a poorer prognosis for patients in the high-risk group.

The amount of immune cell infiltration substantially impacts immunotherapy and chemotherapy, according to previous studies, with the high-risk group having less anti-tumor immune cell infiltration than the low-risk group. To determine if the high-risk group responded less frequently to these four drugs, we examined the IC50 values of cisplatin, vinblastine, paclitaxel, and docetaxel. Compared to the low-risk group, the IC50 values of the high-risk group were much greater (Figure 9A). We employed IPS, IPS-PD1/PD-L1/PD-L2, IPS-CTLA4, and IPS-PD1/PD-L1/PD-L2 + CTLA4 simultaneously to evaluate variations in immune checkpoint inhibitor therapy response rates across various risk groups. IPS, IPS-PD1/PD-L1/PD-L2, and IPS-CTLA4 levels were considerably lower in the high-risk group than in the low-risk group (Figure 9B).

We hypothesize that the distinct immune infiltration microenvironments and drug responsiveness observed among risk groups are based on significantly different signaling pathways and biological differences (Zeng et al., 2021). Consequently, we initially employed the limma R package to find genes differently expressed amongst various risk groups. Then, using cancer hallmarks as a guide, we ran a GSEA analysis on these differentially expressed genes. According to the findings, the high-risk group had significantly higher levels of several signaling pathways and significantly lower levels of others, including IL6-JAK-STAT3, hypoxia, and glycolysis. In summary, the prognostic model based on inflammatory genes exhibits distinct immune infiltration microenvironments, drug responsiveness, and deregulated signaling pathways, which has important implications for clinical decision-making (Figure 10).

We integrated four single-cell datasets and obtained 88,089 ovarian cancer single cells after quality control, dimensionality reduction, and clustering. With a clustering resolution set to 0.4, we identified 20 clusters (Figures 11A, B). Out of these 20 clusters, we identified 12 cell subpopulations by combining manual annotation methods with SingleR automatic annotation. We then analyzed the risk scores of these 12 cell subpopulations (Figures 11C, D). The results showed that the prognostic model genes related to inflammation were primarily expressed in immune-related cells, further validating the findings from traditional transcriptome analysis at the single-cell level.

Following the bioinformatics analysis, the study explored the expression levels of CCL2, IL10, IL6, and TGFβ1 in normal ovarian epithelial cells (HOSE) and ovarian cancer cells (SKOV3). The research team used qPCR to detect the mRNA expression levels in cells. Compared with HOSE, CCL2 (p = 0.003), IL10 (p = 0.003), IL6 (p = 0.002), and TGFβ1 (p = 0.002) were all highly expressed in SKOV3, and the differences were statistically significant (Figure 12A). Subsequently, the study used siRNA IL6 and TGFβ1 to transfect the SKOV3 cell line. qPCR results revealed that following transfection, there was a drop in the mRNA expression levels of TGFβ1 (p = 0.0002) and IL6 (p = 0.002) (Figures 12B, C). Further, CCK8 was used to evaluate the proliferation of SKOV3 after transfection. Following the knockdown of IL6 and TGFβ1, the capacity of SKOV3 cells to proliferate decreased in comparison to the siRNA NC group. There was a statistically significant difference (p < 0.0001, p < 0.0001) (Figure 12D). Finally, we used Western blot to further evaluate the protein expression levels of IL6 and TGFβ1 in HOSE and SKOV3 cells. The SKOV3 group had higher levels of TGFβ1 (p = 0.0001) and IL6 (p = 0.004) protein expression compared to HOSE, and these changes were statistically significant (Figures 12E, F).