Exploring the classification of cancer cell lines from multiple omic views

Cancer cell lines and data pre-processing

Our study involved 1,019 cell lines from 31 previously established cancer types. The mRNA, miRNA, CNV, METHY and RPPA data were downloaded from the CCLE database for all cell lines (https://portals.broadinstitute.org/ccle/data) (Ghandi et al., 2019). The number of cancer cell lines and the cancer types involved were shown in Table 1.

For mRNA sequence data, we used RSEM values in gene level shared by CCLE database. We used miRNA expression data from CCLE for miRNA analysis. For DNA methylation, the promoter CpG data was used for clustering analysis. And reverse phase protein array data was downloaded for protein analysis. In parallel, we downloaded segmented copy number profiles from CCLE database for CNV analysis. This SNP6.0 arrays data was used as the input data for Gistic2.0 software (Mermel et al., 2011). Before pre-processing the data, we mapped segmented copy number to the chromosome arm level using Gistic2.0. This copy number variation by the chromosome arm level was the input data of CNV clustering analysis. Next, the following steps were performed to improve the dataset quality for single omics clustering.

1. For each omics dataset, cell lines with more than 20% features missing, and features with more than 20% cell lines missing were filtered out.

2. For each omics dataset, the missing data points were filled in using average imputations.

3. For mRNA and miRNA data, log2 (x + 1) (x is the value of mRNA and miRNA) transformation were performed before feature selection.

4. For mRNA and METHY data, only features in the top 5,000 in terms of variance were selected. For miRNA, RPPA and CNV data, all features were considered.

For integrated multiple omics clustering, 670 cell lines from 24 cancer types with complete five omics data were used after samples alignment and deletion of cancer types with few number of cell lines.

Single and multiple omics clustering of cell lines

For single omics dataset, we performed hierarchical clustering with different methods and distance measurements (Table 2). We used 30 clustering validity indices to select the optimal clustering number using the R package “NbClust” (version 3.0) (Charrad et al., 2014). The optimal parameters for function “NbClust” were set as follows: min.nc = 10, max.nc = 30, method = “average”.

The Similarity Network Fusion and Consensus clustering algorithm (SNF-CC), a method combining Similarity Network Fusion (SNF) and Consensus clustering (CC) together to take advantage of both for cancer type identification, was applied to integrate multiple omics data (Monti et al., 2003; Wang et al., 2014). The method was implemented by the R package “CancerSubtypes” (version 1.8.0) (Xu et al., 2017). The optimal parameters for function “ExecuteSNF.CC” from “CancerSubtypes” were set as follows: K = 20, alpha = 0.5, t = 20, maxK = 30, pItem = 0.8, reps = 500. The silhouette coefficient, a measurement of consistency of each object within clusters, also derived using the function “silhouette_SimilarityMatrix” from the R package “CancerSubtypes”.

The function “pamr.listgenes” from the “pamr” R package (version 1.56.1) was used to find the most suitable clustering features for the illustration of clustering by heatmap (Tibshirani et al., 2002).

Dominant cancer type and functional enrichment analysis

A hypergeometric distribution was used to calculate the P-value for cancer types in each cluster. Cancer types with a P-value <10⁻³ were chosen as the dominant types of each cluster. The −lg(P-value) represents the enrichment score that cancer types gathered in each cluster.

We explored the differential expressed genes (DEGs) among clusters using the R package “limma” (version 3.38.3) (Ritchie et al., 2015). Genes with adjusted P-value <0.05 were selected, and were further screened according to —log2 fold-change—>1. Finally, the above DEGs were fed into enrichment analysis with GO and KEGG terms using the R package “clusterProfiler” (version 3.10.1) (Yu et al., 2012). The significantly enriched pathways were identified using false discovery rate <0.05.

Feature contribution of integrated multiple omics clustering

We used the normalized mutual information (NMI), which was a measure of the interdependence between two random variables, to measure the contribution of each omics type feature. The function “rankFeaturesByNMI” in the R package “SNFtool” (version 2.3.0) were used to compute NMI (Liu et al., 2018a; Wang et al., 2014). Codes are provided in Script S1.

Tumor maps of cancer cell lines

We used the TumorMap website to create pan-cancer cell lines maps from the above integrated data. TumorMap is an interactive website for assisting in exploring high-dimensional and complicated omics data (https://tumormap.ucsc.edu/) (Newton et al., 2017). In TumorMap, samples are distributed on a hexagonal grid based on their similarity and rendered using Google’s Map technology. The distances were used as input to generate a 2D layout of the samples. We used features that NMI ranks top 20% to calculate Euclidean similarity between each cell line. The Euclidean similarity is equal to 1/(1 + Euclidean distance). All parameters in TumorMap were set default.

Clustering based on single omics data

We initially clustered cell lines based on each type of omics data, which were mRNA, miRNA, CNV, METHY and RPPA data. The optimal clustering numbers were set to 10 (Fig. 1 and Fig. S1).

In the hierarchical clustering result of 901 cell lines by mRNA (Fig. 1A, Table S1 and Fig. S1A), we found that one cluster was mainly formed from a single type of cancer (C7 [SKCM]). Additionally, hematopoietic lymphatic malignancies were separated into two clusters, (C6 [ALL-DLBC-MM] and C9 [LAML-LCML]). Cancer cell lines with histological similarity or proximity tended to group together. These include C2: pan-gastrointestinal [COAD/READ-STAD], C4: nervous system tumors [GBM-LGG and some SARC whose features were the same as them] and C8: pan-gynecological [OV-UCEC and other SARC cell lines]. KEGG enrichment analysis indicated that C1 was enriched in human cytomegalovirus infection, transcriptional misregulation in cancer, proteoglycans in cancer, TNF signaling pathway and NF-kappa B signaling pathway (Figs. S2A and S2B). Meanwhile, the cell lines in C1 were enriched in GO terms including reproductive system development and morphogenesis of embryonic epithelium (Figs. S2C–S2H).

In the clustering result of miRNA data of 879 cell lines (Fig. 1B, Table S2 and Fig. S1B), five clusters predominately contained a single cancer type (C2 [COAD/READ], C5 [STAD], C6 [HNSC], C8 [SKCM] and C10 [NB]). And tumors of the hematopoietic lymphatic system were distributed in two clusters (C4 [ALL-DLBC-LAML-LCML-MM] and C9 [ALL-DLBC-LAML]). The significant signature of these two clusters were high expression of has-miR-142-5p and has-miR-142-3p, which played an important role in lineage differentiation of hematopoietic cells (Sharma, 2017).

CNV data sorted at the chromosome arm-level for 897 cell lines were divided into 10 clusters through hierarchical clustering, four clusters mainly formed from a single cancer type (C3 [GBM], C5 [SKCM], C7 [PAAD] and C8 [SCLC]) (Fig. 1C, Table S3 and Fig. S1C). C8 was characterized by the deletion of chr3p and chr17p and the amplification of chr3q. This characterization had been reported in previous studies (Carter et al., 2017; George et al., 2015; Peifer et al., 2012). C1 and C10 were enriched for ALL, DLBC, LAML and LCML. We observed fewer alterations in C1 but more alterations in C10. For example, C10 was characterized by chr8, chr19 and chr6 copy number increase.

Among the unsupervised clustering result of 755 cell lines using METHY data (Fig. 1D, Table S4 and Fig. S1D), there was one cluster that virtually consisted of one cancer type (C7 [SKCM]). Meanwhile, hematopoietic lymphatic malignancies were still enriched in two clusters (C3 [DLBC-MM] and C6 [ALL-LAML-LCML]). Cancer cell lines originating from same organ often gathered in the same cluster, such as C2 [COAD/READ-STAD-PAAD], a group of digestive system cancers whose common features were the high expression of CNKSR1, FOLH1, ADGRG1, SMAD7, LRATD1 and MVP (Haffner et al., 2009; Ji et al., 2018; Kobayashi et al., 2006; Quadri et al., 2017; Slattery et al., 2010; Teng et al., 2017). Additionally, squamous morphology cancer cell lines aggregated by METHY patterns (C10 [ESCA-HNSC]), particularly in terms of ARHGDIB and SEPTIN9 loss (Bennett et al., 2008).

In hierarchical clustering of RPPA data from 854 cell lines (Fig. 1E, Table S5 and Fig. S1E), C4 [SKCM], C8[LUAD] and C10 [BRCA] mostly contained one cancer type. The characteristics of C10 in this analysis had high level of ER, GATA3, AR, ERBB2, FASN, PREX1, CDH1 and CLDN7, and low level of CAV1 (Barrio-Real et al., 2016; Neve et al., 2006; Taherian-Fard, Srihari & Ragan, 2015). Hematopoietic lymphatic malignancies were enriched in one cluster (C3 [ALL-DLBC-LAML-LCML-MM]). Consistent with METHY analysis, a pan-gastrointestinal carcinoma cluster with COAD/READ-STAD was gathered in C2, which had high level of CDH1, CLDN7 and TYRO3, and a low level of CAV1 (Burgermeister et al., 2007; Di Bartolomeo et al., 2016; Qin & Qian, 2018). In addition, C7 was an enrichment cluster of squamous morphology cancer cell lines, mostly made of HNSC and ESCA and was characterized by high level of CDH1, CLDN7 and CAV1 (Ando et al., 2007; Bello et al., 2008; Shah et al., 2009).

Interestingly, some cancer types, such as SKCM, were individually classified in all five omics data, whereas some cancer types such as SCLC and BRCA were clustered individually only in one or two omics data. This result indicates that exposed information of each omics data is different at molecular level. According to investigation, SKCM is a paradigm of invasive cancer characterized by the highest mutational frequency among all cancer types and a large accumulation of changes in transcriptome (Cancer Genome Atlas N, 2015; Lawrence et al., 2013). Compared to other cancer cell lines, there are a large number of special molecular characteristics in SKCM cell lines. For example, we found that the levels of miR-188-3p and miR-514 were increased significantly, whereas in other cell lines, the levels of the two miRNA were decreased. At the CNV level, amplification of chr7 was found in most SKCM cell lines. It is generally known that several common mutations of SKCM, are on the chromosome 7 (Hayward et al., 2017). Moreover, pan-cancer clusters could be found based on mRNA, METHY and RPPA, but were not individually clustered in miRNA and CNV. This phenomenon was consistent with feature contributions of integrated multiple omics clustering, and was related to the fact that the characteristic information from mRNA, METHY and RPPA dataset were more representative.

Integrated clustering based on multiple omics data

By using SNF-CC, we integrated all five of omics datasets (mRNA, miRNA, CNV, METHY and RPPA) across 670 cell lines and identified 18 clusters (Fig. 2A, Table S6 and Fig. S3).

For these 18 clusters, 12 of them were dominated by a single cancer type (C1 [SCLC], C2[GBM], C4 [ALL], C6 [SARC], C7 [BRCA], C8 [ALL], C10 [MM], C12 [DLBC], C14 [LAML], C15 [SKCM], C16 [NB] and C17 [KIRC]) (Figs. 2A and 2B, Fig. S4). And each clusters also mixed with few amounts of other cancer types. Except SKCM, C15 also contained one glioblastoma multiforme cell line (LN229) with low level of VHL and high expression of has-miR-146a, has-miR-29b and has-miR-188-3p (Aurich, Fleming & Thiele, 2017). It is notable that although SARC is the dominant cancer types in C6, the proportion within the cluster is relatively low.

There were six clusters that dominated by two cancer types (C3 [PAAD-LUAD], C5 [HNSC-ESCA], C9 [LAML-LCML], C11 [COAD/READ-STAD], C13 [LUAD-LIHC] and C18 [OV-UCEC]) (Figs. 2A and 2B, Fig. S4). On the one hand, the proportion of two dominant cancer types were almost equal in C3, C13 and C18. And C3 was characterized by high level of CDH1. On the other hand, in C5, C9 and C11, one of the two dominant cancer type was over 50%. And C9 had high levels of VAV1 and STAT5A and low level of CTNNB1 (Bertagnolo et al., 2011; Harir et al., 2007; Ysebaert et al., 2006).

Three pan-cancer clusters influenced by organ of origin or cell morphology pattern were obtained. These clusters included pan-gastrointestinal cluster (C11 [COAD/READ-STAD]), pan-gynecological cluster (C18 [OV-UCEC]) and pan-squamous morphology carcinoma cluster (C5 [HNSC-ESCA]). For pan-gastrointestinal cluster (C11 [COAD/READ-STAD]), KEGG enrichment analysis results showed that these cell lines shared down-regulated in cytokine-cytokine receptor interaction and TNF signaling pathway (Figs. 3A and 3B and File S1). And cell lines in C11 had high levels of protein binding involved in heterotypic cell–cell adhesion, SNARE binding and G-protein beta-subunit binding in GO terms (Figs. 3C–3H and File S1). Meanwhile, pan-squamous morphology carcinoma cell lines (C5 [HNSC-ESCA]) were characterized by up-regulated of nicotine addiction and cell adhesion molecules pathway. And these cell lines had high levels of CAV1, EGFR and ITGA2 (Ando et al., 2007; Song et al., 2015).

We also observed two clusters with the same cancer type dispersed. For instance, cell lines from ALL were divided into two clusters, C4 and C8, despite the common characteristics such as Human T-cell leukemia virus 1 infection, Th17 cell differentiation, and TNF signaling pathway. The ALL cell lines in C8 were enriched in KEGG terms including up-regulated in ECM-receptor interaction, down-regulated in antigen processing and presentation pathway, while the ALL cell lines in C4 had low level of cellular senescence (Figs. 3A, 3B and File S1). GO enrichment analysis results showed that the ALL cell lines in C8 had low level of cell–cell junction and high levels of calcium channel activity, while the ALL cell lines in C4 were down-regulated in growth factor receptor binding and sulfur compound binding (Figs. 3C–3H and File S1). At other four omics levels, the features of these two clusters were different as well. For example, the levels of PTEN (a tumor suppressor gene), LCK and Syk (two immune-related genes) and has-miR-151-5p (related to tumor invasion and metastasis) was completely inconsistent.

Integrated multiple omics clustering provided a global view of cancer types because it could capture both shared and complementary information from each omics data. Several cancer types which mixed together in one single omics data were divided in other single omics data or integrated omics data. For example, BRCA and SCLC were mixed together based on miRNA data, but they were separated into two distinct molecular clusters based on integrated omics data. Besides, in single omics clustering, three pan-organ system clusters were only found based on mRNA data and the pan-squamous morphology carcinoma cluster was only found based on METHY and RPPA data. But pan-gastrointestinal cluster, pan-gynecological cluster and pan-squamous morphology carcinoma clusters were simultaneously identified by integrated multiple omics clustering.

The relative contribution of each omics data to the integrated clustering was computed based on the NMI value. On the basis of the top 20% statistical features from the five omics data, we found that RPPA and mRNA contributed 32.24% and 29.62% respectively, followed by METHY (16.24%) (Fig. 2C and Table 3). This result demonstrated that mRNA and proteomics data were particular important for cancer molecular classification. Meanwhile, more information was showed based on mRNA and RPPA data than other omics data in single omics clustering. For instance, pan-organ system clusters were identified based on mRNA and RPPA data, but not in miRNA and CNV. This results indicated that mRNA and proteomics data could be preferred if multiple omics data were not able to be measured simultaneously.

The comparison of classification between cancer samples and cell lines

We compared the classification results of 19 cancer types shared by cancer cell lines from CCLE and patient samples from TCGA (Hoadley et al., 2018). Clusters of patient samples and cell lines were divided into three types respectively, namely clusters dominated by single cancer type, pan-cancer clusters and clusters mixed with other cancer types (Table 4).

For hematopoietic lymphatic malignancies, the classification of cancer cell lines is more abundant than patient samples. For example, some LAML cell lines were clustered together in a group, while others were mixed with LCML cell lines into another group (LAML-LCML) in our findings. For patient samples, there is only one LAML group (Hoadley et al., 2018). The classification of DLBC cell lines was consistent with patient samples. Just like hematopoietic malignancies, the SARC patient samples were clustered individually into a group. However, for cell lines, except gathering in a single group, a few other SARC cell lines were mixed with GBM.

For most solid tumors, the classification of patient samples is generally more abundant and diverse than the corresponding cell lines. In general, patient samples with same cancer type can be divided into multiple groups, while cell lines with same cancer type are clustered in one group. For example, the samples of breast cancer were classified into three subgroups (chr8q amp, HER2 amp and Luminal). In addition, there were a large number of BRCA samples gathered with other cancer types in a mixed cluster. Except a few of BRCA cell lines were mixed in a pan-gynecological cluster, whereas almost all BRCA cell lines were clustered in a single group. And there are similar clustering results in gastrointestinal cancer and squamous cell carcinoma. There were two pan-gastrointestinal groups and a single gastric cancer group in patient samples (Hoadley et al., 2018). While most colorectal and gastric cancer cell lines were clustered in a pan-gastrointestinal cluster without forming a STAD group in the classification result of cell lines. Most patient samples of ESCA were divided into a pan-squamous morphology carcinoma group and a pan-gastrointestinal group, while in our study, most ESCA and HNSC cell lines were clustered in a pan-squamous morphology carcinoma cluster. This indicates that not all molecular subtypes of patient samples can be represented by current panels of cancer cell lines.

On the other hand, while the number of clustering groups were same, the types of clusters were different in cell lines and patient samples of some cancer types. For example, there is a pan-kidney group in the classification of patient samples. But this cluster was not in our results because there is only one cell line associated with kidney cancer, renal clear cell carcinoma (KIRC), in the CCLE database (Hoadley et al., 2018; Ricketts et al., 2018). Although there was no pan-kidney cluster in our research, the KIRC cell lines were clustered individually in our clustering results. And the patient samples of OV and UCEC were gathered in different clusters, while cell lines of the two cancer types were clustered together in a pan-gynecological group. GBM patient samples were clustered with LGG to form a pan-cancer group, however due to lacking LGG cell lines, GBM cell lines were clustered individually. The classification of patient samples in some cancer types, such as LIHC, LUSC and THCA, the clusters they formed were dominated by single cancer type (Hoadley et al., 2018). However, the clusters formed by these three types of cancer cell lines were mixed with large number of other cancer types.

In general, due to the larger numbers of patients and the greater heterogeneity in TCGA samples, the classification results of patient samples are more diverse and abundant than cancer cell lines. It is clear that the existing cancer cell lines do not fully represent all the molecular types of corresponding patient samples. However, in some cancer types, the classification of patient samples and cell lines were consistent. This shows that cancer cell lines can represent primary samples to some extent.

The TumorMap landscape of pan-cancer cell lines

We used TumorMap web tool to visualize the landscape of pan-cancer cell lines. The same layout and four different color schemes (SNF-CC cluster, TCGA disease, Pan-organ system and histology) were used to reveal that most cancer cell lines gathered based on organ systems and histopathological similarity (Fig. 4). More nuance within a cancer type were apparent. The SARC was widely distributed and separated into three parts. Most of them enriched in C6 (SARC), while others were fell into C2 (GBM) characterized by amplification of chr7p and C16 (NB) characterized by chr17q amplification (Figs. 4A and 4B). Additionally, major cell lines from STAD were assembled in C11 (COAD/READ-STAD), but few gathered in C18 (OV-UCEC) and C13 (LUAD-LIHC). Pan-organ system clusters and pan-squamous morphology carcinoma cluster reported previously were shown on the map (Fig. 4C) (Berger et al., 2018; Campbell et al., 2018; Liu et al., 2018b). We found that cell lines within C11 (COAD/READ-STAD) and C5 (HNSC-ESCA) were tightly gathered, while cell lines within C18 (OV-UCEC) were relatively dispersed. The TumorMap landscape showed that cancer cell lines with similar histology characterization tended to get together, even though histological information were not used during calculating similarities (Fig. 4D). The hematopoietic lymphatic malignancies were remote from other cancer types on the map. This result underscored that the molecular characteristics of hematopoietic lymphatic malignancies were different from other cancer types (Fig. 4D). Moreover, C15 (SKCM) and C17 (KIRC) were also far away from other solid tumor groups on the map.

We downloaded the drug susceptibility data for 24 anticancer drugs across 504 cell lines in CCLE database. We used TumorMap web tool to analyze the relationship between the drug susceptibility and the pan-cancer clustering. We divided the analysis results into four types and chose the representative drugs as examples (Fig. S5).

1. These anticancer drugs have a strong effect on almost all cancer cell lines. For example, LUAD-LIHC (C13), GBM (C2), SARC (C6), DLBC (C12), SKCM (C15), OV-UCEC (C18) and HNSC-ESCA (C5) cell lines are sensitive to Paclitaxel, a broad-spectrum anticancer drug (Fig. S5A).

2. Almost all cell lines are not sensitive to these drugs, for example, L-685458, a gamma-secretase inhibitor (Fig. S5B).

3. Only one cell line or few cell lines are sensitive to these anticancer drugs. For example, as a BRAF inhibitor, PLX4720 has an obvious effect on some SKCM (C15) cell lines, but has no effect on other cell lines (Fig. S5C).

4. These anticancer drugs have a strong effect on many cancer cell lines, but have a weak effect on others. For example, RAF256 is a dual inhibitor of mutant BRAF and vascular endothelial growth factor receptor 2. This drug can inhibit proliferation of SKCM (C15), GBM (C2), LUAD-LIHC (C13), OV-UCEC (C18) and some COAD/READ-STAD (C11) cell lines (Fig. S5D).