Integrated immunogenomic analyses of high-grade serous ovarian cancer reveal vulnerability to combination immunotherapy

Background: The efficacy of immunotherapies in high-grade serous ovarian cancer (HGSOC) is limited, but clinical trials investigating the potential of combination immunotherapy including poly-ADP-ribose polymerase inhibitors (PARPis) are ongoing. Homologous recombination repair deficiency or BRCAness and the composition of the tumor microenvironment appear to play a critical role in determining the therapeutic response.

Methods: We conducted comprehensive immunogenomic analyses of HGSOC using data from several patient cohorts. Machine learning methods were used to develop a classification model for BRCAness from gene expression data. Integrated analysis of bulk and single-cell RNA sequencing data was used to delineate the tumor immune microenvironment and was validated by immunohistochemistry. The impact of PARPi and BRCA1 mutations on the activation of immune-related pathways was studied using ovarian cancer cell lines, RNA sequencing, and immunofluorescence analysis.

Results: We identified a 24-gene signature that predicts BRCAness. Comprehensive immunogenomic analyses across patient cohorts identified samples with BRCAness and high immune infiltration. Further characterization of these samples revealed increased infiltration of immunosuppressive cells, including tumor-associated macrophages expressing TREM2, C1QA, and LILRB4, as specified by single-cell RNA sequencing data and gene expression analysis of samples from patients receiving combination therapy with PARPi and anti-PD-1. Our findings show also that genomic instability and PARPi activated the cGAS-STING signaling pathway in vitro and the downstream innate immune response in a similar manner to HGSOC patients with BRCAness status. Finally, we have developed a web application (https://ovrseq.icbi.at) and an associated R package OvRSeq, which allow for comprehensive characterization of ovarian cancer patient samples and assessment of a vulnerability score that enables stratification of patients to predict response to the combination immunotherapy.

Conclusions: Genomic instability in HGSOC affects the tumor immune environment, and TAMs play a crucial role in modulating the immune response. Based on various datasets, we have developed a diagnostic application that uses RNA sequencing data not only to comprehensively characterize HGSOC but also to predict vulnerability and response to combination immunotherapy.

1 Introduction

Despite newer therapeutic concepts, ovarian cancer, particularly high-grade serous ovarian cancer, is still the deadliest gynecologic malignancy, with 13,270 expected deaths in 2023 in the U.S (1). While immunotherapy, such as immune checkpoint inhibition monotherapy (e.g., antibodies against PD-1 or PD-L1), has dramatically changed the therapeutic concepts of different cancer types, especially those with mismatch repair deficiency (2), the benefit for ovarian cancer patients with an objective response rate of approximately 10% was found to be rather modest (3–6). However, poly-ADP-ribose polymerase inhibitors (PARPis) and antiangiogenic therapy have improved the survival outcomes of ovarian cancer patients beyond standard care, namely, debulking surgery and platinum-based therapy (7). Furthermore, a number of clinical trials of combination therapies, including immune checkpoint blockade, are underway (8–12). Whereas the recent primary analysis of the double-blind placebo-controlled ENGOT-Ov41/GEICO 69-O/ANITA phase III trial showed that the addition of the anti-PD-L1 antibody (atezolizumab) did not significantly improve the clinical outcome (12), early analysis of the MEDIOLA phase II study adding the PD-L1 inhibitor (durvalumab) and the angiogenesis inhibitor (bevacizumab) to a PARPi (olaparib) was promising, with an objective response rate >90% for a specific patient group with platinum-sensitive relapsed ovarian cancer harboring germline BRCA mutations (11).

PARP is involved in DNA damage and repair, binds to single-strand DNA breaks, and performs posttranslational modifications of histones and DNA-associated proteins by poly-ADP-ribosylation, also known as parylation. PARP inhibitors trap PARP and stall the replication fork, which can subsequently cause DSBs. PARP inhibition is synthetic lethal with deleterious BRCA1 and BRCA2 mutations because homologous recombination repair (HRR) cannot restore these double-strand breaks, introducing genome instability by nonhomologous end joining or leading to tumor cell death (13). In high-grade serous ovarian cancer, approximately 14% harbor a germline and 6% a somatic mutation in the BRCA1 or BRCA2 gene, and approximately 50% are HRR deficient (HRD), indicating that they are favorable for PARPi therapy (14, 15). Sequencing approaches enable researchers to detect mutations in other genes involved in HRR. However, the concept of HRD or BRCAness goes beyond, as it encompasses instabilities and genomic scars, including large-scale transitions, loss of heterozygosity, telomeric allelic imbalance and specific mutational processes with uneven base substitution patterns (mutational signature 3). Several diagnostic assays from commercial providers for the detection of HRD have already been approved (16). However, further efforts are undertaken to identify various biomarkers based on different modalities, such as gene expression or methylation, in the context of different cancer types (17–21). Deleterious BRCA1 mutations and/or PARP inhibition can trigger an immune response at least in part through the cGAS-STING pathway (22–26), suggesting advantages for combined immunotherapies. However, biomarkers or phenotypes to predict the response to therapies, including PARPis and immune checkpoint blockers, are lacking.

In this study, we conducted comprehensive immunogenomic analyses of HGSOC using data from multiple patient cohorts. Integrated gene expression analysis and machine learning on bulk and single-cell RNA sequencing data enabled the 1) development of a 24-gene expression classification model for BRCAness, 2) stratification of patient samples with BRCAness and high immune infiltration, whereby tumor-associated macrophages (TAMs) proved to be an important suppressive component, 3) identification of the activation of immune-related pathways such as the cGAS-STING or JAK-STAT pathway and downstream signaling by PARPi and BRCA1 mutation (BRCAness), and 4) the development of a diagnostic application from RNA sequencing data to comprehensively characterize HGSOC samples and predict vulnerability and response to combination immunotherapy.

2 Methods

2.1 Patient cohorts and datasets

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.

2.2 Cell line experiments

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.

2.3 Immunohistochemistry analyses

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.

2.4 RNA sequencing analyses

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.

2.5 RNA sequencing data analyses

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

2.6 Whole exome sequencing analyses and variant calling

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.

2.7 Functional analysis of gene expression and the tumor immune environment

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

2.8 Determination and classification of BRCAness

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

2.9 Single-cell RNA sequencing analysis

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.

2.10 Gene expression analysis

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.

2.11 Vulnerability score and maps

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

$\begin{array}{ll}CYT to C1QA ratio=0.5*(GZMB+PRF1)/C1QA& (1)\end{array}$

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

$\begin{array}{ll}C2C=1/\left(1+e^{-\frac{CYT to C1QA ratio-0.301}{0.0433}}\right)& (2)\end{array}$

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.

$\begin{array}{ll}Vulnerability score=2.597 * BRCAness probability+1.166 * C2C& (3)\end{array}$

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.

2.12 Statistical analysis

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.

3 Results

We performed immunogenomic characterization and multimodal integrative analyses of data from several patient cohorts, including the TCGA cohort (n=226), the MSK cohort (n=29), the ICON7 cohort (n=327), the CPTAC cohort (n=71), patients from the TOPACIO study receiving combination immunotherapy (n=22), and a new MUI cohort for validation (n=60). The used data modalities and complete analysis workflow is outlined in Supplementary Figure S1.

3.1 A 24-gene signature predicts BRCAness in HGSOC patients

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

Figure 1

Figure 1. BRCAness classification based on the expression of 24 genes. (A) Determination of BRCAness in the TCGA-OV cohort and the development of a gene expression-based BRCAness classifier. (B) Different BRCAness parameters in the TCGA cohort compared between the HRD score and the mutation signature 3 ratio. Samples with mutated homologous recombination repair pathway genes are marked in red, BRCA1/2 promoter methylation in blue and samples with an HRD score > 63 and/or a signature 3 ratio > 0.25 but no mutation or BRCA1/2 promoter methylation are marked in yellow. Samples without BRCAness are marked in white. (C) Mean ROC curve with 10-fold cross-validation of the classifier tested on the TCGA dataset. (D) Confusion matrices with correctly and incorrectly classified instances when the classifier was tested in independent test cohorts of single-cell RNA sequencing and bulk RNA sequencing data. (E) Z scores of log2(TPM+1) normalized expression of the 24 genes of the BRCAness signature in the TCGA cohort as a heatmap clustered by BRCAness and non-BRCAness samples.

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.

3.2 Genome instability is associated with immune-related processes

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

Figure 2

Figure 2. Association between BRCAness and immune parameters. (A) Results of correlation analysis of selected immune signatures and BRCAness parameters in the TCGA-HGSOC cohort (CYT, cytolytic activity; CTL, cytotoxic T lymphocytes; IFNG, interferon gamma signature; HRR mutations, mutations in the homologous recombination repair pathway; NeoAG load, neoantigen load; TMB, tumor mutational burden); white dots indicate significance (FDR<0.1). (B) Direct comparison of selected immune parameters between BRCAness and noBRCAness samples with significant differences, Wilcoxon rank-sum test (FDR<0.1) in the TCGA cohort (n=226). (C) Kaplan−Meier curves according to overall survival (OS) and for 4 patient groups of the TCGA dataset (n=226) based on BRCAness information and median dichotomized estimated CD8 T cell fraction (quanTIseq): patients with 1) BRCAness and high estimated CD8 T cell fraction, 2) BRCAness and low estimated CD8 T cell fraction, 3) noBRCAness and high estimated CD8 T cell fraction, and 4) noBRCAness and high estimated CD8 T cell fraction (p-value is from logrank test). (D) Kaplan−Meier curves according to progression free survival (PFS) for the same groups of patients from the TCGA cohort (n=226) as in (C). (E, F) Waterfall plot of normalized enrichment scores (NES) for the footprint analysis of immune-related pathways with PROGENy between BRCAness and non-BRCAness samples in the MUI (n=60) and TCGA (n=226) cohort.

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.

3.3 PARP inhibition activates the cGAS-STING pathway in vitro

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

Figure 3

Figure 3. Results from cell line experiments with olaparib treatment (A) Top up- and downregulated genes for the cell lines OVCAR3 and UWB1.289 under olaparib treatment when compared to DMSO control. (B) General distribution of up- and downregulated genes after olaparib treatment compared to the DMSO control in both cell lines as volcano plots. Red indicates significantly upregulated genes (FDR<0.1, log2-fold change>1), and blue indicates significantly downregulated genes (FDR<0.1, log2-fold change<-1). (C) Normalized enrichment score of pathways associated with activation of the cGAS STING pathway in BRCA1 mutated (cell lines) and BRCAness samples (cohorts) as well as olaparib-treated cell lines. (D) ClueGO network indicating overrepresented biological processes in the olaparib-treated UWB1.289 cell line. (E) Immunofluorescence staining of the DNA damage marker γH2AX in OVCAR3 and UWB1.289 cell lines with and without olaparib treatment. Comparing the different response to PARPi treatment between BRCA1 mutation and wild type BRCA1 (F) Immunofluorescence staining of cGAS, double stranded DNA (dsDNA) and STING in the OVCAR3 and UWB1.289 cell line comparing the difference between BRCA1 mutation and wild type BRCA1.

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.

3.4 BRCAness and immune subtype stratifies HGSOC patients

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

Figure 4

Figure 4. Profiles of immune parameters in the TCGA HGSOC cohort (n=226) (A) Heatmap of z scores of log2(TPM+1) expression of immune-related genes and fraction of tumor infiltrating immune cells assessed with quanTIseq and patient samples categorized by BRCAness, tumor-immune phenotype, molecular subtype and BRCA1/2 mutation. Furthermore, samples are stratified into four different tumor subtypes 1) BRCAness immune type samples (BRIT), which show an immunoreactive molecular subtype and an infiltrated tumor-immune phenotype, 2) noBRIT samples, which only have BRCAness but do not fulfill the other two requirements, 3) samples with an immune type (IMT) including an immunoreactive molecular subtype and an infiltrated tumor-immune phenotype but noBRCAness, and 4) remaining noIMT samples. (B) Distribution of estimated CD8 T cell fractions and estimated M2 macrophage fraction (from quanTIseq analyses) in the four different tumor subtypes. Benjamini-Hochberg adjusted p-values from pair-wise two-sided Dunn’s posthoc test are indicated.

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.

3.5 Tumor-associated macrophages inform therapy response

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

Figure 5

Figure 5. Single cell analysis of ovarian cancer adnexal samples from the MSK dataset (n=29) (A) UMAP showing the different cell types in of the ovarian cancer samples and which cells and cell types are associated with BRCAness samples. Distribution of major cell types in the BRCAness and noBRCAness are summarized as stacked bar plots. (B) UMAP plots of the myeloid cell compartment showing the association of macrophages with BRCAness cells and the expression of the macrophage marker gene C1QA and the TAM marker gene TREM2 especially in cell clusters associated with BRCAness. (C) Heatmap of expression of macrophage associated marker genes in the different cell types in the myeloid cell compartment. (D) Log₂ fold changes of differentially expressed genes between responder [R] and non-responder [NR] to PARPi-immune checkpoint inhibition combination therapy (niraparib and pembrolizumab) from the TOPACIO clinical trial (n=22) (p<0.05) (E) Distribution of expression visualized by UMAP and violin plots indicating in which (myeloid) cell types LYZ, LILRB, or ITGB2 are expressed (F) Dotplot indicating the distribution of expression and fraction of cells in various cell type for genes up-regulated (red) or down-regulated (blue) in responders vs non-responders to combination therapy as indicated in (D).

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.

3.6 Analyses of an independent cohort indicate vulnerability to combination immunotherapy

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

Figure 6

Figure 6. Expression profiles in the MUI cohort, immunohistochemistry validation, and vulnerability map (A) Heatmap of z-scores log2(TPM+1) expression of immune related genes and fraction of tumor infiltrating immune cells assessed with quanTIseq in all samples (n=60) from the MUI cohort categorized by BRCAness, tumor-immune phenotype, molecular subtype and BRCA1/2 mutation. (B) Immunohistochemistry images stained for CD8, CD4, CD163, γH2AX, and STING for three selected patients from the MUI cohort. Two BRIT samples one with a BRCA1 mutation and one without and one other sample without BRCAness, a deserted tumor-immune phenotype and a differentiated molecular subtype. (C) Vulnerability map showing the ratio between cytolytic activity CYT and C1QA (C2C) on the x-axis and the BRCAness score on the y-axis colored by the vulnerability score. The three selected samples were mapped to the vulnerability map based on their CYT to C1QA ratio (C2C) and BRCAness score.

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

4 Discussion

Here, we described how genomic instability in HGSOC affects the tumor immune environment and the consequences and vulnerabilities of combination immunotherapy combining PARP inhibitors with immune checkpoint inhibitors. A particular status in which patients respond well to PARP inhibitors and platinum-based chemotherapy is given when genes of the HRR pathway such as BRCA1 or BRCA2 are mutated. Genomic scars are consequences of HRD and are used to define an HRD score, often measured by established commercial assays, which allows the assignment of a responsive status beyond BRCA1 and BRCA2 mutations. The applicability and associated cutoff values for different assays and cancer types are under discussion, as the HRD algorithm has been used in clinical studies including different cancer types, such as breast cancer and ovarian cancer (46, 70, 71). Genomic scars are predictive but do not allow direct functional interpretation, whereas gene expression signatures could be an alternative in this regard. Very few approaches have associated gene expression with HRD status (18, 19, 62). Whereas a sixty-gene signature (18) and a two-gene signature (CXCL1, LY9) (19) have focused on microarray data, a recent approach using RNA sequencing data identified a 249-gene signature to predict HRD (62). We observed a number of overlaps with our 24-gene BRCA signature and a high concordance of signature scores in our training (TCGA) and validation (MUI) cohorts, indicating the reliability of our approach. This was underscored by comparison with mutational signature 3 (SigMA) in an independent cohort. The performance of the BRCAness classifier is reasonable, with AUC=0.91 (10-fold cross-validation) and positive predictive value for validation on both bulk RNA sequencing in the validation cohort (MUI) (PPV=0.86) and sample-wise single-cell RNA sequencing data (PPV=0.87).

There is evidence that BRCA1/2-mutated tumors exhibit significantly increased CD8+ TILs (22), although in breast cancer, differential modulation between BRCA1 and BRCA2 mutations in the tumor immune microenvironment has been found (72). We found a significant association between BRCAness and several immune-regulated signatures and evidence that several signaling pathways and processes known to modulate the immune system are activated by BRCA1 mutations or a BRCAness-related phenotype, such as JAK-STAT signaling or an interferon type I response, which are activated by free double-stranded DNA in the cytoplasm of tumor cells via the cGAS-STING pathway and affect dendritic cells (23–25). By expression and immunofluorescence analyses of ovarian cancer cell lines and by treatment with PARPi, we demonstrated that this axis is actually activated. Notably, the STAT3 pathway, which is activated by PARP inhibition, may, however, mediate treatment resistance by promoting the polarization of protumor TAMs, which could be overcome by STING agonism (26). STING, CSF1R, SREBP-1, and VEGFA might also be targets to overcome resistance to PARPi-immunotherapy combinations (73). The upregulation of many chemokines and chemokine receptors indicates that BRCAness tumors are actively involved in immune cell attraction and interaction. For example, CCL5 produced by tumor cells or CXCL9 and CXCL10 also expressed by tumor-resident myeloid cells determine effector T-cell recruitment to the tumor microenvironment (74). We detected significant upregulation of CCL5 and CXCL10 by PARP inhibition, which was also identified as a downstream target of STING (24). Another interesting chemokine that is strongly upregulated in cancer cell lines, particularly by olaparib treatment, is CCL20. CCL20 could be associated with cancer metastasis and progression by interacting with its cognate receptor CCR6 in an ovarian cancer mouse model. However, the higher expression in the myeloid cell compartment, as evident from single-cell analyses, overlies the intrinsic tumor effect.

One of our basic hypotheses was that samples with BRCAness respond better to PARPi therapy and that hot tumors with an activated immune milieu respond better to immune checkpoint inhibition, as has been shown, for example, in melanoma for the activated IFNG pathway (75). However, when we compared the BRCAness immune type (BRIT) with other samples, we observed by using deconvolution approaches that suppressive cell types such as M2 macrophages, MDSCs, and Tregs were more abundant. In particular, TAMs could be a major factor together with low mutational burden, abnormal neovascularization, altered metabolism, and failure to reverse T-cell exhaustion for the limited immunotherapy response in ovarian cancer (76). By using single-cell RNA sequencing data analyses of adnexal cancer tissue, we demonstrated that myeloid cells are the most abundant immune cells, and the majority were characterized as TAMs and rather polarized towards M2-like macrophages compared to classical (M1-like) macrophages, although these better described as continuum of different stages than isolated cell types. A majority of these TAMs are immune suppressive, as indicated by TREM2 expression. TREM2 is a promising therapeutic target for TAM depletion (68). Inhibition of TREM2 has been shown to improve the anti-PD1 response in various mouse models and is currently being investigated in a clinical trial (64, 77). Another recent study underscored the role of TAMs and demonstrated that specifically, the Siglec-9-positive TAM subset is associated with an immune-suppressive phenotype and adverse prognosis in HGSOC patients (78).

Interestingly, a previous work using cyclic immunofluorescence highlighted the role of exhausted T cells in the response to niraparib/pembrolizumab. In responders, particularly in extreme responders, frequent proximity between exhausted T cells and PD-L1+ (CD163+, IBA1+, CD11b+) TAMs was observed (9). Noticeably, based on the selected marker expression, we observed an overlap with the CD169/SIGLEC-1 macrophage cluster (Supplementary Figure S10). In addition, in patients who responded to this combination therapy, we identified a number of downregulated genes that were also highly expressed in TAMs, such as LYZ, LILRB4, and ITGB2. Whereas lysozyme (LYZ) is an antimicrobial ligand and is involved in central macrophage function and is therefore nonspecifically and highly expressed, LILRB4 is an immune checkpoint on myeloid cells, indicating a more regulatory role. High expression of the integrin ITGB2 was previously shown to be associated with poor survival outcome (79), underscoring that high expression in TAMs is crucial. In contrast, ITGB2 is also associated with CD8+ T cells, as it encodes the beta chain of the LFA-1 protein, which has been shown to be essential in the assembly of the immune synapse or to influence lymphocyte extravasation and T-cell recruitment to the tumor and is regulated by GDF-15 (80).

Because stratification of patients based on gene expression in our validation cohort was very similar to the analysis on the TCGA cohort, we set out to adapt our hypothesis and also include elements of the suppressive environment. Already, it was shown that regulatory T cells (Tregs) are an important component of the suppressive milieu and are associated with unfavorable survival outcomes in ovarian cancer (81, 82). We performed immunohistochemistry analyses using FOXP3 and CD163 antibodies in the validation cohort and found very pronounced macrophage infiltration (CD163) but hardly Treg infiltration (FOXP3) into the tumor site in some samples. The results of the single-cell RNA sequencing analyses and the fact that various TAM marker genes were associated with poorer overall survival also suggest that TAMs play a more dominant role in ovarian cancer.

While infiltration of various cell types from the adaptive immune system (83) and other markers, such as tumor mutational burden (TMB) (84) or IFNG signature (75), have been associated with good prognosis and immunotherapy response in various cancer types, the suppressive immune environment with tumor-supportive CD68+CD163+ macrophages is becoming more important (66). Accordingly, a signature of the immune activation ratio of CD8A/C1QA has been found to be prognostic and predictive for immunotherapy response (66). We considered the mean PRF1 and GZMB expression as a proxy for cytolytic activity (45) as predominantly exerted by cytotoxic T lymphocytes. The specific expression pattern of C1QA on TAMs was comparable to that of TREM2 but at a much higher level. Therefore, we also used the member of the complement system C1QA as a surrogate for TAMs and the suppressive tumor immune environment and finally built a ratio of cytolytic activity (CYT) to the expression of C1QA (C2C), indicating the pro- and antitumoral balance of the immune environment. Finally, to build a predictive algorithm for combination therapy response, we included both C2C on the one hand and BRCAness on the other hand into one model. Since HRD measured with companion diagnostic tests is not able to predict all PARPi responders, as shown in several clinical trials, and since PARPi treatment can activate a number of immune-related pathways even in situations with proficient HRR, which is also underlined by our in vitro analyses, this model is considered to be relevant for combination immunotherapy.

Our studies have some limitations in that the training and validation patient cohorts were retrospective studies, and RNA sequencing was performed at a later time point. Additionally, only a limited number of patients who received combination therapy could be included; therefore, the conclusion about the predictive power for the treatment is limited and requires further validation in larger cohorts. One component that was not considered in this study is malignant ascites, which has been shown to contain various cell types, such as macrophages, many soluble factors and cytokines, that influence the protumorigenic phenotype and promote metastatic spread of HGSOC through transcoelomic dissemination (85).

In conclusion, our approach using RNA sequencing data to comprehensively characterize both genome instability and the tumor immune environment enabled us to stratify HGSOC patients. Further analyses indicate that suppressive TAMs in the tumor immune microenvironment may play an essential role in understanding why receiving (combination) immunotherapy shows limited efficacy in ovarian cancer. Based on multiple datasets, we have developed a methodology and corresponding easy-to-use diagnostic application (https://ovrseq.icbi.at) and an R package OvRSeq that uses RNA sequencing data not only to comprehensively characterize newly diagnosed HGSOC patients but also to inform therapy response. Ultimately, this approach will be very useful to obtain comprehensive information about the phenotype of a tumor sample, support clinical decisions, and stimulate further research.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material. RNA sequencing data can be found at Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE237361). The raw RNA sequencing data from the MUI cohort cannot be made publicly available due to data protection restrictions and are available on request from the corresponding author.

Ethics statement

The studies involving humans were approved by Ethics Committee of the Medical University of Innsbruck (reference number: 1189/2019). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

RG: Data curation, Formal Analysis, Methodology, Visualization, Writing – original draft, Writing – review & editing. LM: Investigation, Writing – original draft, Writing – review & editing. PM: Software, Writing – original draft, Writing – review & editing. GF: Formal Analysis, Writing – original draft, Writing – review & editing. SS: Investigation, Validation, Writing – original draft, Writing – review & editing. AZ: Writing – original draft, Writing – review & editing. CM: Resources, Supervision, Writing – original draft, Writing – review & editing. HF: Investigation, Resources, Validation, Writing – original draft, Writing – review & editing. HH: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This research was funded in whole, or in part, by the Anniversary Fund of the National Bank of Austria (OeNB) (grant number 18279 to HH).

Acknowledgments

A version of the manuscript was deposited at the medRxiv preprint server (86). Open Access Funding provided by the Medical University of Innsbruck.

Conflict of interest

SS was employed by Innpath GmbH. AZ reports consulting fees from Amgen, Astra Zeneca, GSK, MSD, Novartis, PharmaMar, Roche, Seagen; honoraria from Amgen, Astra Zeneca, GSK, MSD, Novartis, PharmaMar, Roche, Seagen; travel expenses from Astra Zeneca, Gilead, Roche; participation on advisory boards from Amgen, Astra Zeneca, GSK, MSD, Novartis, Pfizer, PharmaMar, Roche, Seagen. CM reports consulting fees and honoraria from Roche, Novartis, Amgen, MSD, PharmaMar, Astra Zeneca, GSK, Seagen; travel expenses from Roche, Astra Zeneca; participation on advisory boards from Roche, Novartis, Amgen, MSD, Astra Zeneca, Pfizer, PharmaMar, GSK, Seagen. HH has received research funding via Catalym and Secarna.

The remaining authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

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/fimmu.2024.1489235/full#supplementary-material

References

1. Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics, 2023. CA A Cancer J Clin. (2023) 73:17–48. doi: 10.3322/caac.21763

PubMed Abstract | Crossref Full Text | Google Scholar

2. Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, et al. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science. (2017) 357:409–13. doi: 10.1126/science.aan6733

PubMed Abstract | Crossref Full Text | Google Scholar

3. Disis ML, Taylor MH, Kelly K, Beck JT, Gordon M, Moore KM, et al. Efficacy and safety of avelumab for patients with recurrent or refractory ovarian cancer: phase 1b results from the JAVELIN solid tumor trial. JAMA Oncol. (2019) 5:393. doi: 10.1001/jamaoncol.2018.6258

PubMed Abstract | Crossref Full Text | Google Scholar

4. Varga A, Piha-Paul S, Ott PA, Mehnert JM, Berton-Rigaud D, Morosky A, et al. Pembrolizumab in patients with programmed death ligand 1–positive advanced ovarian cancer: Analysis of KEYNOTE-028. Gynecol Oncol. (2019) 152:243–50. doi: 10.1016/j.ygyno.2018.11.017

PubMed Abstract | Crossref Full Text | Google Scholar

5. Konstantinopoulos PA, Cannistra SA. Immune checkpoint inhibitors in ovarian cancer: can we bridge the gap between IMagynation and reality? JCO. (2021) 39:1833–8. doi: 10.1200/JCO.21.00571

PubMed Abstract | Crossref Full Text | Google Scholar

6. Matulonis UA, Shapira-Frommer R, Santin AD, Lisyanskaya AS, Pignata S, Vergote I, et al. Antitumor activity and safety of pembrolizumab in patients with advanced recurrent ovarian cancer: results from the phase II KEYNOTE-100 study. Ann Oncol. (2019) 30:1080–7. doi: 10.1093/annonc/mdz135

PubMed Abstract | Crossref Full Text | Google Scholar

7. Ray-Coquard I, Pautier P, Pignata S, Pérol D, González-Martín A, Berger R, et al. Olaparib plus bevacizumab as first-line maintenance in ovarian cancer. N Engl J Med. (2019) 381:2416–28. doi: 10.1056/NEJMoa1911361

PubMed Abstract | Crossref Full Text | Google Scholar

8. Musacchio L, Cicala CM, Camarda F, Ghizzoni V, Giudice E, Carbone MV, et al. Combining PARP inhibition and immune checkpoint blockade in ovarian cancer patients: a new perspective on the horizon? ESMO Open. (2022) 7:100536. doi: 10.1016/j.esmoop.2022.100536

PubMed Abstract | Crossref Full Text | Google Scholar

9. Färkkilä A, Gulhan DC, Casado J, Jacobson CA, Nguyen H, Kochupurakkal B, et al. Immunogenomic profiling determines responses to combined PARP and PD-1 inhibition in ovarian cancer. Nat Commun. (2020) 11:1459. doi: 10.1038/s41467-020-15315-8

PubMed Abstract | Crossref Full Text | Google Scholar

10. Lampert EJ, Zimmer A, Padget M, Cimino-Mathews A, Nair JR, Liu Y, et al. Combination of PARP inhibitor olaparib, and PD-L1 inhibitor durvalumab, in recurrent ovarian cancer: a proof-of-concept phase II study. Clin Cancer Res. (2020) 26:4268–79. doi: 10.1158/1078-0432.CCR-20-0056

PubMed Abstract | Crossref Full Text | Google Scholar

11. Drew Y, Kim J-W, Penson RT, O’Malley DM, Parkinson C, Roxburgh P, et al. Olaparib plus Durvalumab, with or without Bevacizumab, as Treatment in PARP Inhibitor-Naïve Platinum-Sensitive Relapsed Ovarian Cancer: A Phase II Multi-Cohort Study. Clin Cancer Res. (2023) 30(1):50–62. doi: 10.1158/1078-0432.CCR-23-2249

PubMed Abstract | Crossref Full Text | Google Scholar

12. Gonzalez Martin A, Rubio Perez MJ, Heitz F, Christensen RD, Colombo N, Van Gorp T, et al. LBA37 Atezolizumab (atezo) combined with platinum-based chemotherapy (CT) and maintenance niraparib for recurrent ovarian cancer (rOC) with a platinum-free interval (TFIp) >6 months: Primary analysis of the double-blind placebo (pbo)-controlled ENGOT-Ov41/GEICO 69-O/ANITA phase III trial. Ann Oncol. (2023) 34:S1278–9. doi: 10.1016/j.annonc.2023.10.031

Crossref Full Text | Google Scholar

13. Lord CJ, Ashworth A. PARP inhibitors: Synthetic lethality in the clinic. Science. (2017) 355:1152–8. doi: 10.1126/science.aam7344

PubMed Abstract | Crossref Full Text | Google Scholar

14. Konstantinopoulos PA, Ceccaldi R, Shapiro GI, D’Andrea AD. Homologous recombination deficiency: exploiting the fundamental vulnerability of ovarian cancer. Cancer Discovery. (2015) 5:1137–54. doi: 10.1158/2159-8290.CD-15-0714

PubMed Abstract | Crossref Full Text | Google Scholar

15. Cancer Genome Atlas Research Network. Integrated genomic analyses of ovarian carcinoma. Nature. (2011) 474:609–15. doi: 10.1038/nature10166

PubMed Abstract | Crossref Full Text | Google Scholar

16. Stewart MD, Merino Vega D, Arend RC, Baden JF, Barbash O, Beaubier N, et al. Homologous recombination deficiency: concepts, definitions, and assays. Oncologist. (2022) 27:167–74. doi: 10.1093/oncolo/oyab053

PubMed Abstract | Crossref Full Text | Google Scholar

17. Takamatsu S, Brown JB, Yamaguchi K, Hamanishi J, Yamanoi K, Takaya H, et al. Utility of homologous recombination deficiency biomarkers across cancer types. JCO Precis Oncol. (2022) 6:e2200085. doi: 10.1200/PO.22.00085

PubMed Abstract | Crossref Full Text | Google Scholar

18. Konstantinopoulos PA, Spentzos D, Karlan BY, Taniguchi T, Fountzilas E, Francoeur N, et al. Gene expression profile of BRCA ness that correlates with responsiveness to chemotherapy and with outcome in patients with epithelial ovarian cancer. JCO. (2010) 28:3555–61. doi: 10.1200/JCO.2009.27.5719

PubMed Abstract | Crossref Full Text | Google Scholar

19. Chen T, Yu T, Zhuang S, Geng Y, Xue J, Wang J, et al. Upregulation of CXCL1 and LY9 contributes to BRCAness in ovarian cancer and mediates response to PARPi and immune checkpoint blockade. Br J Cancer. (2022) 127:916–26. doi: 10.1038/s41416-022-01836-0

PubMed Abstract | Crossref Full Text | Google Scholar

20. Kraya AA, Maxwell KN, Wubbenhorst B, Wenz BM, Pluta J, Rech AJ, et al. Genomic signatures predict the immunogenicity of BRCA-deficient breast cancer. Clin Cancer Res. (2019) 25:4363–74. doi: 10.1158/1078-0432.CCR-18-0468

PubMed Abstract | Crossref Full Text | Google Scholar

21. Hoppe MM, Sundar R, Tan DSP, Jeyasekharan AD. Biomarkers for homologous recombination deficiency in cancer. J Natl Cancer Inst. (2018) 110:704–13. doi: 10.1093/jnci/djy085

PubMed Abstract | Crossref Full Text | Google Scholar

22. Strickland KC, Howitt BE, Shukla SA, Rodig S, Ritterhouse LL, Liu JF, et al. Association and prognostic significance of BRCA1/2-mutation status with neoantigen load, number of tumor-infiltrating lymphocytes and expression of PD-1/PD-L1 in high grade serous ovarian cancer. Oncotarget. (2016) 7:13587–98. doi: 10.18632/oncotarget.7277

PubMed Abstract | Crossref Full Text | Google Scholar

23. Bruand M, Barras D, Mina M, Ghisoni E, Morotti M, Lanitis E, et al. Cell-autonomous inflammation of BRCA1-deficient ovarian cancers drives both tumor-intrinsic immunoreactivity and immune resistance via STING. Cell Rep. (2021) 36:109412. doi: 10.1016/j.celrep.2021.109412

PubMed Abstract | Crossref Full Text | Google Scholar

24. Ding L, Kim H-J, Wang Q, Kearns M, Jiang T, Ohlson CE, et al. PARP inhibition elicits STING-dependent antitumor immunity in brca1-deficient ovarian cancer. Cell Rep. (2018) 25:2972–2980.e5. doi: 10.1016/j.celrep.2018.11.054

PubMed Abstract | Crossref Full Text | Google Scholar

25. Zheng J, Mo J, Zhu T, Zhuo W, Yi Y, Hu S, et al. Comprehensive elaboration of the cGAS-STING signaling axis in cancer development and immunotherapy. Mol Cancer. (2020) 19:133. doi: 10.1186/s12943-020-01250-1

PubMed Abstract | Crossref Full Text | Google Scholar

26. Ding L, Wang Q, Martincuks A, Kearns MJ, Jiang T, Lin Z, et al. STING agonism overcomes STAT3-mediated immunosuppression and adaptive resistance to PARP inhibition in ovarian cancer. J Immunother Cancer. (2023) 11:e005627. doi: 10.1136/jitc-2022-005627

PubMed Abstract | Crossref Full Text | Google Scholar

27. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell. (2018) 173:400–416.e11. doi: 10.1016/j.cell.2018.02.052

PubMed Abstract | Crossref Full Text | Google Scholar

28. Desbois M, Udyavar AR, Ryner L, Kozlowski C, Guan Y, Dürrbaum M, et al. Integrated digital pathology and transcriptome analysis identifies molecular mediators of T-cell exclusion in ovarian cancer. Nat Commun. (2020) 11:5583. doi: 10.1038/s41467-020-19408-2

PubMed Abstract | Crossref Full Text | Google Scholar

29. Vázquez-García I, Uhlitz F, Ceglia N, Lim JLP, Wu M, Mohibullah N, et al. Ovarian cancer mutational processes drive site-specific immune evasion. Nature. (2022) 612:778–86. doi: 10.1038/s41586-022-05496-1

PubMed Abstract | Crossref Full Text | Google Scholar

30. McDermott JE, Arshad OA, Petyuk VA, Fu Y, Gritsenko MA, Clauss TR, et al. Proteogenomic characterization of ovarian HGSC implicates mitotic kinases, replication stress in observed chromosomal instability. Cell Rep Med. (2020) 1:100004. doi: 10.1016/j.xcrm.2020.100004

PubMed Abstract | Crossref Full Text | Google Scholar

31. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | Crossref Full Text | Google Scholar

32. Van der Auwera GA, Carneiro MO, Hartl C, Poplin R, Del Angel G, Levy-Moonshine A, et al. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinf. (2013) 43:11.10.1–11.10.33. doi: 10.1002/0471250953.bi1110s43

PubMed Abstract | Crossref Full Text | Google Scholar

33. Larson DE, Harris CC, Chen K, Koboldt DC, Abbott TE, Dooling DJ, et al. SomaticSniper: identification of somatic point mutations in whole genome sequencing data. Bioinformatics. (2012) 28:311–7. doi: 10.1093/bioinformatics/btr665

PubMed Abstract | Crossref Full Text | Google Scholar

34. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. (2012) 22:568–76. doi: 10.1101/gr.129684.111

PubMed Abstract | Crossref Full Text | Google Scholar

35. Kim S, Scheffler K, Halpern AL, Bekritsky MA, Noh E, Källberg M, et al. Strelka2: fast and accurate calling of germline and somatic variants. Nat Methods. (2018) 15:591–4. doi: 10.1038/s41592-018-0051-x

PubMed Abstract | Crossref Full Text | Google Scholar

36. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The ensembl variant effect predictor. Genome Biol. (2016) 17:122. doi: 10.1186/s13059-016-0974-4

PubMed Abstract | Crossref Full Text | Google Scholar

37. Jurtz V, Paul S, Andreatta M, Marcatili P, Peters B, Nielsen M. NetMHCpan-4.0: improved peptide-MHC class I interaction predictions integrating eluted ligand and peptide binding affinity data. J Immunol. (2017) 199:3360–8. doi: 10.4049/jimmunol.1700893

PubMed Abstract | Crossref Full Text | Google Scholar

38. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | Crossref Full Text | Google Scholar

39. 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 USA. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | Crossref Full Text | Google Scholar

40. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. (2009) 25:1091–3. doi: 10.1093/bioinformatics/btp101

PubMed Abstract | Crossref Full Text | Google Scholar

41. Schubert M, Klinger B, Klünemann M, Sieber A, Uhlitz F, Sauer S, et al. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nat Commun. (2018) 9:20. doi: 10.1038/s41467-017-02391-6

PubMed Abstract | Crossref Full Text | Google Scholar

42. 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:34. doi: 10.1186/s13073-019-0638-6

PubMed Abstract | Crossref Full Text | Google Scholar

43. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | Crossref Full Text | Google Scholar

44. Chen GM, Kannan L, Geistlinger L, Kofia V, Safikhani Z, Gendoo DMA, et al. Consensus on molecular subtypes of high-grade serous ovarian carcinoma. Clin Cancer Res. (2018) 24:5037–47. doi: 10.1158/1078-0432.CCR-18-0784

PubMed Abstract | Crossref Full Text | Google Scholar

45. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | Crossref Full Text | Google Scholar

46. Telli ML, Timms KM, Reid J, Hennessy B, Mills GB, Jensen KC, et al. Homologous recombination deficiency (HRD) score predicts response to platinum-containing neoadjuvant chemotherapy in patients with triple-negative breast cancer. Clin Cancer Res. (2016) 22:3764–73. doi: 10.1158/1078-0432.CCR-15-2477

PubMed Abstract | Crossref Full Text | Google Scholar

47. Alexandrov LB, Kim J, Haradhvala NJ, Huang MN, Tian Ng AW, Wu Y, et al. The repertoire of mutational signatures in human cancer. Nature. (2020) 578:94–101. doi: 10.1038/s41586-020-1943-3

PubMed Abstract | Crossref Full Text | Google Scholar

48. Takaya H, Nakai H, Takamatsu S, Mandai M, Matsumura N. Homologous recombination deficiency status-based classification of high-grade serous ovarian carcinoma. Sci Rep. (2020) 10:2757. doi: 10.1038/s41598-020-59671-3

PubMed Abstract | Crossref Full Text | Google Scholar

49. Abkevich V, Timms KM, Hennessy BT, Potter J, Carey MS, Meyer LA, et al. Patterns of genomic loss of heterozygosity predict homologous recombination repair defects in epithelial ovarian cancer. Br J Cancer. (2012) 107:1776–82. doi: 10.1038/bjc.2012.451

PubMed Abstract | Crossref Full Text | Google Scholar

50. Birkbak NJ, Wang ZC, Kim J-Y, Eklund AC, Li Q, Tian R, et al. Telomeric allelic imbalance indicates defective DNA repair and sensitivity to DNA-damaging agents. Cancer Discovery. (2012) 2:366–75. doi: 10.1158/2159-8290.CD-11-0206

PubMed Abstract | Crossref Full Text | Google Scholar

51. Popova T, Manié E, Rieunier G, Caux-Moncoutier V, Tirapo C, Dubois T, et al. Ploidy and large-scale genomic instability consistently identify basal-like breast carcinomas with BRCA1/2 inactivation. Cancer Res. (2012) 72:5454–62. doi: 10.1158/0008-5472.CAN-12-1470

PubMed Abstract | Crossref Full Text | Google Scholar

52. Sztupinszki Z, Diossy M, Krzystanek M, Reiniger L, Csabai I, Favero F, et al. Migrating the SNP array-based homologous recombination deficiency measures to next generation sequencing data of breast cancer. NPJ Breast Cancer. (2018) 4:16. doi: 10.1038/s41523-018-0066-6

PubMed Abstract | Crossref Full Text | Google Scholar

53. Favero F, Joshi T, Marquard AM, Birkbak NJ, Krzystanek M, Li Q, et al. Sequenza: allele-specific copy number and mutation profiles from tumor sequencing data. Ann Oncol. (2015) 26:64–70. doi: 10.1093/annonc/mdu479

PubMed Abstract | Crossref Full Text | Google Scholar

54. Manders F, Brandsma AM, de Kanter J, Verheul M, Oka R, van Roosmalen MJ, et al. MutationalPatterns: the one stop shop for the analysis of mutational processes. BMC Genomics. (2022) 23:134. doi: 10.1186/s12864-022-08357-3

PubMed Abstract | Crossref Full Text | Google Scholar

55. Gulhan DC, Lee JJ-K, Melloni GEM, Cortés-Ciriano I, Park PJ. Detecting the mutational signature of homologous recombination deficiency in clinical samples. Nat Genet. (2019) 51:912–9. doi: 10.1038/s41588-019-0390-2

PubMed Abstract | Crossref Full Text | Google Scholar

56. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. (2018) 19:15. doi: 10.1186/s13059-017-1382-0

PubMed Abstract | Crossref Full Text | Google Scholar

57. Gayoso A, Lopez R, Xing G, Boyeau P, Valiollah Pour Amiri V, Hong J, et al. A Python library for probabilistic analysis of single-cell omics data. Nat Biotechnol. (2022) 40:163–6. doi: 10.1038/s41587-021-01206-w

PubMed Abstract | Crossref Full Text | Google Scholar

58. Franzén O, Gan L-M, Björkegren JLM. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford). (2019) 2019). doi: 10.1093/database/baz046

PubMed Abstract | Crossref Full Text | Google Scholar

59. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. (2020) 15:1484–506. doi: 10.1038/s41596-020-0292-x

PubMed Abstract | Crossref Full Text | Google Scholar

60. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

61. Lord CJ, Ashworth A. BRCAness revisited. Nat Rev Cancer. (2016) 16:110–20. doi: 10.1038/nrc.2015.21

PubMed Abstract | Crossref Full Text | Google Scholar

62. Takamatsu S, Yoshihara K, Baba T, Shimada M, Yoshida H, Kajiyama H, et al. Prognostic relevance of HRDness gene expression signature in ovarian high-grade serous carcinoma; JGOG3025-TR2 study. Br J Cancer. (2023) 128:1095–104. doi: 10.1038/s41416-022-02122-9

PubMed Abstract | Crossref Full Text | Google Scholar

63. Liu Z, Wu H, Deng J, Wang H, Wang Z, Yang A, et al. Molecular classification and immunologic characteristics of immunoreactive high-grade serous ovarian cancer. J Cell Mol Med. (2020) 24:8103–14. doi: 10.1111/jcmm.15441

Crossref Full Text | Google Scholar

64. Binnewies M, Pollack JL, Rudolph J, Dash S, Abushawish M, Lee T, et al. Targeting TREM2 on tumor-associated macrophages enhances immunotherapy. Cell Rep. (2021) 37:109844. doi: 10.1016/j.celrep.2021.109844

PubMed Abstract | Crossref Full Text | Google Scholar

65. Cassetta L, Fragkogianni S, Sims AH, Swierczak A, Forrester LM, Zhang H, et al. Human tumor-associated macrophage and monocyte transcriptional landscapes reveal cancer-specific reprogramming, biomarkers, and therapeutic targets. Cancer Cell. (2019) 35:588–602.e10. doi: 10.1016/j.ccell.2019.02.009

PubMed Abstract | Crossref Full Text | Google Scholar

66. Mezheyeuski A, Backman M, Mattsson J, Martín-Bernabé A, Larsson C, Hrynchyk I, et al. An immune score reflecting pro- and anti-tumoural balance of tumour microenvironment has major prognostic impact and predicts immunotherapy response in solid cancers. EBioMedicine. (2023) 88:104452. doi: 10.1016/j.ebiom.2023.104452

PubMed Abstract | Crossref Full Text | Google Scholar

67. Mullen MM, Lomonosova E, Toboni MD, Oplt A, Cybulla E, Blachut B, et al. GAS6/AXL inhibition enhances ovarian cancer sensitivity to chemotherapy and PARP inhibition through increased DNA damage and enhanced replication stress. Mol Cancer Res. (2022) 20:265–79. doi: 10.1158/1541-7786.MCR-21-0302

PubMed Abstract | Crossref Full Text | Google Scholar

68. Truxova I, Cibula D, Spisek R, Fucikova J. Targeting tumor-associated macrophages for successful immunotherapy of ovarian carcinoma. J Immunother Cancer. (2023) 11:e005968. doi: 10.1136/jitc-2022-005968

PubMed Abstract | Crossref Full Text | Google Scholar

69. Wieser V, Tsibulak I, Reimer DU, Zeimet AG, Fiegl H, Hackl H, et al. An angiogenic tumor phenotype predicts poor prognosis in ovarian cancer. Gynecol Oncol. (2023) 170:290–9. doi: 10.1016/j.ygyno.2023.01.034

PubMed Abstract | Crossref Full Text | Google Scholar

70. Rempel E, Kluck K, Beck S, Ourailidis I, Kazdal D, Neumann O, et al. Pan-cancer analysis of genomic scar patterns caused by homologous repair deficiency (HRD). NPJ Precis Oncol. (2022) 6:36. doi: 10.1038/s41698-022-00276-6

PubMed Abstract | Crossref Full Text | Google Scholar

71. Perez-Villatoro F, Oikkonen J, Casado J, Chernenko A, Gulhan DC, Tumiati M, et al. Optimized detection of homologous recombination deficiency improves the prediction of clinical outcomes in cancer. NPJ Precis Oncol. (2022) 6:96. doi: 10.1038/s41698-022-00339-8

PubMed Abstract | Crossref Full Text | Google Scholar

72. Samstein RM, Krishna C, Ma X, Pei X, Lee K-W, Makarov V, et al. Mutations in BRCA1 and BRCA2 differentially affect the tumor microenvironment and response to checkpoint blockade immunotherapy. Nat Cancer. (2020) 1:1188–203. doi: 10.1038/s43018-020-00139-8

PubMed Abstract | Crossref Full Text | Google Scholar

73. Konstantinopoulos PA, Matulonis UA. Clinical and translational advances in ovarian cancer therapy. Nat Cancer. (2023) 4:1239–57. doi: 10.1038/s43018-023-00617-9

PubMed Abstract | Crossref Full Text | Google Scholar

74. Kandalaft LE, Dangaj Laniti D, Coukos G. Immunobiology of high-grade serous ovarian cancer: lessons for clinical translation. Nat Rev Cancer. (2022) 22:640–56. doi: 10.1038/s41568-022-00503-z

PubMed Abstract | Crossref Full Text | Google Scholar

75. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, et al. IFN-γ–related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. (2017) 127:2930–40. doi: 10.1172/JCI91190

PubMed Abstract | Crossref Full Text | Google Scholar

76. Johnson RL, Cummings M, Thangavelu A, Theophilou G, de Jong D, Orsi NM. Barriers to immunotherapy in ovarian cancer: metabolic, genomic, and immune perturbations in the tumour microenvironment. Cancers (Basel). (2021) 13:6231. doi: 10.3390/cancers13246231

PubMed Abstract | Crossref Full Text | Google Scholar

77. Molgora M, Esaulova E, Vermi W, Hou J, Chen Y, Luo J, et al. TREM2 modulation remodels the tumor myeloid landscape enhancing anti-PD-1 immunotherapy. Cell. (2020) 182:886–900.e17. doi: 10.1016/j.cell.2020.07.013

PubMed Abstract | Crossref Full Text | Google Scholar

78. Wang Y, He M, Zhang C, Cao K, Zhang G, Yang M, et al. Siglec-9+ tumor-associated macrophages delineate an immunosuppressive subset with therapeutic vulnerability in patients with high-grade serous ovarian cancer. J Immunother Cancer. (2023) 11:e007099. doi: 10.1136/jitc-2023-007099

PubMed Abstract | Crossref Full Text | Google Scholar

79. Li C, Deng T, Cao J, Zhou Y, Luo X, Feng Y, et al. Identifying ITGB2 as a potential prognostic biomarker in ovarian cancer. Diagnostics (Basel). (2023) 13:1169. doi: 10.3390/diagnostics13061169

PubMed Abstract | Crossref Full Text | Google Scholar

80. Haake M, Haack B, Schäfer T, Harter PN, Mattavelli G, Eiring P, et al. Tumor-derived GDF-15 blocks LFA-1 dependent T cell recruitment and suppresses responses to anti-PD-1 treatment. Nat Commun. (2023) 14:4253. doi: 10.1038/s41467-023-39817-3

PubMed Abstract | Crossref Full Text | Google Scholar

81. Curiel TJ, Coukos G, Zou L, Alvarez X, Cheng P, Mottram P, et al. Specific recruitment of regulatory T cells in ovarian carcinoma fosters immune privilege and predicts reduced survival. Nat Med. (2004) 10:942–9. doi: 10.1038/nm1093

PubMed Abstract | Crossref Full Text | Google Scholar

82. Wolf D, Wolf AM, Rumpold H, Fiegl H, Zeimet AG, Muller-Holzner E, et al. The expression of the regulatory T cell-specific forkhead box transcription factor FoxP3 is associated with poor prognosis in ovarian cancer. Clin Cancer Res. (2005) 11:8326–31. doi: 10.1158/1078-0432.CCR-05-1244

PubMed Abstract | Crossref Full Text | Google Scholar

83. Galon J, Costes A, Sanchez-Cabo F, Kirilovsky A, Mlecnik B, Lagorce-Pagès C, et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science. (2006) 313:1960–4. doi: 10.1126/science.1129139

PubMed Abstract | Crossref Full Text | Google Scholar

84. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol. (2019) 30:44–56. doi: 10.1093/annonc/mdy495

PubMed Abstract | Crossref Full Text | Google Scholar

85. Almeida-Nunes DL, Mendes-Frias A, Silvestre R, Dinis-Oliveira RJ, Ricardo S. Immune tumor microenvironment in ovarian cancer ascites. Int J Mol Sci. (2022) 23:10692. doi: 10.3390/ijms231810692

PubMed Abstract | Crossref Full Text | Google Scholar

86. Gronauer R, Madersbacher L, Monfort-Lanzas P, Floriani G, Sprung S, Zeimet AG, et al. Integrated immunogenomic analyses of high-grade serous ovarian cancer reveal vulnerability to combination immunotherapy. (2024). doi: 10.1101/2024.01.09.24301038

Crossref Full Text | Google Scholar