Gene expression profiles (RNA-sequencing data) of The Cancer Genome Atlas-ovarian cancer (TCGA-OV) cohort were retrieved from TCGA using the University of California, Santa Cruz, Xena browser (USCS Xena) datahub (https://xenabrowser.net/; accessed on October 12, 2020), which is a high-performance visualization and analysis tool for both large public repositories and private datasets [7]. The data included HTSeq-counts and fragments per kilobase of transcript per million values of 379 cases with patient follow-up information. TCGA-OV cohort was used as the training dataset.

Raw microarray data of the GSE18520 and GSE26712 datasets were accessed from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). The GSE18520 dataset comprised of gene expression profiles from 53 high-grade primary papillary serous ovarian adenocarcinoma specimens and 10 normal ovarian surface epithelium brushing specimens, which were evaluated using the Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Inc., Santa Clara, CA,USA); while the GSE26712 dataset comprised of 185 primary ovarian tumour samples and 10 normal ovarian surface epithelium samples, which were detected using the Affymetrix Human U133A microarray (Affymetrix, Inc.). For the validation studies, we used the GSE140082 dataset, which comprised of data on the gene expression arrays from 380 formalin-fixed paraffin-embedded ovarian cancer samples and overall survival information of the patients, which was also obtained from the GEO database and evaluated using Illumina HumanHT-12 WG-DASL V4.0 R2 Expression BeadChip Array (GEO; NCBI, Bethesda, MD, USA). Details of the primary data used in this study are presented in Table 1.

The quality of the raw data obtained from the GSE18520 and GSE26712 datasets was assessed using the “array-QualityMetrics” (version 3.44.0; Bioconductor) package [8]. A robust multi-array average algorithm was used for data normalization and background correction. Different microarray scanning times were considered as a known batch effect, which was examined using the ComBat function of the “sva” (version 3.36.0; Bioconductor) package [9]. Batch effect control was performed by using a principal component analysis through the “FaceFactoMineR” (version 2.4; Bioconductor) package [10]. The gene annotation file (gencode.v35.annotation.gtf..gz, https://www.gencodegenes.org/; accessed on October 12, 2020) of the human genome was used to extract data on protein-coding genes and was downloaded from GENCODE (https://www.gencodegenes.org/; accessed on October 12, 2020). Quality annotations of TCGA-OV samples were performed according to the merged_sample_quality_annotations.tsv file (https://gdc.cancer.gov/aboutdata/publications/pancanatlas). The workflow of the quality control process of TCGA-OV samples is shown in Supplementary Fig. 1. A variance inflation factor was used to detect the severity of multicollinearity between the gene signature and clinicopathological parameters.

Differential expression analysis of the protein-coding genes between the normal and cancerous ovarian samples obtained from the GEO datasets was performed using the “limma” (version 3.44.3; Bioconductor) package [11]. The batch effects of the scanning time in the microarray datasets were integrated into the limma model design to conduct further differential expression analysis. Genes with |logFC| >1 and adjusted P<0.05, were considered differentially expressed.

The value of protein-coding gene expression on survival was determined using a univariate Cox proportional hazards model analysis (log-rank P<0.05) in the TCGA-OV training cohort.

The intersection data between the potential survival-associated genes and differentially expressed genes were filtered using a least absolute shrinkage and selection operator (LASSO)-penalized Cox regression analysis. Only genes with non-zero coefficients via the minimum criteria (a λ value of 0.066 with log [λ]= −2.718 were selected by a 10-fold cross validation) in the LASSO regression model were chosen to further estimate the prognostic significance using univariate and multivariate Cox regression analyses. Next, Kaplan-Meier survival curves were plotted and the log-rank method was used to estimate the prognostic significance. The signature score of a specific gene was calculated by multiplying the gene expression values by the coefficients. A time-dependent receiver operating characteristic (ROC) curve was plotted to investigate the predictive accuracy using the “survivalROC” (version 1.0.3; CRAN) package [12]. Finally, the validation cohort GSE140082 was used to verify the prognostic significance of the gene signature.

xCell [13], a webtool and gene signature-based method was used to perform cell type enrichment analysis based on the gene expression data of 64 immune and stromal cell types. Data on immunological signature gene sets, which represent the cellular states and perturbations within the immune system, were downloaded from the Molecular Signatures Database v7.2 (MSigDB), Broad Institute (https://www.gseamsigdb.org/; accessed on October 12, 2020).

In the TCGA-OV cohort, the following parameters were selected to construct a nomogram using the “survival” (version 3.2-10; CRAN) [14] and “rms” (version 6.1-1; CRAN) packages: age, stage, grade, and residual disease along with the gene signature score [15]. Calibration curves were plotted to assess the concordance between the actual and predicted survival. A time-dependent ROC curve was plotted to investigate the predictive accuracy of the nomogram.

An adjusted P-value <0.05, with |logFC| >1, was considered significant for differential gene expression analysis. The Kaplan-Meier method and log-rank test were used to estimate the impact of the gene signature on survival. Statistical significance was set at P<0.05. A significant Spearman’s correlation coefficient was estimated with a P-value <0.05. All data were processed using the R software (version 4.0.2 [×64]; R-project) platform.

The clinicopathologic characteristics (including age, stage, grade, residual disease, and/or debulking surgery) of the training and validation cohorts are summarized in Table 2.

Normal ovarian samples and ovarian cancer samples from the GSE18520 and GSE26712 datasets were compared according to the adjusted P-value <0.05, and |logFC| >1. The batch effect, which was the scanning time in the two microarray profiles, was integrated into the differential expression analysis model and visualized using principal component analysis (Fig. 1D, E, respectively). In the GSE18520 dataset, 1,226 and 712 protein-coding genes were significantly upregulated and downregulated, while in the GSE26712 dataset, 552 and 765 protein-coding genes were upregulated and downregulated, respectively. Volcano plots of the differentially expressed genes are shown in Fig. 1A. Based on the data of these genes, an intersection of 502 differentially expressed (227 upregulated and 275 downregulated) protein-coding genes was observed, and these genes were selected for further analysis (Fig. 1B, C).

The correlation between the differentially expressed protein-coding genes and survival in the TCGA-OV cohort was analyzed using a univariate Cox regression analysis (P<0.05). A total of 1,354 protein-coding genes were significantly associated with survival in ovarian cancer, and 32 intersecting protein-coding genes from TCGA-OV, GSE18520, and GSE26712 datasets were screened for prognostic value in ovarian cancer (Fig. 1F). Furthermore, data on these 32 genes were filtered using the LASSO-penalised Cox regression analysis (Fig. 1G, H). Moreover, seven genes, including isocitrate dehydrogenase (IDH2), FA complementation group I (FANCI), C-X-C motif chemokine receptor 4 (CXCR4), PRAME nuclear receptor transcriptional regulator (PRAME), praja ring finger ubiquitin ligase 2 (PJA2), enoyl-CoA delta isomerase 2 (ECI2), and AADAC, were identified for conducting further univariate and multivariate Cox regression validation analyses. Finally, univariate and multivariate Cox proportional hazard regression analyses performed using both the training and validation cohorts verified that AADAC was significantly and independently correlated with prognosis in ovarian cancer (Table 3, Fig. 2).

The samples were classified according to the median AADAC expression score then grouped as high -and low-score signature groups. The Kaplan-Meier analysis revealed a significant difference in the outcome of the patients with ovarian cancer between the high- and low-score signature groups (log-rank test P<0.05; Fig. 3A, D). A better survival outcome was significantly correlated with the high-score AADAC expression signature in the subgroup analysis (Fig. 3C, F). The time-dependent area under the ROC curves (AUC) for 1-, 3-, and 5-year survival in the training and validation cohorts are illustrated in Fig. 3B, E, respectively. After combining the AADAC signature score with the clinicopathological parameters (age, stage, grade, residual disease, and/or debulking surgery) in both cohorts, multivariate Cox regression and subgroup analyses were conducted, which further indicated that the high-score AADAC expression signature was significantly and independently correlated with favorable prognosis in ovarian cancer (Fig. 3C, F). The evaluation of the variance inflation factor suggested that the AADAC expression signature showed less collinearity with the other clinical variables in the model (Supplementary Fig. 2). In order to illustrate the reliability of the signature, we sudivided the groups based on the upper and lower quantiles of the AADAC expression signature score into the following: top 25%, bottom 25%, and 25-75% subgroups, and the results still demonstrated a favorable prognosis in the top 25% high-score subgroup (Supplementary Fig. 3).

Immune infiltration in ovarian cancer was estimated using the xCell algorithm. We found that an increased extent of CD4+ memory T cell infiltration was significantly correlated with the upregulation of AADAC (ρ=0.18; P<0.001; 95% confidence interval [CI], 0.22-0.42; Fig 4A, B), and was associated with a better survival outcome (Fig. 4C). Furthermore, the immunological signature enrichment analysis demonstrated that upregulated AADAC expression was significantly enriched in genes that were overexpressed in the T cell receptor pathway and CD4+ T cell regulation (Fig. 4D).

We combined the AADAC expression signature score with feasible clinicopathologic parameters, including age, stage, grade, residual disease, and signature score (Fig. 5A). Each patient was assigned one point for each prognostic parameter: the higher the total number of points, the poorer the outcome. The calibration plots of the nomogram and an ideal model were compared; the nomogram showed moderate performance (Fig. 5B-D). Moreover, the time-dependent AUC of 1-, 3-, and 5-year overall survival demonstrated that our model had a moderate predictive ability (Fig. 5E).