To spatially map HGSC in the setting of metastatic disease, we applied in situ imaging with high-plex RNA detection at single-cell resolution to 130 HGSC tumors from a total of 94 patients to generate 202 tissue profiles, yielding a total of 2,598,277 high-quality single-cell spatial transcriptomes (Fig. 1a and Supplementary Table 1a,b). Tumor tissues were obtained from the adnexa (ovaries/fallopian tube, n = 84), and/or omentum (n = 46), with 36 patient-matched pairs of adnexal and omental tumors. All tumor tissue profiles were obtained from debulking surgeries in either the treatment-naive or the neoadjuvant chemotherapy-treated setting (Fig. 1b and Supplementary Table 2a,b). Associated patient clinical data including treatment history (for example, PARP inhibitor, bevacizumab and immune checkpoint blockade (ICB)) and survival outcomes are also available (Fig. 1a,b, Supplementary Table 2a,b and Methods). Eight patients in this cohort received ICB treatment, but in all cases ICB treatment was after sample collection. For 40 patients, we also obtained DNA sequencing data spanning a 648-gene panel (Fig. 1b, Supplementary Table 1b and Methods), focused on actionable single-nucleotide variations, somatic CNAs, chromosomal rearrangements and tumor mutational burden (TMB), providing a basis to link tissue structure and somatic genetic aberrations.

a, Cohort and data overview. Created with BioRender.com. b, Summary of clinical history, tumor genetics and ST profiles per patient. Each column represents 1 of 94 patients. NACT, neoadjuvant chemotherapy; Beva, bevacizumab; PARPi, PARP inhibitor. c, Uniform manifold approximation and projection (UMAP) embedding of high-confidence spatial single-cell transcriptomes from the different datasets. n denotes number of cells within each cell-type annotation in the Discovery dataset. d, Representative tumor tissue ST images (11 of 202) with cells plotted in situ and colored based on cell-type annotations. e, Co-embedding of spatial single-cell transcriptomes from this study with six publicly available HGSC scRNA-seq datasets17,24,25,26,27,28,29. f, Cell-type composition (y axis) per tissue profile (x axis) from this study and in six publicly available scRNA-seq HGSC datasets17,24,25,26,27,28,29. g, Pairwise colocalization analysis: the number of tissue profiles (x axis) where each pair of cell types (y axis) shows significantly (BH FDR < 0.05, hypergeometric test) higher (red), lower (blue) or expected (gray) colocalization frequencies compared to those expected by random. h, log2 colocalization quotient (CLQ) of T/NK cells with fibroblasts (CLQT/NK cell→fibroblast, blue) and T/NK cells with malignant cells (CLQT/NK cell→malignant, green, x axis) in ST tissue profiles from the Discovery dataset (n = 87 CLQ pairs, P = 4.31 × 10−10, paired Wilcoxon rank-sum test). Light gray lines connect paired values in each ST tissue core (black dots). In the box plots, the middle line denotes the median, box edges indicate the 25th and 75th percentiles, and whiskers extend to the most extreme points that do not exceed ±1.5 times the interquartile range (IQR); further outliers (minima and maxima) are marked individually as black points beyond the whiskers; ****P < 1 × 10−4, paired Wilcoxon rank-sum test. NS, not significant.

The spatial data were collected using three spatial transcriptomics (ST) platforms, allowing rigorous cross-platform validation of these recently developed technologies (Fig. 1a, Supplementary Fig. 1a–c and Supplementary Table 1b). All three ST platforms used here provide detection of RNA molecules at subcellular resolution. As the spatial molecular imaging (SMI, also known as CosMx)21 platform had the largest gene panels, we used SMI to generate most of the data in this study. Data were generated with the SMI platform in two rounds (Discovery and Test; Fig. 1a). In the first round, we applied SMI to formalin-fixed and paraffin-embedded (FFPE) tissue microarrays (TMAs) to generate a Discovery dataset that spans 960 genes measured across 491,792 cells from 94 tumors. The discovery phase of the study was performed exclusively based on analyses of the Discovery dataset; all the spatiomolecular transcriptional programs identified in this study were derived from this dataset. In the second round, we applied SMI to another TMA from 34 additional (unseen) patients as well 4 whole-tissue sections (from patients included in the Discovery set), spanning 6,175 and 1,000 genes (Supplementary Table 1b and Supplementary Fig. 1a), respectively, and a total number of 1,233,033 cells. The Test datasets were used only in the testing phase of the study to test if the key spatiomolecular findings generalize to unseen patients and larger fields of view (FOVs) within a tumor. For technical comparison and validation, in situ sequencing (ISS; via Xenium22) was applied to profile 280 genes in an FFPE TMA of 32 tissue cores, and MERFISH23 (multiplexed error-robust fluorescence in situ hybridization (ISH)) was applied to profile 140 genes in four fresh-frozen tissue sections (Methods and Supplementary Table 1b).

For robust data processing, we developed an analytical procedure that mitigates segmentation inaccuracies (Extended Data Fig. 1a–f and Methods) and results in robust cell-type annotation through recursive clustering of the spatial single-cell gene expression profiles (Extended Data Fig. 2a,b and Methods). Applying our pipeline to the Discovery dataset identified malignant cells (n = 314,191), T cells and natural killer (NK) cells (n = 28,676), B cells (n = 16,373), monocytes (n = 45,549), mast cells (n = 606), fibroblasts/stromal cells (n = 72,861) and endothelial cells (n = 13,536; Fig. 1c,d and Extended Data Fig. 2a,b). The same procedure resulted in similar annotations of the Validation and Test datasets (Fig. 1c,d). T/NK cells were then further stratified to NK (n = 6,897), CD4+ T (n = 6,040), CD8+ T (n = 8,439) and regulatory T cells (Treg cells; n = 1,905) in the Discovery dataset (Extended Data Fig. 2c–g and Methods), with similar T/NK stratification results obtained in the Test and Validation 1 datasets (Extended Data Fig. 2h,i and Methods).

Cell-type annotations were validated in five ways. First, de novo cell-type signatures identified based on the annotated cells recapitulate well known cell-type markers (Methods, Supplementary Table 3a,b and Extended Data Fig. 2a). Second, cell-type annotations are aligned with matching hematoxylin and eosin (H&E) and immunohistochemical markers (Extended Data Fig. 3a–d). Third, cell-type annotations were aligned both spatially (Extended Data Fig. 3e) and compositionally (Extended Data Fig. 3f) when examining biological and technical replicate-matched tissue samples profiled in the Discovery and Validation 1 datasets. Fourth, we generated a unified HGSC single-cell transcriptomic atlas by integrating the HGSC spatial data with six publicly available single-cell RNA-sequencing (scRNA-seq) datasets17,24,25,26,27,28,29 (Fig. 1e,f and Supplementary Fig. 2a–c). The unified co-embedding corroborates the cell-type annotations obtained independently based on the ST data alone (Fig. 1e, Supplementary Fig. 2a,b and Methods). Fifth, using patient-matched CNA data, we examined for each gene in each cell type whether its expression is significantly associated with its copy number in the tumor (Benjamini–Hochberg false discovery rate (BH FDR) < 0.05, linear mixed-effects model (LMM); Methods). In malignant cells, gene expression is patient specific and the expression of 42% of the genes matches their CNAs (that is, ‘in cis’), compared to only 0–2% in the nonmalignant cell types (Supplementary Fig. 4a). The genes that do not show in cis RNA-to-CNA associations in malignant cells are significantly more copy number stable (median = 2, range = 0–7) compared to those showing the association (median = 3, range = 0–20, P = 5.04 × 10−10, one-sided Wilcoxon rank-sum test). In contrast to CNAs, a relatively small number of genes (0–4%) are associated with clinical and other genetic covariates in malignant and nonmalignant cell types (that is, age at diagnosis, disease stage, pathogenic BRCA1/BRCA2 mutational status and TMB; Supplementary Fig. 4a).

Initial analyses of the data reveal heterogeneous tumor tissue compositions across patients (Fig. 1f), yet a generalizable spatial organization at the macro level, where malignant cells and fibroblasts form spatially distinct compartments (which we refer to as the malignant and stromal compartments; Fig. 1g and Supplementary Fig. 3a,b), such that T/NK cells preferentially localized in the stromal rather than the malignant compartment (P < 1 × 10−4, paired Wilcoxon rank-sum test; Fig. 1h, Supplementary Fig. 3c–g and Methods). This macro-organization principle was first observed in the Discovery dataset (Fig. 1g,h and Supplementary Fig. 3c) and subsequently validated in the Validation and Test datasets (Supplementary Fig. 3a,b,d–g).

Patients with higher T/NK cell abundance had improved clinical outcomes (P = 2.03 × 10−2, univariate Cox regression, P = 2.1 × 10−4, log-rank test), demonstrating that even a relatively small area within the tumor (Supplementary Fig. 1a) is predictive of patient outcomes 5 years and even 8 years later.

Using the Discovery dataset, we set out to map the immune cell-intrinsic and cell-extrinsic factors that mark immune infiltration and exclusion. First, taking an unbiased data-driven approach, we used unsupervised methods to embed and cluster the transcriptomes of cells from each immune cell type without any spatial information. The resulting embedding and clustering shows that immune cells residing in the malignant compartment are transcriptionally distinct from those that reside outside of it (Fig. 2a). Next, we took a spatially supervised approach and identified for each of the five immune cell subtypes robustly represented in the data a tumor infiltration program (TIP), consisting of genes that are significantly (BH FDR < 0.05, LMM; Methods) overexpressed or underexpressed in the pertaining immune cell subtype as a function of proximity to malignant cells (Fig. 2b,c, Extended Data Fig. 4a and Supplementary Table 4a).

a, UMAP embedding of CD8+ T cells (Discovery dataset) derived from gene expression of all genes (top) or only T cell-specific genes (bottom). b, The association (P value and effect size, LMM) of each gene (row) from the CD8 TIP with immune cell infiltration status, when considering CD8+ T cells and other immune cell types in the Discovery dataset (columns). c, Representative ST images from the Validation 1 dataset depicting the CD8 TIP identified in the Discovery dataset. P values denote per tissue core if the expression of the CD8 TIP is significantly higher in CD8+ T cells with a high (above median) versus low (below median) abundance of malignant cells within a radius of 30 μm (one-sided t-test). d, UMAP embedding of CD8+ T cells (Validation 1 dataset) from gene expression alone. e, Average gene expression (z score) in fibroblasts (Discovery dataset) of the top 50 desmoplasia-associated genes (columns) in each tissue profile (rows, n = 87). f, Representative tissue section (HGSC24, adnexa, Discovery dataset, 1 of 100), wherein the desmoplasia-associated genes capture stromal morphology on the per-cell level (n = 1,968 fibroblasts, P value = 7.23 × 10−80, one-sided Wilcoxon rank-sum test). These results were observed repeatedly across samples, as shown in e. g, Ligand–receptor interactions (lines) consisting of genes from the CD8 TIP and their respective ligand/receptor in the malignant compartment and stromal compartment. The arrows connect each gene to the cell type where it was found to mark the malignant or stromal compartment.

CD8+ T cell TIP (CD8 TIP) demonstrates that effector and exhausted CD8+ T cells more frequently colocalize with malignant cells (P = 3.24 × 10−53, LMM; Fig. 2b,c), as also confirmed by annotating CD8+ T cells based on predefined signatures30 (Fig. 2d and Extended Data Fig. 4b–e). CD8 TIP upregulated genes include effector cytotoxicity genes as granzymes (GZMA, GZMB and GZMH) and perforin (PRF1), chemokines (CCL3, CCL4, CCL4L2 and CCL5), interferon gamma (IFN-γ; encoded by IFNG), interferon signaling genes (for example, IFITM1, IFNG, JAK1 and STAT1) and immune checkpoints (CTLA4, HAVCR2, PDCD1, TIGIT and LAG3), as well as KLRB1 (that is, CD161) and CXCR6, which have been previously reported to suppress31 or sustain32,33,34 the cytotoxic function of exhausted CD8+ T cells, respectively. CD8 TIP downregulated genes include naive and memory T cell markers (SELL, IL7R and CD44), the co-stimulatory gene CD28, the granzyme encoded by GZMK and the chemokine receptor encoded by CXCR4 (Fig. 2b). Extending CD8 TIP to the whole-transcriptome level based on scRNA-seq data17 (Methods and Supplementary Table 4b,c) revealed the upregulation of other exhaustion-associated genes30,35 (that is, ENTPD1, BST2, CD63, MIR155HG, MYO7A and NDFIP2), and downregulation of additional naive T cell markers (for example, CCR7, TCF7, SATB1 and KLF2), with MALAT1, KLF2, CCR7, GPR183 and TCF7 being the topmost downregulated genes in the extended CD8 TIP (P < 1 × 10−16, rs > 0.23, Spearman correlation).

Testing these findings in the Test datasets, the CD8 TIP identified in the Discovery dataset was validated as an infiltration marker both in unseen patients and in whole-tissue sections (P = 4.22 × 10−3 and P < 1 × 10−17, LMM, Test datasets 1 and 2, respectively; Extended Data Fig. 4f,g), and exhausted and effector CD8+ T cell subsets were enriched in proximity to malignant cells (P = 1.09 × 10−5 and P < 1 × 10−17, hypergeometric tests, for Test datasets 1 and 2, respectively; Extended Data Fig. 4h).

CD8 TIP is not associated with age at diagnosis, disease stage (III or IV), pathogenic BRCA1/BRCA2 mutations, or TMB, but does show lower expression levels in samples after neoadjuvant treatment (P < 4.42 × 10−3, LMM, also when controlling for malignant cell abundance).

To investigate the role of the stroma in preferentially colocalizing with naive and memory T cells compared to effector and exhausted T cells, we integrated the ST data with sample-matched H&E stains independently annotated by a gynecologic pathologist (Supplementary Fig. 5a–e), revealing two fibroblast subsets, one marking normal stroma and the other marking desmoplasia (that is, a neoplasia-associated alteration in fibroblasts and extracellular matrix with distinct tissue morphology36,37,38,39,40; Fig. 2e,f, Supplementary Fig. 5a–g and Supplementary Table 5a). As expected41,42, desmoplastic fibroblasts not only overexpress collagen fibril organization and extracellular matrix genes (P < 1 × 10−2, permutation test; Supplementary Fig. 5b and Supplementary Table 5b), but also upregulate CXCL12 (the cognate ligand of CXCR4; Fig. 2e) and were associated with niches enriched with T/NK cells (P < 1 × 10−4, LMM).

To link these findings to paracrine signaling, we compiled 2,678 ligand–receptor pairs based on three public resources43,44,45 (Supplementary Table 5c and Methods). Focusing on CD8+ T cells, we identified 24 ligand–receptor pairs that mark the interactions of CD8+ T cells with other cells in the malignant or stromal compartment (Methods). The resulting network (Fig. 2g and Supplementary Table 5d) manifests suppressive ligand–receptor interactions in the malignant compartment (for example, CD80/CD86–CTLA4, CD8+ T cell–monocyte; TIM3–LAGLS9, CD8+ T cell–malignant cell) and CD8+ T cell-mediated chemoattraction of other immune cells via CCL2 and CCL5. Colocalization of CXCR6–CXCL16 (CD8+ T cell–malignant cell) and CXCR4–CXCL12 (CD8+ T cell–fibroblasts) mark chemoattraction of infiltrating and stroma-residing CD8+ T cells, respectively (Fig. 2g; BH FDR < 1 × 10−10, LMM; Supplementary Fig. 5j–l). Of note, TCF1 (encoded by TCF7 and downregulated in the extended CD8 TIP) has been shown to directly repress CXCR6 expression34 (upregulated in the CD8 TIP), suggesting that this central regulator of naive and resting T cells46 represses CD8+ T cell chemoattraction to CXCL16 expressed by the malignant cells.

Mapping the spatial distributions of T/NK cells within the malignant compartment revealed that TILs preferentially colocalize with a subset of malignant cells (Methods, Fig. 3a–c, Extended Data Fig. 5a and Supplementary Table 6a). Although malignant cell states are highly patient specific (Supplementary Fig. 4b) and vary also within patients (Supplementary Fig. 4c–f), we found that the connection between TIL location and malignant cell gene expression appeared repeatedly across the heterogeneous tumors in our cohort and external cohorts (Fig. 3 and Extended Data Figs. 5 and 6).

a, Heat map of genes in the MTIL (malignant transcriptional program that robustly marks the presence of TILs) program (Discovery dataset). Average expression of the top 104 MTIL genes (rows) across spatial frames (columns). b, Top gene sets enriched in MTIL based on Gene Ontology (GO) enrichment analysis. c, MTIL spatial distributions in six representative tumor tissue profiles (6 of 100, Discovery dataset). P values denote if MTIL expression is significantly (one-sided t-test) higher in frames with high-versus-low T/NK abundance, defined based on the median level in the respective tissue section. Matching cumulative analysis is provided in Extended Data Fig. 5e. d, MTIL expression in each malignant cell (Discovery dataset, n = 297,960 cells), stratified based on the relative abundance of T/NK cells in their surroundings (top) and the presence of T/NK cells at different distances (bottom). e, ROC curves obtained for cross-validated SVM classifier using MTIL expression in malignant cells (Discovery dataset) to predict T/NK cell levels, at the sample, spatial frame and single-cell levels. f, MTIL spatial distributions in a representative region from one (of four) whole-tissue section (HGSC1, adnexa, Test 2 dataset; MTIL expression in TIL-high versus TIL-low niches, P = 2.87 × 10−107, one-sided Wilcoxon rank-sum test). A full view of the whole-tissue section is provided in Extended Data Fig. 6g. g, Mean MTIL expression in malignant cells in each FOV (Test 2 dataset, n = 878 FOVs), stratified based on the relative abundance of T/NK cells in each FOV. In d and g, in the box plots, the middle line denotes the median, box edges indicate the 25th and 75th percentiles, and whiskers extend to the most extreme points that do not exceed ±1.5 times the IQR; further outliers are marked individually with circles (minima/maxima). P values of group comparisons are derived from a one-sided Student’s t-test.

Formulating these findings, we used the Discovery dataset to identify a malignant transcriptional program that robustly marks the presence of TILs, abbreviated as the malignant TIL (MTIL) program (Fig. 3a, Extended Data Fig. 5a,b and Supplementary Table 6a). The program consists of 100 upregulated and 100 downregulated genes whose expression in malignant cells is significantly (BH FDR < 0.05, LMM) positively (MTIL-up) and negatively (MTIL-down) correlated with and predictive of T/NK cell infiltration (Fig. 3d,e). MTIL overall expression in malignant cells (Methods) reflects both inter-sample and intra-sample variation in T/NK cell levels (Fig. 3d, e), irrespective of anatomical site (P < 1 × 10−30, LMM; Extended Data Fig. 5c). MTIL continuously increases as a function of T/NK cell abundance and proximity, also when stratifying the T/NK population into its respective cell subtypes (Extended Data Fig. 5d,e). MTIL is associated and predictive of T/NK cell levels both in the Validation (Extended Data Fig. 6a,b,f) and Test datasets, generalizing to unseen patients (Extended Data Fig. 6c,f) and whole-tissue sections (Fig. 3f,g, Extended Data Fig. 6d,f,g and Supplementary Fig. 6). Likewise, an independent scRNA-seq dataset47 demonstrates that MTIL expression in malignant cells is highest in tumors annotated as ‘infiltrated’, moderate in tumors annotated as ‘excluded’ and lowest in tumors annotated as ‘immune desert’ (Extended Data Fig. 6e).

Gene-set enrichment analyses demonstrate the connection between MTIL and immune evasion48,49,50,51,52,53. MTIL-up includes chemokines (for example, CCL5, CXCL10, CXCL9 and CXCL16 the cognate ligand to CXCR6), and oxidative stress genes (for example, GPX3 and SOD2; Fig. 3a,b), and is enriched with multiple immune response genes, including antigen presentation (for example, B2M, CIITA and HLA-A/HLA-B/HLA-C), interferon gamma response genes (for example, IDO1, IFI27, IFIH1, OAS1/OAS2/OAS3, JAK1 and STAT1) and cell adhesion molecules (for example, ICAM1, ITGAV and ITGB2; BH FDR = 1.91 × 10−9, 2.86 × 10−10, 4.59 × 10−2, respectively, hypergeometric test; Fig. 3b and Supplementary Table 6b). MTIL-up also includes immune suppression genes, most notable is LGALS9, encoding galectin 9—the ligand of the immune checkpoint TIM3 (that is, HAVCR2), which is upregulated in the infiltrating T/NK cells (Fig. 2g). However, there is no significant correlation between MTIL and the expression of exhaustion signatures in the nearby T cells (rs < 0.046, P > 0.05, Spearman correlation, Discovery dataset; Supplementary Table 6c). MTIL-down reflects diverse processes including Wnt signaling (for example, CTNNB1, FZD3/FZD4/FZD6, SMO, FGFR2 and WNT7A), epigenetic regulation (DNMT3A and HDAC1/HDAC11/HDAC4/HDAC5), insulin signaling (for example, IGFR1 and IGFBP5) and