Abstract

   Ovarian cancer represents a severe gynecological malignancy with a dire
   prognosis, underscoring the imperative need for dependable biomarkers
   that can accurately predict drug response and guide therapeutic
   choices. In this study, we harnessed online single-cell RNA sequencing
   (scRNAseq) and bulk RNA sequencing (RNAseq) datasets, applying the
   Scissor algorithm to identify cells responsive to paclitaxel. From
   these cells, we derived a gene signature, subsequently used to
   construct a prognostic model that demonstrated high sensitivity and
   specificity in predicting patient outcomes. Moreover, we conducted
   pathway and functional enrichment analyses to uncover potential
   molecular mechanisms driving the prognostic gene signature. This study
   illustrates the critical role of scRNAseq and bulk RNAseq in developing
   precise prognostic models for ovarian cancer, potentially transforming
   clinical decision-making.

   Keywords: Ovarian cancer, Paclitaxel resistance, scRNAseq

1. Introduction

   Ovarian cancer (OV) ranks as the fifth leading cause of cancer-related
   mortality among women, with the incidence increasing significantly
   after the age of 50 [[27]1]. Despite a general five-year survival rate
   of around 50 %, prognosis dramatically varies depending on the stage at
   diagnosis [[28]2]. Often termed the ‘silent killer,’ OV's nondescript
   symptoms frequently lead to late diagnosis, compounding treatment
   challenges [[29]3]. Chemotherapy, incorporating agents such as
   paclitaxel, carboplatin, and doxorubicin, remains the cornerstone of
   post-surgical treatment for OV [[30]4,[31]5]. However, patient
   responses to these chemotherapeutic regimens differ markedly due to
   factors such as cancer stage, tumor cell heterogeneity, and specific
   drug efficacy [[32]6]. Consequently, there is a critical need to
   stratify patients into distinct groups based on their chemotherapy
   responses [[33]7], enabling clinicians to tailor more effective
   treatment strategies for OV.

   Previous studies have leveraged RNAseq to investigate drug responses in
   OV cells, illuminating the genomic landscape of chemoresistance
   [[34]8,[35]9]. For instance, one study employed RNAseq to discern a
   gene set distinctly expressed between drug-sensitive and drug-resistant
   OV cells [[36]10]. This research revealed that drug-resistant cells
   exhibit elevated expression of genes related to drug detoxification and
   cell survival, shedding light on the molecular underpinnings of drug
   resistance in OV. Another study utilized RNAseq to profile the
   transcriptomes of OV cell lines, differentiating those sensitive to
   paclitaxel from those resistant to this commonly used chemotherapeutic
   agent [[37]11]. These investigations underscore the value of RNAseq in
   elucidating drug responsiveness in OV, offering potential pathways for
   targeted therapeutic interventions.

   Over the past decade, scRNAseq has revolutionized our understanding of
   cell biology and disease mechanisms [[38]12]. Unlike traditional bulk
   RNA sequencing, which measures gene expression from a mixed cell
   population potentially obscuring individual cell differences, scRNAseq
   assesses gene expression in each individual cell. This approach
   provides a nuanced view of the cellular diversity and states within
   tissues. Cancer cells exhibit significant heterogeneity, influenced by
   their varied tumor microenvironment (TME) and interactions with
   distinct cell types [[39]13]. This heterogeneity even extends to clonal
   cancer cells, which can display diverse phenotypic features through
   clonal evolution. The TME, comprising immune cells, stromal cells,
   blood vessels, the extracellular matrix, and signaling molecules, plays
   a crucial role in drug resistance [[40]14]. Drug resistance mechanisms
   facilitated by the TME include the secretion of signaling molecules
   such as cytokines and growth factors that activate pathways like
   PI3K/AKT, promoting cell survival and chemotherapy resistance [[41]15].
   Additionally, the TME can physically impede drug delivery to tumor
   cells through increased interstitial fluid pressure, abnormal blood
   vessel structures, or extracellular matrix remodeling. Dynamic changes
   in the TME in response to treatment, such as the induction of growth
   factors post-chemotherapy, further complicate drug resistance by
   promoting tumor cell proliferation and survival [[42]16].

   To delineate the heterogeneity inherent in OV, recent studies have
   employed scRNAseq to dissect subpopulations of cancer and stromal cells
   within the tumor microenvironment (TME), offering novel insights into
   ovarian cancer progression, metastasis, clonal evolution, and drug
   resistance [[43][17], [44][18], [45][19], [46][20], [47][21], [48][22],
   [49][23]]. Despite its profound impact, single-cell technology remains
   prohibitively expensive, limiting the number of samples that can be
   sequenced. Moreover, the practical challenges of treating patients with
   varying drug combinations compared to controlled in vitro experiments
   pose significant hurdles. To bridge these gaps, algorithms have been
   developed that integrate bulk RNAseq data from large cohorts or cell
   line assays with scRNAseq datasets. These algorithms aim to identify
   cell subpopulations driving critical phenotypes such as disease stage,
   tumor metastasis, treatment response, and survival outcomes. One such
   computational method, SCISSOR, is designed to pinpoint subpopulations
   within a single-cell dataset that correlate with specific phenotypes of
   interest [[50]24]. This method holds promise for advancing the
   development of cell type-specific therapies and the discovery of
   prognostic biomarkers.

   In this study, we utilized the Scissor analysis to integrate bulk
   RNAseq data from paclitaxel-resistant OV cells with scRNAseq data from
   OV patient samples. We examined the association between cellular
   composition, cell-cell communication, and clinical outcomes. This
   analysis aims to enhance our understanding of paclitaxel responsiveness
   across different OV patient groups and to establish a foundation for
   personalized management strategies for OV. The insights gained could
   lead to more targeted and effective treatments, tailored to the unique
   cellular landscapes of individual patients.

2. Materials and methods

2.1. Acquisition of datasets and data preprocessing

   Three scRNA-seq datasets [51]GSE154600, [52]GSE158937 and [53]GSE184880
   of total 15 OV samples were first downloaded from the Gene Expression
   Omnibus (GEO, [54]https://www.ncbi.nlm.nih.gov/geo/) database. We then
   used “Seurat” R package (version 4.1.1) to process the datasets
   [[55]25]. Initially, the percentage of mitochondrial genes, ribosomal
   genes and erythrocyte genes were determined using the
   PercentageFeatureSet function. Cells with a gene number <200 and
   >5,000, mitochondrial gene content >10 %, ribosomal gene content <20 %
   and erythrocyte gene content >10 % and were excluded. Subsequently,
   scRNA-seq data were normalized by the NormalizeData function, and the
   top 2000 genes with highly variable features were identified by
   FindVariableFeatures function. We used ScaleData function to
   standardize the count data. The batch effects between different
   datasets were removed by “harmony” R package (version 0.1.1) [[56]26].
   The top 2000 highly variable genes were used for principal component
   analysis (PCA). Thus, the top 30 principal components were manually
   selected for cell clustering analysis using the Uniform Manifold
   Approximation and Projection (UMAP). Canonical marker genes of specific
   cell types in were used to annotate the clusters of cells. We used
   FindAllMarker function to identify the genes upregulated in the
   Scissor-cells. “DoubletFinder” R package (version 2.0.3) was carried
   out to remove the potential doublet cells in the scRNAseq dataset with
   default parameters [[57]27].

   Moreover, the bulk RNA-seq data of TCGA-OV samples were accessed from
   TCGA ([58]https://portal.gdc.cancer.gov/) database. The statistical
   results of somatic mutation in TCGA dataset were visualized with the
   “maftools” R package software (version 2.12.0). Additional OV samples
   were obtained from the [59]GSE14764, [60]GSE23554 and [61]GSE26712
   datasets. In this study, OV samples from patients without a complete
   follow-up information were excluded. Bulk RNAseq data of paclitaxel
   resistance in ovarian cancer cell lines was obtained from the
   [62]GSE172016.

2.2. Cell-cell communication analysis by CellChat

   To perform CellChat (version 1.6.1) analysis, we first removed
   low-quality cells and genes, normalized the expression values, and
   performed clustering to identify cell types using Seurat R package. We
   then used “generate_cell_chat” function to generate the cell-cell
   communication network, which contains information about the
   ligand-receptor interactions between different cell types. To visualize
   the cell-cell communication network, the “plot_cell_chat” function was
   performed.

2.3. Functional and immune cell infiltration analysis

   We used ClusterProfiler R package (version 4.4.4) to perform GO and
   KEGG analysis. To perform GO analysis, we utilized the enrichGO
   function, which performs enrichment analysis using the GO database. To
   perform KEGG pathway analysis, we used the enrichKEGG function, which
   performs enrichment analysis using the KEGG database. To carry out
   immune infiltration analysis in the TCGA-OV datasets, we uploaded our
   data to the Timer website ([63]http://timer.cistrome.org/), which
   allowed us to obtained the immune cell infiltration information of
   TCGA-OV sample based on different algorithms, including CIBERSORT.

2.4. Identification of prognostic genes and construction of scissor genes
prognosis risk model

   The LASSO Cox regression analysis was performed to construct a
   predictive signature. The R package “glmnet” (version 4.1.4) was used
   to achieve the variable selection and shrinkage of the LASSO algorithm.
   The risk scores of each patient were calculated by the following
   formula: Risk score = sum (coeffcients* expression of gene n). The
   “survminer” (version 0.4.9) and “survival” R packages (version 3.4.0)
   were used to determine the Kaplan-Meier survival curve. The nomogram of
   the OV prognosis model was established via univariate and multivariate
   analyses, combined with clinical characteristics. Then, calibration
   curves were constructed for the prediction of survival of different
   time points in OV.

2.5. Statistical analysis

   R software (version 4.0.0) were used for statistical analysis.
   Continuous variables are presented as the mean ± standard deviation.
   Univariate and multivariate Cox regression analyses were used to
   demonstrate their prognostic value. Kaplan–Meier curves with a
   two-sided log-rank test were used to compare the overall survival (OS)
   of the patients in the high and low expression groups. p-values <0.05
   were considered statistically significant.

3. Results

3.1. The cell atlas of OC

   To dissect the heterogeneity of tumor cells and cells in tumor
   microenvironment, we analyzed the transcriptome of 15 OC patients from
   3 public scRNA-seq datasets which contains high quality of sequencing
   data. After a quality screening based on the percentage of the
   expression of the mitochondrial, ribosomal and erythrocyte genes, we
   obtained a total of 52,874 high-quality cells with 29,126 total genes
   for further study ([64]Fig. 1A). We identified twelve clusters
   (clusters 0–11) in these OC sample through dimension reduction and
   cluster analysis ([65]Fig. 1B). Canonical markers for different cell
   type were then used to annotate each cluster of cells and revealed six
   major cell types, including epithelial cells, fibroblasts, endothelial
   cells, T cells, B cells and myeloid cells ([66]Fig. 1C). [67]Fig. 1D
   illustrated typical markers for different cell types.

Fig. 1.

   [68]Fig. 1
   [69]Open in a new tab

   Characterization of major cell populations of OV patient samples.

   (A) UMAP clustering of OV cells integrated from three different
   scRNAseq datasets including 15 OV samples. Each point indicates one
   single cell and cells in the same cluster indicates high similarity in
   transcriptome profile. Cells were color-coded with patient origin. (B)
   UMAP plot of the integrated dataset depicts 12 transcriptional distinct
   cell clusters and cells were color-coded with cluster information. (C)
   UMAP plot of the integrated dataset depicts 6 major cell types and
   cells were color-coded with cell types. (D) Dot plot of top 3 marker
   genes upregulated in each cell type.

3.2. Identifying the subpopulation of tumor cells that are sensitive to
paclitaxel

   To further delineate the transcriptional landscape of OC tumor cells,
   we re-clustered the epithelial cells, identifying seven distinct
   clusters that illustrate pronounced transcriptional heterogeneity among
   OC patients ([70]Fig. 2A–B). Subsequently, we conducted cell cycle and
   epithelial-mesenchymal transition (EMT) score analyses, finding that
   most high-EMT tumor cells were in the G1 phase ([71]Fig. 2C–D).

Fig. 2.

   [72]Fig. 2
   [73]Open in a new tab

   Identification of tumor cells that are sensitive to by Scissor.

   (A) UMAP clustering of epithelial cells. Cells were color-coded with
   cluster information. (B) UMAP clustering of epithelial cells. Cells
   were color-coded with patient origin. (C) UMAP clustering of epithelial
   cells. Cells were color-coded with cell phase. (D) UMAP clustering of
   epithelial cells. Cells were color-coded with EMT score. (E) UMAP plot
   of the integrated dataset depicts Scissor status. Scissor-cells were
   colored in blue. Background cells were colored in grey. (F) Violin
   plots of Scissor score of Scissor-cells and background cells. (G)
   Overall survival analysis of TCGA cohort based on expression of Scissor
   genes by Kaplan-Meier plotter. The patients were grouped into two
   groups based on median Scissor score. (H) Bar plot of the GO enrichment
   analysis for Scissor signature genes. Y-axis represents the enriched GO
   terms; X-axis represents the amount of the Scissor signature genes GO
   terms. GO, Gene Ontology. (I) Bar plot of the KEGG pathway enrichment
   analysis for Scissor signature genes. Y-axis represents pathways;
   X-axis represents the amount of the Scissor signature genes enriched in
   KEGG pathways. KEGG, Kyoto Encyclopedia of Genes and Genomes.

   To identify cells sensitive to paclitaxel treatment, we applied Scissor
   analysis to all tumor cells. The results identified 812
   paclitaxel-sensitive cells (Scissor-) out of 9161 tumor cells ([74]Fig.
   2E). Further investigations revealed 90 significantly upregulated genes
   in the Scissor-cells, termed paclitaxel sensitive genes (PSG) ([75]Fig.
   2F). We then assessed PSG expression in The Cancer Genome Atlas
   (TCGA-OV) samples, categorizing them into two groups based on median
   expression levels ([76]Fig. 2G). Notably, there was a significant
   difference in overall survival (OS) between the high and low Scissor
   groups (p = 0.038, [77]Fig. 2G), suggesting that higher PSG expression
   correlates with improved survival outcomes. To elucidate the biological
   pathways and functions associated with PSG, we conducted Gene Ontology
   (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway
   analyses. The GO analysis showed that PSG were predominantly enriched
   in protein folding and DNA transcription in response to stress
   ([78]Fig. 2H). The KEGG analysis highlighted significant enrichment in
   processes related to protein synthesis in the endoplasmic reticulum,
   lipid metabolism, atherosclerosis, and estrogen signaling pathways
   ([79]Fig. 2I).

3.3. Cell-cell interaction profile of OV with different response to
paclitaxel

   We assessed the PSG score in OC single-cell RNA sequencing samples,
   dividing them into two groups based on their scores: six samples with
   positive scores and nine with negative scores. Cell proportion analysis
   revealed a higher T cell infiltration in the PSG-high group compared to
   the PSG-low group ([80]Fig. 3A). We then conducted a CellChat analysis
   to examine the differences in cell-cell communication between the
   PSG-high and PSG-low groups. This analysis showed a decrease in the
   number of cell-cell interactions but an increase in their strength in
   the PSG-high group relative to the PSG-low group ([81]Fig. 3B). Further
   analysis revealed that the enhanced strength of interactions in the
   PSG-high group was primarily observed between epithelial cells and T
   cells/B cells ([82]Fig. 3C). Additionally, ligand-receptor profiles
   indicated a significant upregulation of MIF/CD74+CXCR4 and
   MIF/CD74 + CD44 interactions between epithelial cells and T cells in
   the PSG-high group ([83]Fig. 3D). These findings suggest that the tumor
   microenvironment (TME) may play a critical role in the response of
   tumor cells to paclitaxel treatment in OC.

Fig. 3.

   [84]Fig. 3
   [85]Open in a new tab

   Analysis of cell-cell communication between the Scissor high and
   Scissor low groups using CellChat

   (A) The proportion of each cell type in the Scissor high and Scissor
   low groups. (B) Total interaction number and strength between cells in
   the Scissor high and Scissor low groups. (C) Overview of selected
   ligand-receptor interactions of tumor cells and other cells in the
   Scissor high and Scissor low groups. The line thickness is proportional
   to the number of ligands when cognate receptors are present in the
   recipient cell type. The loops indicate autocrine circuits. (D) Bubble
   plot showing the selected ligand-receptor interactions between tumor
   cells and other cells in the Scissor high and Scissor low groups. P
   values are indicated by circle size, with the scale to the right
   (permutation test).

3.4. Mutational landscape in TCGA OV sample with different PSG signature

   Given the critical role of genetic mutations in tailoring cancer
   treatment, we analyzed the somatic mutation profiles in TCGA-OV
   samples. Our analysis identified genes with high mutation frequencies,
   including CSMD3, RYR2, and MUC16 ([86]Fig. 4A). Moreover, we observed
   distinct pathway impacts between risk groups: the NRF2 pathway was more
   significantly affected in the high-risk group, while the TP53 pathway
   showed greater alterations in the low-risk group ([87]Fig. 4B–C). We
   identified five genes—ADGRG7, HIVEP1, TECTA, and TMEM131—that exhibited
   significant differences in mutation rates between the high and low-risk
   groups ([88]Fig. 4D). However, there was no observed difference in the
   overall tumor mutation burden (TMB) between these groups ([89]Fig. 4E).

Fig. 4.

   [90]Fig. 4
   [91]Open in a new tab

   Statistics of mutation information in the Scissor high and Scissor low
   TCGA OV samples.

   (A) Statistical results of the top 5 hyper abrupt genes in the scissor
   high samples. (B) Bar plot depicts the distributions of 9 pathways
   affected by mutation in the scissor high samples. (C) Bar plot depicts
   the distributions of 9 pathways affected by mutation in the scissor low
   samples. (D) Bar plot depicts the comparison of 4 driver genes in the
   scissor high and scissor low samples. (E) Comparison of TMB levels
   between the scissor high and scissor low samples. TMB, tumor mutation
   burden. (F) Differences of 22 immune cells infiltration between the
   scissor high and scissor low samples according to CIBERSORT algorithm.

   To further distinguish the immune cell landscape between the high- and
   low-risk groups, we employed CIBERSORT to analyze 22 types of immune
   cells. The results showed that the high-risk group had higher
   proportions of CD4^+ memory resting T cells and activated dendritic
   cells (DCs), whereas the low-risk group was characterized by an
   enrichment of M2 macrophages, resting DCs, and plasma B cells, which
   are typically associated with immunosuppression ([92]Fig. 4F). These
   findings underscore the potential of immune profiles as indicators of
   risk and therapeutic response in ovarian cancer.

3.5. Identification of a five-gene signature

   In our study of the TCGA-OV cohort, all PSGs were evaluated for their
   association with overall survival (OS) using univariate Cox regression
   analysis. This analysis revealed that ten genes were significantly
   correlated with patient prognosis (P < 0.05). Subsequently, a
   paclitaxel-sensitive signature comprising five genes was developed for
   prognosis prediction using the LASSO regression algorithm ([93]Fig. 5A
   and B). The formula for the risk score is as follows: Risk score =
   (−9.842e-07 × expression level of HSP90AA1) + (1.179e-05 × expression
   level of JUN) + (−1.748e-06 × expression level of
   CALM1) + (−1.424e-05 × expression level of
   RBP1) + (−5.665e-05 × expression level of GLRX5). The 342 OC patients
   were categorized into high- and low-risk groups based on the median
   risk score, with the high-risk group demonstrating better survival
   outcomes than the low-risk group ([94]Fig. 5C). The model was also
   validated across different datasets, including [95]GSE14764,
   [96]GSE23554, and [97]GSE26712. Notably, in the [98]GSE23554 and
   [99]GSE26712 datasets, patients with higher signature scores exhibited
   significantly better survival times compared to those with lower scores
   ([100]Fig. 5D–F). Further analysis involved both univariate and
   multivariate Cox regression to evaluate if the risk score could serve
   as an independent prognostic factor relative to other common
   clinicopathological parameters for the TCGA-OV patients. The
   multivariate Cox regression analysis identified both ovarian tumor
   stage and risk score as significant predictors of OS ([101]Fig. 5G).
   Similarly, univariate analysis confirmed the significance of the
   prognostic risk score ([102]Fig. 5H). In the [103]GSE26712 cohort, both
   multivariate and univariate analyses indicated that the risk score was
   a prognostic factor for OS ([104]Fig. 5I and J). The risk graph and
   heatmap were used to detail survival outcomes and display the
   expression differences of the five genes within the model across the
   risk groups for both the TCGA and [105]GSE26712 cohorts ([106]Fig. 6A
   and B). The risk score also showed robust performance in predicting OS
   for individuals in these cohorts ([107]Fig. 6C and D). Based on these
   findings, we developed nomograms that incorporate clinical
   characteristics and the risk score for both the TCGA and [108]GSE26712
   cohorts, as illustrated in [109]Fig. 6E and G. These tools offer a
   quantitative method for predicting survival, enhancing the
   personalization of treatment strategies in ovarian cancer. The
   calibration curve effectively depicts the concordance between the
   predicted and observed survival probabilities at 1-year, 3-year, and
   5-year intervals for ovarian cancer patients. The data points closely
   approach the ideal line, indicating that the prognostic model
   accurately reflects the survival outcomes based on the gene signature.
   This precise calibration suggests that the model captures essential
   biological variables influencing patient survival and highlights its
   potential reliability for clinical application ([110]Fig. 6F and H).
   Taken together, the nomogram may serve as a useful quantitative tool
   for predicting prognoses of OV.

Fig. 5.

   [111]Fig. 5
   [112]Open in a new tab

   LASSO Cox regression and overall survival analysis of Scissor signature
   genes in OV.

   (A) LASSO coefficient profile plots of each independent variable. (B)
   The partial likelihood deviance for the LASSO Cox regression analysis.
   (C) Analysis of the prognostic value of the RiskScore in TCGA-OV
   dataset. (D) Analysis of the prognostic value of the RiskScore in
   [113]GSE14764 dataset. (E) Analysis of the prognostic value of the
   RiskScore in [114]GSE23554 dataset. (F) Analysis of the prognostic
   value of the RiskScore in [115]GSE26712 dataset. (G) Forest plot of the
   multivariate regression analysis in the TCGA cohort. (H) Forest plot of
   the univariate regression analysis in the TCGA cohort. (I) Forest plot
   of the multivariate regression analysis in the [116]GSE26712 cohort.
   (J) Forest plot of the univariate regression analysis in the
   [117]GSE26712 cohort.

Fig. 6.

   [118]Fig. 6
   [119]Open in a new tab

   Validation of a 6 scissor genes-based OS model in OV.

   (A) Risk score and survival status of each patient in the TCGA cohort.
   Heatmap of 6 scissor genes in the TCGA cohort. (B) Risk score and
   survival status of each patient in the [120]GSE26712 cohort. Heatmap of
   6 scissor genes in the [121]GSE26712 cohort. (C) The ROC analysis of
   TCGA-OV dataset for prognosis prediction by riskscore. (D) The ROC
   analysis of [122]GSE26712 dataset for prognosis prediction by
   riskscore. (E) Columnar plots to predict 1-year, 3-year, and 5-year
   overall survival of OV patients in TCGA cohort. (F) Calibration curves
   of the overall survival line plot model in the TCGA group. (G) Columnar
   plots to predict 1-year, 2-year, and 3-year overall survival of OV
   patients in [123]GSE26712 cohort. (H) Calibration curves of the overall
   survival line plot model in the [124]GSE26712 group.

4. Discussion

   Ovarian cancer drug resistance is a common problem in the treatment of
   ovarian cancer, resulting in the cancer continuing to grow and spread,
   making treatment more difficult [[125]28]. There are several factors
   that can contribute to drug response in ovarian cancer, including
   genetic mutations, changes in tumor microenvironment, and the
   activation of certain cellular pathways [[126]29]. Platinum-based
   chemotherapeutics remain the principal agents for treating both newly
   diagnosed ovarian cancer and platinum-sensitive recurrent forms of the
   disease. The integration of poly (ADP-ribose) polymerase inhibitors
   (PARP inhibitors) into the ovarian cancer treatment regimen represents
   a significant advancement, offering substantial benefits to patients
   exhibiting defects in DNA repair mechanisms. Recent findings underscore
   the utility of PARP inhibitors in patients with advanced-stage ovarian
   cancer, extending beyond those harboring BRCA mutations. However, the
   emergence of resistance to PARP inhibitors among certain patient
   cohorts presents a critical challenge, tempering the optimism generated
   by these breakthroughs [[127]30]. To overcome this, effort has been
   made to develop new drugs or combination therapies that target multiple
   pathways or mechanisms of resistance [[128]31]. However, it is also
   important to identify biomarkers that can predict drug response,
   allowing for the selection of more personalized treatment options for
   each patient.

   Previous research has identified multiple biomarkers capable of
   predicting drug response in OC. For instance, patients with mutations
   in the TP53 gene or elevated HE4 protein levels exhibit reduced
   sensitivity to chemotherapy, suggesting a need for alternative
   treatment approaches [[129]32,[130]33]. Conversely, mutations in genes
   associated with the DNA damage response (DDR) pathway, such as ATM,
   ATR, and CHK1, or low ERCC1 expression levels may enhance
   responsiveness to platinum-based chemotherapy [[131]34]. However, it is
   crucial to recognize that biomarkers do not always guarantee accurate
   drug response predictions; thus, treatment strategies should consider a
   variety of factors. In this context, RNA sequencing (RNAseq) has been
   employed to refine drug response prediction accuracy. One notable study
   utilized TCGA online data and OV GEO datasets to develop a predictive
   model based on 10 genes [[132]35]. Our findings indicate that the NRF2
   pathway is notably active in paclitaxel-sensitive patients, whereas
   alterations in the TP53 pathway are more prevalent in other patient
   groups.

   Drug treatments, notably chemotherapy, can significantly impact T cell
   functionality, playing a pivotal role in cultivating an
   immunosuppressive tumor microenvironment. Our study highlights that in
   OV, chemotherapy enhances the interaction between epithelial cells and
   T cells/B cells. Specifically, interactions involving MIF/CD74+CXCR4
   and MIF/CD74 + CD44 were found to be substantially upregulated,
   suggesting their crucial role in the tumor's response to paclitaxel.
   Chemotherapy is known to prompt the secretion of immunosuppressive
   cytokines such as TGF-β and IL-10 from both tumor and stromal cells
   [[133]36]. These cytokines can inhibit T cell activation and
   proliferation, facilitating the proliferation of regulatory T cells
   (Tregs) and myeloid-derived suppressor cells (MDSCs), which further
   diminish T cell functionality [[134]37]. Additionally, chemotherapy can
   alter the expression of major histocompatibility complex (MHC)
   molecules and components of the antigen-processing machinery in tumor
   cells, impairing the recognition and activation of T cells [[135]38].
   It also influences the expression of co-stimulatory molecules, such as
   CD80 and CD86, which are essential for robust T cell activation
   [[136]39]. These effects collectively highlight the complex interplay
   between chemotherapy and the immune landscape within tumors,
   underscoring the need for integrated therapeutic strategies that
   consider both tumor cell killing and the modulation of the immune
   response.

   Our research significantly contributes to advancing ovarian cancer
   prognostics, yet it is essential to recognize the inherent limitations
   that could moderate our conclusions. One major challenge is the
   validation of our prognostic model using publicly available datasets
   such as TCGA and [137]GSE26712. These datasets often lack comprehensive
   clinical parameters, especially detailed treatment histories, which are
   crucial for enhancing the predictive accuracy of the model [[138][40],
   [139][41], [140][42]]. Therefore, robust validation in a broader range
   of prospective cohorts, enriched with extensive longitudinal data on
   treatments and outcomes, is imperative. Such data will not only
   validate the reliability of the model but also expand its clinical
   applicability [[141]43].

   Furthermore, despite the proven robustness of the Scissor algorithm,
   which has shown superior performance in multiple studies, it does not
   capture all sources of biological variability [[142]44,[143]45]. The
   heterogeneity of the tumor microenvironment, a pivotal factor in
   resistance and prognosis, is still not fully represented in the
   datasets we analyzed. This limitation may introduce biases in our
   findings, potentially affecting their applicability across diverse
   patient cohorts or in different clinical settings. This underscores the
   need for ongoing refinement of prognostic models to incorporate a more
   comprehensive view of the biological landscape, ensuring that they
   remain robust and predictive in a real-world clinical context.

   To mitigate the potential limitations identified in our current study,
   it is imperative for future research efforts to focus on integrating
   diverse datasets, particularly emphasizing the inclusion of data from
   multi-institutional and multicentric cohorts. Enriching these datasets
   with a broad range of clinical variables and comprehensive demographic
   data is crucial. These enhancements will not only strengthen the
   external validity of our prognostic models but will also significantly
   boost their applicability in the field of precision oncology. By
   ensuring that these models are robust and adaptable, we can better
   accommodate variations in clinical practices and patient populations
   across different settings. This approach will enable more personalized
   and effective treatment strategies, ultimately improving patient
   outcomes in ovarian cancer care. The strategic integration of varied
   data sources and detailed clinical information will provide a more
   holistic understanding of the factors influencing disease progression
   and treatment response, facilitating the development of universally
   applicable and highly accurate prognostic tools.

   In summary, combining the public scRNAseq dataset and bulk RNAseq
   databases, we selected the feature DEGs to construct a model to
   determine the patients' responses to platinum-based chemotherapy.
   Meanwhile, a prognostic risk model was established to help predict
   patients’ prognoses.

Data availability

   GEO data analyzed during this study can be downloaded from
   [144]https://www.ncbi.nlm.nih.gov/geo website, including
   [145]GSE154600, [146]GSE158937, [147]GSE184880, [148]GSE14764,
   [149]GSE23554, [150]GSE26712 and [151]GSE172016 datasets. TCGA data can
   be downloaded from [152]https://portal.gdc.cancer.gov/database.

CRediT authorship contribution statement

   ZhenWei Zhang: Writing – review & editing, Writing – original draft,
   Visualization, Methodology, Formal analysis, Data curation,
   Conceptualization. MianMian Chen: Writing – review & editing, Software,
   Investigation. XiaoLian Peng: Writing – review & editing,
   Visualization, Validation, Supervision, Methodology, Formal analysis,
   Data curation, Conceptualization.

Declaration of competing interest

   The authors declare that they have no known competing financial
   interests or personal relationships that could have appeared to
   influence the work reported in this paper.

References