The cell lines, IC50 of which are smaller than maximal tested concentrations, are defined as sensitive, otherwise resistant. Our method iteratively splits the cancer cells into two groups in a way that gains the best separation of drug sensitivity between two groups until reaching two terminal conditions: The p value of the difference of two groups’ drug sensitivity values is smaller than 0.05 or all cancer cells in a node are sensitive or resistant.

Suppose the set of all cells is S total , and each cell has its drug sensitivity value. The procedure to grow a pharmacological subtype tree is as follows:

1) Define the root node. Set S as the set of cancer cells in this node, and assign S total to S .

We used a heuristic algorithm to search for a reasonable number of divisions to split S into S l and S r . Suppose the number of the cells set is n s . We first sorted the cells by IC50 ascending and divided the cells with higher drug sensitivity into the S l group and the cells with lower drug sensitivity into the S r group. Since each node at least 6 cells, a total of n s − 5 divisions were considered.

For each division, the Mann-Whitney U test was applied to calculate the p value of drug sensitivity between S l and S r . The division with the smallest p value was selected as the optimized one to split the cancer cells in the current node.

The similarity of cells is calculated based on their response to the different drugs. For each drug k , a pharmacological subtypes tree S k is generated, consisting of n subtypes: S k = ( S 1 k , S 2 k , … ,   S n k ) = ( S t k ) t = 1 → n

Based on this tree, a matrix that defines a similarity of cancer cells, regarding whether they are in the same subtype and their drug response to drug k , is calculated. We assume there are m cell lines tested for drug k . The cell lines are labelled as sensitive, resistant, or NA (not available). A sensitive cell line has its IC50 value smaller than maximum tested concentration; a resistant cell line has its IC50 value greater than maximal tested concentration. For two given cell lines i and j , their similarity matrix ( D i j k ) is defined as: D i j k = { 1 ,   if ( i , j ) ∈ S t k   and   both   are   labelled   as   sensitive   ones ; − 1 ,   if ( i , j ) ∈ S t k   and   both   are   labelled   as   resistant   ones ; 0 ,   otherwise .

If two cell lines i and j are in the same subtype and both are sensitive to drug k , their score is 1, and −1 for both cell lines with a resistant response, while if they are not in the same subtype or show opposite response to the drug, their score is 0. A m × m similarity matrix of cell lines is then generated.

Then, the similarity of cell lines i and j can be further calculated as: Cell sim = ∑ k = 1 q | D i j k | q ,   i ≠ j Wherein i and j denote the cell lines tested for q drugs.

Furthermore, for two given drugs a and b , their similarity is calculated as the summarized similarity across the same cell lines which are tested for both drugs. Dru g sim = ∑ i = 1 m ∑ j = 1 m D i j a · D i j b ∑ i = 1 m ∑ j = 1 m D i j a · D i j a + ∑ i = 1 m ∑ j = 1 m D i j b · D i j b − ∑ i = 1 m ∑ j = 1 m D i j a · D i j b ,   i ≠ j Wherein i and j denote the cell lines tested for both drugs, and m is the total number of the same cell lines.

The maximum concentration of the drug was used to distinguish between the resistant pharmacological subtypes and sensitive pharmacological subtypes. The chi-squared test was used to calculate the connection between genomic alterations, including mutation and translocation, and pharmacological subtypes. Genomic alterations were determined with p values corrected with the Benjamini-Hochberg method for controlling the false discovery rate (Benjamini and Hochberg, 1995). Here we used the corrected p value <0.1 as the threshold. We defined two types of connections: sensitive and resistant. If pharmacological subtypes with higher drug sensitivity had more frequent alteration occurrences, we defined this connection as sensitive. If the opposite was true, we defined the connection as resistant.

The connection between mRNA expression and pharmacological subtypes was determined with p values from the student t test corrected with the Benjamini-Hochberg method. If pharmacological subtypes with higher drug sensitivity had higher expression, we defined this connection as sensitive. If the opposite was true, we defined the connection as resistant. Here we used the corrected p value <0.1 as the threshold. To obtain more confident connections, the Pearson correlation coefficient (PCC) was used to measure the correlation between drug sensitivity (IC50) and gene expression in the sensitive subtypes, and only genes with PCC greater than 0.3 were considered. The individual subtypes with sensitive and resistant cells mixed were excluded in the above analyses.

Gene set enrichment analysis (GSEA) (Subramanian et al., 2005) was employed to determine the 333 gene sets from KEGG, enriched by a pre-ranked list of all genes, which were sorted by the statistical significance of differential expression defined by DESeq2 analysis (Love et al., 2014). Gene sets with FDR <0.05 were statistically significant.

The Fisher’s exact test was used to respectively determine whether there is a significant association between pathways and clusters of drugs determined by the number of drug pharmacological subtypes across 6 cancer types or not, whether any difference of phenotypes (cancers/pathways) in clusters of cell lines/drugs derived from the similarity matrix was significant, and whether cancer censor genes are enriched in a list of genes with the most connections to pharmacological subtypes. A hypothesis test based on hypergeometric distribution was used to determine whether a histological cancer is enriched in the most sensitive pharmacological subtype of a drug. The p value <0.05 was regarded as statistically significant. Kolmogorov-Smirnov test (K-S test) was employed to compare cumulative distribution function (CDF) of the number of pharmacological subtypes across the drugs between two histological cancers.

We built a pharmacological tree for each drug based on the divisibility of the drug sensitivity of cancer cells. Taking the tree of the MEK1/2 inhibitor Refametinib as an example (Figure 1A), all cells in the root node were first divided into left and right child nodes with relatively high and low sensitivity, respectively, and further cells in these two nodes were capable of being divided into five final subgroups, C1, C2, C3, C4, and C5. Each subgroup had varied drug sensitivity and reached maximum divisibility, thus representing distinct pharmacological subtypes. C1 denotes the most sensitive and C5 the most resistant subtype to the drug Refametinib. Except for C3 where both sensitive and resistant cells were mixed, the other subgroups contained the homologous sensitive or resistant cells (Figure 1B). There were much more sensitive cells in C1 and C2 than resistant cells in C4 and C5 (Figure 1C). We then investigated how histological subtypes were distributed in the pharmacological subtypes (Figure 1D). C1 was composed of 25 histological subtypes, with top 2 cancers, SKCM (19.7%) and COREAD (12.4%). C5 was composed of 22 histological subtypes, with leading cancers including SCLC (31.8%) and BRCA (17.0%).

Our analysis established 367 pharmacological trees. The number of pharmacological subtypes with distinct drug sensitivities varied from 2 to 9 across 367 trees (Figure 1E). This suggested that heterogeneity of treatment effects across histological subtypes existed widely, and one pharmacological subtype presenting concordant drug response indeed comprised histological disparate cells to some extent. A high degree of treatment heterogeneity, meaning more than seven pharmacological subtypes, was observed for 54.0% (16.9% of the drugs had 7; 15.8% of the drugs had 8; and 21.3% of the drugs had 9 pharmacological subtypes) of the drugs. These drugs included Dabrafenib (BRAF inhibitor), Selumetinib (MEK1/2 inhibitor), Erlotinib (EGFR inhibitor), and Alectinib (ALK inhibitor). 44.1% of the drugs had 5 or 6 subtypes, showing a moderate degree of treatment heterogeneity, while 1.9% of the drugs had subtypes lower than 4, showing a low degree of heterogeneity related to drug therapy (Figure 1E).

We also particularly examined how pharmacological subtypes were constituted within the same histological tumor subtypes, including LUAD, BRCA, SKCM, SCLC, COREAD, and HNSC, in which the average number of cell lines per drug was more than 35 (Figures 1F–H). We built approximately 310 pharmacological trees for each type of the above cancers. The density distribution of the number of pharmacological subtypes indicated the varied treatment heterogeneity across 6 histological subtypes (Figure 1F). The shape of the distribution curve in LUAD was characterized by double kurtosis and significantly positively biased, and oppositely the curve in HNSC was negatively biased, suggesting the highest treatment heterogeneity in LUAD and lowest treatment heterogeneity in HNSC. Furthermore, K-S test was applied to evaluate if two different histological subtypes had the same level of treatment heterogeneity by comparing their cumulative distributions. The results showed the cumulative distribution curve of LUAD is significantly different from that of HNSC (p = 1.11e-16, K-S test) (Figure 1G). The pairwise comparisons between the curves of any two histological subtypes showed BRCA, SCLC, and SKCM had the statistically similar distributions (p > 0.05, K-S test), and other cancers had not (p < 0.05, K-S test) (Supplementary Figure S1).

To observe the varied treatment heterogeneity of the different histological subtypes in detail, the number of pharmacological subtypes of drugs across 6 histological cancers were shown in the heatmap (Figure 1H). Six histological tumors were arranged in columns according to their treatment heterogeneity, from high to low. Interestingly we observed that EGFR inhibitors (Afatinib, Cetuximab, and Gefitinib) targeting in EGFR signalling pathway (green drugs in Figure 1H) were clustered together (p < 0.001, Fisher’s exact test). And in LUAD and HNSC there were more pharmacological subtypes of EGFR inhibitors than other drugs, suggesting high treatment heterogeneity for EGFR inhibitors in LUAD and HNSC. In addition, drugs targeting ERK/MAPK signalling (red drugs in Figure 1H) were clustered together (p < 0.001, Fisher’s exact test). There was high treatment heterogeneity for these drugs in SKCM.

To identify the molecular alterations associated with varied drug sensitivity, we examined whether the genomic alterations, including mutations and translocations, significantly changed across the pharmacological subtypes. The statistical significance was determined by p values calculated with the Chi-squared test and corrected with the Benjamini-Hochberg method for controlling the false discovery rate. If the pharmacological subtypes with higher drug sensitivity had more frequent mutation or translocation events, we defined their association as a sensitive connection. If the opposite was true, we defined their association as resistant. Consequently, genes with specific genomic alterations were associated with the respective drugs. Since one gene could have multiple connected drugs and one drug could have multiple connected genes, these gene-drug connections ultimately composed a network that allowed us to investigate the contribution of either a drug or gene to the holistic understanding of how altered genes connect to drug responses (Figure 2A). The hub genes with the most connections to the drugs, including known cancer genes KRAS, RB1, BRAF, TP53, NOTCH2, and other genes LTB, TNFRSF9, DERA et al. were highlighted in the network.

We first investigated the distribution of genes on their connected drugs (Figure 2B). Our analysis identified 624 genes whose genomic alterations were connected with at least two drugs’ pharmacological subtypes. Most of the genes (70.67%) had fewer than 3 connections with drugs, but 40 (6.41%) genes whose genomic alterations were associated with more than 5 drugs. By ranking genes by their number of connections to drugs, we found that the top ranked genes were enriched in cancer genes (Figure 2C; Supplementary Figure S2), a catalogue of genes with mutations or fusions that are causally implicated in cancer provided by the COSMIC database (Tate et al., 2019). The top 10 genes connected to more than 12 drugs, 30% of which were remarkably annotated as cancer genes (Figure 2C). These cancerous percentages for the top 10 (p = 3.62e-03, Fisher’s exact test), top 20 (p = 1.27e-07, Fisher’s exact test), and top 30 genes (p = 1.62e-09, Fisher’s exact test) were significantly higher than the background when using genome-wide 21,972 genes as a reference (Figure 2C). Among the top ranked oncogenes, ABL1_BCR ranked first with connections to 28 drugs, followed by BCR_ABL1 and KRAS with 20 and 12 drugs (Figure 2D), implying that the highly drug-connected genes that are not currently identified as cancer genes could actually be causally implicated in cancer therapy, such as KCTD1_PCAT18 fusion, TNFRSF9, LTB, and NFATC4, which showed connections to over 18 drugs and ranked as the top five genes (Supplementary Figure S2).

We then investigated the distribution of drugs on their connected genes. In this analysis, we focused only on cancer censor genes. For the majority of drugs (86.36%), the number of connected genes ranged from 2 to 5 (Figure 2E). Among the pharmacological subtypes, genes whose alteration was sensitively connected with drug resistance were significantly more abundant than those that were connected in a resistant manner (p = 0.0102, student t test). (Figure 2F). We ranked the drugs with more than 5 genes’ connection with their pharmacological subtypes (Figure 2G). Drugs with gene connections greater than 5 included GSK690693, affecting the PI3K/MTOR signaling pathway, TAK−715, affecting the JNK and p38 signaling pathway, VX−11e, affecting the ERK MAPK signaling pathway, and IGFR_3801, affecting the IGF1R signaling pathway.

Next, we observed how the mutation rate of genes associated with the variability of drug sensitivity changed among drug-sensitive and drug-resistant subtypes. Taking the MEK1/2 inhibitor Refametinib as an example, three genomic alterations were found to connect to the drug and be associated with pharmacological subtypes. As shown in Figure 2H, BRAF mutations and KRAS mutations occurred significantly more frequently (over 25%) in the sensitive subtype C1 and then decreased gradually in the subtypes C2, C3, C4, and C5. Here one subtype C3 with sensitive and resistant cells mixed was excluded in the analysis. Conversely, RB1 mutation, occurred more frequently in the resistant subtypes. Therefore, C1 group, consisting mainly of SKCM (19.7%) and COREAD (12.4%) (Figure 1D), was significantly sensitive to Refametinib due to its high occurrence of BRAF and KRAS mutations and low occurrence of RB1 mutations (Figure 2H). Other types of MEK inhibitors, including PD0325901, Trametinib, Selumetinib, and CI-1040, further confirmed that BRAF, KRAS, and RB1 mutations were robustly connected to pharmacological subtypes, contributing to the drug sensitivity of MEK inhibitors (Supplementary Figure S3). Our analysis revealed that the varied distribution of genomic alterations across pharmacological subtypes could lead to their differences in response to anticancer therapies.

We also associated mRNA expression with the pharmacological subtypes. The connection between gene expression and pharmacological subtypes was determined by p values calculated with the student t test and corrected with the Benjamini-Hochberg method for controlling the false discovery rate. If the subtypes with higher drug sensitivity had higher gene expression, we defined this connection as sensitive. If the opposite was true, we defined the connection as resistant. The connections between genes and drugs constituted a network (Figure 3A). To narrow down the genes whose expression levels were directly associated with drug sensitivity, we only kept genes which act in core cancer pathways (Supplementary Table S1), and whose Pearson correlation between drug sensitivity and gene expression were greater than 0.3 in the sensitive subtypes. The hub genes with the most connections to the drugs, including known cancer genes WWTR1, JUN, BCL9L and other genes NQO1, CAPN2, and PERP et al., were highlighted in the network.

We first investigated the distribution of genes on their connected drugs (Figure 3B). Our analysis identified 684 genes whose expression was connected to at least two drugs’ pharmacological subtype. Among them, 4.97%, 14.04%, and 24.85% of genes respectively had connections with greater than 20, 10, and 5 drugs. The top 30 genes connected to more than 20 drugs, 25% of which were annotated as cancer genes (Figure 3C; Supplementary Figure S4). These cancerous percentages for the top 20 (p = 4.37e-06, Fisher’s exact test), and top 30 genes (p = 5.84e-52, Fisher’s exact test) were significantly higher than the background. The top genes with their expression alterations associated with subtypes were different from those with genomic alterations associated with subtypes (Figures 2D, 3D). Among them, WWTR1, BCL9L, and JUN ranked the top three with connections to 26 drugs and 25 drugs (Figure 3D).

We then investigated the distribution of drugs on their connected genes. The number of genes whose expression was associated with the subtypes varied from 0 to 223 across the drugs, with 4.65% drugs had more than 100 gene connections (Figure 3E). Comparing with genomic alterations, there were more genes whose expression was connected to the drugs. Moreover, genes whose expression was sensitively connected with drug resistance across pharmacological subtypes were significantly more abundant than genes that were resistant connected (p = 0.0462, student t test) (Figure 3F). The top 20 drugs with the most connections to genes included Nilotinib, affecting the ABL signaling, PHA−793887 and NPK76−II−72−1, affecting cell cycle, Lapatinib and KIN001−236, affecting RTK signaling pathway, Navitoclax, affecting Apoptosis regulation, and some compounds such as Ara−G, XMD14−99 et al., affecting multiple or unknown pathways (Figure 3G).

Methotrexate, as an example, is a chemotherapy that specifically acts during DNA and RNA synthesis, and cancer cells were classified into seven subtypes (Figure 3H) in terms of drug sensitivity to methotrexate. The first subtype, C1, was most sensitive to the drug, and the seventh subtype, C7, was the most resistant to the drug (Figure 3H). C1 included 12 histological subtypes, with the top 2 shown in the bottom of Figure 3H, ALL and DLBC. C2 and C3 were composed of 18 and 19 histological subtypes, respectively, with SCLC both included in the top 2 (Figure 3H). A total of 38 genes with high expression in sensitive group (C1, C2, and C3) and low expression in resistant group (C5, C6, and C7), which had sensitive connections to pharmacological subtypes of methotrexate, were selected (Figure 3I). Pathway enrichment analysis showed that these genes were enriched in pathways including cell cycle and DNA replication et al. (Figure 3J), suggesting that the ectopic activation of cell cycle led to the sensitive response of ALL, DLBC, LAML, and SLCL to Methotrexate.

Pharmacological subtypes provide us a new point of view to redefine the similarity of cells purely from therapeutic concordance. Thus, we calculated the similarity of pairwise cells based on pharmacological subtypes of drugs. If there are quite a few cases that two cells are in the same sensitive or resistant pharmacological subtypes of drugs, these cells gain a high similarity. Otherwise, they gain a low similarity. The similarity of cells defined a hierarchical cluster as shown in Figure 4A. We further applied the Fisher’s exact test to observe if the cells in the same clusters were from the same histological subtypes (Figures 4B–D). The points with black border in the figure represent statistical significance (p < 0.05), with the smaller p values shown in red.

The results showed that the cancer cell lines were respectively divided into 2, 4, and 15 clusters at the different hierarchical cluster levels (Figure 4A). Two clusters at the high level were reflected by two distinct patterns of similarity (Figure 4A). One cluster (Cluster H1) was significantly correlated with solid tumors including BRCA, COREAD, LUAD, MESO, OV, and PAAD, and the other (Cluster H2) hematological tumors including ALL, DLBC, LAML, LCML, and MM (Figure 4B). Here, the point with black border represents statistical significance (p < 0.05, Fisher’s exact test) of the association between the cluster and cancer type. So in a holistic view, the different types of hematological tumors shared a similar drug sensitivity profile, which was quite different from the profile of solid tumors. However, for a given cancer type (the right panel of Figure 4B), such as ACC, BLCA, CESC, CLL, or ESCA et al., its cells were unbiasedly distributed in two clusters (p > 0.05, Fisher’s exact test), showing the treatment heterogeneity for the cells from the same tissue origin. At the median hierarchical levels, 4 clusters were identified (Figure 4C). Among them, Cluster M3 was specifically associated with haematological tumors (p < 0.001, Fisher’s exact test) including ALL, DLBC, LAML, LCML, and MM. The treatment heterogeneity was still observed in cancers.

When cell lines were divided into 15 clusters at the lower hierarchical level (Figure 4D), we observed that the cells from ACC, CLL, STAD, or UCEC were unbiasedly distributed in the different clusters, still showing their treatment heterogeneity. Interestingly, the cells from some cancers such as HNSC, KIRC, and SCLC were enriched in more than 2 clusters with distinctly different drug response. For instance, HNSC cells were enriched in Clusters L12, L4 and L5; KIRC cells were enriched in Clusters L15 and L4; SCLC cells were enriched in Clusters L13, L3, and L7. Interestingly, the SCLC cells in Cluster L13 shared a more similar profile with hematological tumors including DLBC, LAML, and LCML. To investigate the underlying mechanism of such treatment heterogeneity in SCLC, we applied DESeq2 analysis to call the differentially expressed genes between Cluster L13 and Cluster L3, and further identified the differential functions via gene set enrichment analysis (GSEA) (Subramanian et al., 2005). The results showed that the functions like DNA replication, RNA transport et al. were significantly upregulated in Cluster L13 cells, while in Cluster L3 cells the functions like immune response, antigen processing and presentation et al. were mainly upregulated (Supplementary Figure S5A). Moreover, MYC, as an oncogenic transcription factor of cell growth and proliferation via enhancing the cell cycle regulated genes (Oster et al., 2002), was significantly upregulated in Cluster L13 than in L3, with a log2-fold change of 4.28 (Supplementary Figure S5B). These suggested that SCLC cell lines in Cluster L13 were higher Myc-dependent than those in Cluster L3, characterized with increased cell proliferation and evasion of immune response.

The similarity of pairwise drugs was also defined as the summarized similarity of drug response across the same cells which were tested for both drugs. For two drugs, if there are quite a few of cell pairs which are in the sensitive (or resistant) pharmacological subtypes of one drug are also in the sensitive (or resistant) subtypes of the other drug, they gain a high similarity. Otherwise, they gain a low similarity. We then generated a drug similarity profile across 367 drugs, and similarly, a hierarchical clustering algorithm was used to group drugs based on drug similarity (Figure 4E). The results showed that at the high level, drugs were divided into two clusters, demonstrating two different similarity patterns. The first cluster (Cluster H1) mostly included chemotherapy drugs such as those effecting in DNA replication and mitosis, and the second cluster was mixed with target and chemotherapy drugs, specifically enriched by the drugs on RTK signaling, other kinases, and cytoskeleton (Figure 4F). When the drugs were divided into 4 clusters (Figures 4E,G), we observed that the drugs enriched in L2 included chemotherapy drugs effecting in cell cycle, chromatin histone accelylation, and DNA replication, the drugs enriched in L3 were functioning in PI3K/mTOR signaling pathway, and the drugs enriched in L4 were hormone related ones. Moreover, L4, accounting for a massive part of drugs, also included the drugs targeting the different types of signaling pathways such as ABL, EGFR, Jak, and P38, RTK, and WNT signalling. These drugs did not reach the significant level but provided the obvious association.

To enable knowledge sharing and reuse, we developed a database (http://www.hywanglab.cn/dtdb/, Supplementary Figure S6), which provided the comprehensive knowledge of pharmacological subtypes of drugs in this study, for supporting the further exploration of therapeutic response-based classification of tumors. The web interface of the database is user-friendly and allows users to search, browse and download data. A record of the database includes the following information: the pharmacological subtypes of a drug, the drug sensitivity values of each subtype, the distribution of histological subtypes in each pharmacological subtype, genomic and expression alterations connected to subtypes. In addition, the users can acquire the other information of a drug, a gene, or a specific cancer type in terms of our analysis. For example, given a drug, the other drugs with the most similar drug response across cancer cells are provided. Similarly, given a cancer cell line, the other cell lines with the most similar drug response across drugs are provided. Given a histological cancer type, such as SCLC, the pharmacological subtypes with the distinct treatment response were depicted, as well as the underlying mechanism of such treatment heterogeneity. The resources are expected to support the researches of precision medicine from a purely therapy-oriented perspective.