- 1Department of Pharmacy, First Affiliated Hospital of Zhengzhou University, Zhengzhou, Henan, China
- 2Department of Obstetrics and Gynecology, Zhongshan Hospital, Fudan University, Shanghai, China
Background: Ovarian cancer (OC) is a highly lethal and aggressive gynecologic cancer, with an overall survival rate that has shown little improvement over the decades. Robust models are urgently needed to distinguish high-risk cases and predict reliable treatment options for OC. Although anoikis-related genes (ARGs) have been reported to contribute to tumor growth and metastasis, their prognostic value in OC remains unknown. The purpose of this study was to construct an ARG pair (ARGP)-based prognostic signature for patients with OC and elucidate the potential mechanism underlying the involvement of ARGs in OC progression.
Methods: The RNA-sequencing and clinical information data of OC patients were obtained from The Center Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. A novel algorithm based on pairwise comparison was utilized to select ARGPs, followed by the Least Absolute Shrinkage and Selection Operator Cox analysis to construct a prognostic signature. The predictive ability of the model was validated using an external dataset, a receiver operating characteristic curve, and stratification analysis. The immune microenvironment and the proportion of immune cells were analyzed in high- and low-risk OC cases using seven algorithms. Gene set enrichment analysis and weighted gene co-expression network analysis were performed to investigate the potential mechanisms of ARGs in OC occurrence and prognosis.
Results: The 19-ARGP signature was identified as an important prognostic predictor for 1-, 2-, and 3-year overall survival of patients with OC. Gene function enrichment analysis showed that the high-risk group was characterized by the infiltration of immunosuppressive cells and the enrichment of adherence-related signaling pathway, suggesting that ARGs were involved in OC progression by mediating immune escape and tumor metastasis.
Conclusion: We constructed a reliable ARGP prognostic signature of OC, and our findings suggested that ARGs exerted a vital interplay in OC immune microenvironment and therapeutic response. These insights provided valuable information regarding the molecular mechanisms underlying this disease and potential targeted therapies.
1 Introduction
Ovarian cancer (OC) is the most lethal of all female gynecologic cancers worldwide (1). OC is often referred to as the “silent killer” owing to atypical symptoms, such as abdominal bloating and ascites, and a lack of definitive screening tools, leading to its diagnosis in the advanced stages (2). High-grade serous OC accounts for up to 70% of all cases, with many patients experiencing relapse and drug resistance after initial treatment (3). The conventional management of OC includes optimal debulking surgery and platinum-based chemotherapy. With advances in neoadjuvant chemotherapy, cytoreductive surgery, and molecular targeted therapy, the clinical outcome of OC patients has improved significantly (4). However, overall survival (OS) in OC has improved minimally over the last decades. Therefore, developing more effective prognostic models and exploring innovative therapeutic targets for predicting and enhancing outcomes are urgent in OC.
Anoikis is a form of programmed cell death that occurs due to detachment from the extracellular matrix. This regulatory process is essential for inhibiting tumor invasion and metastasis in cancer (5). A critical step in tumor metastasis is the development of anoikis resistance (6). Numerous studies have shown that anoikis plays a vital role in OC. In vitro experiments suggest that LRRC15, NOTCH3, and PCMT1 promote OC cells’ migration and adhesion by hindering the anoikis-induced cell death (7–9). CPT1A and CBX2 are overexpressed in high-grade serous OC and act as drivers of anoikis resistance and tumor dissemination (10, 11). In contrast, angiopoietin-like protein 2, an inflammatory factor, inhibits peritoneal metastasis of OC cells by suppressing anoikis resistance (12). Recently, some prognostic models based on anoikis-related genes (ARGs) using RNA-sequencing data have been developed for predicting the clinical outcomes of several cancers, including hepatocellular carcinoma (13), cutaneous melanoma (14), gastric cancer (15), lung adenocarcinoma (16), glioblastoma (17), and colorectal cancer (18). However, few studies have systematically evaluated the prognostic value of ARGs in OC. Therefore, a more in-depth and detailed analysis is urgent to explore the prognostic value of ARGs in OC.
In the present study, we adopted a novel algorithm that employed pairwise comparison to identify valid anoikis-related gene pairs (ARGPs) and generated a prognostic signature of 19 ARGPs. The 19-ARGP signature was shown to be reliable and valid in predicting the OS of OC patients and was identified as an independent predictive factor in OC. Furthermore, we revealed potential relationships between ARGs and the immune microenvironment, as well as their potential value in therapy.
2 Materials and methods
2.1 Data collection and preprocessing
We obtained the HTSeq-FPKM raw data and relevant clinical information data of OC cases from The Center Genome Atlas (TCGA) database ((https://portal.gdc.cancer.gov/) as the training set. After excluding patients without complete survival data, a total of 365 OC patients with complete follow-up data and their follow-up time greater than 30 days were included. In the process of further verification, we adopted the same inclusion criteria and downloaded the GSE9891 (n = 229) dataset of OC cases from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). All datasets were freely available and publicly accessible from the TCGA database or the GEO database. This study was strictly complied with the data extraction policies of the database, and the ethics committee approval was not required to conduct the current study. The raw data were preprocessed to carry out batch normalization to offset the deviations in datasets using “Combat” function in R’s “sva” package (19). ARGs were obtained from the GeneCards database (https://www.genecards.org/Search/Keyword?queryString=anoikis). A total of 434 ARGs were acquired from the GeneCards database (relevance score>0.4).
2.2 Identification of ARGPs in patients with OC
Univariate Cox regression analysis of datasets was performed using R’s “survival” package for identifying prognostic ARGs, and those ARGs with p< 0.05 were selected (n = 52). We pairwise compared the relative expression levels of prognostic ARGs in each sample and scored each ARG (20). ARGs were paired to constitute ARG pairs (ARGPs), including ARG1 and ARG2. If the expression level of ARG1 was lower than that of ARG2 in a particular ARGP, the score of this pair was recorded as 1; otherwise, the score was recorded as 0. This approach based on pairwise comparison could remove the technical heterogeneity and does not require additional standardization because the score entirely depends on gene expression levels in samples. ARGPs generally received the same score (0 or 1) in more than 80% of the samples, possibly attributable to (1) the dependence on platform priority measurements, which could lead to bias and may not be reproducible across platforms, and (2) the biological priority transcription, which could not provide discriminatory information about patients’ survival. Therefore, we deleted these ARGPs with constant value in each dataset and obtained 431 ARGPs. Finally, univariate Cox regression analysis of ARGP was performed using R’s “survival” package to identify prognostic AGRPs (n = 50).
2.3 Construction, validation, and assessment of an ARGP prognostic signature
The data of OC cases from the TCGA database were utilized as the training cohort. The Least Absolute Shrinkage and Selection Operator (LASSO) penalized Cox regression analysis was performed to further screen and remove collinearity for 50 ARGPs. Finally, a risk score model with β (Coef) value multiplied by the ARGP score was established. Risk score = (β1*ARGP1+β2*ARGP2+β3*ARGP3+…+βn*ARGPn), where xi and β (Coef) are the relative expression level and coefficient of ARGP (21, 22), respectively. The risk score formula was as follows:
Finally, a total of 19 ARGPs were identified for prognostic signature building. A total of 366 patients were divided into high- and low-risk groups with the cutoff thresholds of median risk score. Kaplan–Meier survival analysis and log-rank test were performed to generate survival curves for OS and compare the differences between high- and low-risk cases using R’s “survival” and “survminer“ packages, respectively. The 1-, 2-, and 3-year receiver operating characteristic (ROC) curves and their area under the ROC curves (AUC) were evaluated. The prognostic accuracy of the risk score and other clinical information, such as age, were assessed using univariate and multivariate Cox regression analyses. The predictive power was assessed using a nomogram and time-dependent ROC curves. The discriminative ability of the nomogram was evaluated using calibration curves.
2.4 Weighted gene co-expression network analysis and functional enrichment analysis
Weighted gene co-expression network analysis (WGCNA) could cluster genes with similar expression patterns and find the association between modules and specific phonotypes. The signature has well quantified the risk in patients with OC. To elucidate the reason for the significant heterogeneity between the high- and low-risk groups, we constructed a weighted gene co-expression network using R’s “WGCNA” package with approximately scale-free characteristics (23). An adjacency matrix was determined by these genes’ expression correlation. The gene modules were produced using the method of topological overlap measure (TOM) (24). Co-expression genes modules and clustering graphs were performed using the dynamic tree cutting algorithm (25). Afterwards, the modules with related genes were merged into new modules. Gene significance (GS) and module significance (MS; the mean value of all GS values) were calculated to measure the correlation between genes and modules. We further identified the gene module (blue module) that was significantly co-expressed in the high-risk group and performed functional enrichment analysis for this module. Furthermore, we used the dataset of C2–C8 from the Molecular Signature Database (MSigDB) database for gene set enrichment analysis (GSEA) to comprehensively analyze the activity differences in the pathway between high- and low-risk groups, and verify the results of the previous functional enrichment analysis (26).
2.5 Immune infiltration analysis
We analyzed the difference in immune cell infiltration between high- and low-risk groups, and the association between signature and immuno-infiltration scores using the following R packages: “MCPcounter” (27), “CIBERSORT” (28), “xCell” (29), “TIMER” (30), “EPIC” (31), “CIBERSORT-ABS” (32), and “QUANTISEQ” (33). Stromal scores in malignant tissue were estimated using the “Estimation of STromal and Immune cells in Malignant Tumors using Expression data” (ESTIMATE) algorithm (34). The Wilcox test was utilized to compare the differences of stromal scores between high- and low-risk groups. Fisher’s exact test was utilized to evaluate the correlation between risk scores and stromal scores.
2.6 Mutation profile analysis
We obtained the somatic mutation data of patients with OC from the TCGA database (https://portal.gdc.cancer.gov/). Waterfall graphs were utilized to show the differences in mutation between high- and low-risk groups. R’s “maftools” package was utilized to analyze, annotate, and visualize the mutation profile of this signature. The “plotmafSummary” function was applied to show the mutation landscape in high- and low-risk groups.
2.7 Prediction of drug therapy response
Drug IC50 values were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). We calculated the differences of drug IC50 values between the high- and low-risk cases using R’s “pRRophetic” package to guide clinical medication (35).
2.8 Statistical analysis
R software was utilized to perform all statistical analyses in this study, and statistical significance was set at p< 0.05.
3 Results
3.1 Identification of prognostic ARGs and ARGPs
The flowchart of this study is shown in Figure 1. OC samples from the TCGA database (n = 365, excluding the patients without complete survival data) were analyzed as the training cohort, and the dataset GSE9891 (n = 229) from the GEO database was analyzed as part of the validation cohort. A total of 434 ARGs (relevance score > 0.4) were obtained from the GeneCards database. Univariable Cox regression analysis was performed to identify prognostic ARGs (n = 52, p value <0.05) and prognostic ARGPs (n = 50) using R’s “survival” package (Figure 2A; Supplementary Table 1).
Figure 2 Construction of the ARGP signature in OC patients. (A) The forest plot of 52 candidate prognostic ARGs using univariable Cox regression analysis. ARGs with a hazard ratio > 1 are considered as prognostic risk factors. (B) Coefficient of the 19 ARGPs. (C) Parameter selection by the LASSO algorithm. (D) Trend graph of LASSO coefficients. ARGPs, anoikis-related gene pairs; OC, ovarian cancer; TCGA, The Cancer Genome Atlas; LASSO, Least Absolute Shrinkage and Selection Operator.
3.2 Construction and validation of an ARGP prognostic signature
After screening and removing collinearity via the LASSO penalized Cox regression analysis, a risk score model with the β (Coef) value multiplied by the ARGP score was established. A total of 19 ARGPs had significant associations with OS (p< 0.001) in OC and were identified for prognostic signature construction (Table 1; Figures 2B–D). A total of 366 patients were divided into high- and low-risk groups with the cutoff thresholds of median risk score. We predicted the OS in the training and validation datasets. One-, 2-, and 3-year time-dependent ROC curves were plotted for the training cohort dataset. The AUCs at 1, 2, and 3 years were 0.683, 0.734, and 0.715 in the training cohort, respectively (Figure 3A). The Kaplan–Meier curves showed that OS in the high-risk group was worse than that in the low-risk group (p = 1.6062 × 10−11) (Figure 3B). An independent validation set, GSE9891 (n = 229), was utilized for confirming the consistency of the prognostic value of the novel ARGP signature. The same approaches were performed to calculate risk scores and divide them into high- and low-risk groups with the cutoff thresholds of median risk score. In this validation cohort, the AUCs at 1, 2, and 3 years were 0.745, 0.752, and 0.766, respectively (Figure 3C). The patients in the high-risk group had worse OS than the patients in the low-risk group (Figure 3D). Next, univariable and multivariable Cox regression analyses were performed for evaluating the effects of age, tumor grade, and ARGP risk score on the time-independent ROC curve of the risk score in the training cohort. Univariable Cox regression analysis exhibited a significant association of the 19-ARGPs prognostic signature with OS (hazard ratio [HR] =3.852, 95% confidence interval [CI]: 2.770–5.282; p< 0.001) (Figure 3E). Multivariable Cox regression analysis showed that the 19-ARGPs prognostic signature was a prognostic factor in patients with OC, independently of age and tumor grade (HR = 3.605, 95% CI: 2.606–4.987; p< 0.001) (Figure 3F). Then, we established a nomogram including age and risk scores to evaluate OC prognosis (Figure 3G). ROC curves (Figure 3H) and calibration curves (Figure 3I) of the model proved that the nomogram had a great predictive performance.
Figure 3 Validation of the ARGP model. OC cases were classified into high- and low-risk groups by median risk score. (A, B) Time-dependent ROC curves for predicting OS at 1, 2, and 3 years in the training and validation cohort. (C, D) Kaplan–Meier curves of OS in the training and validation cohort. (E) The forest plot of univariate Cox regression analysis for the prognosis of OC patients. (F) The forest plot of multivariate Cox regression analysis for the prognosis of OC patients. (G) The nomogram was constructed based on age and risk scores. (H) ROC curves of the nomogram predicting OS at 1, 2, and 3 years. (I) Calibration plots of OS at 1, 2, and 3 years. OS, overall survival; ROC, receiver operating characteristic. **p<0.01; ***p<0.001.
3.3 WGCNA and functional enrichment analysis
R’s “pickSoftThreshold” package was utilized to calculate the soft thresholding power β (set at 3), in which the scale independence reached 0.9 (Figure 4A) and had a relatively high-average connectivity (Figure 4B). We constructed a weighted gene co-expression network using R’s “WGCNA” package with approximately scale-free characteristics. Co-expression gene modules and clustering graphs were performed using the dynamic tree cutting algorithm. Five gene co-expression modules were finally constructed (Figure 4C). We used the method of TOM and mapped the relationships between the identified gene modules (Figure 4D). The light color represented a low overlap, and the dark red color represented a high overlap. These results revealed that the gene expression level was relatively independent between modules. Afterwards, the correlated modules were merged into new modules. The results showed that five modules could be clustered into blue, turquoise, brown, yellow, and gray modules, respectively (Figure 4E). We further calculated the GS and MS, and identified the gene module blue, which was mostly significantly correlated with high-risk cases (p< 0.001). Then, we performed GO and KEGG analyses in the blue module. The GO analysis showed that (Figure 5A), regarding the biological process, the genes were mainly enriched in extracellular structure and matrix organization, as well as external encapsulating organization. As for the cellular component, the genes were mainly enriched in extracellular matrix structural constituent, collagen binding, and extracellular matrix binding. As for molecular function, the genes were mainly enriched in basement membrane and collagen-containing extracellular matrix. The KEGG analysis showed that the blue module pathways include the PI3K-Akt signaling pathway, human papillomavirus infection, focal adhesion, and ECM–receptor interaction (Figure 5B). Finally, the dataset of C2–C8 from the MSigDB database was used for GSEA to comprehensively analyze the differences in signaling pathway activation between high- and low-risk groups, and confirm the results of the previous functional enrichment analysis. The GSEA revealed that the genes were mainly enriched in pathways in cancer, ECM–receptor interaction, focal adhesion, extracellular structure organization, immune response, and collagen-containing extracellular matrix (Figure 5C). Consequently, the GSEA results were almost in agreement with the conclusions obtained from GO and KEGG analyses.
Figure 4 The co-expression network established using WGCNA. (A) The x-axis suggests the soft-thresholding power value. The y-axis represents the scale-free fit index. (B) The y-axis represents the mean connectivity. (C) Clustering dendrogram of genes on the basis of a dissimilarity measure. Different colors reflect gene modules. (D) The heatmap shows the visualization of WGCNA network. (E) Module–trait relationships. WGCNA, weighted gene co-expression network analysis.
Figure 5 Functional enrichment analysis of 19-ARGPs. (A) Significant enriched GO analysis. (B) Significantly enriched KEGG pathways. (C) The GSEA reveals these most significant enrichment pathways. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis.
3.4 Assessment of immune infiltration
A major direction is to investigate the infiltration of immune cells in the immune microenvironment of OC. We compared the differences in immune cell infiltration between high- and low-risk groups, and the association between signature and immuno-infiltration scores using the following R packages: “MCPcounter”, “CIBERSORT”, “xCell”, “TIMER”, “EPIC”, “CIBERSORT-ABS”, and “QUANTISEQ”. Stromal scores were calculated by the “ESTIMATE” algorithm. Fourteen immune cell types showed significantly different proportions between the high- and low-risk groups (p< 0.05) (Figure 6A). In particular, M2 macrophages, neutrophils, endothelial cells, and cancer-associated fibroblasts (CAFs) had elevated infiltration levels in high-risk cases compared with those in low-risk cases, whereas M1 macrophages, CD4+ memory T cells, CD4+ Th1 T cells, CD8+ T cells, CD8+ central memory T cells, class-switched memory B cells, common lymphoid progenitors, hematopoietic stem cells, and plasmacytoid dendritic cells had decreased proportions (all p< 0.01). Among them, the infiltration of M2 macrophages, neutrophils, endothelial cells, and CAFs was positively correlated with the risk score of OC patients (all p< 0.001). The infiltration of M1 macrophages, CD4+ memory T cells, CD4+ Th1 T cells, CD8+ central memory T cells, class-switched memory B cells, common lymphoid progenitors, hematopoietic stem cells, and plasmacytoid dendritic cells was negatively correlated with the risk score of OC patients (all p< 0.001) (Figure 6B). Furthermore, the Wilcox test showed that stromal scores in high-risk cases were higher than those in low-risk cases (p< 0.001) (Figure 6C). Fisher’s exact test showed that stromal scores were positively correlated with risk scores (p< 0.001) (Figure 6D). In high-risk patients with OC, high stromal scores reflected the presence of mesenchymal cells that had been reported previously (36). These results suggested that OC patients in the high-risk group could be more prone to immune escape due to the existence of immunosuppressive microenvironment, resulting in tumor recurrence and metastasis.
Figure 6 Assessment of immune cell infiltration between the two risk groups. (A) The heatmap shows that differences of immune cell infiltration between high- and low-risk groups using seven algorithms. (B) The correlation analysis of immune cells using bubble chart. (C) The differences of stromal scores between high- and low-risk cases. (D) The correlation analysis between stromal scores and risk scores.
3.5 Identification of mutation profiles
We used R’s “maftools” package to analyze, annotate, and visualize the mutation profiles of this signature. The “plotmafSummary” function was performed to show the somatic mutation profiles in high- and low-risk groups of OC patients. In the high-risk group, a waterfall plot revealed that approximately 116 of 117 (99.15%) samples showed somatic mutations (Figure 7A). In the low-risk group, approximately 123 of 127 (96.85%) samples showed somatic mutations (Figure 8A). Among them, the top 10 mutated genes in the high- and low-risk groups are shown in Figures 7G, 8G, respectively. TP53 was the most frequent mutated gene in high-risk (85%) and low-risk (87%) cases, followed by TTN (27%) in high-risk (27%) and low-risk (26%) cases. Other mutated genes were different in these two groups. In the high-risk group, CSMD3 (17%) and USH2A (10%) were common mutated genes. In contrast, MUC16 (9%) and RYR2 (9%) were common mutated genes in the low-risk group. The overall levels of tumor mutation burden (TMB) in high-risk cases (192 mutations) were lower than those in low-risk cases (248 mutations). Next, gene mutations were classified into eight types. Missense mutation was the most frequent variation type in both groups, followed by nonsense mutation, deletion, and insertion (Figures 7B, 8B). Single-nucleotide polymorphism mutated more frequently than insertion and deletion, and C>T was the most frequent single-nucleotide variant class (Figures 7D, E, 8D, E). Furthermore, variant number in each sample was calculated in high-risk (median number: 53) and low-risk groups (median number: 58) (Figures 7C, 8C). Boxplot revealed the variant classification with different colors (Figures 7F, 8F). The differences of co-occurrence and mutually exclusive expression in mutation profiles between the two groups are shown in Figures 7H, 8H.
Figure 7 Mutation profiles in the high-risk group of OC patients. (A) Waterfall plot reflects the top 30 frequent mutation genes. (B) Eight common variant classifications in all samples. (C) Variations per sample. (D) Three variant types in all samples. (E) Six classes of SNV. (F) Summary of variant classification. (G) Top 10 mutated genes. (H) The correlation analysis of mutated genes. *p<0.001; ·p<0.05.
Figure 8 Mutation profiles in the low-risk group of OC patients. (A) Waterfall plot reflects the top 30 frequent mutation genes. (B) Eight common variant classifications in all samples. (C) Variations per sample. (D) Three variant types in all samples. (E) Six classes of SNV. (F) Summary of variant classification. (G) Top 10 mutated genes. (H) The correlation analysis of mutated genes. *p<0.001; ·p<0.05.
3.6 Prediction of drug therapy response
Targeted therapy plays a critical role in many aspects of OC management, and the relationship between risk score and targeted therapeutic response was analyzed. We downloaded the half-maximal inhibitory concentration (IC50) values of drugs from the GDSC website and calculated the differences in IC50 values between the high- and low-risk groups using R’s “pRRophetic” package. High-risk cases were more sensitive to AP.24534 (ponatinib), AZD.0530 (saracatinib), Bexarotene, Embelin, Dasatinib, Imatinib, Midostaurin, Pazopanib, and Pyrimethamine (p< 0.001) (Figures 9A–I, 10A–I). These results confirmed the validity of the ARGP signature in predicting targeted therapy response. It is a valuable guidance for the medication of OC patients.
Figure 9 Prediction of drug therapy response. Boxplot of IC50 values for drugs including AP.24534 (A), AZD.0530 (B), Bexarotene (C), Embelin (D), Dasatinib (E), Imatinib (F), Midostaurin (G), Pazopanib (H), and Pyrimethamine (I) in high- and low-risk groups.
Figure 10 The correlation analysis between the IC50 value of drugs and the risk score. These drugs include AP.24534 (A), AZD.0530 (B), Bexarotene (C), Embelin (D), Dasatinib (E), Imatinib (F), Midostaurin (G), Pazopanib (H), and Pyrimethamine (I).
4 Discussion
OC is known for its insidious onset and poor prognosis. Even with decades of advancements in OC treatment, the 5-year OS rate has not shown significant improvement (37). Anoikis, a type of programmed cell death, plays an important role in tumor progression and metastasis. ARGs have been confirmed to be involved in the occurrence and development of OC. However, the prognostic values of ARGs in OC have not been elucidated. Considering the important role of ARGs in OC, we constructed an ARGPs-based model to predict the prognosis of OC cases. In this study, we identified 52 prognostic ARGs and created 431 ARGPs. Among them, 50 ARGPs with prognostic values were found (Supplementary Table 1). We found that BCL2L11|PAK1, BCL2L11|BAG1, BCL2L11|CASP2, EGFR|CD44, EGFR|CEACAM1, ITGA5|CD44, ITGA5|PIN1, FN1|IFI27, CXCL12|CD4, CXCL12|XIAP, CXCL12|RANBP9, SIK2|CASP2, CRYAB|BAG1, CRYAB|ENDOG, PDGFRB|PIN1, PDGFRB|SERPINB1, PAK4|SERPINB1, LRP1|BAG1, LRP1|APOBEC3G, LRP1|CASP2, LRP1|ENDOG, and ITGB5|SERPINB1 were prognostic risk factors (HR > 1). Other gene pairs were prognostic protection factors in OC. In particular, BCL2L11, CD44, and PAK1 were reported as the independent prognosis markers of OC (38–40). Expect for the three genes, other genes were identified for the first time as prognostic factors in OC. Genes EGFR, XIAP, PIN1, and PDGFRB play important roles in the occurrence and progression of OC and may act as the therapeutic targets in OC (41–43). The expression levels of ITGA5, CRYAB, IFI27, and CEACAM6 genes were significantly increased in high-grade OC, which promoted OC metastasis (44–46). IFI27 may be related to platinum resistance of OC (47). Endonuclease G (ENDOG) and BAG1 promote OC cell proliferation (48) by regulating EGF signaling pathway and enhancing anoikis resistance (49). Furthermore, the roles of LRP1, RANBP9, SERPINB1, and CASP2 genes in OC have not been reported in previous studies, which need to be explored.
We pairwise compared the expression levels of prognostic ARGs in each sample and innovatively calculated the ARGP score. A total of 19 ARGPs were used for the construction of the prognostic model. Among them, CASP8|SIK2, PAK1|PTPN1, CCR7|SFRP1, RANBP9|ITGB5, ELK1|BIN1, ELK1|CCDC80, and BAG1|CCDC80 are prognostic protection factors (HR< 1). The downregulated expression of caspase-8 in OC may be linked to high aggressiveness with chronic inflammation and immune resistance (50). As for chemotherapy, SIK2 decreases sensitivity of OC cells to paclitaxel and promotes migration of OC cells (51, 52). PTPN1 is overexpressed in high-grade serous OC, and may act as a marker of better chemotherapy response (53). Furthermore, SIK2 and CCR7 both participate in epithelial–mesenchymal transition and promote OC cell migration and invasion (54), whereas SFRP1 and BIN1 both inhibit the epithelial OC through inhibiting Wnt/β-catenin signaling (55, 56). Steroid hormones also regulate ARGs and are involved in OC occurrence. For example, follicle-stimulating hormone can stimulate phosphorylated Elk-1 in OC cells (57). High expression levels of APOBEC3G and ITGB5 may be correlated with improved outcomes in high-grade OC patients (58, 59). These ARGs were found for the first time as prognostic protection factors in OC. More in-depth clinical research is necessary to explore their prognostic protection function.
In this study, we found that OC patients in the high-risk group had worse OS than the low-risk group. External validation dataset and internal hierarchical verification revealed the great prediction power of this signature, providing a reliable and effective tool for clinicians to evaluate the survival of OC patients. Furthermore, we explored the potential function of a 19-ARGP signature in OC using immune cell infiltration analysis and WGCNA combined with functional enriched analysis. High proportions of immunosuppressive cells, such as M2 macrophages, and CAFs were found in the high-risk cases. M1 macrophages, CD4+ memory T cells, CD4+ Th1 T cells, CD8+ T cells, CD8+ central memory T cells, class-switched memory B cells, common lymphoid progenitors, hematopoietic stem cells, and plasmacytoid dendritic cells had high proportions in low-risk cases. These results suggest that immunosuppressive microenvironment probably promotes OC progression and metastasis by mediating immune scape. CAFs could recruit ITGA5high ascitic tumor cells to form metastatic units, which further sustain ascitic OC cells’ ITGA5 expression by EGF secretion (44). CAFs also induce epithelial–mesenchymal transition in OC by CXCL12/CXCR4 and PAK1/β-catenin signaling pathways (60). The regulation of ARGs in the immune microenvironment needs further investigation.
At present, many clinical trials focus on targeted therapeutics. We further explored the drug response of patients with OC in two groups. Except for the nine drugs mentioned above, other compounds were also identified, for which the high-risk group showed higher sensitivity than the low-risk group (Supplementary Figures 1, 2). Somatic mutations are considered as the key of immunotherapy in tumors (61). In this study, the proportion of somatic mutations in high-risk cases (99.15%) were higher than that in low-risk cases (96.85%). High proportions of somatic mutations will lead to an increase in neoantigen production, thus improving the immune therapy response. These results require clinical trials to further elucidate and verify.
There are some limitations in our study. This is a retrospective study with data from the TCGA and GEO databases, lacking some clinical information, such as therapeutic and prognostic data. Our results require further investigation about patients with OC in clinical settings. At the same time, a portion of the prognostic genes related to the anoikis that we have obtained still lack sufficient in vitro and in vivo experiment verification in OC.
5 Conclusion
In conclusion, we established a novel and reliable 19-ARGP prognostic signature for predicting the OS of OC patients and analyzed the association between the immune microenvironment and OC. Furthermore, this study indicated that the risk score was relevant to immune cells’ infiltration and somatic mutation. Finally, the ARGP signature may help identify patients with OC suitable for targeted therapy and act as a promising predictive factor to offer insights into therapeutic strategies.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Author contributions
YD and XX performed the study and analyzed the data. XX wrote the manuscript. YD participated in the study design and helped revised the manuscript. All authors contributed to the article and approved the submitted version.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2023.1193622/full#supplementary-material
Supplementary Table 1 | Identification of 50 prognostic ARGPs using the R’s “survival” package.
Supplementary Figure 1 | Differences in drug therapy response between high- and low-risk groups.
Supplementary Figure 2 | The correlation analysis between IC50 value of compounds and the risk score.
References
1. Slatnik CL, Duff E. Ovarian cancer: ensuring early diagnosis. Nurse Pract (2015) 40(9):47–54. doi: 10.1097/01.NPR.0000450742.00077.a2
2. Jayson GC, Kohn EC, Kitchener HC, Ledermann JA. Ovarian cancer. Lancet (2014) 384(9951):1376–88. doi: 10.1016/S0140-6736(13)62146-7
3. Bowtell DD, Böhm S, Ahmed AA, Aspuria PJ, Bast RC Jr., Beral V, et al. Rethinking ovarian cancer II: reducing mortality from high-grade serous ovarian cancer. Nat Rev Cancer (2015) 15(11):668–79. doi: 10.1038/nrc4019
4. Liu J, Matulonis UA. New strategies in ovarian cancer: translating the molecular complexity of ovarian cancer into treatment advances. Clin Cancer Res (2014) 20(20):5150–6. doi: 10.1158/1078-0432.CCR-14-1312
5. Adeshakin FO, Adeshakin AO, Afolabi LO, Yan D, Zhang G, Wan X. Mechanisms for modulating anoikis resistance in cancer and the relevance of metabolic reprogramming. Front Oncol (2021) 11:626577. doi: 10.3389/fonc.2021.626577
6. Kim YN, Koo KH, Sung JY, Yun UJ, Kim H. Anoikis resistance: an essential prerequisite for tumor metastasis. Int J Cell Biol (2012) 2012:306879. doi: 10.1155/2012/306879
7. Ray U, Jung DB, Jin L, Xiao Y, Dasari S, Sarkar Bhattacharya S, et al. Targeting LRRC15 inhibits metastatic dissemination of ovarian cancer. Cancer Res (2022) 82(6):1038–54. doi: 10.1158/0008-5472.CAN-21-0622
8. Zhang J, Li Y, Liu H, Zhang J, Wang J, Xia J, et al. Genome-wide CRISPR/Cas9 library screen identifies PCMT1 as a critical driver of ovarian cancer metastasis. J Exp Clin Cancer Res (2022) 41(1):24. doi: 10.1186/s13046-022-02242-3
9. Brown CW, Brodsky AS, Freiman RN. Notch3 overexpression promotes anoikis resistance in epithelial ovarian cancer via upregulation of COL4A2. Mol Cancer Res (2015) 13(1):78–85. doi: 10.1158/1541-7786.MCR-14-0334
10. Sawyer BT, Qamar L, Yamamoto TM, McMellen A, Watson ZL, Richer JK, et al. Targeting fatty acid oxidation to promote anoikis and inhibit ovarian cancer progression. Mol Cancer Res (2020) 18(7):1088–98. doi: 10.1158/1541-7786.MCR-19-1057
11. Wheeler LJ, Watson ZL, Qamar L, Yamamoto TM, Post MD, Berning AA, et al. CBX2 identified as driver of anoikis escape and dissemination in high grade serous ovarian cancer. Oncogenesis (2018) 7(11):92. doi: 10.1038/s41389-018-0103-1
12. Takeshita Y, Motohara T, Kadomatsu T, Doi T, Obayashi K, Oike Y, et al. Angiopoietin-like protein 2 decreases peritoneal metastasis of ovarian cancer cells by suppressing anoikis resistance. Biochem Biophys Res Commun (2021) 561:26–32. doi: 10.1016/j.bbrc.2021.05.008
13. Chen Y, Huang W, Ouyang J, Wang J, Xie Z. Identification of anoikis-related subgroups and prognosis model in liver hepatocellular carcinoma. Int J Mol Sci (2023) 24(3):2862. doi: 10.3390/ijms24032862
14. Zhou Y, Wang C, Chen Y, Zhang W, Fu Z, Li J, et al. A novel risk model based on anoikis: predicting prognosis and immune infiltration in cutaneous melanoma. Front Pharmacol (2022) 13:1090857. doi: 10.3389/fphar.2022.1090857
15. Zhao Z, Li C, Peng Y, Liu R, Li Q. Construction of an original anoikis-related prognostic model closely related to immune infiltration in gastric cancer. Front Genet (2022) 13:1087201. doi: 10.3389/fgene.2022.1087201
16. Diao X, Guo C, Li S. Identification of a novel anoikis-related gene signature to predict prognosis and tumor microenvironment in lung adenocarcinoma. Thorac Cancer (2023) 14(3):320–30. doi: 10.1111/1759-7714.14766
17. Sun Z, Zhao Y, Wei Y, Ding X, Tan C, Wang C. Identification and validation of an anoikis-associated gene signature to predict clinical character, stemness, IDH mutation, and immune filtration in glioblastoma. Front Immunol (2022) 13:939523. doi: 10.3389/fimmu.2022.939523
18. Cai Z, Zhou F. A novel anoikis and immune-related genes marked prognostic signature for colorectal cancer. Med (Baltimore) (2022) 101(46):e31127. doi: 10.1097/MD.0000000000031127
19. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical bayes methods. Biostatistics (2007) 8(1):118–27. doi: 10.1093/biostatistics/kxj037
20. Li X, Wang R, Wang S, Wang L, Yu J. Construction of a b cell-related gene pairs signature for predicting prognosis and immunotherapeutic response in non-small cell lung cancer. Front Immunol (2022) 13:989968. doi: 10.3389/fimmu.2022.989968
21. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01
22. Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for cox's proportional hazards model via coordinate descent. J Stat Softw (2011) 39(5):1–13. doi: 10.18637/jss.v039.i05
23. Langfelder P, Horvath S. Fast r functions for robust correlations and hierarchical clustering. J Stat Softw (2012) 46(11):i11. doi: 10.18637/jss.v046.i11
24. Li A, Horvath S. Network module detection: affinity search technique with the multi-node topological overlap measure. BMC Res Notes (2009) 2:142. doi: 10.1186/1756-0500-2-142
25. Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for r. Bioinformatics (2008) 24(5):719–20. doi: 10.1093/bioinformatics/btm563
26. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U.S.A. (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102
27. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol (2016) 17(1):218. doi: 10.1186/s13059-016-1113-y
28. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337
29. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18(1):220. doi: 10.1186/s13059-017-1349-1
30. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res (2017) 77(21):e108–e10. doi: 10.1158/0008-5472.CAN-17-0307
31. van Veldhoven CM, Khan AE, Teucher B, Rohrmann S, Raaschou-Nielsen O, Tjønneland A, et al. Physical activity and lymphoid neoplasms in the European prospective investigation into cancer and nutrition (EPIC). Eur J Cancer (2011) 47(5):748–60. doi: 10.1016/j.ejca.2010.11.010
32. Tamminga M, Hiltermann TJN, Schuuring E, Timens W, Fehrmann RS, Groen HJ. Immune microenvironment composition in non-small cell lung cancer and its association with survival. Clin Transl Immunol (2020) 9(6):e1142. doi: 10.1002/cti2.1142
33. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med (2019) 11(1):34.
34. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
35. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468
36. Zhang Q, Wang C, Cliby WA. Cancer-associated stroma significantly contributes to the mesenchymal subtype signature of serous ovarian cancer. Gynecol Oncol (2019) 152(2):368–74. doi: 10.1016/j.ygyno.2018.11.014
37. Morand S, Devanaboyina M, Staats H, Stanbery L, Nemunaitis J. Ovarian cancer immunotherapy and personalized medicine. Int J Mol Sci (2021) 22(12):6532. doi: 10.3390/ijms22126532
38. Cai J, Qiu J, Wang H, Sun J, Ji Y. Identification of potential biomarkers in ovarian carcinoma and an evaluation of their prognostic value. Ann Transl Med (2021) 9(18):1472. doi: 10.21037/atm-21-4606
39. Bartakova A, Michalova K, Presl J, Vlasak P, Kostun J, Bouda J. CD44 as a cancer stem cell marker and its prognostic value in patients with ovarian carcinoma. J Obstet Gynaecol (2018) 38(1):110–4. doi: 10.1080/01443615.2017.1336753
40. Siu MK, Wong ES, Chan HY, Kong DS, Woo NW, Tam KF, et al. Differential expression and phosphorylation of Pak1 and Pak2 in ovarian cancer: effects on prognosis and cell invasion. Int J Cancer (2010) 127(1):21–31. doi: 10.1002/ijc.25005
41. Cremona M, Vandenberg CJ, Farrelly AM, Madden SF, Morgan C, Kalachand R, et al. BRCA mutations lead to XIAP overexpression and sensitise ovarian cancer to inhibitor of apoptosis (IAP) family inhibitors. Br J Cancer (2022) 127(3):488–99. doi: 10.1038/s41416-022-01823-5
42. Li X, Wang H, Ding J, Nie S, Wang L, Zhang L, et al. Celastrol strongly inhibits proliferation, migration and cancer stem cell properties through suppression of Pin1 in ovarian cancer cells. Eur J Pharmacol (2019) 842:146–56. doi: 10.1016/j.ejphar.2018.10.043
43. Sun T, Bi F, Liu Z, Yang Q. TMEM119 facilitates ovarian cancer cell proliferation, invasion, and migration via the PDGFRB/PI3K/AKT signaling pathway. J Transl Med (2021) 19(1):111. doi: 10.1186/s12967-021-02781-x
44. Gao Q, Yang Z, Xu S, Li X, Yang X, Jin P, et al. Heterotypic CAF-tumor spheroids promote early peritoneal metastatis of ovarian cancer. J Exp Med (2019) 216(3):688–703. doi: 10.1084/jem.20180765
45. du Manoir S, Delpech H, Orsetti B, Jacot W, Pirot N, Noel J, et al. In high-grade ovarian carcinoma, platinum-sensitive tumor recurrence and acquired-resistance derive from quiescent residual cancer cells that overexpress CRYAB, CEACAM6, and SOX2. J Pathol (2022) 257(3):367–78. doi: 10.1002/path.5896
46. Kim YS, Hwan JD, Bae S, Bae DH, Shick WA. Identification of differentially expressed genes using an annealing control primer system in stage III serous ovarian carcinoma. BMC Cancer (2010) 10:576. doi: 10.1186/1471-2407-10-576
47. Guo K, Li L. Prediction of key candidate genes for platinum resistance in ovarian cancer. Int J Gen Med (2021) 14:8237–48. doi: 10.2147/IJGM.S338044
48. Choi YN, Seo TW, Lee YT, Jeong DH, Yoo SJ. Nuclear endonuclease G controls cell proliferation in ovarian cancer. FEBS Open Bio (2023) 13(4):655–69. doi: 10.1002/2211-5463.13572
49. Pagé V, Côté M, Rancourt C, Sullivan M, Piché A. BAG-1 p29 protein prevents drug-induced cell death in the presence of EGF and enhances resistance to anoikis in SKOV3 human ovarian cancer cells. Biochem Biophys Res Commun (2005) 328(4):874–84. doi: 10.1016/j.bbrc.2004.12.193
50. Kostova I, Mandal R, Becker S, Strebhardt K. The role of caspase-8 in the tumor microenvironment of ovarian cancer. Cancer Metastasis Rev (2021) 40(1):303–18. doi: 10.1007/s10555-020-09935-1
51. Zhao J, Zhang X, Gao T, Wang S, Hou Y, Yuan P, et al. SIK2 enhances synthesis of fatty acid and cholesterol in ovarian cancer cells and tumor growth through PI3K/Akt signaling pathway. Cell Death Dis (2020) 11(1):25. doi: 10.1038/s41419-019-2221-x
52. Lu Z, Mao W, Yang H, Santiago-O'Farrill JM, Rask PJ, Mondal J, et al. SIK2 inhibition enhances PARP inhibitor activity synergistically in ovarian and triple-negative breast cancers. J Clin Invest (2022) 132(11):e146471. doi: 10.1172/JCI146471
53. Davidson B, Bock AJ, Holth A, Nymoen DA. The phosphatase PTPN1/PTP1B is a candidate marker of better chemotherapy response in metastatic high-grade serous carcinoma. Cytopathology (2021) 32(2):161–8. doi: 10.1111/cyt.12921
54. Cheng S, Guo J, Yang Q, Yang X. Crk-like adapter protein regulates CCL19/CCR7-mediated epithelial-to-mesenchymal transition via ERK signaling pathway in epithelial ovarian carcinomas. Med Oncol (2015) 32(3):47. doi: 10.1007/s12032-015-0494-1
55. Zhang H, Sun D, Qiu J, Yao L. SFRP1 inhibited the epithelial ovarian cancer through inhibiting wnt/β-catenin signaling. Acta Biochim Pol (2019) 66(4):393–400. doi: 10.18388/abp.2019_2757
56. Lv W, Jia Y, Wang J, Duan Y, Wang X, Liu T, et al. Long non-coding RNA SNHG10 upregulates BIN1 to suppress the tumorigenesis and epithelial-mesenchymal transition of epithelial ovarian cancer via sponging miR-200a-3p. Cell Death Discovery (2022) 8(1):60. doi: 10.1038/s41420-022-00825-9
57. Choi KC, Kang SK, Tai CJ, Auersperg N, Leung PC. Follicle-stimulating hormone activates mitogen-activated protein kinase in preneoplastic and neoplastic ovarian surface epithelial cells. J Clin Endocrinol Metab (2002) 87(5):2245–53. doi: 10.1210/jcem.87.5.8506
58. Leonard B, Starrett GJ, Maurer MJ, Oberg AL, Van Bockstal M, Van Dorpe J, et al. APOBEC3G expression correlates with T-cell infiltration and improved clinical outcomes in high-grade serous ovarian carcinoma. Clin Cancer Res (2016) 22(18):4746–55. doi: 10.1158/1078-0432.CCR-15-2910
59. Zhu T, Chen R, Wang J, Yue H, Lu X, Li J. The prognostic value of ITGA and ITGB superfamily members in patients with high grade serous ovarian cancer. Cancer Cell Int (2020) 20:257. doi: 10.1186/s12935-020-01344-2
60. Zhang F, Cui JY, Gao HF, Yu H, Gao FF, Chen JL, et al. Cancer-associated fibroblasts induce epithelial-mesenchymal transition and cisplatin resistance in ovarian cancer via CXCL12/CXCR4 axis. Future Oncol (2020) 16(32):2619–33. doi: 10.2217/fon-2020-0095
Keywords: anoikis-related gene pairs, ovarian cancer, prognostic signature, immunoinfiltration, therapeutic response
Citation: Duan Y and Xu X (2023) A signature based on anoikis-related genes for the evaluation of prognosis, immunoinfiltration, mutation, and therapeutic response in ovarian cancer. Front. Endocrinol. 14:1193622. doi: 10.3389/fendo.2023.1193622
Received: 25 March 2023; Accepted: 23 May 2023;
Published: 13 June 2023.
Edited by:
Zhenyu Song, Shanghai Jiao Tong University, ChinaCopyright © 2023 Duan and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yiqi Duan, a2FyZW5keXFAMTYzLmNvbQ==
†These authors share first authorship