The analysis workflow and used datasets from various cohorts are summarized in Supplementary Figure S1 . Patient characteristics for the TCGA-OV cohort (n=226) (15) and the new HGSOC cohort from the Medical University in Innsbruck (MUI) (n=60) are listed in Supplementary Tables S1 and S2 . RNA sequencing data and clinical data for the validation cohort (Medical University of Innsbruck; MUI) were deposited at https://doi.org/10.5281/zenodo.10251467. Controlled access data for whole exome sequencing and RNA sequencing data for the TCGA-OV cohort were obtained through dbGaP access permission (phs000178). Processed data (including methylation beta values) and clinical data were downloaded from Firebrowse (firebrowse.org, BROAD Institute). Additional clinical data were retrieved from the supplementary data of another resource (27). Raw RNA sequencing data and clinical annotations for the ICON7 cohort (28) were downloaded from the EGA archive (EGAS00001003487). Single-cell RNA-seq data (29) were downloaded from the Gene Expression Omnibus (GEO) (GSE180661) as an annotated count matrix (anndata-object) in h5ad-format. Data files from the TOPACIO clinical trial (9) were retrieved from Synapse (https://doi.org/10.7303/syn21569629). Data from the Clinical Proteomic Tumor Analysis Consortium ovarian cancer cohort (CPTAC-OV) (30) were downloaded from (https://proteomics.cancer.gov) (n=71). Data from RNA sequencing analysis of OVCAR 3 and UWB1.289 cancer cell lines performed in this study were deposited in GEO (GSE237361). Only complete data sets were used, and observations (rows) with missing values were deleted before specific analyses were performed.

Two epithelial ovarian carcinoma cell lines, UWB1.289 harboring a deleterious BRCA1 and OVCAR3 with intact BRCA1, were obtained from ATCC. OVCAR3 cells were grown in RPMI 1640 with 0.01 mg/ml bovine insulin and 20% FBS, whereas UWB1.289 cells were grown in a mixture of 48.5% MEGM Bullet Kit medium (Lonza) and 48.5% RPMI 1640 with 3% FBS. Viability assays were used to determine the IC50 for olaparib. Both cell lines were treated with olaparib or DMSO for 96 hours in four replicates. Treated and untreated UWB1.289 and OVCAR3 cells were stained with indirect immunofluorescent antibodies to detect γH2AX as an indicator of double-strand breaks. To determine activated STING signaling, double-stranded DNA and its presence in the cytosol, cGAS, STING, and phosphorylated STING were detected. The antibodies used are listed in Supplementary Table S3 .

Slices of 10 selected tumor blocks were subjected to immunohistochemistry analyses performed on the BenchMark ULTRA automated staining device (Ventana, Oro Valley, AZ/Roche, Vienna, Austria). The examined markers were CD163 for macrophages and CD8, PD-1, CD4, and FOXP3 for T cells. Furthermore, the markers γH2AX and STING were analyzed. All antibodies used are listed in Supplementary Table S4 .

RNA from cancer cell line samples was isolated from 2x10⁶ cells each using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s protocols. RNA quantity and quality were assessed using NanoDrop™ 2000c and Bioanalyzer 2100 with Agilent 6000 Nano Kit and cDNA libraries were generated using the QuantSeq 3’ mRNA-Seq Library Prep Kit (Lexogen) according to the manufacturer’s instructions. Paired-end sequencing (150 bp) was performed on a NovaSeq 6000 sequencing device at GENEWIZ/Azenta. RNA isolation from 60 fresh frozen tumor samples from the MUI HGSOC validation cohort was conducted in a similar manner resulting in sufficient quality (RIN factors from 6.4 to 9.9), and sequencing was performed at Novogene (Cambridge, UK) for paired end sequencing (PE150) on an Illumina NovaSeq 6000 sequencing device using TrueSeq (Illumina) strand-specific total RNA libraries.

Raw reads were quality checked using FastQC. Reads were mapped to the human reference genome version hg38 (GRch38) using STAR (version 2.7.1) in 2pass mode (31). Gene level expression quantification was performed with featureCounts (version 2.0.0) using GENCODE annotations (v36). Raw counts were normalized using TPM (transcripts per million). RNA sequencing raw data from the MUI HGSOC cohort and the ICON7 cohort were analyzed in the same way. For sequencing data of the cell lines, single-end reads were processed by trimming adapter and low-quality sequences using BBDuk with the parameters specified by Lexogen. The trimmed reads were mapped to the human reference genome version hg38 (GRch38) using STAR (version 2.7.9a) in 2-pass mode. Gene level expression quantification was performed with featureCounts (version 2.0.0) and GENCODE annotations (v38).

Raw exome sequencing reads in fastq format were quality checked using FastQC. Reads of paired tumor and normal samples were mapped against the human reference genome version hg38 (GRch38) using BWA. For the detection of germline variants the HaplotypeCaller was used. To assess somatic variants in the tumor samples, four different variant callers, Mutect2 (32), SomaticSniper (33), Varscan2 (34), and Strelka2 (35) were used. If a variant was called by two of four variant callers and the variant allele frequency was ≥ 0.05 in the tumor sample and <0.05 in the normal sample, the variant passed filtering. Variants were annotated using VEP (36) with the ClinVar extension. Only pathogenic (class V) and likely pathogenic (class IV) variants were considered to affect the function of HRR genes such as BRCA1 or BRCA2. Tumor mutational burden was calculated based on the number of nonsynonymous single nucleotide variants per megabase for each tumor sample. For neoantigen prediction, from somatic mutation derived peptide sequences with lengths between 8-11 amino acids - taking phasing into account - were generated and tested for the respective HLA alleles with NetMHCpan-4.0 (37), whereby %rank<2 was considered a weak binder and %rank<0.5 was considered a strong binder. Dissimilarity to the normal human proteome (hg38) was identified by the antigen.garnish package. Neoantigen load was calculated for each tumor based on predicted weak and strong binding neoantigens – irrespective of their peptide length and taking all HLA alleles into account – per megabase.

Differential gene expression analysis was conducted using the R package DESeq2 (38). P values were adjusted for multiple testing based on the false discovery rate (FDR) according to the Benjamini−Hochberg method. Genes with more than a twofold change at an FDR<0.1 and average expression across all samples (baseMean>10) were considered significantly differentially expressed. To identify functional annotation and affected biological processes, log2-fold change preranked gene set enrichment analyses (GSEA) (39) using hallmark and selected immune-related gene sets from MSigDB were performed. ClueGO was used to build a network and group significantly overrepresented pathways, which share genes (40). The STRING database (https://string-db.org/) was used to identify an interaction network within the differentially expressed genes, and subnetworks were found by MCL clustering with inflation parameter=3. Footprint analyses of response genes of perturbed cancer signaling pathways were performed using PROGENy (41). To assess tumor infiltration of immune cells quanTIseq (42) using the immunedeconv R package was applied to bulk RNA sequencing data. To characterize the immune-related processes, well-described immune signatures, such as T-cell inflammation, IFN gamma signature, cytolytic activity, cytotoxic T lymphocyte function, and T-cell exhaustion ( Supplementary Table S5 ), were analyzed. Based on log2(TPM+1) normalized expression data, single sample gene set enrichment using GSVA (43) was performed for signatures with more than 10 genes or otherwise average expression was calculated. The tumor-immune phenotype (infiltrated, excluded, desert) was determined based on a previously developed classification model based on 157 genes using digital pathology describing the presence and position of CD8+ T cells relative to the center or margin of the tumor (28) and a random forest model was used to characterize samples from the TCGA and the MUI cohorts. To classify ovarian cancer samples into molecular subtypes, the consensusOV R package (44) was used. The immunophenoscore (IPS) was determined as described previously (45).

BRCAness was determined based on HRD scores (46), mutational signature 3 (47), mutations in homologous recombination repair pathway genes and methylation of promoter regions of BRCA1. All HGSOC samples of the TCGA cohort for which paired tumor and normal exome sequencing and matched RNA sequencing data were available (n=226) were used. Samples were classified with a BRCAness phenotype when they had either a deleterious mutation in the homologous recombination pathway, an ovarian cancer-specific HRD score of ≥ 63 (48), a mutational signature 3 ratio > 0.25 or a methylation level beta value >0.7 of the BRCA1 promoter. HRD scores were calculated as the unweighted sum of the three genomic scar values, loss of heterozygosity (LOH) (49), telomeric allelic imbalance (TAI) (50), and large-scale state transitions (LST) (51). To compute the genomic scar values, scarHRD (52) was used on genome segmentation files generated with Sequenza (53). The mutational signature 3 score was computed using MutationalPatterns (54) and we calculated the ratio between mutational signature 3 supporting mutations and all detected mutations. To classify BRCAness, genes expressed in at least one ovarian cancer cell were identified using single-cell RNAseq data. Normalized gene expression values (log2 (TPM+1)) of these genes in the TCGA dataset were then subjected to recursive feature elimination in a balanced design with three different machine learning models (random forest, AdaBoost and gradient boosting) to identify the 50 most important features for each model. Since ensemble methods and random subsampling (bootstrap) were included, we did not use nested cross-validation to avoid overfitting. Genes that were among the top 50 in at least two of the three models (24 genes) were then used subsequently to train a random forest classification model to discriminate between BRCAness and noBRCAness samples based on gene expression data. The performance of the classifier was evaluated by analysis of the receiver operating characteristic curve with 10-fold cross-validation. The area under curve (AUC) was used as a performance measure. A cutoff for BRCAness (P>0.5266) was selected using the Youden index. Furthermore, the classifier was tested in 29 patients of the independent validation cohort (MUI) with HRD information based on SNP arrays and further validation using Myriad MyChoice CDx. Samplewise BRCA classification in single-cell RNA sequencing data from 29 patients was performed with an optimized cutoff (P>0.45) and based on the majority of classified tumor cells. A method to detect mutational signature 3 (Sig3), termed SigMA, from clinical panel sequencing data have been previously developed and associated with HRD (55). To further compare our BRCAness classifier, which is based on a BRCAness definition more completed as compared to the other methods described in literature such as the SigMA score, using available RNA sequencing and mutation data from the Clinical Proteomic Tumor Analysis Consortium ovarian cancer (CPTAC-OV) cohort (30).

All analyses of single-cell data were performed in Python using scanpy (56) and scvi-tools (57). Since the samples were sequenced separately for sorted CD45+ and CD45- cells, the raw read counts were integrated using scvi-tools with a batch effect correction. Counts were normalized to counts per million (CPM) and log2 transformed, adding a pseudocount of 1. Quality metrics were determined using scanpy and filtered for genes that are expressed in at least one cell. The dataset was filtered for samples from the primary tumor (adnexal tumor tissue). Principal component analysis and nearest neighbor analyses were calculated with default settings, and clustering was performed with the Leiden algorithm. Super cell types were annotated as previously defined. Subtypes of T-cell and myeloid cell clusters were assigned based on the expression of marker genes using published marker genes for different cell types and the PanglaoDB (58). Differentially expressed genes between clusters were calculated using the Wilcoxon ranked sum test. For visualization, we used uniform manifold approximation and projection (UMAP) dimensional reduction. Gene expression between cell types was compared by heatmaps, violin plots, and bubble plots. Distribution of cell fractions for each sample or combined for BRCAness and noBRCAness group were compared by stacked barplots and two-sided Wilcoxon rank sum test. Pseudobulk analyses and DESeq2 analyses was performed for selected immune response markers to test effect of BRCAness versus noBRCA samples on the immune response. To assess ligand−receptor interactions between cancer cells and cells from the TME, CellPhoneDB (59) analysis was used.

Gene expression analysis in the TOPACIO cohort was performed on the NanoString platform. We used the nSolver software from NanoString (Seattle, US) to obtain normalized data. Differential expression analysis was performed using the R package limma (60), and genes with p<0.05 were considered differentially expressed.

Vulnerability maps consist of three variables: the vulnerability score, the BRCAness score and the cytolytic activity (CYT) to C1QA ratio. For the BRCAness score, the prediction probability from the random forest classifier was used. The CYT to C1QA ratio was calculated from the log2 (TPM+1) values of GZMB, PRF1, and C1QA (Equation 1).

The CYT to C1QA ratio was transformed to values between 0 and 1 using a sigmoid function with softmax transformation and parameters derived from the TCGA cohort and termed C2C (Equation 2).

The vulnerability score was defined as the weighted sum of BRCAness probability and C2C (Equation 3), whereby the weights were identified using a logistic regression model on the CYT to C1QA ratio using log2 intensity expression values and SigMA status (mutational signature 3) as proxy for BRCAness from the TOPACIO cohort and the treatment response as a binary dependent variable.

For visualization of the vulnerability map, a two-dimensional map was created with C2C as one coordinate, BRCA probability as the other coordinate, and the color-coded vulnerability score.

Survival analyses were performed for both HGSOC cohorts (TCGA, MUI) for selected genes, immune parameters, or immune cell fractions by dichotomization of patients based on the median or maximum log-rank statistics using the R package survival. For the TCGA cohort, overall survival (OS) and progression free survival (PFS) survival status were derived from a clinical data resource for TCGA (27) and for the cohort from Medical University Innsbruck from the clinical data as provided by the Department of Obstetrics and Gynecology. Univariate and multivariable Cox regression analyses taking clinical parameters into account (age, FIGO stage, residual tumor) were performed, Kaplan-Meier survival curves were generated, and compared by log rank test. To determine the association between continuous or binary variables, point biserial correlation analysis was used. For the correlation between binary variables, the Phi coefficient and chi-square test were used, and for the correlation between continuous variables, Pearson’s correlation coefficient was used. To compare parameters between two groups, the Wilcoxon rank-sum test was used. For multiple group comparisons the non-parametric Kruskal-Wallis test followed by pairwise two-tailed Dunn posthoc tests with p-value adjustment based on the false discovery rate (FDR) were used. Where indicated, p values were adjusted for multiple testing based on the FDR according to the Benjamini−Hochberg method. P<0.05 or FDR<0.1 were considered significant.

Because the response to platinum-based chemotherapies or therapy with PARP inhibitors in ovarian cancer is not limited to patients with tumors harboring BRCA1 or BRCA2 mutations, we expanded the group of patients by using a genomic characterization termed BRCAness, which has very much in common with homologous recombination repair deficiency (HRD) status (61). BRCAness status includes mutations of genes in the homologous recombination DNA repair pathway (HRR), genomic scars, loss of heterozygosity, telomeric allelic imbalance, or large-scale transitions, mutational signature 3, or promoter methylation of the BRCA1 or BRCA2 gene. We assessed these parameters based on whole exome sequencing data and methylation data from the TCGA OV cohort ( Figure 1A ). Very few patients harboring HRR mutations or BRCA1 promoter methylation fell below the combination of the HRD cutoff (HRD>63) and the MutSig3 ratio cutoff (0.25), indicating a reasonable selection of the cutoff values ( Figure 1B ). To identify BRCAness solely based on gene expression data, we developed a machine learning classifier that can discriminate between BRCAness and non-BRCAness samples using bulk and single-cell RNA sequencing data ( Figure 1A ). Recursive feature elimination based on multiple models resulted in a BRCAness gene expression signature with 24 genes, which was used to train a random forest model discriminating between BRCAness and noBRCAness. The receiver operating characteristics with 10-fold cross-validation on the training dataset showed an area under the curve (AUC) of 0.91 ± 0.04 ( Figure 1C ). To validate our BRCAness classifier we tested its performance in two independent ovarian cancer cohorts. We could classify BRCAness in a new HGSOC cohort from the Medical University Innsbruck (MUI) (n=60) based on RNA sequencing data with an accuracy of 0.79, an F1-score of 0.86, and a positive prediction value of 0.86 ( Figure 1D ). Furthermore, we demonstrated that in addition to classifying bulk RNAseq samples the classifier is also capable of classifying samples from single cell RNA-seq data at the sample level in the HGSOC MSK cohort (n=29) with an accuracy of 0.86, an F1 score of 0.87 and a positive prediction value of 0.87 ( Figure 1D ).

There was also good agreement with a recently defined gene expression-based HRDness signature including 173 up- and 76 downregulated genes (62) using a single sample gene set enrichment derived score in the TCGA cohort as well as the MUI validation cohort with Spearman’s rank correlation of ρ=0.72 (P<0.001) and ρ=0.63 (P<0.001), respectively ( Supplementary Figures S2 , S3 ). Interestingly, six genes from our 24-gene signature ( Figure 1E ) to classify BRCAness (CCDC90B, CRABP2, FZD4, GPAA1, PRCP, SNRP1) were also among the upregulated and two genes (RAD17, LTA4H) among the downregulated genes. The 24-gene BRCAness signature was further compared to the SigMA score (mutational signature 3) of the CPTAC-OV cohort (n=71). Although our BRCAness signature is based on a more complete BRCAness definition than the mutational signature, a significant correlation was observed with a Spearman’s rank correlation of ρ=0.43 (P<0.001) ( Supplementary Figure S4 ). We used this BRCAness classification for all remaining analyses.

In summary, we developed a 24-gene-based BRCAness model validated in several single-cell and bulk RNAseq datasets with reasonable classification performance.

To identify the relationship between genomic instability and the activation of the immune system, we performed correlation analyses between the BRCAness status and various immune-related signatures in the TCGA cohort (n=226). BRCAness could be significantly positively associated with the enrichment of immune-related signatures, such as those for IFNG response (rho=0.38, p=0.004) and T-cell inflamed tumor microenvironment (rho=0.46, p=0.0014), even to a larger extent with high tumor mutational burden (p<0.001) and high neoantigen load (p<0.001) ( Figures 2A, B ). However, compared to other cancer types with defective DNA mismatch repair the TMB or neoantigen load in ovarian cancer is rather low. Thus, increased immune activity is more indicative of deficient HRR. It is known that BRCAness is associated with longer overall survival, indicating that those patients are more responsive to platinum-based chemotherapy. In order to determine to which extend this could be explained by higher immune cell infiltration we estimated CD8+ T-cell infiltration from RNA sequencing data using quanTIseq and divided patients into 4 groups: BRCAness patients with high CD8 T cell fraction, BRCAness patients with low CD8 T cell fraction, noBRCAness patients with high CD8 T cell fraction, and noBRCAness patients with low CD8 T cell fraction. We observed a significant association with overall survival (p<0.0001; log rank test) ( Figure 2C ) and progression-free survival (p=0.0016; log rank test) ( Figure 2D ), with the BRCAness group with a high proportion of positive CD8 T cells being associated with the longest survival. Multivariable Cox regression analysis showed a significant effect of BRCAness vs. no BRCAness on overall survival (p=0.00062, HR=0.51with 95%-CI 0.35-0.75), while the impact of CD8 T cell was not significant ( Supplementary Table S6 ) indicating a more pronounced role of the BRCAness status. Nevertheless, analyses of signaling pathways by downstream target expression using PROGENy indicated for the TCGA cohort (n=226) as well as the MUI validation cohort (n=60) that immune-related pathways, including TNFa, NFkB, and JAK-STAT, were activated in the BRCAness samples ( Figures 2E, F ). Using STRING analyses in the MUI cohort (n=60), we also identified a highly connected network including various chemokines and interleukins and their respective receptors (CCL7, CCL11, CXCL5, CXCL9, CXCL13, CCR2, CCR3, CCR4, CCR8, CXCR3, and IL6), which were significantly higher expressed in BRCAness tumors than in non-BRCAness tumors, indicating attraction and interaction with various immune cells ( Supplementary Figure S5 ).

Essentially, we observed a correlation of BRCAness with various immune-related processes and, in particular, a group of patients with BRCAness and high CD8 T-cell proportion associated with longer survival times in HGSOC patients.

To study the effect of PARPis on immune activation, we performed in vitro analyses. As tumor models, an ovarian cancer cell line with a proficient BRCA1 gene (OVCAR3) and a cell line with a mutation in the BRCA1 gene (UWB1.289) were used. We performed RNA sequencing analyses to identify differentially expressed genes between olaparib (PARPi)-treated and control (DMSO)-treated cell lines. Significantly upregulated genes ( Figures 3A, B ) indicate activation of various processes ( Figures 3C, D ), including pattern recognition receptor activation, response to cytokine signaling, interferon alpha response (type I), NFkB pathway, and cGAS-STING signaling. To further validate the results at the protein level, we performed immunofluorescence analyses indicating effects on γH2AX –a surrogate marker of DNA damage - by mutation in the BRCA1 gene and an even stronger effect by olaparib (PARPi) treatment ( Figure 3E ). Similarly, we observed a different activation of cGAS and STING – and indicating activation of the (innate) immune system – in the BRCA1-deficient versus the BRCA1-proficient cell model ( Figure 3F ). Furthermore, using gene set enrichment analyses, a significant interferon alpha response was also observed in BRCAness samples of both the TCGA cohort and the MUI validation cohort ( Figure 3C ).

In summary, we observed cGAS-STING activation by olaparib treatment in vitro and an interferon type I response as well as chemokine expression in HGSOC patient cohorts with BRCAness status.

We next focused on characterizing the presence of cytotoxic T lymphocytes and their spatial distribution in the tumor, following a recent approach in which digital pathology could be linked to gene expression (28). With the reported list of 157 genes and using random forest analysis, we were able to divide the patients of the TCGA cohort into a group with infiltrated, excluded, or desert tumor-immune phenotypes. Interestingly, the excluded phenotype was associated with upregulation of TGFβ and high expression of markers for cancer-associated fibroblasts, such as FAP or PDPN, which could form a barrier to prevent T-cell infiltration ( Figure 4A ). The immunoreactive molecular subtype (IMR) is very relevant to delineate immunoreactivity because many of the immunity genes, including cytotoxic effectors, factors involved in antigen processing and presentation, or immune checkpoints, are highly expressed in this condition ( Figure 4A , Supplementary Table S7 ). We have selected a group of patients with tumor BRCAness, an infiltrated tumor immune phenotype and an immunoreactive molecular subtype called BRCAness immune type (BRIT), which we expect to respond well to combination immunotherapy. When comparing the estimated immune cell infiltrates in these cancer samples with BRCAness cancers without immune type (noBRIT), we found that not only cytotoxic T lymphocytes such as CD8+ T cells were significantly more abundant (p=0.001) but also a number of suppressive immune cells (M2 macrophages (p<0.001), regulatory T cells (p=0.005), myeloid-derived suppressor cells; MDSCs (p<0.001) ( Figure 4B , Supplementary Figure S6 ). This is in line with previous observations (63) and in order to identify an effect by BRCAness we additionally defined an immune type (IMT) with noBRCAness, infiltrated tumor immune phenotype (INF), and an immunoreactive molecular subtype (IMR). However no significant difference between BRIT and IMT as well as noBRIT and noIMT could be observed for the analyzed cell types ( Figure 4B , Supplementary Figure S6 ) indicating a more pronounced effect by immune infiltration than by BRCAness. Furthermore, we did not observe a significant difference in overall survival between BRIT versus noBRIT patients (p=0.56, HR =0.81, 95% CI 0.42-1.60).

These observations underscore the importance of the suppressive immune environment and suggest that suppressive immune cells may be an important factor, which is why ovarian cancer patients have a limited response to immunotherapy.

Single-cell RNA sequencing analyses allow a more comprehensive characterization of the tumor environment and evaluation of the cell interplay. Analyses of more than 300 thousand cells of adnexal ovarian tumor tissue from 29 patients of the MSK cohort (29) allowed a clear separation between major cell type populations by clustering and nonlinear projection (UMAP) ( Figure 5A , Supplementary Figure S7 ). In contrast to cell types from the tumor microenvironment, tumor cells showed a clear separation between BRCAness and noBRCness samples ( Figure 5A ). To further investigate the effect of BRCAness on the tumor microenvironment we compared the distribution of major cell types between BRCAness and noBRCAness samples. We did not observe any significant differences for the major cell types between these patient groups ( Supplementary Figure S8 ) but in a summarized analysis ( Figure 5A ) the proportion of myeloid cell was slightly higher and the proportion T/NK cells was slightly lower in BRCAness compared to noBRCAness. Furthermore, we investigated the differential expression of immune response marker genes ( Supplementary Table S7 ) in the BRCAness versus the noBRCAness group by pseudobulk analysis with DESeq2. Interestingly, 21 genes have been found significantly upregulated in BRCAness compared to noBRCAness including immune checkpoints such as CD274 and a number of antigen processing and presentation genes ( Supplementary Figure S9 ). Because cells from the suppressive environment have a major impact, we focused on the myeloid cell compartment and demonstrated that the majority of these cells were macrophages, and we identified subpopulations based on most dominant marker genes, including CD169 (SIGLEC1) macrophages, CX3CR1 macrophages, and MARCO macrophages ( Figure 5B ). One described hallmark marker of TAMs is TREM2, which has been identified as an attractive target for cell depletion therapy and is being tested in an ongoing clinical trial (64). Notably, the expression patterns of TREM2 and BRCAness are very similar, showing high expression in all macrophage subtypes and, to a lesser extent, in monocytes ( Figure 5B ). To search for further genes with similar expression patterns in myeloid subpopulations, we analyzed known tumor-associated macrophage and monocyte marker genes (65). As indicated by this analysis, C1QA showed a similar but even more pronounced expression pattern than TREM2 ( Figures 5B, C ). C1QA was also recently described as a surrogate marker for the CD68+CD163+ macrophage subset (66). We found that several genes are highly expressed in macrophages ( Supplementary Figure S10 ) and that, based on marker genes, a polarization towards M2 macrophages occurs ( Supplementary Figure S11 ). The immune suppressive effect of tumor associated macrophages were underscored by a number of myeloid immune checkpoint genes, which show a worse effect on overall survival (hazard ratio>1) in the TCGA cohort ( Supplementary Figure S12 ).

Since our main goal was to predict vulnerability to combination immunotherapy we took advantage of the availability of gene expression data from a clinical trial (TOPACIO) and could identify 7 upregulated genes and 22 down regulated genes between responder and non-responder to PARPi-immune checkpoint inhibition combination therapy (niraparib and pembrolizumab) ( Figure 5D ). We next analyzed in which cell types these genes are generally expressed using the results from our single cell RNA-sequencing data analyses from the MSK cohort and found a number of downregulated genes in responders, such as LYZ, LILRB4, and ITGB2, were most highly expressed in myeloid cells (macrophages), LILRB4 in dendritic cells, and integrin subunit beta 2 (ITGB2) in other cell types, such as T/NK cells ( Figures 5E, F ). ITGB2 was also found correlated in the TCGA cohort with estimated M2 macrophages and CD8 T-cell infiltration ( Supplementary Figure S13 ).

Interestingly, we identified various ligand−receptor interactions with expressed ligands in tumor cells and respective receptors expressed in tumor-associated macrophage subsets using CellPhoneDB (59) ( Supplementary Figure S14 ). The growth arrest-specific protein 6 (GAS6) – AXL tyrosine kinase (AXL) interaction, for example, which are both associated with poor outcome, have already been evaluated in clinical trials in ovarian cancer by inhibiting their interaction (67). LILRB1 and LILRB2 expressed in macrophage subsets were found to interact with the nonclassical human leukocyte antigen HLA-F expressed in cancer cells. Blocking macrophage colony-stimulating factor CSF1 and its receptor CSF1R axis and several drugs that target these factors have been under investigation (68).

These observations summarized together suggest that TAMs may not only play a role in immunotherapy alone but are also essential in informing about therapy response when combined with PARP inhibitors.

To validate the results, we performed RNA sequencing analyses of an HGSOC cohort of patients from Medical University Innsbruck (n=60). Stratification of these patients resulted in very similar expression patterns evident from a number of immune marker genes, which were highly expressed in the BRCAness immune type patient group (BRIT) ( Figure 6A ). To further characterize immune infiltrates in different patient groups, we performed immunohistochemistry analyses on ten selected samples for various markers and highlight the results from three patient samples. One BRIT tumor sample showed high γH2AX activity, STING activation, CD8+ T-cell infiltration, CD4+ T-cell infiltration, and strong CD163+ tumor-associated macrophage populations ( Figure 6B , left panel). These effects were even more pronounced in one sample with no detected BRCA1 or BRCA2 mutation, underscoring the importance and validity of predicted BRCAness ( Figure 6B , middle panel). Another tumor sample with no BRCAness, a desert tumor-immune phenotype, and a differentiated molecular subtype was used as a negative control, and in fact, no activity for any of the tested markers was observed ( Figure 6B , right panel). To better address the potential for combination immunotherapy response, we again took advantage of data from the TOPACIO trial and, based on the clinical response, trained a logistic regression model and learned weights for three surrogate variables: MutSig3 as an indicator for BRCAness, average expression of PRF1 and GZMB as indicators for cytolytic activity, and expression of C1QA as an indicator for tumor-associated suppressive macrophages. Based on the HGSOC samples from TCGA, we developed a two-dimensional vulnerability map, with the ratio of cytolytic activity and C1QA expression as one variable (C2C) and the BRCAness prediction probability as the other variable. The vulnerability score is indicated by color ( Figure 6C ). When applied to the selected examples from the MUI validation cohort, these differed significantly for areas with high vulnerability scores (indicating response to combination immunotherapy) compared to the negative control with low vulnerability scores ( Figure 6C ). Furthermore, we observed a significant difference in overall survival between patients with high and low vulnerability scores (p<0.001, HR = 0.47, 95% CI 0.33-0.66) in the TCGA HGSOC cohort ( Supplementary Figure S15 ), indicating a positive association of a high vulnerability score with longer overall survival. For patients in the MUI validation cohort, no significant difference in overall survival (p=0.368, HR = 0.78, 95% CI 0.45-1.34) could be revealed ( Supplementary Figure S16 ). To enable the characterization of newly diagnosed HGSOC samples based on RNA sequencing data, we developed an easy-to-use R package (OvRSeq), which allows us to not only estimate the parameters to determine the vulnerability score (and generate the vulnerability maps) but also comprehensively annotate the sample for BRCAness, tumor-immune phenotype, molecular subtype, estimate immune infiltrates, enrichment of immune-related signatures, and individual marker genes. This also includes other clinically relevant parameters, such as the angiogenesis score we previously defined, which might be useful for the prediction of anti-VEGF therapy (69). The web application (https://ovrseq.icbi.at) allows the generation of summary information as a report of individual samples ( Supplementary Figure S17 ).

The developed application should ultimately be useful to identify vulnerabilities and support clinical therapy decisions for HGSOC patients.