Integrated single-dose kinome profiling data is predictive of cancer cell line sensitivity to kinase inhibitors

Data sources

The primary data sources we used can be split into two categories: the integrated kinome profiling data set and the cancer cell line set:

The following were downloaded from the respective Supplementary Materials to create the integrated set of kinome profiling data:

1. Kinome profiling data from the kinobeads assay

   1. Klaeger et al. (2017)

2. Kinome profiling data from KINOMEscan assay

   1. LINCS: kinome profiling datasets for individual compounds downloaded programmatically from http://lincs.hms.harvard.edu/db/datasets/ (Koleti et al., 2018)

   2. Kinome profiling data for the PKIS drug set was downloaded from the Supplementary Data of https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0181585 (Drewry et al., 2017)

   3. Kinome profiling data for the KCGS drug set was downloaded from the Supplementary Data of https://www.mdpi.com/1422-0067/22/2/566 (Wells et al., 2021) and from internal data provided by SGC-UNC.

The following were downloaded from the DepMap portal (https://depmap.org/portal/download/all/) to create the set of cancer cell line sensitivities and their gene expression characteristics:

1. DepMap secondary repurposing screen (“secondary-screen-dose–response-curve-parameters.csv”)

2. CCLE gene expression set (“CCLE_expression.csv”)

The high-throughput screening data gathered in PDAC patient-derived cell lines was gathered from Lipner et al. (2020) with methods as described in Berginski et al. (2021a).

Data preprocessing

The scripts implementing these descriptions are all available through github.

Klaeger et al. Kinobead Kinase Inhibition Profiles: As previously described (Berginski et al., 2022), we read the values from the Supplemental Data table into R and produced a filtered list of kinase and kinase interactor relative intensity values. We imputed missing values with the default “no interaction” value of 1, and truncated likely outlier values to the 99.99 percentile (3.43).

KINOMEscan Inhibition Profiles: We read in the three datasets mentioned above into R and concatenated them into a single combined data set. All the individual data sets contain identical protein lists because of the same assay type. Values are reported as “Percent Control”, a ratio of protein pulled down in experimental condition (with inhibitor) vs control condition (without inhibitor). These were divided by 100 to convert the scale to 0–1 to match the Kinobeads relative intensity data.

Creating the Combined Kinome Inhibition Profiling Set: We took the kinobeads dataset and the KINOMEscan dataset and concatenated them into one large set containing inhibitor-kinase interaction states for ∼800 total kinases and kinase interactors. We left out assays that included recombinantly mutated kinases but left those with naturally occurring post-translational modifications. The vast majority (99.95%) of the inhibitor-kinase pairs represented was unique for either assay type, but for the 0.05% inhibitor-kinase pairs, we took the mean value of the measurements across the two assay types. Additionally, any missing values were imputed with the default “no interaction” value of 1. In the end we were left with kinome inhibition states for ∼1,000 kinase inhibitors.

Dataset of Cancer Cell Line Sensitivity to Drugs from DepMap: The DepMap repurposing dataset contains cell viability measurements across multiple doses, but since our dataset of kinome states is restricted to single-dose measurements, we extracted single summary statistics of cell line sensitivity to kinase inhibitors: Dose–response area under the curve (AUC) and half-maximal inhibitory concentration (IC₅₀). We extracted these by reading in the “secondary-screen-dose–response-curve-parameters” dataset into R, which contains curve parameters for a log–logisitic curve fit to the cell viability dose response curve and filtered it down to cell line name, IC₅₀, AUC and other associated metadata.

Matching of Kinase Inhibitors between Profiling Dataset and Cell-Line Sensitivity Dataset: The compound names from each dataset were read into R, and the package Webchem (Szöcs et al., 2020) was used to retrieve PubChem compound IDs. The two sets of compound names were then matched based on these reference IDs. There were 252 matches between the two sets, forming a final set of ∼70,000 inhibitor-cell line combinations.

Baseline Gene Expression: As described before (Berginski et al., 2022) the RNAseq data provided in the “CCLE_expression.csv” file needed no modifications while preprocessing. Our only modification was to add identifiers to each gene label (“exp_”), to ensure that kinome inhibition data and expression data related to the same gene weren’t accidentally combined.

String: The STRING database (Szklarczyk et al., 2021) was processed as described previously (Berginski et al., 2022) to annotate kinases and kinase interacting genes.

Modeling techniques

To assess our models we used a random 10-fold cross validation strategy. The number of features was varied as specified by the feature selection scheme described in the results section. We compared the performance of three model types using this strategy: LASSO (Least Absolute Shrinkage and Selection Operator) regression using the glmnet engine (Friedman, Hastie & Tibshirani, 2010), random forest using the ranger engine (Wright & Ziegler, 2017) and gradient boosting using the XGBoost (eXtreme Gradient Boosting) engine (Chen & Guestrin, 2016). Model performance was assessed by the R² value between predicted and actual outcome within the cross-validation scheme. For each model type, we tuned sets of 30 hyperparameters to find the best possible performer as follows:

1. LASSO

   1. Penalty (1E−10–0.9)

2. Random Forest

   1. Trees (100–2000)

3. XGBoost

   1. Trees (100–1000)

   2. Tree Depth (4–30)

After final model selection, we fit the model on the entire dataset and then made predictions on inhibitor-cell line pairs not found in the original DepMap screening data.

Compound testing

BT-474, HCC1806, SUM-159 and SKBR-3 cells were assayed as described previously (Berginski et al., 2022). Briefly, cells were grown in ATCC recommended media and seeded in triplicate at 4000, 2000, 4000 and 500 cells per well respectively. 24 h after seeding, cells were treated with inhibitors at 30 µM, 3 µM, 1 µM, 300 nM, 100 nM, 30 nM, 10 nM, and 3 nM, along with the appropriate DMSO controls. seeded at, in white flat-bottom 96-well plates (Corning). Seventy-two hours post-treatment, cells were lysed with CellTiter-Glo (Promega) and luminescence was read using the PHERAstar FS microplate reader (BMG Labtech) and gain adjustments were conducted for each cell line. Notably, we used a 3-day treatment period instead of the 5-day treatment period used in the highly multiplexed PRISM assay, since we saw over-confluence in our 96-well plate assay in times longer than 3 days. However, it is likely that summary measures of trends in response like IC50 and AUC are minimally affected by the small difference in treatment times. Data were averaged over replicates normalized row-wise to the DMSO-only (0.1% on cells) control samples on each plate to calculate relative viability. Quality checks were performed to look at the data distribution and the presence of spatial bias on a plate. A quality control metric of <120% of DMSO was applied to all rows analyzed. The full file of cell viability curves is available in Table S1.

The functions “ComputeAUC” and “ComputeIC50” from the R package dr4pl (Gadagkar & Call, 2015) was used to fit a four-parameter log–logistic curve to the cell viability data, and extract AUC and IC₅₀ values from the 9-point cell viability curves.

Creating an integrated set of kinome profiling data across a wide chemical space

Kinase inhibitors have been profiled using a number of assays, but for this study we have focused on a specific subset of kinase inhibitors that have been assayed using the kinobead/MS-based method (Klaeger et al., 2017) or the KINOMEscan® (DiscoverX) method. These methods assess their specific kinase targets as well as the magnitude of inhibition of each kinase in response to different inhibitor concentrations (Klaeger et al., 2017). We combined kinome profiling datasets from Klaeger et al. (2017) (Kinobeads), LINCS (Koleti et al., 2018) (KINOMEscan), and UNC (Drewry et al., 2017) (KINOMEscan), filtering down to profiles measured only at 1uM. For the small amount of overlap between datasets, the mean inhibition value was taken across drug-kinase combinations. Given that both assays measure the engagement of inhibitors to kinases, most of the proteins that appear in assay results are either known kinases or closely associated proteins. Specifically, kinase inhibitor profiles include measurements on all wild-type and phosphorylated kinases (∼500), along with a set of associated proteins (∼300). As such, we will refer to this data as “kinase inhibition states”, and the profile of each individual drug as its “kinome inhibition state” (Fig. 1A).

After integration, we were left with a final set of ∼1,000 compounds with corresponding information on their inhibitor-induced kinome states, describing changes in ∼800 kinases and kinase interactors. To summarize the relationship between inhibition states induced by all inhibitors, we performed a UMAP dimensionality reduction (McInnes et al., 2018) on the dataset, only using kinases profiled in both assays for visualization (Fig. 1B). The UMAP coordinates represent the aggregate effect of each inhibitor on the kinome, i.e., it is a representation of the uniqueness of its kinome inhibition state. Inhibitors that have similar effects on the kinome will have similar coordinates, while disparate inhibitors will have coordinates that are far apart. Using this, we can examine the diversity of kinome space targeting in our dataset, based on the origin of kinome profiling data. Our analysis shows that integrating KINOMEscan data for 800 inhibitors vastly increases the kinome space targeted (Fig. 1B), compared to kinobeads alone.

However, it is important to note that despite the similarity in the assay readouts, the KINOMEscan and kinobeads assays are very different molecularly, and care needs to be taken while integrating these assay types. To assess the degree of difference between the two assay types, we retained all kinases captured by both assays (∼250) and measured Pearson’s correlation of inhibition states for the same compound between both assays. We found a reasonable degree of correlation (Pearson’s R ∼0.5) across all data. In addition, we additionally found a much higher degree of agreement between data types (90%) when both assays were cut off at 80% inhibition (a “strong” hit in either assay) (Fig. S1). As a result, keeping in mind the advantages of integrating both assay types, we moved forward with connecting these integrated kinome inhibition states to cell line responses.

Connecting inhibited kinome states to cancer cell line sensitivities from the DepMap repurposing screen

To connect these kinase inhibitors and their induced inhibition states with their corresponding phenotypes in cancer cell lines, we make use of the DepMap repurposing screen, which uses the PRISM assay (Yu et al., 2016) to run highly multiplexed cell viability assays. This dataset contains cell viability measurements for over 1,500 drugs profiled in 450 cell lines. From within this data, we found ∼200 drugs for which we also have corresponding profiling data as described above.

The DepMap repurposing dataset provides cell viability measurements across multiple drug doses. However, since our dataset of kinome states is restricted to single-dose measurements, we extracted two single summary statistics for describing cell line sensitivity to kinase inhibitors (Fig. 1C): Dose–response area under the curve (AUC) and half-maximal inhibitory concentration (IC₅₀). These properties are highly correlated with each other, having a Pearson’s correlation coefficient ∼0.9 (Fig. 1D). We extracted these properties from DepMap and matched them to our kinome states. The final integrated dataset has ∼250 drugs tested across ∼450 cell lines, representing ∼70,000 inhibitor-cell line combinations representing nearly all cancer types.

Examining bivariate association of features to cell line sensitivities provides a means for feature selection

The 450 cell lines tested in the DepMap dataset also have accompanying baseline RNAseq gene expression data, so we integrated the ∼20,000 TPM values for each cell line into the kinome-state and cell line sensitivity dataset. This adds baseline cell line-specific gene expression information to our cell line agnostic inhibitor-induced kinome states.

To identify which features were most correlated with drug sensitivity, we examined bivariate associations of each of the ∼21,000 features (kinome inhibition states, baseline gene expression, gene essentiality, protein expression and copy number variation) were compared individually to the outcome of cell line sensitivities (dose response AUC and IC₅₀), and their Pearson’s R correlation was calculated. We additionally calculated Spearman’s correlation coefficients, and ensured all correlated features were significant at the 0.5 level. The full table of feature correlations with significance values is attached in Table S2. We found that the most correlated feature is the drug-induced kinase inhibition state of TP53RK with a correlation coefficient R ∼0.3, while the most correlated baseline gene expression value was OGFRL1 with a correlation coefficient R ∼0.05. Overall, inhibitor-induced kinome states showed stronger correlation with cell line sensitivity metrics (Fig. 2A) despite there being 40× more baseline gene expression features than kinome states. However, it is important to note that this large difference in overall correlation may be partly due to the sample imbalance in abundance of kinome inhibition states (∼200 per cell line) compared to baseline gene expression (one per cell line).

After exploring the relationship between features and cell line sensitivity, we sought to use machine learning models to integrate these features into predictive models of cell line sensitivity to kinase inhibitor treatment. To create a feature selection metric, all features were ranked by the absolute value of their Pearson’s correlation to the cell sensitivity outcomes. Features were incorporated into models based on this ranking, i.e., including only a top-ranked subset of best-correlated features. Seeking models with the best balance of predictive power and feature number, we tested model performance with different numbers of features, again selected starting from the most informative feature. We tested models with 100, 200, 300, 400, 500, 1,000, 2,000, 3,000, 4,000, and 5,000 features to see how adding more and more lower ranked features affected model performance. Figure 2B shows the ranking for all the feature sets, and the order they are added into the model.

Machine learning models predict cancer cell line sensitivity from a combination of kinome inhibition states and baseline transcriptomics

To build machine learning models to predict cancer cell line AUC and IC₅₀ in response to treatment with kinase inhibitors, the highest ranked 100-5000 features were selected from the dataset linking drug-induced kinome states to cancer cell line responses (Fig. 2B). We compared three model types: LASSO regression, random forest and XGBoost. All models were trained with 10-fold random cross validation to minimize overfitting on the training data. Here, 10% of the dataset was randomly sampled and held out for testing, with the remaining 90% used for training. This was repeated ten times, and model performance was averaged across all the ten folds. This ensures comparable accuracy of the model predictions on new kinase inhibitors and cell lines. For each “batch” of feature numbers modelled from 100–5,000 we tuned sets of 30 hyperparameters for all model types (Fig. 3A). The R² value between predicted and actual value was utilized as the metric for model comparison. Overall, the 5000 feature XGBoost model performed the best with a cross-validation R² of ∼0.7 (Fig. 3B).

Since tree-based machine learning models like XGBoost offer in-built explainability, it is possible to interrogate and explain which features were most important in predicting the outcome of cell line sensitivities. These importances generated via Shapley values (Rozemberczki et al., 2022) show kinase inhibition states to be overwhelmingly more important for predicting cell line responses when compared to baseline gene expression. Kinases involved in cell cycle and proliferation are overrepresented in the top 25 features (Fig. 3C) (MAP2K, MEK2, CDKL5 etc.), and we further annotated each kinase with whether they are considered by current NIH guidelines as understudied (“Dark”) or well-characterized (“Light”) (Berginski et al., 2021b; Oprea et al., 2018). Interestingly four kinase interactor proteins are included as well, suggesting that interactions between inhibitors and non-kinases (off-target effects) have important consequences for cell viability. Baseline gene expression features show much lower model importances (Fig. 3D), but 40% of the top 25 genes have known interactions with kinases whose inhibition states are used in the model.

Inclusion of various multi-omics data with kinome inhibition states and gene expression did not improve model predictive performance

In addition to the baseline gene expression data, all the cell lines in the DepMap database have three other profiling data types available: copy number variation, gene essentiality from CRISPR/KO, and baseline proteomics. To see if inclusion of these data into models would improve predictions, we integrated these with the modeling dataset of kinome inhibition states and gene expression and used identical modeling strategies described above to select correlated features, build, and evaluate LASSO, random forest, and XGBoost models (Fig. S2). We found that adding in the various multi-omic data types did not significantly outperform the models limited to kinase inhibition states combined with baseline gene expression (R² of ∼0.69 for predicting IC₅₀).

Experimental validation of model predictions were successful in characterized and novel cell lines

After fitting the model on 70,000 cell line-drug combinations, predictions were made on 1.2 million unseen (not previously seen by the model) drug-cell line combinations (Fig. 4A). Approximately 90% of the untested inhibitors were associated with KINOMEscan datasets. As an initial validation, we tested a subset of the predictions in well-characterized breast cancer cell lines, two unseen by the models and two previously seen (HER2 positive: SK-BR-3, BT-474 and two triple negative: SUM159, HCC1806). We analyzed the performance of the model on experimental data for unseen drug-cell line combinations, arriving at an R value ∼0.6 for all but one (SKBR3) breast cancer cell line (Fig. 4B). Notably, all the drugs tested had kinome profiling data from the Kinobeads assay.

We then further validated the model by predicting inhibitor effects from collected RNAseq data in tumor (two samples) and stroma (one sample) derived cell lines from PDAC patients (Lipner et al., 2020; Berginski et al., 2021a). Importantly, these patient-derived cell lines were profiled for baseline gene expression in-house and represent a challenging and highly heterogeneous transcriptional landscape which the model has not seen before. Dose–response AUC predictions were made by the model for 58 drugs with kinome profiling data from the Kinobeads assay and 18 drugs with kinome profiling data from the KINOMEscan assay. The model predicted AUCs were compared to experimentally generated AUCs, revealing an average R ∼0.5 for drugs with kinome profiling data from the Kinobeads assay tested in patient stroma-derived cell lines (Fig. 4C), and R ∼0.4 for drugs with kinome profiling data from the KINOMEscan assay. On the other hand, in patient tumor-derived cell lines, drugs from kinobeads had model accuracy R ∼0.49, and drugs from KINOMEscan had R ∼0.3 (Fig. 4C).