Twenty-two ADO metabolism related genes (18 ADO metabolism enzymes and 4 ADO receptor) were collected from a previous study (8) and the description of their biological function is shown in Supplementary Table 1 . A total of 372 ovarian cancer patients from TCGA was used as training set and 1,137 patients from 6 GEO datasets were as validation set. The RNAseq, whole-exon sequencing (WES), and corresponding clinical data of TCGA were downloaded from portal.gdc.cancer.gov, and the six GEO validation datasets (GSE14764, GSE23554, GSE26712, GSE32062, GSE49997, GSE140082) were downloaded from www.ncbi.nlm.nih.gov/gds. The basic information of these cohorts is summarized in Supplementary Table 2 . The gene expression and drug response inhibitory concentration (IC50) values for 198 drugs in ovarian cancer cell lines were obtained from Genomics of Drug Sensitivity in Cancer (GDSC, www.cancerrxgene.org) (14).

We used the non-negative matrix factorization (NMF) algorithm to obtain adenosine related Meta genes based on the expression data from TCGA. In detail, the R package “NMF” (15) with the default methods ‘brunet’ was applied. The rank was set as 2 to 10 and the number of runs for each rank set was 200. The cophenetic correlation coefficient was used to determine optimal clustering number. The factorization rank was selected once cophenetic correlation coefficient started to decrease (16). After that, we extracted the top 15 genes for each basis from the results of NMF which were used for the downstream analysis. A total of 20 genes was included in the signature construction. Two gene clusters were identified by consensus clustering algorithm with a correlation matrix of genes expression (17). Further, we calculated the enrichment score by ssGSVA method implemented in R package ‘GSVA’ (18) with the above two gene clusters respectively and defined the ADO score as the difference between the enrichment score of two gene clusters.

Patients were classified into high and low ADO groups by the median of ADO score in each cohort. Kaplan-Meier survival curve and univariable/multivariable Cox regression model were used to analyze the association between ADO score and overall survival (OS). In order to explore the interaction effects of BRCA1/BRCA2 mutation and Homologous Recombination Deficiency (HRD) with ADO score, hierarchical analysis was applied and interaction P value was evaluated by Wald test. The Number of Telomeric Allelic Imbalances (NtAI) count, Large-scale State Transitions (LST) count, Homologous Recombination Deficiency (HRD-LOH) score and HRD score for TCGA samples were obtained from Andrea M Marquards’ article (19). In the validation set, meta-analysis with random effect was used to pool the hazard ratio (HR) of the six GEO datasets.

To explore the molecular mechanism, the gene set enrichment analysis (GSEA) was performed in TCGA patients to analyze the significantly enriched pathways in high and low ADO score groups. In addition, the association between ADO score and ten cancer hallmark signatures related Gene Ontology (GO) terms were also analyzed (20). In order to explore the relationship between ADO score and immune microenvironment, immune infiltration signature calculated by CIBERSORT (21), ESTIMATE (22), MCPcounter (23), single-sample gene set enrichment analysis (ssGSEA) (24), TIMER (25) and TIDE (26) were compared between two ADO groups using Wilcoxon rank-sum test. The gene expression of immune checkpoint was similarly analyzed. For drug sensitivity, drugs with missing value in more than 80% of the samples were removed from the analysis. Finally, 180 drugs involving 23 pathways were used for Pearson correlation analysis between ADO score and IC50. Benjamini-Hochberg was used for P value correction for multiple comparisons.

The flowchart of this study is demonstrated in Figure 1A . All analyses were performed with R version 4.0.3. The Wilcoxon rank-sum test was used to test difference for all the calculated score or the gene expression between two groups if not specifically described. The P value less than 0.05 was regard as statistically significant and all the tests were two-side.

For the 22 ADO related genes, only 38 of 434 (8.76%) ovarian cancer samples had non-silent mutation in these genes ( Supplementary Figure 1A ). Among them, mutation frequency of ENPP1 was the highest, which occurred in 5 of 434 samples (1.46%). There was no mutation event in ADA2, ENTP2, ENTPD8, and NME1 ( Supplementary Figure 1A ). For the copy number variants, it showed that the copy number variants were ubiquitous for ADO genes which indicated the mRNA expression of ADO genes might be discrepant among samples ( Supplementary Figure 1B ). According to the expression data of ADO genes, patients could be classified into four clusters in consensus cluster analysis ( Figure 1B ). Through NMF analysis, the optimal cluster number was set as 4 according to cophenetic correlation coefficient ( Supplementary Figure 2A ). The four Metagenes of ADO genes had different expression patterns and showed discrepant expression in all samples ( Figure 1C , Supplementary Figure 2B ). Top 15 genes in each metagene were extracted and formed a gene list for the following analysis. According to the correlation among these genes, we got two ADO gene groups ( Figure 1D ). In gene set 1, 8 of 12 genes were involved in the generation of ADP, AMP or ADO. However, in gene set 2, expect for ADA gene, the other 7 genes were enriched in the pathways related to ADP generation or ADO degradation ( Supplementary Table 1 ). Furthermore, the enrichment scores by ssGSVA were calculated for the two gene sets, respectively. The ADO score was defined as enrichment score of gene set 1 minus gene set 2, which was significantly different among the four sample clusters ( Figure 1E ).

In TCGA cohort, patients with high ADO score had shorter OS than those with low ADO score (HR, 1.33; 95% CI, 1.02-1.72, Figure 2A ). After adjusted for age, grade, stage, BRCA1/2 mutation and debulking surgery, the association was still significant (HR, 1.42; 95% CI, 1.06-1.88, Supplementary Table 3 ). Among the six validation cohorts, similar results were observed in GSE49997 and GSE140082 dataset (Log-rank P = 0.04 for GSE49997 and 0.01 for GSE140082, Figures 2B , C ). In the meta-analysis of the six validation cohorts, the pooled HR was still significant in both univariable (HR, 1.27; 95% CI, 1.00-1.62; Figure 2D ) and multivariable Cox models adjusted for age, stage, debulking surgery, histology and histological grade (HR, 1.24; 95% CI, 1.02-1.51, Figure 2E ). The detailed results were reported in Supplementary Table 3 . Besides, the hierarchical meta-analysis was applied stratified by serous and non-serous cancers, different grade and different stage, and the high ADO score was still significantly associated with worse OS in subgroups of stage III/IV (HR, 1.24; 95% CI, 1.01-1.52; Supplementary Table 4 ) and serous histology (HR, 1.23; 95% CI, 1.00-1.50; Supplementary Table 4 ) in the multivariable model. In grade G3/4 subgroup, the association were similar and borderline significant (HR, 1.28; 95% CI, 0.96-1.70; Supplementary Table 4 ).

There is no significant difference in BRCA1/2 mutation frequency, NtAI count, LST count, HRD-LOH score and HRD between low and high ADO score groups ( Supplementary Figures 3A, B ). For BRCA1/2 mutation, the association of ADO score and OS was significant in the patients with BRCA1/2 mutation, but not significant in other patients (P interaction = 0.018) ( Figures 3A , B , Supplementary Table 5 ). Besides, the patients were classified in two subgroups based on the median of HRD score. There was a significant association between ADO score and OS in high HRD subgroup, but not in low HRD subgroup. However, the P interaction was not significant (P interaction = 0.35, Figures 3C , D , Supplementary Table 5 ). The results showed that ADO score tended to be significantly associated with prolonged survival only in BRCA1/2 mutant or HRD patients.

According to GSEA analysis, we found that genes in ATP synthesis related pathways were significantly enriched in the low ADO score group (adjusted P value = 0.02, Figure 4A ). The associations of ADO score with ten cancer hallmarks signatures were analyzed ( Figure 4B ). Patients with high ADO scores had higher tissue invasion and metastasis, sustained angiogenesis, self-sufficiency of growth signals, but lower genome instability than those with low ADO scores (all adjusted P values < 0.01, Figure 4B , Supplementary Table 6 ). For ovarian cancer, patients with high genome instability tend to prolong OS (27, 28). Consistently, we found that patients with high HRD score had better prognosis (Log-rank P < 0.001, Supplementary Figure 4 ).

In order to explore the association between ADO scores and immune microenvironment, five immune signatures algorithms were applied. The heatmap showed that the high ADO group had enriched immunosuppressive cells, including stromal cells and CAF, while the low ADO group had enriched immune activating cells, including cytotoxic cells and M1 macrophage cells (adjusted P value = 0.01 for stromal cells, <0.001 for CAF, <0.001 for cytotoxic cells, 0.03 for M1 macrophage cells, Figure 4C , Supplementary Figure 5 , Supplementary Table 6 ). Besides, for the immune checkpoint genes, compared with low ADO group, the high ADO group had higher expression of CD276, but lower expression of other genes, including CD40, IDO1, LAG3 and TIGIT (adjusted P value = 0.001, 0.02, <0.001, 0.05 and 0.02 respectively, Figure 4D , Supplementary Table 6 ). Further, based on the TIDE algorithm which could estimate the overall immune status, similar results were found. In detail, IFNG, representative of immune activating status was significantly higher in low ADO group than in high ADO group, while the signatures of T cell dysfunction, T cell exclusion and CAF were significantly lower in low ADO group (all adjusted P values < 0.001, Figure 4E , Supplementary Table 6 ).

We constructed ADO score based on the expression data of DGSC ovarian cancer cell lines and calculated the Pearson correlation between ADO score and IC50 of 180 drugs. A total of 26 cell lines was used to obtain the IC50 ( Supplementary Table 7 ). IC50 of 65 drugs was significantly associated with ADO score after Benjamini-Hochberg correction (adjusted P value < 0.05, Figure 5A ). The detailed adjusted P values are shown in Supplementary Table 8 . High ADO score group had significant higher IC50 than low ADO score group for conventional drugs, including Talazoparib, Paclitaxel and Olaparib (adjusted P value = 0.04, 0.04 and 0.05 respectively, Figure 5B ), while not for Cisplatin and Docetaxel (adjusted P value = 0.10 and 0.10 respectively, Figure 5B ). For MK-1775 (WEE1 Inhibitor) which was researched in a phase II clinical trials (29), high ADO score group also had significant higher IC50 than low ADO score group (adjusted P value = 0.007 and 0.05, Figure 5B ).