The Impact of Stroma Admixture on Molecular Subtypes and Prognostic Gene Signatures in Serous Ovarian Cancer

There have been several attempts to identify molecular subtypes and prognostic genomic signatures for advanced-stage serous ovarian cancers (1–4). The first large study to report molecular subtypes studied mostly serous tumors and identified a subtype that was associated with poor overall survival, C1 (1). Subsequently, The Cancer Genome Atlas (TCGA) study reported four similar subtypes in high-grade serous ovarian cancer (HGSOC) but found no difference in clinical outcome between the subtypes (2), which led researchers to question the clinical utility of these classifications (5). Subtype predictions using the TCGA classification are not consistently robust, sometimes producing overlapping clusters, with most tumors assigned to more than one subtype (6, 7). A recent comprehensive meta-analysis could not robustly classify subtypes across data sets (8), leading to alternative subtype classifications that were associated with survival, tumor purity, age, and lymphocyte infiltration (8).

Molecular subtype and biomarker discovery in HGSOC has mostly been performed on gene-expression profiles of bulk tumors that may include fibroblasts, immune cells, and endothelial cells (9). The gene-expression profile of a bulk tumor is the mean of the admixture of malignant epithelial, stromal, and other cells, but gene-expression studies of bulk tumors rarely record or adjust for tumor purity, which can be highly variable within and between studies. Stromal content is variable in two of the largest gene-expression profiling studies (1, 2). In the Australian Ovarian Cancer Study (AOCS), 40% of tumors in the molecular subtype C1 had low tumor percentage (≤50%) compared with 9% in the other molecular subtypes combined (1). In the TCGA study, the proportion of stromal cells exceeded 50% in about 1% of tumors, and the median proportion of stromal cells was 10% (2).

The relative proportions of tumor and stromal cells may also be an independent prognostic indicator in several epithelial cancers, including breast (10), colon (11), and possibly ovarian cancer (12). Cells in the tumor microenvironment are important for chemosensitivity (13), and genes expressed by cells in the tumor microenvironment are prognostic in ovarian cancer (4, 14–16). For example, activation of Smad signaling in cancer-associated fibroblasts in the tumor microenvironment is associated poor patient survival (16).

Additionally, cancer type – or subtype-specific biomarkers obtained in bulk gene-expression profiling studies may include genes expressed by stromal or immune cells. Recent studies report discovery of six immune subtypes in 33 cancer types (17) and two stromal subtypes in ovarian cancer (16). The potentially complex relationships between nonmalignant cells in the HGSOC tumor microenvironment and malignant epithelial cell gene expression are unknown, but sampling variability in the cell admixture may affect subtype discovery (Fig. 1). Single-cell gene-expression profiling may resolve many of these issues, but the extent to which nonmalignant cell admixture has influenced the field to date is unknown.

We systematically explored the extent to which stromal cell admixture in tumor samples influences molecular subtype and prognostic gene signature discovery (Fig. 1). We hypothesized that, if unadjusted in statistical analyses, variance in the percent of stromal cells may lead to false discovery of tumor molecular subtypes, low precision of tumor classification, and reduced power to detect true ovarian epithelial carcinoma molecular subtypes and prognostic signatures. We tested this by identifying genes that were differentially expressed between microdissected tumor and stromal cells and examining whether these genes were enriched in published serous ovarian cancer subtypes and prognostic gene-expression signatures. We show that stromal gene expression and gene signatures of stromal cell admixture are associated with overall survival in a meta-analysis of HGSOC gene expression data sets. We also investigate the number of published prognostic ovarian cancer gene signatures that add prognostic information to tumor–stroma proportion.

Details and sources of all data sets are listed in Supplementary Table S1. Raw gene-expression data for two AOCS data sets (2) were downloaded from Gene-Expression Omnibus (GEO). These AOCS data sets were (i) high-grade (G2 and G3) ovarian adenocarcinomas (n = 215; 11 endometrioid, 203 serous and 1 undifferentiated) that were previously classified as high-grade molecular subtypes C1, C2, C4, or C5 (the “AOCS data set,” GEO GSE9891) and (ii) a subset of four C1 tumors that had matched microdissected tumor and stroma (the “AOCS microdissected dataset,” n = 8, GEO GSE9890). We excluded low-grade and low malignant potential AOCS subtypes (C3 and C6 subtypes, n = 36) and unclassified tumors (n = 34). The Massachusetts General Hospital (MGH) ovarian cancer data set was received directly from the authors and included Affymetrix gene-expression profiles of pairs (n = 38) of microdissected tumor and stroma of HGSOC from women hospitalized at the MGH (GSE18520 and GSE40595, the “MGH dataset”; refs. 18, 19). The MGH data set also included unmatched normal ovary stroma (n = 10). In TCGA ovarian cancer data set (n = 518, downloaded November 24, 2010), all but one tumor was high-grade serous (2). Public gene-expression data from 16 studies were extracted from the curatedOvarianData (20) database. Nonovarian validation data sets were downloaded from GEO (Supplementary Table S1) and included microarray gene-expression profiles of microdissected stroma and tumor from noninflammatory breast cancer (n = 34 pairs; GSE5847, Boersma and colleagues, ref. 21) and invasive breast cancer (n = 28 pairs; GSE10797, Casey and colleagues, ref. 22). These also included two prostate cancer data sets with known percent stromal content (GSE17951, n = 136; GSE8218, n = 109; ref. 9).

To compare bulk and microdissected AOCS molecular profiles, raw data (cel files) were normalized using frozen robust multiarray analysis (fRMA; ref. 23), which uses an external reference normalization approach, allowing one to combine data from different batches (when these data were generated using the same platform). Other GEO raw data sets were normalized using RMA (24, 25). Data in curatedOvarianData were used as provided because these have been extensively manually curated, uniformly preprocessed, and normalized.

Analyses were performed using R and Bioconductor, versions 2.15 and 2.10, respectively, or later.

Genome-wide differential gene-expression analyses were performed using a moderated t test, and P values were adjusted for multiple testing using the Benjamini–Hochberg method (26) from the Bioconductor package limma (27). For paired microdissected stroma and tumor samples, a paired moderated t test with a P value cutoff of 0.01 was used. A separate analysis of the AOCS microdissected data set with a less stringent cutoff (P < 0.05) was also performed to identify a broader set of tumor and stromal genes. Genes that were differentially expressed between each AOCS or TCGA molecular subtypes were detected using an unpaired test. In differential gene-expression analysis of the MGH data set, a more stringent P value cutoff was applied to generate gene lists of similar length. Hierarchical clustering was performed using Pearson correlation with average linkage clustering (28).

Tissue-specific tumor or stroma marker genes were identified by transforming gene expression to binary on/off calls (29) and ranking genes that were expressed in one tissue and not expressed in the other microdissected tissues. Genes were considered tissue-specific if they had high sensitivity and specificity (area under a receiver operating curve ≥0.9).

To generate a single gene signature covariate in statistical models, a gene signature score was computed using a weighted sum approach as previously described (30). Positive weights (+1) were applied to genes upregulated in poor-prognosis tumors, BRCA mutation–like tumors, metastases, an angiogenic subtype, or malignant tumors. Negative weights (−1) were applied to the other genes in the signature. Linear regression was used to test for association between continuous variables and percent stromal content. Fisher exact test was used to test for overlap between gene sets, with the background number of genes equal to the number of probes or genes on the platform (e.g., the AOCS n = 18,769).

The gene lists reported to predict AOCS and TCGA molecular subtypes were downloaded as online publication supplements (1, 2), curated, and submitted to GeneSigDB (31). To reproduce the AOCS and TCGA molecular subtype classifiers, gene lists were used to build subtype classifiers using the single sample predictor (SSP) approach (32). The centroids of gene-expression profiles of the four molecular subtypes in the AOCS study (C1, C2, C4, and C5; ref. 1) and TCGA study (mesenchymal, immunoreactive, differentiated, proliferative; ref. 2) were calculated (Supplementary Fig. S1), and cases (tumor or stroma) were assigned to the molecular subtype with highest Pearson correlation coefficient. A case was not classified to any molecular subtype (“unclassified”) if its correlation to each one of the subtypes was less than 0.7. There was high concordance in classification between our implementation of the subtype classifiers and their original studies. The AOCS SSP molecular subtype classifier had a self-validation accuracy of 95%, sensitivity 0.93–0.97, and specificity 0.96–1.0 (Supplementary Fig. S1C). The TCGA molecular subtype SSP classifier had an overall self-validation accuracy of 89%, 0.84–0.94 sensitivity, and 0.95–0.98 specificity for each subtype (Supplementary Fig. S1D).

Paired microdissected tumor and stroma (n = 38) were computationally mixed in increasing increments of 10%, as described previously (33). The gene-expression values of each tumor (G_t) were mixed with gene-expression values of its adjacent stroma (G_s) using the linear combinations mG_t + (1 − m)G_s. Mixing parameter m is the proportion of tumor and 1 − m is the proportion of stroma, incremented by 0.1 from 0 to 1. This was performed for all genes in the molecular subtype classifier to produce a computationally simulated mixture of tumor and stroma. The assumption of linearity has been previously supported by comparing the computational values to those observed in gene expression of cocultured cells for the breast molecular subtype PAM50 gene signature (33).

Cox proportional hazards regression was used to test for association between continuous variables and survival variables [relapse-free survival, overall survival (OS)]. In survival analyses of gene signatures using independent data sets, the P values were not corrected for multiple testing. Meta-analysis of the tumor–stroma gene signatures was performed using the metafor R package and summary hazard ratios and confidence intervals of the 16 studies (n = 2,527 ovarian tumors) were calculated with a fixed-effects model.

Additional clinical data, stromal fraction, leukocyte fraction, five core immune signatures (wound-healing, macrophage regulation, lymphocyte infiltration, interferon gamma response, and TGFβ response), CIBERSORT estimates of >25 immune cell types, and mutation load for each TCGA HGSOC tumor were downloaded from the TCGA immune landscape study (17) article supplement (https://isb-cgc.shinyapps.io/shiny-iatlas/). The Spearman rank correlation coefficient was used to summarize the relationship between stromal gene-expression score and immune features in TCGA HGSOC tumors.

We investigated the extent to which non-tumor cell gene expression distinguishes HGSOC molecular subtypes using three different approaches: (i) cluster analysis to explore correlations between gene-expression profiles of microdissected tumor, microdissected stroma, and HGSOC tumors of different molecular subtypes, (ii) differential gene-expression analysis to determine if genes that are differentially expressed between molecular subtypes are also differentially expressed in microdissected tumor or stroma, and (iii) in silico cell admixture analysis to discover the robustness of molecular subtype classifications to increasing proportion of stromal cells.

We defined gene-expression signatures that characterize HGSOC and its adjacent stroma using paired laser microdissected tissue from two independent studies. Using the AOCS microdissected data set, we defined a 688-gene signature using C1/stromal tumors (n = 8; P < 0.01, paired t test with modified variance (27) and FDR correction (ref. 26; Supplementary Table S2). This AOCS tumor–stroma gene signature contained 461 and 227 unique genes with increased expression in C1 stroma and tumor, respectively, which are referred to as the “stromal gene set” and “tumor gene set.” A list with a larger number of genes extracted using lower differential expression stringency (P < 0.05) is available in Supplementary Table S2. In a second, independent analysis of paired microdissected tumor and stroma from HGSOC tumors of 38 patients in an MGH study (18), we identified a similar gene signature. This MGH tumor–stroma gene signature of 519 unique genes contained 429 and 90 genes that had higher expression in stroma and tumor, respectively (Supplementary Table S2). Despite the small sample size from which the AOCS tumor–stroma signature was derived, the MGH signature had significant overlap with it: 24/90 genes overexpressed in tumor (Fisher exact test P < 10⁻¹⁰) and 122/429 genes overexpressed in stroma (P < 10⁻²⁵). The concordance of these signatures was supported by the observation that the AOCS signature correctly discriminated all but one of the microdissected tumor and stroma MGH specimens using unsupervised hierarchical cluster analysis (Supplementary Fig. S2).

Using these tumor–stroma gene signatures, we investigated the contribution of stromal cell gene expression to the HGSOC molecular subtype using unsupervised cluster analysis. The 688-gene AOCS tumor–stroma gene signature partitioned the AOCS high-grade serous tumors (n = 215) into two clusters (Fig. 2): those with low versus high stromal gene expression. These two clusters could be further subdivided to form four clusters, which largely reflected the four AOCS HGSOC molecular subtypes: C4/CA125, C5/MYCN deregulation (34), C2/immunoreactive, and C1/stromal. We observed that matched microdissected tumor samples did not cluster with the corresponding bulk tumor. Instead, microdissected C1 tumors clustered with C4/CA125 molecular subtype tumors, whereas microdissected C1 stroma clustered with the C1 bulk tumors.

The MGH tumor–stroma gene signature similarly partitioned the AOCS data set into clusters (Supplementary Fig. S2C) which reflected AOCS HGSOC molecular subtypes, particularly the split of C1 compared with others. The observation that C1 tumors had higher stromal content is known, but the high concordance of tumor and stroma clusters with published HGSOC molecular subtypes suggested that genes that are differentially expressed between tumor and stroma might be sufficient to discriminate the HGSOC molecular subtypes.

These molecular subtypes were then further examined using gene-expression analysis, in which we quantified the relative contribution of stroma and tumor gene expression to each molecular subtype. We performed global differential gene expression analysis between each pair of AOCS tumor molecular subtypes [C1, C2, C4, and C5; limma (ref. 27), 2-fold change and P < 0.001 after FDR correction; Supplementary Table S3A]. The number of genes that discriminated between each pair of subtypes varied in size from 343 (C1 vs. C2) to 1,267 (C1 vs. C5; Fig. 3A). We tested the overlap between these gene lists and the AOCS stromal and tumor gene sets. The stromal and tumor gene sets derived from the AOCS microdissected data set's C1 tumors were again used so that the gene sets were most comparable. We observed that any contrast involving C1 had strong overlap with the stromal gene set (Fisher test P < 10⁻⁶¹), as did the contrast of C2 versus C4 (Fisher test P < 10⁻⁷³). Specifically, 54% and 41% of genes significantly different in contrasts of C1 versus C2 or C1 versus C4, respectively, were upregulated in stromal cells (Fig. 3A). This suggests that at least C1 is primarily driven by different amounts of stromal gene expression. In contrast, only 9% of genes differentially expressed between C4 and C5 were stromal genes (P = 1), consistent with the subtypes' reported low amount of desmoplasia (1). The tumor gene set and “other” genes dominated C5 differential gene expression, consistent with C5's distinct molecular properties (2, 34).

Based on these results, we also investigated the influence of stromal gene expression in TCGA HGSOC molecular subtypes (Fig. 3B; Supplementary Table S3B). The TCGA mesenchymal, immunoreactive, differentiated, and proliferative subtypes have highest global gene-expression correlation with the AOCS subtypes C1, C2, C4, and C5, respectively (Fig. 3C), and the mesenchymal subtype has significantly higher percent stromal content than other subtypes (t test P = 1.35 × 10⁻⁷). Results from analysis of the TCGA study confirmed findings in the AOCS subtype analysis. Genes differentially expressed (P < 0.001, FDR correction) in the TCGA mesenchymal subtype, similar to C1, were mostly in the stromal gene set (Fig. 3D, 59%), and overlapped strongly (35%) with genes differentially expressed in C1 (Fig. 3D). Similar to the AOCS C4 and C5 subtypes, few genes in the stromal gene set were differentially expressed between TCGA's differentiated and proliferative molecular subtypes (<5%, P = 1). Lastly, genes that were differentially expressed in the proliferative subtype, similar to C5, were relatively enriched in tumor-specific and “other” genes.

Similar results were observed when these analyses were repeated with the MGH tumor and stroma signatures (Supplementary Table S3C). Similar to the AOCS tumor signature, the MGH tumor signature did not significantly overlap with gene lists from pairwise comparisons of any two AOCS or TCGA molecular subtypes (P = 1 for overlap with any comparison). Consistent with the AOCS stromal signature, the MGH stromal signature strongly overlapped with comparisons between the AOCS C1 subtype and either C2 or C4 and similar comparisons between the TCGA mesenchymal subtype and either the immune or differentiated subtypes. We observed that 24% and 19% of genes differentially expressed between C1 versus C2 (P < 10⁻¹³) and C1 vs. C4 (P < 10⁻¹⁵), respectively, were MGH stromal genes. There was no significant overlap between MGH stromal genes and any pairwise comparison involving the AOCS C5 subtype (P = 1). Therefore, in two independent studies, stroma-associated expression explains a large portion of the differential gene expression of select subtypes.

We next investigated if the gene signature classifiers for these molecular subtypes were stable when the proportion of stromal cells changed. We trained an SSP classifier using the gene lists reported to predict each of the AOCS C1, C2, C4, and C5 molecular subtypes (ref. 1; Supplementary Fig. S1A and S1C). When applying the AOCS classifier to epithelial carcinoma cells microdissected from HGSOC tumors in an independent cohort (MGH data set n = 38), most microdissected tumors (20/38) were assigned to subtype C4 (Supplementary Fig. S1E and S1I). Matched microdissected stroma was all classified as C1 (n = 23) or unclassified (n = 15; Supplementary Fig. S1G and S1J). However, computational mixing of only 10% stromal gene expression with paired microdissected tumor caused 6/38 tumors to be reclassified (Supplementary Fig. S3A), including C4 tumors that were reclassified as C1 (n = 2) or C2 (n = 2). Molecular subtype C5 tumors remained stable, reflecting the lack of stromal gene expression determining this subtype. Five additional tumors were predicted to have a different molecular subtype when the stromal content was increased to 20%, and at 30% stroma, over one third (15/38) of tumors' subtypes were reassigned.

We observed similar results regarding the robustness of TCGA molecular subtypes with increasing proportions of stroma. We built an SSP classifier of TCGA molecular subtypes using their 100-gene signature (ref. 7; Supplementary Fig. S1B and S1D). When it was applied to classify microdissected tumors in the MGH data set, most were assigned either to the differentiated subtype (23/38), which is similar to the AOCS C4 subtype, or to the proliferative subtype (8/38), which is similar to C5 (Supplementary Fig. S1F and S1I). Microdissected stroma was classified as mesenchymal (20/38) or was unclassified (10/38; Supplementary Fig. S1H and S1J). With 10%, 20%, and 30% stroma, 0/38, 5/38, and 8/38 tumors, respectively, were classified to a different molecular subtype (Supplementary Fig. S3B). Thus, the percentage of stromal content influences subtype classification.

The AOCS C1/stromal subtype is reported to have a worse prognosis, and therefore we explored whether percent stroma itself was prognostic. We confirmed that the AOCS C1/stromal subtype, which has high stromal proportion, was associated with poorer OS. C1 had a hazard ratio (HR) of 1.72 (95% confidence interval 1.15–2.58, P = 0.0083). The AOCS C1/stromal molecular subtype had a HR of 1.49 (95% CI, 0.99–2.26, P = 0.056) within high stage (III or IV) tumors.

In the TCGA data, pathologist scores of percent stromal content had a weak but significant (HR 1.01 for each percent increase in stromal content, CI 1.0–1.02, P = 0.024) association with OS in univariate Cox proportional hazards regression analyses of high stage tumors, but not across all tumors (HR 1.0, CI 0.99–1.02, P = 0.088; Fig. 4C for survival association with stromal content > 30%). In other studies, no estimate of the percentage of stroma was available, but we assessed the prognostic power of a tumor–stroma gene signature in an aggregated cohort of 2,527 HGSOC patients from 16 studies, which included the TCGA and AOCS data sets (20). The 688-gene AOCS tumor–stroma signature had a HR of 1.17 (95% confidence interval 1.11–1.23; Fig. 4A).

The AOCS and MGH tumor–stroma signatures were generated using differential gene-expression analysis and might include genes that are expressed in both tumor and stroma cells, albeit at different levels. Therefore, we derived a list of markers specific to tumor or stroma that better reflect cell admixture. Using expression profiles of paired microdissected epithelial tumor (MGH data set n = 38), adjacent tumor stroma (n = 38), and unrelated normal ovarian stroma (n = 10), we applied a barcode approach (29) to convert gene expression to binary expressed or not expressed calls. We then extracted a list of 28, 11, and 38 tissue markers with high specificity (Supplementary Fig. S4) for tumor, stroma adjacent to tumor, or normal ovary stroma (Supplementary Table S2). The 11 cancer-associated stromal tissue markers had similar prognostic power (HR, 1.17; CI, 1.11–1.23; Fig. 4B), but the tumor markers were not prognostic (HR, 0.96; CI, 0.96–1.01) and performed comparably to a randomly selected 20-gene signature (HR, 0.96; CI, 0.96–1.02). These results support that a higher stroma proportion in HGSOC tumors is associated with worse OS.

Given that stromal gene expression was associated with survival in HGSOC tumors, we hypothesized that some published prognostic ovarian cancer gene signatures might incidentally be proxies for the presence of stromal tissue. Published gene signatures of ovarian cancer (n = 61) were downloaded from the database GeneSigDB (31) or curated from the literature (Supplementary Table S4). Gene signatures varied in size from 4 to 1,903 genes, with a median of 47 genes. Most gene signatures, including the TCGA study's prognostic signature and Konstantinopoulos and colleagues BRCAness signature (35), were not significantly enriched in stromal genes (Supplementary Table S4) and not associated with percent stromal content in TCGA samples (Table 1).

Over a third of published prognostic gene signatures (n = 21/61) were enriched in AOCS stromal genes, and 9/61 were enriched in MGH stroma genes (P < 0.05 after FDR correction). Eight signatures were further analyzed (Table 1); seven had strong enrichment (P < 0.001) in both the AOCS and MGH stromal genes (Supplementary Table S4), and the Biade and colleagues benign tumor signature (36) nearly met these P value cutoffs. All of these eight signatures (AOCS signatures; Biade and colleagues benign tumor signature, ref. 36; Bentink and colleagues, angiogenic signatures, ref. 14; and prognostic signatures from Bignotti and colleagues, ref. 37; Bonome and colleagues, ref. 38; and Spentzos and colleagues, ref. 39) were positively associated with pathologic scores of percent stromal content in TCGA tumors (linear regression, P < 10⁻⁸; Table 1B). Genes reported to be associated with poor prognosis overlapped with the stromal gene set in 7/8 gene signatures (P < 0.001), whereas the good prognosis genes were often enriched in the tumor gene sets (Table 1A). The Biade and colleagues signature (36) was an exception, which had increased stromal cell proportions in benign tumor samples (36).

Most of the eight gene signatures that were enriched in stromal genes were no longer associated with poor OS when adjusted for pathologist's estimates of percent stromal content in TCGA tumors in Cox proportional hazard regression analysis (Table 2A). The AOCS prognostic gene signatures remained associated with OS and progression-free survival (PFS) when adjusted for stromal content (P < 0.01 and P < 0.05, respectively; Table 2A) but were no longer significant after adjusting for both stromal content and stage. Gene signatures not enriched in stromal genes, such as the BRCAness signature, remained associated with OS in multivariate models that adjusted for both stromal content and stage (P < 0.05; Table 2A).

Several types of carcinoma cells have complex interactions with cells in the tumor microenvironment. Therefore, we asked whether stromal genes derived from ovarian cancer are specific to HGSOC or might more broadly characterize epithelial cancers. Each of the eight ovarian cancer signatures enriched in stromal genes were strongly associated with stroma versus epithelial components of breast cancer (Boersma and colleagues, ref. 21; Casey and colleagues, ref. 22) or percent stromal content in prostate cancer (Wang and colleagues, ref. 9; Table 1, linear regression or paired t test P < 10⁻⁴ in at least 3 of 4 data sets).

We observed that the strength of association between greater stromal proportion and survival was variable across data sets. We hypothesized that this could be related to sampling heterogeneity. In the AOCS C1 subtype, 63% of tumor samples were from extra-ovarian sites (most commonly peritoneum, as defined by the AOCS study), which often reflects a higher stage and thus worse prognosis, compared with 19% extra-ovarian sampling in other AOCS subtypes (1). In the TCGA data set, every tumor except for four were reported to be sampled at the ovary/pelvis (2). After adjusting for whether the sampling location was extra-ovarian, C1's association with OS (adjusted P = 0.065) was substantially attenuated. Thus, C1's prognostic power is partially related to sampling location, which may in turn reflect tumor stage.

We show that the proportion of stroma admixture in bulk tumors is a source of differential gene expression that affects molecular subtyping and prognostic gene signature discovery in HGSOC.

Although adjacent stromal cross-talk can affect nearby epithelial cells' gene-expression patterns (16), this is a distinct phenomenon from changes in observed gene expression caused by the admixture of stromal and other cells in bulk tumors (tumor purity). In bulk gene-expression profiling, the “observed” measure is an average of all cells in the admixture, which varies if the proportions of cell types change. We show that stromal genes differentiate the AOCS C1 molecular subtype or the TCGA mesenchymal subtype from other HGSOC molecular subtypes. This finding is supported by a recent IHC study of TCGA mesenchymal subtype tumors showing that several subtype marker genes were expressed exclusively or primarily in stroma (40). The dependence of the C1 subtype classifier on stromal genes may also explain why the TCGA mesenchymal subtype lacks robustness when applied to other data sets that may have variable tumor purity (8) and why few reproducible mutations or copy-number alternations have been described for the subtype. In miRNA analysis, the mesenchymal TCGA molecular subtype clusters with the proliferative subtype (2). In contrast, the AOCS C5 and TCGA proliferative subtypes are reported to have genetic alterations, including copy-number variation, in several genes (2, 34), and consistently, these subtypes' gene-expression patterns were not well explained by stromal genes.

Small increases in stromal gene expression caused subtype misclassification in our computational model that mixed stroma and tumor gene expression profiles. The influence of non-cancer cell gene expression on subtype classification is potentially a source of systematic error in other solid tumors. The accuracy of signature classifiers that predict breast tumor molecular subtype (PAM50) or risk of recurrence (Oncotype DX) decreases considerably with decreasing proportion of tumor cells and increasing proportion of normal cells (33). Molecular subtype definitions developed from gene-expression profiles of bulk ovarian, breast, and other carcinomas may need to be redefined such that molecular subtypes more precisely model the contributions of tumor, stroma, and other cell types. Thus far, at least two HGSOC stroma subtypes have been proposed (16) and their relationships to epithelial HGSOC molecular subtypes are yet to be described.

Deconvolution of cell mixtures in bulk gene-expression studies is an active area of research (41), and we identified tissue-specific marker genes of HGSOC epithelial carcinoma and stroma (Supplementary Table S2) that may be applied to computational prediction of cell proportions and in the study of tumor/stroma interactions. These tissue-specific genes highlight the importance of studying both epithelial and stromal cell biology in cancer. CEP55 had the highest specificity and sensitivity of the tumor-specific genes, and it codes for a protein that interacts with BRCA2 and plays a critical role in regulating the final step of cytokinesis (42). Some stroma-specific were possibly related either to macrophage activity or to the recently described TGFβ-dependent and TGFβ-independent Smad signaling (16) in stroma of HGSOC. Two of the 11 stroma-specific genes (IFFO1 and GLIPR1) were strong predictors both of pathologist's stroma proportion scores in TCGA tumors (Pearson correlation 0.7 and 0.67, respectively, P < 0.0001) and of a gene signature that predicts macrophage/monocyte activity (Supplementary Fig. S4, colony-stimulating factor-1, CSF1 signature; refs. 17, 43). Three other stroma-specific genes were strongly correlated with TGCA immune landscape scores for TGFβ response (Supplementary Fig. S4) in TCGA tumors (17). These included thrombospondin 1 (THBS1), which is associated with tumor growth and metastasis in gastric carcinoma (44), and biglycan (BGN), which has been shown to form a complex with either TGFβ₁ or TGFβ₁ type I receptor to intensify the phosphorylation of Smad2/3 in cultured endothelial cells (45).

Further understanding of stromal cell biology in HGSOC would be valuable, particularly because stromal gene expression was associated with a worse prognosis (Table 2; ref. 12). Few studies provide pathologist estimates of stroma proportion, and we developed gene signature models for stroma proportion, which were weakly associated with worse survival in a meta-analysis of 2,527 tumors (HR 1.17; ref. 20). We performed a systematic analysis of published prognostic signatures in HGSOC and discovered that most, including the BRCAness signature (35), do not predict stroma proportion. Several signatures were highly enriched in stromal genes and did not provide additional prognostic information beyond the stromal cell proportion, suggesting that adjusting for stromal proportions maybe useful in subtype and signature discovery.

In ovarian cancer, the link between stroma and survival may sometimes be related to other variables, such the sampling methods. For example, the prognostic performance of the AOCS stroma-associated subtype, C1, is partly explained by stage or sampling location. This association with sampling location also suggests that stromal content may be variable within a single tumor or between studies. In Biade and colleagues (36), the benign rather than poor-prognosis tumors were enriched in stroma, and a recent analysis of tumor purity across several cancers showed that tumor purity's relationship with survival is often confounded with other variables (46).

These issues of sampling heterogeneity in bulk tissue were recognized as early as 1999 when Golub and colleagues (47) remarked: “Studies will require careful experimental design to avoid potential experimental artifacts – especially in the case of solid tumors. Biopsy specimens, for example, might have gross differences in the proportion of stromal cells. Furthermore, accumulating evidence suggests that the composition and gene expression of stroma may vary depending upon the specific anatomic location of the specimens. Blind application of class discovery could result in identifying classes reflecting the proportion of stromal contamination in the samples, rather than underlying tumor biology. Such ‘classes’ would be real and reproducible but would not be of biological or clinical interest.”

Our results strongly support single-cell analysis or microdissection of tumor samples in gene-expression studies so that the features of the epithelial tumor cells, stromal, and immune components can be each analyzed. Retrospective analysis of bulk gene-expression patterns in light of emerging single-cell sequencing data is needed to refine molecular subtypes in ovarian cancer and other carcinomas and enable a more comprehensive search for molecular targets amenable to therapeutic intervention.

C.S. Mitsiades is an employee at Takeda, reports receiving honoraria for one-time consulting/advisory board from Ionis Pharmaceuticals and Fate Therapeutics, and reports receiving commercial research grants from Janssen/Johnson & Johnson, TEVA, EMD Serono, AbbVie, Karyopharm, Sanofi, and Arch Oncology. W. Wei, G. Parmigiani, and M.J. Birrer have patent ownership in WO2014153442A2, METHODS AND SYSTEMS FOR TREATMENT OF OVARIAN CANCER. G. Parmigiani has ownership interest in CRA Health and Phaeno Biotech. J. Quackenbush is a scientific advisory board member for Caris Life Sciences. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M. Schwede, S.C. Mok, J. Quackenbush, A.C. Culhane

Development of methodology: M. Schwede, D.P. Harrington, A.C. Culhane

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Schwede, S.C. Mok, W. Wei, M.J. Birrer

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Schwede, A. Basunia, G. Parmigiani, D.P. Harrington, M.J. Birrer, A.C. Culhane

Writing, review, and/or revision of the manuscript: M. Schwede, L. Waldron, A. Basunia, M.A. Merritt, G. Parmigiani, D.P. Harrington, J. Quackenbush, M.J. Birrer, A.C. Culhane

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.J. Birrer, A.C. Culhane

Study supervision: J. Quackenbush, M.J. Birrer, A.C. Culhane

This study was supported in part by grants from the NCI at the NIH (RC4CA156551 to M.J. Birrer, P30CA006516 and R01CA174206 to G. Parmigiani, R01CA133057 to S.C. Mok, R03CA191447 to L. Waldron, U19CA148065 and P01CA087969 to J. Quackenbush) and the National Institute on Minority Health and Health Disparities at the NIH (G12MD007599 to L. Waldron) and by the Assistant Secretary of Defense Health Program, through the Breast Cancer Research Program (W81XWH-15-1-0013 to A.C. Culhane). Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the U.S. Department of Defense.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.