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).

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).

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.

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).

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.

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.

Drug IC₅₀ values were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). We calculated the differences of drug IC₅₀ values between the high- and low-risk cases using R’s “pRRophetic” package to guide clinical medication (35).

R software was utilized to perform all statistical analyses in this study, and statistical significance was set at p< 0.05.

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 ).

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⁻¹¹) ( 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.

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.

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.

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 .

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 (IC₅₀) values of drugs from the GDSC website and calculated the differences in IC₅₀ 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.