Ovarian cancer cell line RNA sequencing (RNA-seq) data was obtained from three sources: Klijn et al. (2015), the Cancer Cell Line Encyclopedia (CCLE, Ghandi et al., 2019) and Reyes et al. (2019). Raw sequence files from Klijn et al. (2015) were obtained with permission from the Genentech Data Access Committee (dataset ID EGAD00001000725) and downloaded from the EMBL-EGA servers using the pyEGA3 download client. For the CCLE and Reyes datasets, raw sequence files in FASTQ format were obtained from the European Nucleotide Archive (accessions PRJNA523380 and PRJNA470980). Klijn et al. (2015), CCLE and Reyes et al. (2019) provide RNA-seq data for 48, 47 and 12 OC cell lines respectively. In total, RNA-seq reads were available for 69 unique OC cell lines. An overview of datasets used and associated metadata are detailed in Supplementary Table S1.

OC cell lines that were misclassified, engineered or non-epithelial in origin were removed from analysis (n = 11, Supplementary Table S2). Two cell lines (DOV13 and OVCAR433) could not be processed due to errors with the EGA download client pyEGA3. 56 EOC cell lines remained for further analysis. A literature search on site of origin, original subtype annotation, TP53 mutational status and treatment for each cell line and is detailed in Supplementary Table S3.

Ovarian primary tumor tissue data was obtained from three independent datasets (Corona et al., 2020; Kang et al., 2021; Akasu-Nagayoshi et al., 2022). Corona et al. (2020) provides both RNA-seq and chromatin immunoprecipitation sequencing (ChIP-seq) data on EOC tumors of various subtypes: CCOC (n = 5), HGSOC (n = 5), ENOC (n = 4) and MOC (n = 5). RNA-seq data for these 19 primary tumors were analyzed in this study. Kang et al. (2021) generated transcriptomic profiles for 51 HGSOC tumors of various stages and sensitivities to chemotherapy. Finally Akasu-Nagayoshi et al. (2022) conducted RNA-seq on stage III and IV primary tumor tissue of various subtypes: CCOC (n = 11), HGSOC (n = 8) and ENOC (n = 4). There is no publicly available RNA-seq data for LGSOC primary tumors (although a number of array-based sequencing datasets are available). In total, transcriptomic data for 16 CCOCs, 64 HGSOCs, 8 ENOCs and 5 MOCs were used.

Quality control checks were performed on all FASTQ files using the FastQC package (Andrews, 2010). Forward and reverse reads for each cell line were mapped using the Spliced Transcripts Alignment to a Reference (STAR) software, version 2.7.9a (Dobin et al., 2013). The reads generated by Akasu-Nagayoshi et al. (2022), were single-end, which was taken into account during the alignment step. Reads were mapped to the human reference genome sequence (GRCh38) from Gencode (Frankish et al., 2021). Genome indices were generated using the comprehensive Gencode v39 gene annotation. Separate indices were required for each read length (Supplementary Table S1). Transcript level counts were generated using htseq-count (Anders et al., 2015). Strandedness of each dataset was taken into account during this step. Counts were compiled into a single matrix using R studio (v4.1.1). ComBat (Leek et al., 2012) was used to correct for batch effects between the three independent sources of data.

NMF was carried out on the processed count matrix, filtered to retain only cell lines, in R studio (v4.1.1). A variance stabilizing transform was applied to read counts using the DESeq2 package (Love et al., 2014). In order to retain the most variable genes, transcripts were filtered to exclude those with a median absolute deviation of less than 1.5. There were 2233 transcripts with median absolute deviation of ≥1.5; these were used as input for the NMF R package (Gaujoux and Seoighe, 2010). Estimation of factorization rank (r) was carried out by running NMF for values of r from 2 to 8, using 50 random initiation points. Quality measures including the cophenetic correlation coefficient, residual sum of squares and dispersion were used to select the most suitable value of r. In the presented data, clustering for r values of 2 and 5 displayed high quality metrics (Figure 2). After selection of a suitable r value, 200 runs of NMF were carried out with a factorization rank of 5 and a fixed random seed. Samples were assigned to one of five clusters based on maximum coefficient matrix values. The most likely subtype of each NMF cluster was inferred by comparison with previous studies and literature annotations.

For cell lines present in the CCLE study, mutational and copy number landscape was originally determined by Ghandi et al. (2018) and visualized using cBioportal (Cerami et al., 2012; Gao et al., 2013). Mutational profiles of cell lines not present in the CCLE were determined from multiple sources (Anglesio et al., 2013; Beaufort et al., 2014; Tate et al., 2019). Copy number data was not available for 59M, UWB1289, OVCA420, OVCA432, PE01, HEY, RMGII, TYKNUCPR, COV413A, COV413B, GTFR230, VOA1056 or OVCA429. Mutational data was not available for a subset of genes in the RMGII, TYKNUCPR, COV413A, OVCA431, COV413B, GTFR230, VOA1056, and OVCA429 cell lines. A literature search was also performed to investigate the mutation frequency of 26 cancer driver/tumor suppressor genes in each of the five subtypes of EOC (Supplementary Table S4)

The uniformly processed EOC cell line and tumor tissue counts were upper quartile normalized. ComBat (Leek et al., 2012) was then used to correct for batch effects between different sources of data. The normalized and batch corrected counts were filtered to keep the 2233 most variable transcripts in EOC cell lines (as identified earlier). The Spearman correlation was then calculated between cell lines and tumor tissues. Results were plotted as a correlation matrix using the corrplot package (Wei et al., 2021) and as boxplots using the ggplot2 package (Wickham, 2016), ordered by median correlation.

Non-negative matrix factorization is a method of unsupervised learning often applied to gene expression data to extract biologically relevant information (Brunet et al., 2004; Gaujoux and Seoighe, 2010; Barnes et al., 2021). It functions by transforming a large, non-negative matrix (such as read counts) into a lower dimensional space: two smaller matrices W and H. The basis matrix W delineates the contribution of a small subset of genes to what are termed ‘metagenes’. This is essentially a decomposition of genes into those whose co-expression influences cluster assignment (Brunet et al., 2004). The coefficient matrix H details the co-expression pattern of metagenes in each sample, and can be used to cluster samples into a defined number of groups (r). Seeding is used to initialize the starting point for the NMF algorithm (i.e., to provide starting values for the basis and consensus matrices). The NMF package contains a number of built in seeding methods. Here the ‘random’ method is used, whereby these initial values are obtained from a uniform distribution with the same range as the input matrix. For reproducibility, a numerical value is passed into the seed function to seed the random number generator. In order to achieve a stable result from this random seeding method, multiple runs of NMF are required.

The number of metagenes and sample clusters are defined by the factorization rank (r), a critical parameter that is selected by the user-ideally, r should be small enough to reduce noise but large enough to preserve biologically relevant information. Brunet et al. (2004) suggests that the factorization rank should be chosen as the smallest value of r for which the cophenetic correlation coefficient starts decreasing. The cophenetic correlation coefficient indicates the robustness of clusters for a given choice of r. Frigyesi and Höglund (2008) suggest investigating methods other than the cophenetic correlation coefficient. It is argued that an increase in r is only relevant if the information captured by factorization is greater than that derived from permuted data, otherwise this increase will lead to overfitting. They suggest using the smallest value of r where the decrease in residual sum of squares (RSS) is lower than the decrease observed in randomized data. By comparing the residual error of NMF from the original data to that of permuted data, a solution is identified that is more information than noise, while preventing overfitting.

The consensus matrix also provides insight into the stability of the clusters. The entries of the consensus matrix reflect the probability two samples belong to the same cluster. The dispersion of the consensus matrix measures the reproducibility of clusters obtained from NMF; 1 for a perfect consensus matrix, and between 0 and 1 for a scattered consensus matrix. A consensus map can be plotted, which is a consensus matrix computed over multiple independent NMF runs, which is the average of the connectivity matrices of each separate run (Figure 2A). Quality of clustering can also be assessed using silhouette scores. This takes into account how similar a particular sample is to the other samples in its assigned cluster (mean intra-cluster distance), as well as the similarity to samples in other clusters (mean neighboring-cluster distance). A silhouette score close to 0 indicates the sample is on or near the decision boundary of two clusters whereas a silhouette score of 1 shows that the sample is distinct from clusters to which it does not belong. A negative silhouette score indicates a given sample has likely been erroneously assigned to a given cluster.

r was estimated by performing NMF with 50 runs of each value of r from 2 to 8. This was also completed for variance stabilized counts that had been permuted using the randomize function from the NMF package (Gaujoux and Seoighe, 2010). In the presented data, clustering for r values of 2 and 5 displayed high quality metrics (Figures 2A, B). Cophenetic correlation coefficients for r = 2 and 5 show robustness, and a drop in cophenetic correlation coefficient for r = 3 and r = 6 indicates less stability. For r = 6 the decrease in RSS in the observed data is less than the decrease in RSS in the randomized data. Therefore, a factorization rank of 5 was chosen for NMF analysis. The most likely subtype of each of the 5 NMF clusters was inferred by comparison with previous studies and literature annotations. In order to investigate the intra-cluster similarity, principal component analysis was carried out on the variance stabilized counts (Figure 2C). Each cluster assigned by NMF largely grouped together on the PCA plot. The ENOC and CCOC lines clustered away from the other samples, whereas there was some overlap between certain MOC, LGSOC and HGSOC cell lines. This suggested that some cell lines may have been misannotated by NMF or may display transcriptional characteristics of multiple EOC subtypes.

The results obtained in using NMF analysis predominantly validated subtypes assigned in other studies (Supplementary Table S3). This study analyzed 43, 42, 22, and 20 EOC cell lines in common with Barnes et al. (2021), Domcke et al. (2013), Beaufort et al. (2014) and Anglesio et al. (2013), with subtype assignment agreeing in 39, 36, 15, and 14 cases respectively. Disagreements between our classifications and those of other studies added clarity to certain subtype assignments. We also utilized NMF to assign subtypes to five previously unannotated cell lines.

There were a number of differences between the subtypes assigned to certain cell lines in this study compared to other previous works, often adding clarity to stratifications. OAW42, while classified as serous by Beaufort et al. (2014), was found to cluster with CCOC lines both here and by Barnes et al., (2021). Domcke et al. (2013) also classed this line as unlikely to be of HGSOC origin. The fact that OAW42 is TP53 wild type, and possesses characteristic CCOC mutations (ARID1A and PIK3CA, Figure 3), suggests that this line is likely a model of CCOC. OVCAR8 was classified as LGSOC both here and in Barnes et al. (2021), unlikely to be HGSOC by Domcke et al. (2013) and as HGSOC by Anglesio et al. (2013). Although this line has a TP53 mutation and MYC amplification (both of which are common to HGSOC), it also possesses KRAS and ERRB2 mutations (Figure 3), suggesting it is more likely a model of LGSOC, in agreement with Barnes et al. (2021) and Domcke et al. (2013). HEY was classified here as LGSOC, which directly contradicted with Anglesio et al. (2013), where this line was determined to be of HGSOC origin. HEYA8 is another cell line analyzed, which was derived from the peritoneal cavity of mice injected with HEY cells (Supplementary Table S3). HEYA8 was categorized as a LGSOC line both here and by Barnes et al. (2021). It is likely that HEY is of LGSOC origin, as it possesses characteristic LGSOC mutations (KRAS and BRAF, Figure 3), and it is the clonal ancestor of the HEYA8 cell line. Finally, OV56 was classified as LGSOC in Barnes et al. (2021), unlikely HGSOC in Domcke et al. (2013) and ENOC/CC in Beaufort et al. (2014). It is important to note that although classified as LGSOC by Barnes et al. (2021), the authors suggest that OV56 more likely represents CCOC, due to the presence of KRAS, ARID1A, PIK3R1 and PTEN mutations (Figure 3). Indeed, we found this to be true in our analysis, with OV56 clustering with the CCOC cell lines.

Five cell lines (GTFR230, OVCA420, OVCA429, OVCA432, and TYKNUCPR) have not been investigated in terms of EOC subtype since their creation, with original annotation data unavailable for four of these lines (OVCA420, OVCA429, OVCA432, and TYKNUCPR). Additionally, these cell lines are rarely cited in literature and as such, mutational and copy number data was largely unavailable (Figure 3). TYKNUCPR is from a population of TYKNU cells exposed to cisplatin chemotherapy, so it is likely that these cell lines share a common subtype. Indeed, both lines were shown to cluster together, with putative ENOC cell lines. GTFR230 was originally derived from a stage IC MOC, although we show here that it clustered with LGSOC lines with a relatively high silhouette score (Figure 2A). As no mutational or copy number data was available for this line, an investigation into the similarity of GTFR230 to LGSOC and MOC tumors could not be carried out. This highlights the importance of other methods to determine subtype molecular similarity. OVCA432 was assigned as HGSOC in this study, while OVCA420 and OVCA429 both clustered with MOC lines. Mutational data was available for OVCA420, which showed that this line may possess characteristics of both HGSOC (TP53, BRCA1, CDK12, RB1 mutations) and ENOC/CCOC/MOC (ARID1A) (Figure 3). Ultimately, we assign potential subtypes to these five cell lines based on NMF cluster membership. However, we note that due to the lack of additional mutational and copy number data, further investigation is needed to confirm classification of not only these lines, but all lines analyzed in this study.

As discussed previously, subtype assignment of certain cell lines can prove difficult to resolve using transcriptional data, especially when mutational or copy number data is unavailable. Comparison of cell lines to tumor tissue is required to determine the suitability of these in vitro models. Some progress has been made in identifying suitable cell lines to utilize as HGSOC models. Both Domcke et al. (2013) and Yu et al. (2019) compared various molecular characteristics of cell lines to primary tumors of HGSOC origin, providing researchers with recommendations of the most suitable HGSOC lines. As of yet however, no study has extended analysis of this nature to multiple subtypes, to compile a reference set of cell lines representative of other EOC subtypes.

Publicly available transcriptomic data was available for primary tumors of HGSOC, CCOC, ENOC and MOC origins (93 samples in total, Supplementary Table S1). No RNA-seq data was available for LGSOC primary tumors, therefore this subtype could not be included in correlation analysis. Following normalization and batch corrections, we calculated correlation of gene expression profiles between EOC cell lines and primary EOC tumors (Figures 4, 5). Median correlations ranged from 0.28 to 0.59, with a similar range observed in Yu et al. (2019). In general, we observed that cell lines we assigned as HGSOC, CCOC and MOC were highly correlated to their respective primary tumor subtypes, whereas ENOC lines were poorly correlated to EOC overall (Figures 4, 5). We also note that most cell lines classified as LGSOC correlated poorly with tumors of HGSOC, CCOC, MOC and ENOC origin. Our data suggests that these lines may be useful models of LGSOC, although further analysis is needed to compare similarity to actual LGSOC tumor samples.

All HGSOC cell lines (apart from Fuov1) were within the top 20 most correlated cell lines to HGSOC primary tumors, with Caov4, COV413A, OVCA432, OVCAR3 and OVKATE representing the most highly correlated cell lines (Figure 5A). This study represents the first instance of subtype assignment of OVCA432, demonstrating its high similarity to HGSOC tumors. To validate that our ranking of HGSOC lines was consistent with both Domcke et al. (2013) and Yu et al. (2019), we first calculated the mean correlation of all cell lines with the 56 primary HGSOC tumors. We filtered this to include only the 39 cell lines analyzed by both Domcke et al. (2013) and Yu et al. (2019), and ranked cell lines in order of mean correlation with HGSOC tumors. Our suitability ranking was highly correlated with that observed by both Yu et al. (2019) (Spearman’s rho = 0.89, p-value <2.2e-16) and Domcke et al. (2013) order (Spearman’s rho = 0.58, p-value = 9.592e-05). A lower rho value was observed between our ranking, produced using gene expression data only, compared to that produced by Domcke et al. (2013), which took into account copy number alteration and mutational landscape. A higher rho value was observed between our ranking and that detailed in Yu et al. (2019), which was expected, as both studies also analyze transcriptomic data. However, in this study we analyze a completely independent set of 56 HGSOC tumors yet produce a similar suitability ranking.

Similar to HGSOC, cell lines classified here as CCOC possessed the highest ranking median correlations to CCOC tumors (Figure 5B). The 12 CCOC lines fell within the top 13 most correlated lines to CCOC tumors, with OVCAR4 (a HGSOC cell line) ranking 11th. The top 5 cell lines with highest mean correlation to CCOC tumors are RMGII, RMGI, OVISE, IGROV1 and OAW42. Although no mutational data is available for RMGII, its high correlation to CCOC tumors warrants its use as a CCOC model; however, mutational and copy number analysis on this line would be extremely useful to confirm this similarity. TOV21G and JHOC5 are also highly correlated to CCOC tumors, supporting recommendations for their use as CCOC models from Anglesio et al. (2013). It is important to note that Domcke et al. (2013) observed hypermutated genomes in TOV21G and IGROV1. Therefore, although these cell lines are highly correlated with CCOC tumors, caution should be taken when using these in experiments due to their hypermutated genomes. High correlation of OVCAR4 with CCOC tumors could be explained by the fact that this line harbors mutations that are seen in both HGOSC and CCOC (BRCA2, MYC amplifications), as well as those that are generally not observed in HGSOC (MET and PTEN amplifications, mutations in CDKN2A) (Figure 3). This high mutational overlap with CCOC may influence the high position of OVCAR4 in this ranking.

None of the six lines classified as ENOC ranked within the top 20 most correlated lines to ENOC tumors (Figure 5C). In fact, these cell lines were poorly correlated to tumors of HGSOC, CCOC and MOC. The most striking observation relates to A2780, a cell line reported to represent over 90% of citations in EOC studies (Domcke et al., 2013). A2780 has been previously reported to poorly model HGSOC tumors, and our results show that is a poor model of ENOC, the EOC tumor type it is expected to represent. Expression based clustering by Domcke et al. (2013) also showed that A2780 clusters closer to cancers of non-ovarian origin (such as lung, liver, stomach and small intestine) than to those of ovarian origin. Our analysis bolsters this observation, demonstrating that not only does A2780 display poor correlation with HGSOC tumors, but of EOC tumors of multiple subtypes. Based on our evidence, we recommend that use of A2780 for in vitro EOC studies should be avoided, and importantly, there is a striking need to develop additional ENOC cell lines. We observe a similar pattern for other lines classified as ENOC. Cell lines classified as HSGOC and CCOC may represent the best available models of this subtype, as they possess the highest median correlations to ENOC primary tumors (Figure 5C). Indeed HGSOC tumors have been noted to be highly similar to ENOC tumors of higher grades (Gilks et al., 2008; Karnezis et al., 2013; Lim et al., 2016) and CCOC tumors display a high mutational overlap with characteristic ENOC genomic alterations (Supplementary Table S4).

All lines classified as MOC (apart from JHOM2B) fell within the top 20 most highly correlated lines to MOC primary tumors (Figure 5D). COV644, MCAS, OVCA420, OVCA429, and RMUGS were amongst the lines with highest similarity. Again, this is the first reported instance of subtype assignment for OVCA420 and OVCA429, showing that these lines possess high molecular similarity to MOC tumors, despite a lack of mutational and copy number data. We also note here that GTFR230, while classified in this study as LGSOC, originated from a stage IC MOC tumor (Supplementary Table S2). We observed that this line ranks 48th (out of 56 lines) in terms of median correlation to MOC tumors, much lower than any putative MOC cell line (Figure 5D). We therefore suggest that this cell line is more representative of LGSOC, due to clustering with LGSOC lines. Further studies are needed to confirm this.

Most cell lines classified as LGSOC (HEY, HEYA8, JHOM1, OVCAR8, 59M, OV7, and ES2) were poorly correlated to primary tumors of HGSOC, CCOC, ENOC, and MOC origin (Figures 3–5). Additionally, these lines possess the major genomic alterations of LGSOC (Figure 3), and have been suggested to represent the LGSOC subtype by Barnes et al. (2021). We therefore hypothesize that these lines may potentially represent models of LGSOC. However, a comparison of these lines to primary LGSOC tumors is needed to confirm this, highlighting the need for additional datasets to be generated for this rarer EOC subtype.

Overall, we show that correlation of gene expression patterns between EOC cell lines and primary tumors can identify models with high molecular similarity to specific subtypes and aid in subtype identification, with the use of techniques such as NMF. In summary, we generated a reference dataset of cell lines most representative of HGSOC, CCOC, and MOC (Figure 6). We also highlight potential LGSOC models and those that are unsuitable for use in subtype specific studies.