
==== Front
Commun Biol
Commun Biol
Communications Biology
2399-3642
Nature Publishing Group UK London

6812
10.1038/s42003-024-06812-3
Article
Identifying cell lines across pan-cancer to be used in preclinical research as a proxy for patient tumor samples
http://orcid.org/0000-0003-0842-8768
Bose Banabithi banabithi.bose@northwestern.edu

12
http://orcid.org/0000-0002-4813-4310
Bozdag Serdar serdar.bozdag@unt.edu

3456
1 https://ror.org/000e0be47 grid.16753.36 0000 0001 2299 3507 Center for Genetic Medicine, Feinberg School of Medicine, Northwestern University, Chicago, IL USA
2 grid.16753.36 0000 0001 2299 3507 Department of Pharmacology, Feinberg School of Medicine, Northwestern University, Chicago, IL USA
3 https://ror.org/00v97ad02 grid.266869.5 0000 0001 1008 957X Department of Computer Science and Engineering, University of North Texas, Denton, TX USA
4 https://ror.org/00v97ad02 grid.266869.5 0000 0001 1008 957X Department of Mathematics, University of North Texas, Denton, TX USA
5 https://ror.org/00v97ad02 grid.266869.5 0000 0001 1008 957X BioDiscovery Institute, University of North Texas, Denton, TX USA
6 https://ror.org/00v97ad02 grid.266869.5 0000 0001 1008 957X Center for Computational Life Sciences, University of North Texas, Denton, TX USA
7 9 2024
7 9 2024
2024
7 11015 1 2023
30 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
In pre-clinical trials of anti-cancer drugs, cell lines are utilized as a model for patient tumor samples to understand the response of drugs. However, in vitro culture of cell lines, in general, alters the biology of the cell lines and likely gives rise to systematic differences from the tumor samples’ genomic profiles; hence the drug response of cell lines may deviate from actual patients’ drug response. In this study, we computed a similarity score for the selection of cell lines depicting the close and far resemblance to patient tumor samples in twenty-two different cancer types at genetic, genomic, and epigenetic levels integrating multi-omics datasets. We also considered the presence of immune cells in tumor samples and cancer-related biological pathways in this score which aids personalized medicine research in cancer. We showed that based on these similarity scores, cell lines were able to recapitulate the drug response of patient tumor samples for several FDA-approved cancer drugs in multiple cancer types. Based on these scores, several of the high-rank cell lines were shown to have a close likeness to the corresponding tumor type in previously reported in vitro experiments.

A computational tool has been developed to compute similarity scores, determining how well different cell lines resemble patient tumor samples in various cancer types by integrating multi-omics datasets at genetic, genomic, and epigenetic levels.

Subject terms

Data integration
Cancer genomics
https://doi.org/10.13039/100000057 U.S. Department of Health & Human Services | NIH | National Institute of General Medical Sciences (NIGMS) R35GM133657 Bozdag Serdar issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Drug testing on human subjects in cancer research creates major logistical and ethical concerns. Multiple anti-cancer drug trials on a real cancer patient are extremely dangerous to the patient’s health and survival. As a result, scientists commonly use human tumor-derived cell lines as models for cancer tissues in pre-clinical trials to better understand drug efficacy1. The fundamental benefit of employing a cell line for research is that it is eternal. Cell lines are inexpensive, and they can be cloned and cultured for long periods of time, allowing an experiment to be repeated several times. In vitro culture of cancer cell lines, on the other hand, frequently results in genetic alterations, resulting in systematic phenotypic variations from patient tumor samples2,3. Hence, pre-clinical drug response models based on cell lines must be carefully evaluated before being translated into a clinical setting1,2.

Owing to the rapid technological advancements in high throughput technology, massive quantities of genome, proteome, transcriptome, and epigenome datasets, known as multiple “omes” or multi-omics, are now available to researchers all over the world. Several studies generated multi-omics datasets, such as gene expression, mutation, copy number aberration (CNA), DNA methylation, and microRNA (miRNA) expression from the large set of biospecimens3–5. Among these, multi-omics profiles of tumor samples and cancer cell lines in the Cancer Genome Atlas (TCGA)3 and Cancer Cell Line Encyclopedia (CCLE)4 databases pave the way for a comparative study of the patient tumor samples with cancer cell lines. To this end, using gene expression data, Liu et al.6 used a Transcriptome Correlation (TC) analysis to correlate cell lines and tumor samples. Sandberg et al.7 computed similarity scores, namely, Tissue Similarity Index (TSI) between a cell line and a patient cohort using Singular Value Decomposition (SVD) with gene expression of tumor samples and cell lines. Recently, Warren et al.8 developed an unsupervised alignment method called Celligner that mapped the gene expression of tumor samples to the gene expression of the cell lines. Their study found an alignment of the majority of cell lines with tumor samples of the same cancer type while revealing systematic differences in others.

All these studies relied on a genome-wide correlation-based approach between cell lines and bulk tumor samples without considering the heterogeneous immune cell population present in bulk tumors9,10, that play roles in targeted therapies in cancer11–13. Also, these studies considered single omics profiles, i.e., only gene expression data to relate patient tumor samples and cell lines, whereas several studies have shown that cell lines mirror many aspects of the multi-omics signatures identified in patient tumors14,15. To date, studies leveraging the multi-omics profiles to find systematic differences and similarities between patient tumor samples and cancer cell lines to compute a personalized similarity score between an individual patient tumor sample with an individual cell line have not yet been devised. Furthermore, no previous computational approach was developed to consider the critical role of the biological pathway activities in computing cell line-tumor similarity scores16 despite the knowledge that pathway activation status in cancer cell lines might connect tumor samples17,18. In a conference, we presented a pipeline (hereafter, CTDPathSim1.0) to compute the similarity index between breast and ovarian patient tumor samples and cell lines using DNA methylation and gene expression data19.

Recent studies indicate that CNAs are valuable as prognostic and predictive biomarkers for drug response in cancer patients20. Given the variability in copy number profiles between cell lines and tumor samples20,21, it is crucial to incorporate the effect of CNAs when calculating similarity indices between patient tumor samples and cell lines.

Previously, we introduced CTDPathSim1.0 as a conceptual tool focused on DNA methylation and gene expression data19. Building on this foundation, we now introduce CTDPathSim2.0, a significant enhancement over its predecessor. This upgraded version broadens the analytical scope to include DNA methylation, gene expression, and now CNAs, offering a more holistic omics data analysis. Further, CTDPathSim2.0 has been rigorously tested across 22 cancer types, showcasing its increased utility and reliability in cancer research.

To our knowledge, no large-scale study has yet adopted a personalized approach to systematically map the relationship between each patient tumor sample and cell line across a comprehensive array of genomic, genetic, and biological data on a pan-cancer scale. In this study, we integrate CNA data with DNA methylation and gene expression to identify suitable cell lines as proxies for primary tumor samples in vitro experiments. CTDPathSim2.0, our computational pipeline, leverages these three omics datasets to calculate tumor sample-cell line similarity scores across various cancers. This method offers a more dependable means of identifying cell lines that closely resemble actual tumor samples compared to existing approaches.

We performed an extensive analysis of our pipeline on 22 different cancer types and computed 7,358,104 similarity scores between 7228 patient tumor samples from TCGA and 1018 cancer cell lines from CCLE capturing a wide range of tumor types. Based on these scores, several of the high-rank cell lines were shown to closely resemble the correct tumor type in previously published in vitro tests. Our results show that CTDPathSim2.0 outperformed the state-of-the-art gene expression-based methods in capturing the drug response concordance between patient tumor samples and cell lines for several FDA-approved cancer drugs in several cancer types. Also, CTDPathSim2.0 outperformed all these tools in the selection of cell lines belonging to the same tissue type in several cancer types. Our study provides a way to assess the proper translation of the findings of pre-clinical cell line-based drug screening assays to clinical settings in cancer. Furthermore, we presented similarity scores between cell lines and patient tumor samples that could capture known cancer types and subtypes in patient tumor samples. For the use and replication of our technology, we presented the pipeline of CTDPathSim2.0 as an R software package. In this study, our approach was tested on cancers, however, the proposed methodology could be used to uncover cell lines related to other diseases.

Results

Computing cell line-tumor similarity score using multi-omics data

Considering the current gap in the understanding of similarities and differences between primary tumors i.e., patient tumor samples and cancer cell lines, in this study, we integrated DNA methylation, gene expression, and CNA datasets to compute multi-omics-based similarity scores between patient tumor samples and cell lines. Specifically, we utilized 7228 patient tumor samples in 22 different cancer types from TCGA and 1018 cancer cell lines from CCLE and computed their pairwise similarity score. For TCGA, we provided the cohort sizes of each data modality of the 22 cancer types in Table 1. For CCLE, we retrieved gene expression of 1019 cell lines, DNA methylation of 831 cell lines, and CNA of 916 cell lines belonging to these 22 cancer types. We followed a systematic data processing approach to ensure the quality of the TCGA multi-omics datasets (see “Methods” section).Table 1  Cohort sizes in different data modalities and cell types after deconvolution in 22 TCGA cancer types in the study

TCGA cancer type name (Abbreviation)	Expr	Meth	CNA	Cell types	
Adrenocortical carcinoma (ACC)	79	79	77	MC, VE, AC	
Bladder Urothelial Carcinoma (BLCA)	411	411	408	CD4+ T, VE, AC	
Breast invasive carcinoma (BRCA)	1073	1073	1064	CD4+ T, VE, MC	
Cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC)	302	302	290	CD4+ T, CN, VE, AC	
Lymphoid Neoplasm Diffuse Large B cell Lymphoma (DLBC)	47	47	46	CD4+ T, NK, B	
Esophageal carcinoma (ESCA)	160	160	159	VE, CN, MC, AC, CD4+ T	
Glioblastoma multiforme (GBM)	120	120	144	CN, AC, MC	
Head and Neck squamous cell carcinoma (HNSC)	495	495	490	AC, CD4+ T, VE	
Kidney renal clear cell carcinoma (KIRC)	490	293	486	CD4+ T, VE, AC	
Acute Myeloid Leukemia (LAML)	133	133	131	NK, NP, MC, B	
Brain Lower Grade Glioma (LGG)	492	492	491	VE, CN, MC, AC	
Liver hepatocellular carcinoma (LIHC)	371	371	369	CD4+ T, CN, AC, VE	
Lung adenocarcinoma (LUAD)	508	508	509	AC, VE, C CD4+ T, CN	
Lung squamous cell carcinoma (LUSC)	491	491	490	AC, CD4T, CN, VE	
Mesothelioma (MESO)	85	85	85	VE, AC, CD4+ T	
Ovarian serous cystadenocarcinoma (OV)	368	594	358	CN, VE, AC, MC	
Pancreatic adenocarcinoma (PAAD)	146	146	146	VE, CD4+ T, CN	
Prostate adenocarcinoma (PRAD)	400	400	397	CD4+ T, VE, AC	
Skin Cutaneous Melanoma (SKCM)	103	103	103	CD4+ T, AC, CE	
Stomach adenocarcinoma (STAD)	369	369	368	VE, AC, CD4+ T	
Thyroid carcinoma (THCA)	468	468	466	CN, CD4+ T, AC	
Uterine Corpus Endometrial Carcinoma (UCEC)	546	546	538	CD4+ T, CN, AC	
Expr Gene expression, Meth DNA methylation, CNA Copy Number Aberration.

‘Expr’ stands for gene expression; ‘Meth’ stands for DNA methylation; ‘CNA’ stands for copy number alteration. Cell types are the different cell populations such as B cell (B), NK natural killer, CD4+ T, MC monocytes, AC adipocytes, CN cortical neurons, and VE vascular endothelial found to be present in bulk tumor DNA methylation data after applying deconvolution algorithm in step 1 of CTDPathSim2.0.

To address the fact that the gene expression and DNA methylation signal from bulk tumors could be confounded by the presence of heterogeneous cell population9,10, in the first two steps of CTDPathSim2.0 (Fig.1), we applied a deconvolution algorithm using a quadratic programming approach (see “Methods” section) utilizing different immune cell types, namely, B cell (B), natural killer (NK), CD4+ T, CD8+ T, monocytes (MC), adipocytes (AC), cortical neurons (CN), and vascular endothelial (VE) and computed deconvoluted DNA methylation and expression profiles of each patient tumor sample in each cancer type (Table 1). Our focus on immune cell types stems from the pivotal role that immune cells play in cancer, particularly in the context of targeted therapies11–13,20–23. Several studies have highlighted the significant impact of chemotherapeutic drugs on various immune cells. For example, cyclophosphamide has been shown to notably reduce T cells in cancer24. A treatment regimen comprising Docetaxel, Doxorubicin, and Cyclophosphamide significantly lowered CD4+ T cells in breast cancer25. Furthermore, a combination of cisplatin, bleomycin, etoposide, and granulocyte-macrophage colony-stimulating factor (GM-CSF) was observed to affect NK cell levels in testicular cancer, with these levels not fully recovering by the end of the first treatment cycle26. Inspired by this, our pipeline aimed to establish a connection between tumor samples and cancer cell lines through immune-related signals. We anticipated that immune cell-specific deconvoluted omics profiles of patient tumor samples would be able to capture the immunobiological signals of cancer patients.Fig. 1 Flowchart of the CTDPathSim2.0 pipeline.

To evaluate the impact of the deconvolution steps, we conducted a comprehensive analysis (see Supplementary Note S1) to evaluate the correlation between the original data and the data after deconvolution for both DNA methylation and gene expression. We observed strong correlations between the original data and the deconvoluted data signifying that the deconvolution method effectively preserves the immune cell-type-specific signals present in the bulk tumor data while efficiently reducing undesired noise. The deconvolution step here acted as a filtering step of the bulk tumor data without altering it significantly. Furthermore, using a simulated breast cancer dataset that includes experimentally validated proportions of different cell types, including immune cells, we also presented an evaluation of the deconvolution method used in this study (see Supplementary Note S1).

In the third and fourth steps of CTDPathSim2.0 (Fig. 1), considering the critical role of the biological pathway activities in cancer16 and the previous studies that suggested that pathway activation status in cancer cell lines could connect tumor samples17,18, we computed enriched biological pathways for the patient tumor samples and cancer cell lines utilizing patient-specific and cell line-specific differentially expressed (DE), differentially methylated (DM) and differentially aberrated (DA) genes (see “Methods” section). We excluded one cell line that had no enriched pathway from this analysis. In the final step (Fig. 1), utilizing DE, DM, and DA genes, we computed Spearman correlation coefficients to get gene expression-, DNA methylation-, and CNA-based similarity scores for each tumor sample-cell line pair (see “Methods” section). We scaled these three similarity measures to the range of 0–1 using min-max normalization and computed an average similarity score for each sample-cell line pair in each cancer type. Since we did not have DNA methylation and copy number data for all the cell lines, there were sample-cell line pairs that had only expression-based scores. We also computed the similarity scores only using DNA methylation and gene expression data (hereafter, we will call this approach CTDPathSim1.0) in 22 cancer types to compare with the similarity scores from CTDPathSim2.0. For this, we considered expression-based similarity and DNA methylation-based similarity to compute an average similarity score for each cell line-sample pair. We applied z-normalization to all the computed final scores to set the mean score to 0 and the standard deviation to 1 in each TCGA cancer type. For running each of the steps of our pipeline, we provided our tool as an R software package. For the detailed pipeline with algorithms and R functions, please follow the Methods section of the manuscript.

CTDPathSim2.0 outperformed the existing methods for computing similarity scores that recapitulate the known drug responses between the tumor samples and the cell lines

To examine if the similarity scores are in concordance with the drug response of TCGA samples and CCLE cell lines, we computed a concordance score. The scores computed by CTDPathSim2.0 were compared to the scores computed by CTDPathSim1.0 and three state-of-the-art methods, namely TSI, TC analysis, and Celligner (see “Methods” section). We chose conditional density plots to assess the concordance of drug response with similarity scores obtained from five different approaches. We examined only those drugs and similarity scores that fulfilled the four filtering criteria given in the “Methods” section, using z-normalized scores for all the methods (Fig. 2). Except for Celligner (which had no filtered drug for ACC), we retrieved 24 total drugs in 19 distinct cancer types common to TCGA samples and CCLE cell lines (see Supplementary Table S1) using these filters on the generated scores by the various approaches. Figure 3 displays the conditional density plots for known drug responses vs. similarity scores of five different drugs from five methodologies in five cancer types. The conditional density plots visually represent the percentage of tumor sample-cell line pairs with matching (i.e., same drug response) and mismatching (i.e., opposite drug response) drug responses based on the similarity scores from five different methods. “1” represents the dark gray area indicating the matching percentage, while “0” represents the light gray area indicating the mismatching percentage.Fig. 2 Computed similarity scores in different cancer types by different methods.

Z-normalized scores between CCLE cell lines and TCGA samples of 22 different cancer types computed by A CTDPathSim2.0, B CTDPathSim1.0, C TSI, D TC analysis, and E Celligner. The number of sample-cell line pairs for ACC: 80,422; BLCA: 418,398; BRCA: 1,092,314; CESC: 307,436; DLBC: 47,846; ESCA: 162,880; GBM: 122,160; HNSC: 503,910; KIRC: 498,820; LAML: 135,394; LGG: 500,856; LIHC: 377,678; LUAD: 517,144; LUSC: 499,838; MESO: 86,530; OV: 374,624; PAAD: 148,628; PRAD: 407,200; SKCM: 104,854; STAD: 375,642; THCA: 476,424; UCEC: 555,828.

Fig. 3 Drug response concordance with sample-cell line similarity scores.

Conditional density plots of match and mismatch rate of drug response between TCGA samples and CCLE cell lines with computed z-normalized similarity scores by five different methods: A Gemcitabine in BLCA for 87 samples with 355 cell lines, B Doxorubicin in GBM 20 samples with 356 cell lines, C 5-Fluoruoracil in ESCA 12 samples with 361 cell lines, D Paclitaxel in LUSC 28 samples with 106 cell lines and E Erlotinib in LUAD 20 samples with 94 cell lines. The matching percentage of the response of a drug increases as the similarity score increases in CTDPathSim2.0 in all these cancer types.

For visual assessment, we categorized the drugs based on the trend observed in the plots. If the matching percentage of the response increased as the similarity score increased, we labeled it as “concordant”. Conversely, if the matching percentage decreased as the similarity score increased, we labeled it as “discordant”. A linear trend, where there was no increase or decrease, was labeled as “linear” while all other trends were labeled as “undecided”. In Fig. 3, CTDPathSim2.0 shows the drug response concordance for all these five drugs in all five cancer types. The conditional density plots for all other drugs in 19 cancer types for the five methods with a detailed drug response trend summary from visual inspection are in Supplementary Note S2. There were 86 different conditional density plots for 24 different drugs in 19 different cancer types, with the concordance trend obvious in 42, 32, 23, 19, and 22 plots for CTDPathSim2.0, CTDPathSim1.0, TSI method, TC analysis, and Celligner, respectively. These results indicate that similarity scores computed by our method align with known drug responses between the tumor samples and the cell lines.

To provide a numerical evaluation, we calculated concordance scores (both micro and macro) for each drug by examining drug response concordance using computed similarity scores in high and low similarities. We defined high and low similarities as the top 20 and bottom 20 percentiles of the scores, respectively. We selected these percentiles as they were found to be the most suitable for depicting conditional density plots of drug response matching and mismatching. We presented a detailed analysis behind this choice of threshold in Supplementary Note S3. Because of the skewed distribution of TSI scores (see Fig. 2C), to allow enough data points for our comparison, we used high and low-similarity thresholds of above 0 and below 0, respectively. The drugs were labeled as “concordant” if the difference between response matching percentage in the top and bottom percentiles was ≥0.10, “discordant” if that was ≤−0.10, and “undecided” otherwise. In each cancer type, we computed two different concordance scores (i.e., micro concordance scores): (1) C1 = DConDTot and (2) C2 = DConDCon+DDis; where DTot denotes the total number of drugs, DCon denotes the number of concordant drugs and DDis denotes the number of discordant drugs. Table 2 displays the micro concordant scores for five different methods in 19 different cancer types. For CTDPathSim2.0, CTDPathSim1.0, TSI method, TC analysis, and Celligner, the number of cancer types with the highest micro concordant scores with C1 was eight, six, three, five, and four, respectively. Similarly, for CTDPathSim2.0, CTDPathSim1.0, TSI method, TC analysis, and Celligner, the number of cancer types with the highest micro concordant scores with C2 was eight, five, three, five, and four, respectively. We computed macro-averaged and weighted macro-averaged concordance scores using C1 and C2 independently to measure the overall effectiveness of the five different methods on pan-cancer (Table 3). We considered eighteen distinct cancer types omitting ACC in these scorings because there were no filtered drugs with computed Celligner scores in ACC. For each method, the macro-averaged concordance score = ∑cConcCanTot, where c denotes the cancer type, Conc denotes the micro concordance score (i.e., C1 or C2) in cancer type c and CanTot denotes the number of cancer types; the weighted macro-averaged concordance score = ∑cDcConcDTot, where c denotes the cancer type, Conc denotes the micro concordance score (i.e., C1 or C2) in cancer type c, DTot denotes the total number of drugs across all cancer types and Dc denotes the number of drugs in cancer type c. In comparison to all other methods, CTDPathSim2.0 achieved the greatest macro-averaged and weighted macro-averaged concordant scores with both C1 and C2. CTDPathSim1.0 achieved the second-highest weighted macro-averaged concordance score with C1 and C2. CTDPathSim1.0 outperformed TSI and TC analysis for macro-averaged with C1 and C2 as well. However, in terms of macro-averaged scores, Celligner slightly outscored CTDPathSim1.0. Overall, CTDPathSim2.0 outperformed CTDPathSim1.0 in terms of reproducing known drug response concordance between sample-cell line pairs; CTDPathSim1.0 outperformed the TSI technique, TC analysis, and Celligner, which were all in the same ballpark.Table 2  The micro concordant scores for five different methods in 19 different cancer types with C1 outside of the parentheses and C2 inside the parentheses

Cancer type	CTDPathSim2.0	CTDPathSim1.0	TSI method	TC analysis	Celligner	
ACC	1 (1)	0 (0)	0 (NA)	0 (0)	NA (NA)	
BLCA	0.3 (0.4)	0.6 (0.7)	0.2 (0.5)	0.3 (0.4)	0.2 (0.4)	
BRCA	0.4 (0.7)	0.4 (0.5)	0.4 (0.5)	0.3 (0.3)	0.2 (0.2)	
CESC	0.6 (0.8)	0.5 (0.8)	0.4 (0.4)	0.2 (0.2)	0.7 (0.8)	
DLBC	0.8 (0.8)	0.4 (0.7)	0.2 (0.2)	0.8 (0.8)	0.3 (0.5)	
ESCA	0.7 (1)	0.3 (0.5)	0.7 (1)	0.3 (0.3)	0 (0)	
GBM	0 (NA)	0 (0)	0 (NA)	0 (NA)	1 (1)	
HNSC	0 (0)	0.6 (0.8)	0.5 (0.8)	0.3 (0.7)	0.2 (0.2)	
KIRC	0.2 (0.2)	0.5 (0.5)	0.2 (0.2)	0.5 (0.6)	0.2 (0.2)	
LGG	0 (0)	0 (0)	0 (0)	0 (0)	1 (1)	
LIHC	0.8 (0.8)	0.7 (1)	0.2 (0.3)	0.2 (0.2)	0 (0)	
LUAD	0.2 (0.3)	0.1 (0.1)	0.1 (0.2)	0.5 (0.8)	0.4 (0.4)	
LUSC	0.6 (0.6)	0.6 (0.7)	0.3 (0.5)	0.1 (0.3)	0.3 (0.5)	
MESO	0.2 (0.3)	0 (0)	0 (0)	0.2 (0.3)	0.8 (1)	
OV	0.5 (1)	0 (0)	0 (0)	0.5 (0.5)	0 (0)	
PAAD	0.3 (0.3)	0.7 (1)	0 (0)	0 (0)	0.5 (1)	
SKCM	0.5 (0.5)	0.5 (0.7)	0.5 (0.7)	0.2 (0.2)	0.8 (1)	
STAD	0.3 (0.3)	0.3 (0.3)	0.4 (0.5)	0.5 (0.6)	0.1 (0.2)	
THCA	1 (1)	0 (NA)	1 (1)	0 (0)	0 (0)	
The highest values for each cancer type are shown in bold.

NA stands for not applicable. If there was no filtered drug for the test then NA appears in the case of both C1 and C2; if there is no concordant and discordant drug then NA appears only in the case of C2.

Table 3  The macro-averaged and weighted macro-averaged concordant scores with C1 and C2

Score	CTDPathSim2.0	CTDPathSim1.0	TSI method	TC analysis	Celligner	
Macro (C1)	0.44	0.33	0.27	0.26	0.37	
Weighted macro (C1)	0.50	0.40	0.26	0.29	0.29	
Macro (C2)	0.56	0.46	0.40	0.34	0.47	
Weighted macro (C2)	0.63	0.52	0.38	0.39	0.38	
The highest values in each score are shown in bold.

CTDPathSim2.0 outperformed state-of-the-art methods in the selection of significant cell lines belonging to the same tissue type in several cancers

We expected that the cell lines that have high similarity to a primary tumor sample would be significantly enriched in the cell lines derived from that tissue type, whereas low-similarity cell lines would not have such enrichment. To test which similarity calculation method met this expectation, we picked CCLE cell lines that had high and low similarity to a primary tumor sample and performed enrichment of these CCLE cell lines in the cell lines derived from that tissue type in 21 different TCGA cancer types using all five methods. Since there were no Adrenocortical Carcinoma-derived cell lines in CCLE, we removed ACC samples from this analysis. We performed this evaluation using a high-similarity threshold ≥1 and a low-similarity threshold ≤−1 for the sample-cell line pairs for all the methods, except for TSI method. These thresholds were chosen to include a sufficient number of sample-cell line pairs for all different methods to run this analysis. Because of the skewed distribution of TSI scores (see Fig. 2C), to allow enough data points for our comparison, we used high and low-similarity thresholds of ≥0 and ≤0, respectively. In these thresholds, we computed highly and lowly similar cell lines for each patient tumor sample and checked if the cell lines for each patient tumor sample had highly overlapping tissue-specific cell lines by running hypergeometric test for each patient tumor sample. For instance, we checked if the highly similar cell lines to LUSC samples had significantly overlapping lung cancer cell lines. We found that highly similar cell lines to primary tissue samples computed by CTDPathSim2.0 had significantly lower hypergeometric p-values than hypergeometric p-values for lowly similar cell lines. Figure 4 shows the distribution of hypergeometric p-values (i.e., enrichment density plot) for corresponding tissue-specific cancer cell lines in high and low-similarity thresholds for all five methods in GBM, UCEC, KIRC, LUAD, and LUSC samples in panels A–E, respectively. In CTDPathSim2.0, for all five cancer types (Fig. 4; Column 1), the tissue-specific cell lines were picked with significant hypergeometric p-values (<0.05) in high-similarity thresholds compared to low-similarity thresholds and hypergeometric p-values were separated within two groups. CTDPathSim1.0 was also able to separate the hypergeometric p-values within two groups in four cancer types (Fig. 4; Column 2) except GBM. The TSI method (Fig. 4; Column 3) produced highly overlapping p-values in all five cancer types. For the Celligner (Fig. 4; Column 5) significant p-values appeared to be present in low-similarity groups with higher frequencies than in high-similarity groups. TC analysis (Fig. 4; Column 4) performed well in separating the two groups in four cancer types except GBM, however, there are several insignificant hypergeometric p-values (>0.05) in the high-similarity group in KIRC. The enrichment density plots for all cancer types with different thresholds can be accessed from Supplementary Note S4.Fig. 4 Density plots of hypergeometric p-values of enrichment for known tissue-specific cancer cell lines for the five methods.

Enrichment of highly (cyan) and lowly (red) similar cell lines in A 33 GBM-specific cell lines; B 28 UCEC-specific cell lines; C 33 KIRC-specific cell lines; D 76 LUAD-specific cell lines and E 23 LUSC-specific cell lines. The t-test (alternative = less”) p-values comparing each pair of the density plots in each cancer type are shown in red. NA stands for not applicable in cases where a t-test could not be run due to insufficient data points between two groups.

We compared the density plots of hypergeometric p-values for the enrichment of known tissue-specific cancer cell lines in 21 cancer types for five different methods using a one-sided t-test (Table 4). The significant t-test (alternative = “less”) p-value (p-value < 0.05) shows that the enrichment density plot in the high-similarity group has a lower mean compared to the low-similarity group. There were 20, 17, 14, 16, and 16 distinct cancer types with significant t-test p-values for CTDPathSim2.0, CTDPathSim1.0, TSI method, TC analysis, and Celligner, respectively. Overall, CTDPathSim2.0 outperformed all other methods in this test, with CTDPathSim1.0 coming in second, the TC analysis and Celligner being comparable, and the TSI approach coming in last.Table 4 One-sided t-test p-values for comparing the density plots of hypergeometric p-values for the enrichment of known tissue-specific cancer cell lines in 21 cancer types for five different methods

Cancer type	CTDPathSim2.0	CTDPathSim1.0	TSI method	TC analysis	Celligner	
BLCA	0	0	0	0	1	
BRCA	0	0	0	0	1	
CESC	NA	NA	6.7e-75	NA	NA	
DLBC	6.2e-202	0	1	NA	NA	
ESCA	8.9e-248	0	1	2.3e-02	NA	
GBM	0	1	1	NA	NA	
HNSC	0	0	0	6.1e-06	NA	
KIRC	0	0	1	0	1	
LAML	3.5e-18	NA	1	NA	NA	
LGG	0	0	9.9e-01	NA	1	
LIHC	1.1e-123	0	1	0	1	
LUAD	0	0	2.6e-28	0	1	
LUSC	0	0	0	0	1	
MESO	4.8e-81	7.3e-166	3.8e-21	4.1e-17	NA	
OV	0	0	0	0	1	
PAAD	0	0	0	0	0.98	
PRAD	4.1e-52	0	0	1.8e-258	1	
SKCM	2.3e-115	1	3.1e-01	1.2e-86	1	
STAD	1.5e-225	0	6.4e-116	0	1	
THCA	1.3e-243	2.6e-06	3.9e-19	2.9e-22	1	
UCEC	0	0	5.6e-47	0	1	
The significant t-test (alternative = “less”) p-value (p-value < 0.05) shows that the enrichment density plots (Fig. 4 and Supplementary Note S4) for the high-similarity groups have a lower mean compared to the low-similarity groups. NA stands for not applicable in cases where a t-test could not be run due to insufficient data points between two groups.

CTDPathSim2.0 tended to assign a high rank to tissue-specific cell lines among all the other cell lines in CCLE

In this analysis, for CTDPathSim2.0, we checked if the tissue-specific cell lines in each cancer type were in higher ranks among 1018 cell lines from CCLE. We utilized the similarity scores computed in 21 different cancer types for which there was at least one tissue-specific cell line among 1018 CCLE cell lines and computed the rank of 1018 cell lines based on the median similarity scores in each cancer type. In each cancer type, we additionally assigned an empirical p-value to these median rankings. For example, in BLCA, we had 25 BLCA-specific cell lines with a median rank of 710, hence, we picked 25 cell lines at random and computed their median rank to see if they had a better (i.e., higher) median rank than 710. We repeated this procedure 1000 times and computed the empirical p-value as the percentage of cases where random cell lines had a better median rank. The number of tissue-specific cancer cell lines and their median rank in each cancer type are shown in Table 5 (here, the higher the rank better the result). Tissue-specific cell lines were shown to have a high rank (median of median rank in 21 cancer types was 738) among 1018 CCLE cell lines for all of these cancer types. There were ten cancer types with an empirical p-value < 0.05 and sixteen cancer types with an empirical p-value ≤ 0.1. Figure 5 displays several of the tissue-specific cell lines to be in higher rank among all the other cell lines in eight cancer types. The plots for all other cancer types can be accessed from Supplementary Fig. S1. The rank of each tissue-specific cell line in 21 cancer types can be accessed from Supplementary Note S5.Table 5 Number of tissue-specific cancer cell lines in 1018 CCLE cell lines and their median rank with an empirical p-value in the parentheses based on the median similarity values to the samples in each cancer type

Cancer type	Cell line	Median Rank (p-value)	Cancer type	Cell line	Median Rank (p-value)	
BLCA	25	710 (0.02)	LIHC	25	601 (0.16)	
BRCA	51	608 (0.08)	LUAD	76	728 (0.00)	
CESC	3	827 (0.08)	LUSC	23	663 (0.06)	
DLBC	39	985 (0.00)	MESO	9	749 (0.06)	
ESCA	27	629 (0.11)	OV	44	506 (0.53)	
GBM	33	655 (0.03)	PAAD	40	684 (0.00)	
HNSC	33	744 (0.00)	PRAD	8	680 (0.11)	
KIRC	33	785 (0.00)	SKCM	54	677 (0.00)	
LAML	35	993 (0.00)	STAD	37	573 (0.21)	
LGG	10	646 (0.17)	THCA	11	577 (0.33)	
	UCEC	28	807 (0.00)	
Here, the higher the rank better the result.

Fig. 5 The rank of all 1018 CCLE cell lines based on median similarity scores with all samples in the corresponding cancer type computed by CTDPathSim2.0 in eight cancer types.

The rank of all cancer cell lines in A BRCA; B DLBC; C GBM; D LUAD; E LIHC; F LUSC; G OV; H UCEC. Each dot is a CCLE cell line and the few representatives of highly-ranked tissue-specific cell lines are marked by red.

CTDPathSim2.0 computed similarity scores between cell lines and tumor samples were able to capture known cancer types and subtypes

To further evaluate CTDPathSim2.0, we tested whether the computed similarity scores between CCLE cell lines and primary tumor samples could detect known cancer types and subtypes present in the tumor samples. We used the cancer type and subtype annotation data from8 focusing on situations where there were matched cancer type and subtype annotations across cell lines and tumor samples. We conducted k-means clustering of patient tumor samples based on their cell line-specific similarity scores (see the “Methods” section). We also conducted k-means clustering of the same samples using multi-omics datasets such as gene expression, DNA methylation, and CNA. Then we compared these clustering results with real annotations (Table 6). We also compared the resulting clusters using similarity scores with those derived from the original multi-omics datasets (Table 7). The comparisons were evaluated using the Rand Index (RI) in Table 6, a measure of the similarity between two data clustering groups, with values closer to 1 indicating higher concordance.Table 6 Comparative analysis of clustering concordance using patient-cell line similarity scores and three types of omics data with the known cancer type and subtype labels

Data type	Analysis level	Number of clusters	Number of patients	Genes as features	Cell lines as features	RI	
Similarity scores	Cancer types	19	7228	NA	1018	0.89	
Gene expression	Cancer types	19	7228	32,807	NA	0.79	
DNA methylation	Cancer types	19	7228	9490	NA	0.91	
CNA	Cancer types	19	7228	32,807	NA	0.85	
Similarity scores	Subtypes	28	7228	NA	1018	0.92	
Gene expression	Subtypes	28	7228	32,807	NA	0.89	
DNA methylation	Subtypes	28	7228	9490	NA	0.95	
CNA	Subtypes	28	7228	32,807	NA	0.88	
Data Type: Type of analysis data; Analysis Level: Analysis focus, cancer types or subtypes; Number of Clusters: Clusters identified; Number of Patients: Patient samples analyzed; Genes as Features: Gene count for omics data that used for clustering; Cell Lines as Features: Cell line count for similarity scores that used for clustering (see Table 9); RI (Rand Index): Concordance measure, with values closer to 1 indicating higher similarity.

NA not applicable.

Table 7 Concordance between the clusters from patient-cell line similarity scores and original multi-omics data

Data type	Analysis level	Number of clusters	RI	
Similarity vs. Gene expression	Cancer types	19	0.79	
Similarity vs. DNA methylation	Cancer types	19	0.89	
Similarity vs. CNA	Cancer types	19	0.85	
Similarity vs. Gene expression	Subtypes	28	0.87	
Similarity vs. DNA methylation	Subtypes	28	0.93	
Similarity vs. CNA	Subtypes	28	0.87	
Data Type: Type of analysis data; Analysis Level: Analysis focus, cancer types or subtypes; Number of Clusters: Clusters identified; RI (Rand Index): Concordance measure, with values closer to 1 indicating higher similarity.

In Table 6, the RI values indicate a high concordance between our cell line-specific similarity scores and the known annotated labels for both cancer types and subtypes. Notably, the similarity scores outperform or are comparable to traditional omics data analyses in several instances, particularly for DNA methylation in subtype clustering, where the similarity scores approach closely matches the highest RI value of 0.95. Similarly, while comparing the clusters from cell line-specific similarity scores vs. original omics datasets, we got high RI values indicating their robust concordance (Table 7).

Furthermore, we ran the Random Forest classifier using cell lines as features with the computed similarity matrix to classify the patient tumor samples and assessed this classification based on the known cancer type and subtype labels. This analysis yielded a high classification accuracy of 0.97 for cancer types and 0.90 for subtypes. In our analysis, we also computed the macro-averaged F1 scores to evaluate the performance of our classification model for both cancer types and subtypes. The F1 score is a harmonic mean of precision and recall, providing a balance between the model’s ability to correctly label positive cases and its capacity to not mislabel negative cases. For cancer types, we achieved a macro-averaged F1 score of 0.96, and for subtypes, the score was 0.88 (Supplementary Note S6). These high scores indicate that our model was highly effective in accurately classifying tumor samples according to their respective cancer types and subtypes.

A notable aspect of our results pertains to the HER2-enriched subtype. We observed that 67% of samples in this category were classified as “luminal A” and 22% as “basal,” resulting in an F1 score of “NA” (not applicable) for this subtype. This outcome aligns with recent research27, which suggests that HER2-enriched breast cancer encompasses a diverse range of transcriptional subtypes. This diversity includes luminal A and basal-like subtypes, with the HER2-enriched subtype demonstrating a unique transcriptional landscape, independent of HER2 amplification. The “NA” score in our study reflects this inherent heterogeneity within the HER2-enriched subtype, underscoring the complexity and varied nature of these cancer categories.

Similarly, for the “melanocytic” subtype, 70% of the samples were classified as “transitory,” leading to an F1 score of “NA”. This finding is consistent with existing literature28 that suggests a close relationship between the “transitory” and “melanocytic” subtypes, both of which represent a nuanced differentiation within the proliferative phenotype of cancer cells.

These results established the efficacy of CTDPathSim2.0 in computing similarity scores that accurately reflect the phenotypic characteristics of patient tumor samples. By correlating these samples with specific cell lines, our tool not only demonstrates its precision in classification but also provides valuable insights into the complex nature of cancer subtypes. The cancer type and subtype classification evaluation metrics such as precision, recall, and F1 scores can be accessed from Supplementary Note S6. The details of k-means and Random Forest analysis are provided in the “Methods” section.

Pathway-specific DE/DM/DA gene analysis reveals cancer-associated biological processes

To assess the relevance of the pathway-specific DE/DM/DA genes with cancer-related biological processes, we further evaluated the pathways from the driver DEs whether they are involved in important biological and cancer-related processes or not. Given the variability in patient pathways (i.e., differential pathways), we consolidated the pathways across four cancer types, namely ACC, CESC, GBM, and LUAD for this analysis (see Supplementary Note S7). Among the most common pathways (top ten from each of these cancers) from which driver DE genes were identified, the “Signal Transduction” pathway is associated with regulating cell growth, division, death, fate, and motility, fitting within the broader context of cancer progression through disturbances in signaling networks, including tumor microenvironment alterations, angiogenesis, and inflammation29,30. The “Immune System” pathway highlights the role of inflammatory immune cells in tumors and the pivotal contribution of such cells to cancer-related inflammation over the past two decades31. Additionally, the “Gene Expression (Transcription)” pathway is recognized for its significance in cancer32.

Discussion

We developed a computational pipeline, CTDPathSim2.0, that utilizes a pathway activity-based approach in computing similarity between primary tumors and cell lines at genetic, genomic, and epigenetic levels. We applied CTDPathSim2.0 on patient tumor samples from TCGA and cancer cell lines from CCLE for 22 different cancer types. We rigorously compared the scores computed by CTDPathSim2.0 to those generated by CTDPathSim1.0 and state-of-the-art methods: TSI, TC analysis, and Celligner.

One comparison criterion was based on the drug response concordance between TCGA tumor samples and CCLE cell lines. Our findings revealed that as the similarity score from CTDPathSim2.0 increased, the drug response concordance also increased consistently for all five drugs across five cancer types (Fig. 3). This indicates an alignment between tumor-cell line similarities computed by our method and actual drug responses.

After evaluating the performance across the board, CTDPathSim2.0 achieved the highest macro-averaged and weighted macro-averaged concordant scores (Table 3). CTDPathSim1.0 achieved the second-best performance in weighted macro-averaged concordant scores.

Our hypothesis was that cell lines with a high similarity to primary tumor samples should predominantly belong to the same tissue type. To test this, we assessed the enrichment of such cell lines across 21 TCGA cancer types using all methods. The results, using a one-sided t-test, showed that CTDPathSim2.0 consistently outperformed all other methods, followed by CTDPathSim1.0 (Fig. 4 and Table 4).

For analyses involving only gene expression and DNA methylation data, both CTDPathSim2.0 and its predecessor, CTDPathSim1.0, perform equivalently. The CTDPathSim2.0 R package facilitates using either version based on the data scenario. Our analysis suggests that pipelines including CNA data alongside gene expression and DNA methylation likely produce more accurate similarity indices, especially in multi-omics analyses of various cancer types, demonstrating the value of including CNA data. Furthermore, the tumor samples are inherently heterogeneous due to the presence of multiple cell types. This complex heterogeneity presents both challenges and opportunities when trying to establish similarities with cancer cell lines. Our method is distinctly designed to address this heterogeneity. Unlike other methods that generalize bulk tumor samples, our tool actively decodes the intricacies of this heterogeneity when determining similarity scores with cancer cell lines.

We also evaluated the link between DE, DM, and DA genes and cancer-specific pathways. By analyzing DE genes within patient-specific biological pathways as driver genes, and extending this approach to cell lines, we established a methodology that emphasizes the biological relevance of genes in cancer progression. Our analysis covered multiple cancer types, identifying critical pathways such as Signal Transduction, Immune System, and Gene Expression, and employed a pathway activity-based selection criterion to ensure the biological significance of genes used in computing similarity scores. All these DE/DM/DA genes along with patient-specific and cell-line-specific pathways from 22 different cancer types can be accessed from the accompanied datasets. Utilizing the CTDPathSim2.0 R package, our users will be able to access and analyze these differential genes and pathways for their specific research needs. The package vignette included detailed instructions on how to access and interpret this pathway information.

Using a median similarity score computed by CTDPathSim2.0 across the same tissue-specific tumor samples, we calculated the rank of tissue-specific CCLE cell lines. For example, each lung tissue-derived cell line with a median similarity score was taken for TCGA lung cancer (LUSC) samples. Table 8 shows the top five tissue-specific cell lines in fifteen different cancer types having more than ten tissue-specific cell lines. Several of the high-rank cell lines were shown to have a close likeness to the corresponding tumor type in previous in vitro experiments. For instance, the HCC1599 cell line was reported among highly resembling cell lines for the breast cancer samples with estrogen-positive, progesterone-positive, and human epidermal growth factor receptor 2-positive (HER2+) group21. A study reported CAOV3 as one of the most resembling cell lines for ovarian cancer33. DKMG cell line was found to have EGFR gene amplification, rearrangement, and expression, similar to GBM cancer samples34. In several studies, RAJI was a popular choice as a proxy for DLBC tumors35,36. We also checked the tissue-specific cell lines that were low-ranked in the corresponding tumor samples. For DLBC, the REC-1 cell line had the lowest rank. An in vitro study found that the CD44 gene was hypermethylated and transcriptionally silenced in all DLBC cell lines, however, the same gene was found to be unmethylated and expressed in REC-1 cell lines37. We found that for BRCA, the UACC-893 cell line had the lowest rank. This cell line was reported to resemble the brain “metastasis” than the primary tumor in HER2+ breast cancer38. A study analyzed 57 breast cancer cell lines and found that the UACC-893 cell line was the only cell line that had the NAT2 gene’s expression higher than the NAT1 gene’s expression compared to other cell lines39. The same study reported the KPL-1 breast cancer cell line to be an MCF-7 derivative. Interestingly, our method ranks MCF-7 and KPL-1 as 18 and 17, respectively among all 51 breast cancer cell lines used in the study. These findings point toward the ability of CTDPathSim2.0 for ranking the cell lines based on their genomic profiles concerning primary tumors. CTDPathSim2.0 also generated similarity scores between cell lines and patient tumor samples, which were able to capture known cancer types and subtypes, and these similarity scores can be used to refine current cancer and subtype classifications in tumor samples. Our tool enables the discovery of cell line models that best replicate the multi-omics properties of a tumor type, or even a specific tumor sample of interest, using biological pathway activation status. Our tool also provides separate scores for gene expression, DNA methylation, and CNA, in addition to the combined similarity score. Users can easily access these individual scores to assess the contribution of each data type to the overall similarity. This feature allows for a more tailored analysis, enabling researchers to focus on the most relevant aspect of similarity for their specific research question.

We anticipate that a deeper knowledge of the commonalities between cancer cell lines and tumors will enable better model selection and translation of the understanding of drug response from pre-clinical models to clinical samples.Table 8 Five top-ranked tissue-specific cell lines for each cancer type

Cancer type	Top five cell lines	
BLCA	5637, BC3C, CAL29, KU1919, RT4	
BRCA	BT474, HCC1143, HCC1599, HCC2157, HCC38	
DLBC	BL41, DOHH2, NUDHL1, OCILY3, RAJI	
ESCA	COLO680N, ECGI10, TE10, TE1, TE5	
GBM	DKMG, KG1C, KS1, LN18, SF268	
HNSC	BICR56, FADU, PECAPJ15, PECAPJ34CLONEC12, SCC15	
KIRC	A704, KMRC2, KMRC3, SLR23, SNU1272	
LAML	HL60, HNT34, KG1, KO52, MUTZ3	
LIHC	HUH1, HUH7, JHH5, JHH7, PLCPRF5	
LUAD	HCC4006, NCIH222, NCIH2291, NCIH3122, RERFLCKJ	
LUSC	KNS62, LC1, LUDLU1, NCIH1869, NCIH2170	
OV	59 M, CAOV3, JHOS2, JHOS4, OVTOKO	
PAAD	CAPAN1, MIAPACA2, PANC0213, SNU410, SU8686	
SKCM	HS839T, HS934T, MALME3M, SKMEL3, UACC257	
STAD	NCIN87, NUGC4, SNU5, SNU601, SNU668	

The methodology developed in our cancer research, designed to link patient tumors with laboratory-grown cell lines, has broader applications, including in other diseases like diabetes. For instance, using our computational pipeline, we can connect diabetic patients to relevant diabetes-specific cell lines, like the MIN6 pancreatic beta-cell line used for insulin production studies40. This approach allows for an assessment of how well these cell lines represent a patient’s multi-omics profile, which is crucial in testing and developing effective, targeted therapies. This not only enhances drug relevance and efficacy for diabetes but also illustrates the potential of our pipeline in creating personalized treatment strategies across various diseases.

Finally, we provided CTDPathSim2.0 as an R package to run the program for ease of application and reproducibility.

Methods

All the experiments were conducted in accordance with relevant guidelines and regulations.

Datasets

Data preprocessing for running CTDPathSim2.0

In CTDPathSim2.0, we integrated DNA methylation, gene expression, and CNA datasets from 22 different cancer types. We downloaded the genomic data of cancer patient samples from the TCGA database. Using TCGAbiolinks41, we downloaded RNA sequencing (RNA-Seq) data in HTSeq-FPKM format. We downloaded DNA methylation data from Infinium HumanMethylation27 Bead-Chip (27 K) and Infinium HumanMethylation450 Bead-Chip (450 K) platforms. We retrieved masked copy number variation (Affymetrix SNP Array 6.0) and computed the gene-centric copy number value compatible with hg38 genome assembly using R Bioconductor package CNTools42. From the COSMIC (Catalogue Of Somatic Mutation In Cancer)43 cancer gene census project, we downloaded a list of 702 cancer-driving genes that are found to be frequently mutated in different cancer types.

We downloaded gene-centric copy number data for 1040 cell lines from the CCLE database. We downloaded the CCLE_RNAseq_gene_rpkm _10180929.gct.gz file for gene expression for 1019 cell lines and CCLE_RRBS_cgi_CpG_clusters_20181119.txt.gz file for DNA methylation data of 831 cell lines. Using the R package BioMethyl44, we removed CpG sites that had missing values in more than half of the samples and imputed the rest of the missing values by employing the Enmix45 R package with default parameters.

To prepare the reference DNA methylation profiles for the deconvolution step, we processed raw methylation probes provided in FlowSorted.Blood.450k46 R package that consists of DNA methylation probes from peripheral blood samples with five different cell types, namely, B cell (B), natural killer (NK), CD4+ T, CD8+ T, and monocytes (MC), generated from adult men with replicates. We also processed raw DNA methylation data available at NCBI Gene Expression Omnibus (GEO) database repository with the dataset identifier GSE122126 for three more cell types, namely, adipocytes (AC) (for 450 K DNA methylation only), cortical neurons (CN), and vascular endothelial (VE) cells with replicates. We processed these raw methylation files using the R package minfi47 to prepare our reference DNA methylation profiles of different cell types. Two different reference files, one for 450 K probes and the other for 27 K probes, with eight different cell types were prepared. Please see Table 1 to find the sample sizes which were used after preprocessing in this work.

Data preprocessing for evaluating CTDPathSim2.0’s results

To compare the drug responses of patient tumor sample-cell line pairs, we compiled our ground truth drug response data for the patients from the TCGA database using Overall Survival (OS) and Progression-Free Survival (PFI) as clinical endpoints. These survival endpoints were considered representative of the drug response status of the patients. In OS, the patients who were dead from any cause were considered dead, otherwise censored. In PFI, the patients having a new tumor event whether it was a progression of the disease, local recurrence, distant metastasis, new primary tumor event, or died with cancer without a new tumor event, including cases with a new tumor event whose type is N/A were considered as the event occurred and all other patients were censored. We considered survival status 0 (i.e., censored), representative of sensitive drug response whereas survival status 1 (i.e., dead/event occurred), representative of resistant drug response as we mapped the patients’ survival status to drug response. In order to determine if PFI or OS status can be used as a drug response metric, we looked at the time intervals between drug applications. The following five rules were used to map the patients’ survival status to drug response (the details of these rules can be accessed from Supplementary Note S3):Drug response = OS status	; if PFI status = 0 and PFI.time = OS.time	
Drug response = PFI status	; if PFI status = 1 and PFI.time = OS.time	
Drug response = PFI status	; if drug application end days ≤ PFI days and PFI status = 1 with PFI. time < OS.time	
Drug response = OS status	; if drug application end days ≤ PFI days and PFI status = 0 with PFI.time < OS.time	
Drug response = OS status	; if drug application start days ≥ PFI days	

For cell lines’ drug response data, we considered IC50 (the concentration of a drug that is required for 50% inhibition in vitro) values as drug response quantification. Furthermore, ln(IC50)<−2.0 was used to define positive (i.e., sensitive) anti-cancer drug response, and other ln(IC50) values were considered as negative (i.e., resistant) drug responses. Since common drugs between TCGA and CCLE were few (only twelve), we also added cancer cell lines drug response data with IC50 values from the Genomics of Drug Sensitivity in Cancer (GDSC)48 database in our ground truth drug response data. We found 989 overlapping cell lines between the CCLE and the GDSC databases. We processed IC50 values of the GDSC drug similarly as we did for the cell lines in the CCLE. We retrieved 24 common drugs between TCGA patients of 19 different cancer types and cell lines from CCLE and GDSC after using the five filtering criteria mentioned below. A list of these drugs with the number of patient-cell line pairs, number of resistant and sensitive patients, and cell lines can be accessed from Supplementary Table S1.

To compare the drug responses between sample-cell line pairs according to the computed similarity scores, we used z-normalized scores and considered only those drugs and similarity scores that passed the following filtering criteria to be eligible for the drug response concordance test between sample-cell line pairs (shown in Table 9 and Table 10 as an example):Table 9 Selection of the drugs D1 and D2 based on match-mismatch ratio

Drug	No. of matching sample-cell line pairs	No. of mismatching sample-cell line pairs	Match-mismatch ratio	Drug status	
D1	5	30	0.17	Eliminated	
D2	10	20	0.5	Selected	

Table 10 Selection of the similarity scores for drug D based on hypergeometric test and match-mismatch difference

Drug	Sim	Universe	Set 1	Set 2	Overlap of Set 1 & Set 2	Hypergeometric p-value	Matching percentage	Mismatch percentage	
D	0.305	418398	5504	117	30	0.01964841	40	60	
D	0.506	418398	5504	140	20	0.03838658	60	40	
In this example, both similarity scores would be considered in drug response evaluation as they pass both the hypergeometric test (i.e., p-value < 0.05) and match-mismatch difference (difference ≥10%).

1. Filtering the drugs based on match-mismatch ratio: To check if a drug has a consistent response between sample-cell line pairs in the ground truth data, we computed a match-mismatch ratio, i.e., the ratio between the number of matched (i.e., same drug response) sample-cell line pairs to the number of mismatches (i.e., different drug response) sample-cell line pairs. The drugs that had a match-mismatch ratio < 0.2 were eliminated as even in the ground truth data, the majority of sample-cell line pairs have discordant drug responses. Table 9 shows an example of this filtering step.

2. Filtering the drugs based on similarity thresholds: We considered the drugs for which the sample-cell line pairs had similarity thresholds at least above or below 1 and −1, respectively to allow enough data points in high and low-similarity sample-cell line groups.

3. Filtering the computed similarity scores based on hypergeometric p-value: Since we computed millions of sample-cell line similarity scores, we had discrete similarity scores, each with multiple sample-cell line pairs after rounding the similarity scores to three digits. This rounding precision allowed us to generate a substantial number of data points for our evaluation. Decreasing the rounding threshold by two digits or increasing the rounding threshold by more than five digits would have led to the elimination of several “discrete scores.” Therefore, after careful consideration, we chose to use three digits for rounding. Furthermore, to assess the significance of each discrete similarity score, we implemented a stringent filter as described below:

For a given cancer type, we calculated the number of sample-cell line pairs for all similarity scores that had drug response data for drug D, forming Set 1. We also computed the number of sample-cell line pairs for a specific discrete similarity score, denoted as K, forming Set 2. To determine if K is a significant discrete similarity score, we conducted a hypergeometric test to assess the overlap between Set 1 and Set 2 (hypergeometric p-value < 0.05). The universe for this hypergeometric test encompassed all sample-cell line pairs of the cancer type under study.

This approach ensured that we considered only those discrete similarity scores that exhibited a significant overlap between Set 1 and Set 2, thereby enabling us to evaluate drug response match and mismatch between different sample-cell line pairs for drug D. Table 10 shows an example of this filtering step.

4. Filtering the computed similarity scores based on match-mismatch difference: We filtered the similarity scores that showed at least a ten percent difference between match and mismatch (i.e., the difference between the last two columns in Table 10) drug response within sample-cell line pairs for considering these as confident scores for the drug response concordance test (Table 10).

Running CTDPathSim2.0

CTDPathSim2.0 has five main computational steps. In the following paragraphs, we describe the CTDPathSim2.0 R functions to run these steps. We performed these steps for each of 22 different cancer types separately. The entire pipeline of CTDPathSim2.0 is illustrated in Fig. 1.

Step 1: Computing sample-specific deconvoluted DNA methylation profile

To capture the accurate DNA methylation signal of tumor samples, in the first step, we utilized a deconvolution-based algorithm49 to infer the sample’s deconvoluted DNA methylation profile with proportions of different cell types in the sample’s bulk tumor tissue. The deconvolution algorithm requires a reference DNA methylation profile of different cell types. Hence, we compiled a list of reference DNA methylation profiles with eight different cell types, namely, B, NK, CD4+ T, CD8+ T, MC, AC, CN, and VE with replicates. We considered these cell types since a blood-based DNA methylation profile can depict the distribution of immune cell types50,51, and these cell types play roles in targeted therapies in cancer11–13. This algorithm computes a matrix with DNA methylation profiles of each of the constituent cell types (i.e., W = (wij) where W ∈ Rl×c ; l is the number of loci, c is the number of cell types, i = 1, 2,…., l and j = 1, 2,…., c.) and a matrix with the proportions of constituent cell types in each input tumor sample (i.e., H = (hij) where H ∈ Rc×s ; c is the number of cell types, s is the number of tumor samples, i = 1, 2,…., c and j = 1, 2,…., s.)) minimizing the Euclidean distance between their linear combination (i.e., WH) and the original matrix of samples’ DNA methylation profiles (i.e., V = (vij) where V ∈ Rl×s ; l is the number of loci, s is the number of tumor samples, i = 1, 2,…., l and j = 1, 2,…., s.) (Fig. 1, Step 1). We designed a function RunDeconvMethylFun() that ran a minimization algorithm in an iterative procedure that, in each round, alternates between estimating constituent cell-type proportions and DNA methylation profiles by solving constrained least-squares problems of ∣∣V−WH∣∣F2 using quadratic programming method49 with the constraints that 0 ≤ vij ≤ 1, 0 ≤ wij ≤ 1 and ∑ihij=1. These three constraints for this problem were from the fact that DNA methylation measurements and cell type proportions were numbers in the range of 0 to 1 (hence, 0 ≤ vij ≤ 1 and 0 ≤ wij ≤ 1), and cell type proportions within a sample added up to 1 (hence, ∑ihij=1) (Fig. 1, Step 1). These constraints restrict the space of possible solutions, thus making it possible for the local iterative search to find a global minimum and an accurate solution reproducibly. To initiate the algorithm, we utilized the reference profiles as W with a set of selected 500 marker loci (probes) based on comparisons of each class of reference against all other samples using t-test. These loci were likely to exhibit variation in DNA methylation levels across different constituent cell types.

We designed a function, SampleMethFun() that utilized the estimated proportion of constituent cell types (i.e., HE) and the estimated probe-centric DNA methylation profile of constituent cell types (i.e., WE) from the output of RunDeconvMethylFun() to compute the sample-specific deconvoluted DNA methylation profile for each tumor sample (see Fig. 1, Step 1). We ran this step for 450 K DNA methylation probes (i.e., loci) and 27 K DNA methylation probes separately. For each sample-specific deconvoluted DNA methylation profile, we computed gene-centric DNA methylation beta values from the probe-centric DNA methylation values with ProbeToGeneFun() function in our R package. For the 450 K platform, the average beta value for promoter-specific probes was being considered due to their role in transcriptional silencing52. Given lower coverage in the 27 K platform, the average beta value of all the probes of a gene was considered as the gene’s DNA methylation level.

Step 2: Computing sample-specific deconvoluted expression profile

In this step, we computed the deconvoluted expression profile for each tumor sample. We computed cell-type-specific expression profiles for the reference cell types utilized in Step 1 minimizing the objective function ∣∣Y−ZH∣∣F2 where Y is the original matrix of tumor samples’ expression profiles (i.e., Y = (yij) where Y ∈ Rg×s, g is the number of genes, s is the number of samples, i = 1, 2,…., g and j = 1, 2,…., s.)), and Z is the cell-type-specific expression profiles (i.e., Z = (zij) where Z ∈ Rg×c, g is the number of genes, c is the number of cell types, i = 1, 2,…., g and j = 1, 2,…., c.) (Fig. 1, Step 2) to be estimated. We provided RunDeconvExpr() function in the package that utilized the estimated cell proportions (HE) from Step 1 as a fixed input (i.e., H = HE) and computed average gene expression profiles of constituent cell types through a constrained least-squares of ∣∣Y-ZH∣∣F2 using quadratic programming with solutions constrained to [0,∞) for gene expression measurements49. Then, using estimated cell-type-specific proportions in each sample (HE) and estimated cell-type-specific expression (ZE) in SampleExprFun() function, we computed the deconvoluted gene-centric expression profile of each tumor sample (Fig. 1, Step 2).

Step 3: Computing sample- and cell line-specific DE genes and biological pathways

In this step, we computed enriched biological pathways listed in REACTOME53 database for each sample and cell line using sample-specific DE genes. We designed GetSampleDE() function to compute the median expression of a gene across all samples and further computed the fold-change of that gene to the computed median expression. A gene was considered up if the fold-change ≥4 and down if the fold-change ≤−4. Likewise, using the GetCellLineDEFun() function and similar thresholds, we computed cell line-specific DE genes. Considering the critical role of the biological pathway activities in cancer16 and the previous studies that suggested that pathway activation status in cancer cell lines could connect tumor samples17,18, we identified enriched biological pathways for each sample and cell line using GetPathFun() function in our R package utilizing these DE genes. GetPathFun() found enriched biological pathways listed in the REACTOME database via a pathway enrichment tool Pathfinder54. To prioritize cancer-related biological pathways, specifically, we used the known frequently mutated genes that were considered cancer-driving genes from the DE gene list for this analysis. We obtained the frequently mutated cancer-driving genes from the COSMIC43 database. GetPath() considered the pathways that were significantly enriched with FDR-corrected hypergeometric p-values < 0.05 in this DE cancer gene list.

Step 4: Computing sample- and cell line-specific differentially methylated (DM) and differentially aberrated (DA) genes

In this step, we computed sample-specific and cell line-specific DM and DA genes. For sample-specific DM genes, we computed M-values55 from gene-centric beta values and utilized GetSampleDMFun() function in the R package to compute the median of gene-centric M-values across all the samples in each cancer type. This function considered a gene hypermethylated if the M-value fold-change ≥4 and hypomethylated if fold-change ≤−4. Likewise, using GetCellLineDMFun() with the same thresholds, we computed cell line-specific DM genes among all available cell lines in the study.

To compute sample-specific DA genes, we computed highly amplified and deleted genes utilizing chromosomal copy number profiles of each cancer type cohort as inputs in the GISTIC 2.056 tool in the GenePattern57 web server using a confidence interval of 0.90. We selected the genes that exceeded the high-level GISTIC thresholds for amplification and deletions as 2 and −2, respectively as DA genes.

For cell lines, we selected the highly amplified and deleted genes for each cell line compared with all other cell lines from various cancer types to capture the cancer-specific variation in copy number profiles in computing the similarity score between a patient sample and cell line. For this purpose, we utilized the GetCellLineDA() function, which selects cell line-specific DA genes based on Tukey’s mean-difference58 curve based on gene-centric copy number data of cell lines. To determine whether a gene is DA for a specific cell line, we calculated the mean of that gene’s copy number values across all cell lines and the difference between the copy number value of that gene in the cell line under consideration and the mean of the copy number values of that gene in the other cell lines. For each cell line, the highly amplified genes were chosen based on mean ≥0.5 and difference ≥1, and highly deleted genes were chosen based on mean ≥0.5 and difference ≤−1.

Step 5: Computing sample-cell line pathway activity-based similarity score

In this step, we computed Spearman rank correlation between each sample-cell line pair using the sample-specific deconvoluted expression, sample-specific deconvoluted DNA methylation, and sample’s gene-centric copy number values. All DM genes that occurred in an enriched pathway of a sample or cell line were used to compute the Spearman rank correlation between each sample-cell line pair via FindSim() function in the R package (Fig. 1, Step 5). Similarly, using DE and DA genes that were present in the gene list of the union of enriched pathways, we compute the Spearman rank correlation between each sample-cell line pair. Using DM and DE genes, we used the sample’s deconvoluted profiles to compute Spearman rank correlation, whereas for DA genes, we used the sample’s gene-centric copy number profile to compute Spearman rank correlation (Fig. 1, Step 5). We scaled the Spearman similarities to the range of 0 to 1 using min-max normalization and computed an average similarity score taking a mean of expression-based correlation, DNA methylation-based correlation, and copy number-based correlation for each sample-cell line pair. Since we did not have DNA methylation data and copy number data for all the cell lines, there were sample-cell line pairs that had only expression-based scores.

We presented the CTDPathSim2.0 pipeline in this section, but our R package can be used to compute similarity scores without using CNA data i.e., running CTDPathSim1.0. In that case, the FindSim() function of the package can be used to run the pipeline without DA genes (please check the vignette of the CTDPathSim2.0 in Supplementary Software 1 for these details).

Running state-of-the-art-methods

We compared CTDPathSim2.0 with our previous method CTDPathSim1.0 and three other state-of-the-art methods, namely, TSI, TC analysis, and Celligner by running them on each of 22 different cancer type cohorts from TCGA and cancer cell line datasets from CCLE. For CTDPathSim1.0, we computed only DNA methylation-based and gene expression-based similarity scores and averaged them to get a unified similarity score for each sample-cell line pair in a cohort.

Running TSI method

For the TSI method, we performed z-normalization of RNA-Seq data for each gene in the samples and cell lines and performed all the required steps to compute tumor sample-cell line similarity scores based on the TSI paper7. We applied SVD on the normalized gene expression matrix of samples. Then each sample was projected into SVD space by measuring its correlation to the first 16 Eigenarrays (i.e., left singular vectors) that explain the most variance of the data. Similarly, we projected each cell line into SVD space. We computed a Pearson correlation score between each patient and cell line as a similarity score, namely, the TSI score in the SVD space. We applied z-normalization on TSI scores computed between each cancer sample and cell line.

Running TC analysis

For TC analysis, we used 30,681 common genes between CCLE and TCGA RNA-Seq gene expression datasets with log2 transformation. Then, we rank-transformed gene RPKM (reads per kilobase of transcript per million reads mapped) values for each CCLE cell line and ranked all the genes according to their expression variation across all CCLE cell lines. The 1000 most variable genes were kept as “marker genes” to compute Spearman rank correlation, namely, TC scores, between each cell line and sample pair using the respective expression profiles.

Running Celligner

Considering Celligner’s global aligning approach of gene expression data between cell lines and tumor samples that depends on various cancer types and the composition of the dataset, we ran this pipeline as described in their paper utilizing the RNA-Seq data in TPM (transcript per million reads) of 12,236 tumor samples and 1249 cell lines as log2(TPM + 1) transformation from Cancer Dependency Map [https://depmap.org] data portal for 37 different cancer types. We used the gene expression data without ‘non-coding RNA’ and ‘pseudogene’ as input in the Celligner pipeline. Celligner first performed unsupervised clustering of each dataset and removed the expression signatures with excess intra-cluster variance in the sample compared to cell line data. Then, using DE genes between clusters in the data it ran contrastive Principal Component Analysis (cPCA) followed by a batch correction method, MNN (mutual nearest neighbors) that aligned similar sample-cell line pairs to produce corrected gene expression data. We further filtered this gene expression matrix for 22 different TCGA cancer types utilized in the study and 1015 CCLE cell lines that were common to CTDPathSim2.0. Then, we performed PCA on the Celligner-aligned data and took the Euclidean distance between each cell line and patient sample in PCA space (using the first 70 principal components following59) as a similarity measure.

Performing k-means clustering and Random Forest analysis

We filtered annotation labels with disease type (i.e., cancer type) and subtype for which at least ten tumor samples were available in TCGA. To display the clusters of patient tumor samples based on their similarity scores with 1018 cell lines, we used 7228 tumor samples with 19 distinct cancer type labels and 28 distinct cancer subtype labels. For clustering analysis, we utilized the patient-cell line similarity matrix. This matrix was constructed with similarity scores as entries, where each row represents a patient, and each column represents a cell line.

Using this matrix, we performed k-means clustering on patient tumor samples based on their cell line-specific similarity scores. We then compared these clusters with annotated cancer type groups i.e., 19 disease labels and 28 subtype labels. We used R base package Stats with k = 19 and k = 28 for getting similarity matrix-based clusters of two groups while setting all other parameters as default. We also compared the resulting clusters with those derived from the original omics data such as gene expression, DNA methylation, and CNA. Using the RI, a measure ranging from 0 (no agreement) to 1 (complete agreement), we quantified the similarity between two clustering groups.

For the Random Forest analysis, we treated each cell line as a feature and used the patient-cell line similarity scores as feature value. We partitioned our data, using 75% of the tumor samples (selected randomly) for training and the remaining 25% for testing. To train the model and tune the hyperparameters, we performed 10-fold cross-validation using the training data. Then the trained modeled was tested on the held-out test data. This approach allowed us to evaluate the model’s performance and adjust parameters as needed to optimize accuracy. The Random Forest classifier was implemented using the ‘randomForest’ package in R.

Statistics and reproducibility

The statistical analyses were conducted using publicly available datasets, as detailed in the datasets section. To ensure reproducibility of our results, we have released CTDpathsim2.0 as open-source R software and made all processed datasets60, and scripts available.

Supplementary information

Peer Review File

Supplementary Information

Description of Additional Supplementary Files

Supplementary Data 1

Supplementary Software 1

Supplementary information

The online version contains supplementary material available at 10.1038/s42003-024-06812-3.

Acknowledgements

This work was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number R35GM133657.

Author contributions

B.B. and S.B. conceived the study, B.B. conducted the study, S.B. supervised the study, B.B. developed the software, B.B. wrote the manuscript, B.B. and S.B. reviewed and edited the manuscript.

Peer review

Peer review information

Communications Biology thanks Yuri Kotliarov and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Toril Holien, Eve Rogers, Anam Akhtar, and Johannes Stortz. A peer review file is available.

Data availability

The datasets of this study can be accessed from Science Data Bank (10.57760/sciencedb.01713).

Code availability

The CTDPathSim2.0 pipeline was developed as an R package. For installation, use remotes::install_github(“boseb/CTDPathSim2.0”). The source codes of the package are available at https://github.com/bozdaglab/CTDPathSim2.0 under Creative Commons Attribution Non-Commercial 4.0 International Public License. The vignette of the R package and the scripts for running the evaluation results can be accessed from Supplementary Software 1 and Supplementary Data 1, respectively.

Competing interests

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Ben-David U Genetic and transcriptional evolution alters cancer cell line drug response Nature 2018 560 325 330 10.1038/s41586-018-0409-3 30089904
Ben-David, U. et al. Genetic and transcriptional evolution alters cancer cell line drug response. Nature 560, 325–330 (2018).30089904 10.1038/s41586-018-0409-3
2. Gillet J-P Varma S Gottesman MM The clinical relevance of cancer cell lines J. Natl Cancer Inst. 2013 105 452 458 10.1093/jnci/djt007 23434901
Gillet, J.-P., Varma, S. & Gottesman, M. M. The clinical relevance of cancer cell lines. J. Natl Cancer Inst. 105, 452–458 (2013).23434901 10.1093/jnci/djt007
3. Tomczak K Czerwińska P Wiznerowicz M The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge Contemp. Oncol. 2015 19 A68 A77
Tomczak, K., Czerwińska, P. & Wiznerowicz, M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp. Oncol. 19, A68–A77 (2015).
4. Barretina J The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity Nature 2012 483 603 607 10.1038/nature11003 22460905
Barretina, J. et al. The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature 483, 603–607 (2012).22460905 10.1038/nature11003
5. Yang W Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells Nucleic Acids Res. 2013 41 D955 D961 10.1093/nar/gks1111 23180760
Yang, W. et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961 (2013).23180760 10.1093/nar/gks1111
6. Liu K Evaluating cell lines as models for metastatic breast cancer through integrative analysis of genomic data Nat. Commun. 2019 10 1 12 30602773
Liu, K. et al. Evaluating cell lines as models for metastatic breast cancer through integrative analysis of genomic data. Nat. Commun. 10, 1–12 (2019).30602773
7. Sandberg R Ernberg I Assessment of tumor characteristic gene expression in cell lines using a tissue similarity index (TSI) Proc. Natl Acad. Sci. USA 2005 102 2052 10.1073/pnas.0408105102 15671165
Sandberg, R. & Ernberg, I. Assessment of tumor characteristic gene expression in cell lines using a tissue similarity index (TSI). Proc. Natl Acad. Sci. USA 102, 2052 (2005).15671165 10.1073/pnas.0408105102
8. Warren, A. et al. Global computational alignment of tumor and cell line transcriptional profiles. Nat. Commun. 10.1038/s41467-020-20294-x (2021).
9. Decamps, C. et al. Guidelines for cell-type heterogeneity quantification based on a comparative analysis of reference-free DNA methylation deconvolution software. BMC Bioinformatics 21, 16 (2020).
10. Grün, D. Revealing dynamics of gene expression variability in cell state space, Nat. Methods10.1038/s41592-019-0632-3 (2019).
11. Gun SY Lee SWL Sieow JL Wong SC Targeting immune cells for cancer therapy Redox Biol. 2019 25 101174 10.1016/j.redox.2019.101174 30917934
Gun, S. Y., Lee, S. W. L., Sieow, J. L. & Wong, S. C. Targeting immune cells for cancer therapy. Redox Biol. 25, 101174 (2019).30917934 10.1016/j.redox.2019.101174
12. Neal, J. T. et al. Organoid modeling of the tumor immune microenvironment. Cell 175, 1972–1988 (2018).
13. Tsou P Katayama H Ostrin EJ Hanash SM The emerging role of B cells in tumor immunity Cancer Res. 2016 76 5597 5601 10.1158/0008-5472.CAN-16-0431 27634765
Tsou, P., Katayama, H., Ostrin, E. J. & Hanash, S. M. The emerging role of B cells in tumor immunity. Cancer Res. 76, 5597–5601 (2016).27634765 10.1158/0008-5472.CAN-16-0431
14. Kao J Molecular profiling of breast cancer cell lines defines relevant tumor models and provides a resource for cancer gene discovery PLoS ONE 2009 4 e6146 10.1371/journal.pone.0006146 19582160
Kao, J. et al. Molecular profiling of breast cancer cell lines defines relevant tumor models and provides a resource for cancer gene discovery. PLoS ONE 4, e6146 (2009).19582160 10.1371/journal.pone.0006146
15. Neve RM A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes Cancer Cell 2006 10 515 527 10.1016/j.ccr.2006.10.008 17157791
Neve, R. M. et al. A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes. Cancer Cell 10, 515–527 (2006).17157791 10.1016/j.ccr.2006.10.008
16. Vogelstein B Kinzler KW Cancer genes and the pathways they control Nat. Med. 2004 10 789 799 10.1038/nm1087 15286780
Vogelstein, B. & Kinzler, K. W. Cancer genes and the pathways they control. Nat. Med. 10, 789–799 (2004).15286780 10.1038/nm1087
17. Mattes MJ Apoptosis assays with lymphoma cell lines: problems and pitfalls Br. J. Cancer 2007 96 928 936 10.1038/sj.bjc.6603663 17342089
Mattes, M. J. Apoptosis assays with lymphoma cell lines: problems and pitfalls. Br. J. Cancer 96, 928–936 (2007).17342089 10.1038/sj.bjc.6603663
18. Feng H T-lymphoblastic lymphoma cells express high levels of BCL2, S1P1 and ICAM1 Leading to a Blockade of Tumor Cell Intravasation Cancer Cell 2010 18 353 366 10.1016/j.ccr.2010.09.009 20951945
Feng, H. et al. T-lymphoblastic lymphoma cells express high levels of BCL2, S1P1 and ICAM1 Leading to a Blockade of Tumor Cell Intravasation. Cancer Cell 18, 353–366 (2010).20951945 10.1016/j.ccr.2010.09.009
19. Bose, B. & Bozdag, S. CTDPathSim: cell line-tumor deconvoluted pathway-based similarity in the context of precision medicine in cancer. In Proc. 11th ACM International Conference on Bioinformatics, Computational Biology and Health Informatics, in BCB ’20 1–10 (Association for Computing Machinery, 2020).
20. Cheng L Integration of genomic copy number variations and chemotherapy-response biomarkers in pediatric sarcoma BMC Med. Genom. 2019 12 23 10.1186/s12920-018-0456-5
Cheng, L. et al. Integration of genomic copy number variations and chemotherapy-response biomarkers in pediatric sarcoma. BMC Med. Genom. 12, 23 (2019).10.1186/s12920-018-0456-5
21. Sun Y Liu Q Deciphering the correlation between breast tumor samples and cell lines by integrating copy number changes and gene expression profiles BioMed. Res. Int. 2015 2015 901303 10.1155/2015/901303 26273658
Sun, Y. & Liu, Q. Deciphering the correlation between breast tumor samples and cell lines by integrating copy number changes and gene expression profiles. BioMed. Res. Int. 2015, 901303 (2015).26273658 10.1155/2015/901303
22. Wu Q Cancer-associated adipocytes: key players in breast cancer progression J. Hematol. Oncol. 2019 12 95 10.1186/s13045-019-0778-6 31500658
Wu, Q. et al. Cancer-associated adipocytes: key players in breast cancer progression. J. Hematol. Oncol. 12, 95 (2019).31500658 10.1186/s13045-019-0778-6
23. Mukherjee O Rakshit S Shanmugam G Sarkar K Role of chemotherapeutic drugs in immunomodulation of cancer Curr. Res Immunol. 2023 4 100068 10.1016/j.crimmu.2023.100068 37692091
Mukherjee, O., Rakshit, S., Shanmugam, G. & Sarkar, K. Role of chemotherapeutic drugs in immunomodulation of cancer. Curr. Res Immunol. 4, 100068 (2023).37692091 10.1016/j.crimmu.2023.100068
24. Hughes E T‐cell modulation by cyclophosphamide for tumour therapy Immunology 2018 154 62 68 10.1111/imm.12913 29460448
Hughes, E. et al. T‐cell modulation by cyclophosphamide for tumour therapy. Immunology 154, 62–68 (2018).29460448 10.1111/imm.12913
25. Verma R Lymphocyte depletion and repopulation after chemotherapy for primary breast cancer Breast Cancer Res. 2016 18 10 10.1186/s13058-015-0669-x 26810608
Verma, R. et al. Lymphocyte depletion and repopulation after chemotherapy for primary breast cancer. Breast Cancer Res. 18, 10 (2016).26810608 10.1186/s13058-015-0669-x
26. Kubota Y Ohji H Itoh K Sasagawa I Nakada T Changes in cellular immunity during chemotherapy for testicular cancer Int. J. Urol. 2001 8 604 608 10.1046/j.1442-2042.2001.00392.x 11903686
Kubota, Y., Ohji, H., Itoh, K., Sasagawa, I. & Nakada, T. Changes in cellular immunity during chemotherapy for testicular cancer. Int. J. Urol. 8, 604–608 (2001).11903686 10.1046/j.1442-2042.2001.00392.x
27. Godoy-Ortiz, A. et al. Deciphering HER2 breast cancer disease: biological and clinical implications. Front. Oncol.10.3389/fonc.2019.01124 (2019).
28. Tsoi J Multi-stage differentiation defines melanoma subtypes with differential vulnerability to drug-induced iron-dependent oxidative stress Cancer Cell 2018 33 890 904.e5 10.1016/j.ccell.2018.03.017 29657129
Tsoi, J. et al. Multi-stage differentiation defines melanoma subtypes with differential vulnerability to drug-induced iron-dependent oxidative stress. Cancer Cell 33, 890–904.e5 (2018).29657129 10.1016/j.ccell.2018.03.017
29. Hanahan D Weinberg RA The hallmarks of cancer Cell 2000 100 57 70 10.1016/S0092-8674(00)81683-9 10647931
Hanahan, D. & Weinberg, R. A. The hallmarks of cancer. Cell 100, 57–70 (2000).10647931 10.1016/S0092-8674(00)81683-9
30. Sever, R. & Brugge, J. S. Signal transduction in cancer. Cold Spring Harb. Perspect. Med. 10.1101/cshperspect.a006098 (2015).
31. Gonzalez H Hagerling C Werb Z Roles of the immune system in cancer: from tumor initiation to metastatic progression Genes Dev. 2018 32 1267 1284 10.1101/gad.314617.118 30275043
Gonzalez, H., Hagerling, C. & Werb, Z. Roles of the immune system in cancer: from tumor initiation to metastatic progression. Genes Dev. 32, 1267–1284 (2018).30275043 10.1101/gad.314617.118
32. Vishnoi K Viswakarma N Rana A Rana B Transcription factors in cancer development and therapy Cancers 2020 12 2296 10.3390/cancers12082296 32824207
Vishnoi, K., Viswakarma, N., Rana, A. & Rana, B. Transcription factors in cancer development and therapy. Cancers 12, 2296 (2020).32824207 10.3390/cancers12082296
33. Hernandez L Characterization of ovarian cancer cell lines as in vivo models for preclinical studies Gynecol. Oncol. 2016 142 332 340 10.1016/j.ygyno.2016.05.028 27235858
Hernandez, L. et al. Characterization of ovarian cancer cell lines as in vivo models for preclinical studies. Gynecol. Oncol. 142, 332–340 (2016).27235858 10.1016/j.ygyno.2016.05.028
34. Del Vecchio, C. A. et al. EGFRvIII gene rearrangement is an early event in glioblastoma tumorigenesis and expression defines a hierarchy modulated by epigenetic mechanisms. Oncogene10.1038/onc.2012.280 (2013).
35. Donnou S Murine models of B-cell lymphomas: promising tools for designing cancer therapies Adv. Hematol. 2012 2012 701704 10.1155/2012/701704 22400032
Donnou, S. et al. Murine models of B-cell lymphomas: promising tools for designing cancer therapies. Adv. Hematol. 2012, 701704 (2012).22400032 10.1155/2012/701704
36. Chen Z Dual effect of DLBCL-derived EXOs in lymphoma to improve DC vaccine efficacy in vitro while favor tumorgenesis in vivo J. Exp. Clin. Cancer Res. 2018 37 190 10.1186/s13046-018-0863-7 30103789
Chen, Z. et al. Dual effect of DLBCL-derived EXOs in lymphoma to improve DC vaccine efficacy in vitro while favor tumorgenesis in vivo. J. Exp. Clin. Cancer Res. 37, 190 (2018).30103789 10.1186/s13046-018-0863-7
37. Eberth S Epigenetic regulation of CD44 in Hodgkin and non-Hodgkin lymphoma BMC Cancer 2010 10 517 10.1186/1471-2407-10-517 20920234
Eberth, S. et al. Epigenetic regulation of CD44 in Hodgkin and non-Hodgkin lymphoma. BMC Cancer 10, 517 (2010).20920234 10.1186/1471-2407-10-517
38. Kuroiwa, Y. et al. Proliferative classification of intracranially injected HER2-positive breast cancer cell lines. Cancers10.3390/cancers12071811 (2020).
39. Carlisle SM Hein DW Retrospective analysis of estrogen receptor 1 and N‑acetyltransferase gene expression in normal breast tissue, primary breast tumors, and established breast cancer cell lines Int. J. Oncol. 2018 53 694 702 29901116
Carlisle, S. M. & Hein, D. W. Retrospective analysis of estrogen receptor 1 and N‑acetyltransferase gene expression in normal breast tissue, primary breast tumors, and established breast cancer cell lines. Int. J. Oncol. 53, 694–702 (2018).29901116
40. Kobayashi M Functional analysis of novel candidate regulators of insulin secretion in the MIN6 mouse pancreatic β cell line PLoS ONE 2016 11 e0151927 10.1371/journal.pone.0151927 26986842
Kobayashi, M. et al. Functional analysis of novel candidate regulators of insulin secretion in the MIN6 mouse pancreatic β cell line. PLoS ONE 11, e0151927 (2016).26986842 10.1371/journal.pone.0151927
41. Colaprico A TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data Nucleic Acids Res. 2016 44 e71 10.1093/nar/gkv1507 26704973
Colaprico, A. et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 44, e71 (2016).26704973 10.1093/nar/gkv1507
42. Zhang, J. CNTools: convert segment data into a region by sample matrix to allow for other high level computational analyses. In: R Package Version 1.40.0 https://bioconductor.org/packages/release/bioc/html/CNTools.html (2019).
43. Tate JG COSMIC: the catalogue of somatic mutations in cancer Nucleic Acids Res. 2019 47 D941 D947 10.1093/nar/gky1015 30371878
Tate, J. G. et al. COSMIC: the catalogue of somatic mutations in cancer. Nucleic Acids Res. 47, D941–D947 (2019).30371878 10.1093/nar/gky1015
44. Wang Y Franks JM Whitfield ML Cheng C BioMethyl: an R package for biological interpretation of DNA methylation data Bioinformatics 2019 35 3635 3641 10.1093/bioinformatics/btz137 30799505
Wang, Y., Franks, J. M., Whitfield, M. L. & Cheng, C. BioMethyl: an R package for biological interpretation of DNA methylation data. Bioinformatics 35, 3635–3641 (2019).30799505 10.1093/bioinformatics/btz137
45. Xu Z Niu L Li L Taylor JA ENmix: a novel background correction method for Illumina HumanMethylation450 BeadChip Nucleic Acids Res. 2016 44 e20 10.1093/nar/gkv907 26384415
Xu, Z., Niu, L., Li, L. & Taylor, J. A. ENmix: a novel background correction method for Illumina HumanMethylation450 BeadChip. Nucleic Acids Res. 44, e20 (2016).26384415 10.1093/nar/gkv907
46. Jaffe, A. E. FlowSorted.Blood.450k: Illumina HumanMethylation data on sorted blood cell populations. In: Bioconductor R Package Version https://bioconductor.org/packages/release/data/experiment/html/FlowSorted.Blood.450k.html (2020).
47. Hansen, K. D. et al. minfi: Analyze Illumina Infinium DNA methylation arrays. In: Bioconductor Version: Release (3.11)10.18129/B9.bioc.minfi (2020).
48. Garnett MJ Systematic identification of genomic markers of drug sensitivity in cancer cells Nature 2012 483 570 575 10.1038/nature11005 22460902
Garnett, M. J. et al. Systematic identification of genomic markers of drug sensitivity in cancer cells. Nature 483, 570–575 (2012).22460902 10.1038/nature11005
49. Onuchic V Epigenomic deconvolution of breast tumors reveals metabolic coupling between constituent cell types Cell Rep. 2016 17 2075 2086 10.1016/j.celrep.2016.10.057 27851969
Onuchic, V. et al. Epigenomic deconvolution of breast tumors reveals metabolic coupling between constituent cell types. Cell Rep. 17, 2075–2086 (2016).27851969 10.1016/j.celrep.2016.10.057
50. Koestler DC Blood-based profiles of DNA methylation predict the underlying distribution of cell types: a validation analysis Epigenetics 2013 8 816 826 10.4161/epi.25430 23903776
Koestler, D. C. et al. Blood-based profiles of DNA methylation predict the underlying distribution of cell types: a validation analysis. Epigenetics 8, 816–826 (2013).23903776 10.4161/epi.25430
51. Titus AJ Gallimore RM Salas LA Christensen BC Cell-type deconvolution from DNA methylation: a review of recent applications Hum. Mol. Genet. 2017 26 R216 R224 10.1093/hmg/ddx275 28977446
Titus, A. J., Gallimore, R. M., Salas, L. A. & Christensen, B. C. Cell-type deconvolution from DNA methylation: a review of recent applications. Hum. Mol. Genet. 26, R216–R224 (2017).28977446 10.1093/hmg/ddx275
52. Maunakea AK Conserved role of intragenic DNA methylation in regulating alternative promoters Nature 2010 466 253 257 10.1038/nature09165 20613842
Maunakea, A. K. et al. Conserved role of intragenic DNA methylation in regulating alternative promoters. Nature 466, 253–257 (2010).20613842 10.1038/nature09165
53. Jassal B The reactome pathway knowledgebase Nucleic Acids Res. 2020 48 D498 D503 31691815
Jassal, B. et al. The reactome pathway knowledgebase. Nucleic Acids Res. 48, D498–D503 (2020).31691815
54. Ulgen, E. egeulgen/pathfindR). R https://github.com/egeulgen/pathfindR (2020).
55. Du P Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis BMC Bioinform. 2010 11 587 10.1186/1471-2105-11-587
Du, P. et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinform. 11, 587 (2010).10.1186/1471-2105-11-587
56. Mermel CH GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers Genome Biol. 2011 12 R41 10.1186/gb-2011-12-4-r41 21527027
Mermel, C. H. et al. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 12, R41 (2011).21527027 10.1186/gb-2011-12-4-r41
57. Reich M Liefeld T Tamayo P Mesirov J GenePattern 2.0 Nat. Genet. 2006 38 500 501 10.1038/ng0506-500 16642009
Reich, M., Liefeld, T., Tamayo, P. & Mesirov, J. GenePattern 2.0. Nat. Genet. 38, 500–501 (2006).16642009 10.1038/ng0506-500
58. Tukey’s range test, Wikipedia (2020). Accessed: Jan. 09, 2021. [Online]. Available: https://en.wikipedia.org/w/index.php?title=Tukey%27s_range_test&oldid=979274615
59. Warren, A. et al. Global computational alignment of tumor and cell line transcriptional profiles. Nat Commun 12, 22 (2021).
60. Serdar, B. & Bose, B. CTDPathSim2.0 Dataset [DS/OL]. V2. Science Data Bank 10.57760/sciencedb.01713 (2022).
