Microsatellite instability

We first tested AKLIMATE on predicting microsatellite instability in the colon and rectum adenocarcinoma (COADREAD) and uterine corpus endometrial carcinoma (UCEC) TCGA cohorts. Microsatellite instability (MSI) arises as a result of defects in the mismatch repair machinery of the cell. Tumors with MSI (often accompanied by higher mutation rates) represent a clinically relevant disease subtype that is associated with better prognosis. MSI is also an immunotherapy indicator as such tumors produce more neoantigens. MSI can be predicted with high accuracy from expression alone [7], providing a straightforward benchmark for AKLIMATE performance with a single feature type in a binary classification setting.

We used expression data (upper-quantile normalized RSEM counts of RNA-Seq) and MSI annotations as described in [8] for the COADREAD and UCEC TCGA cohorts. The UCEC cohort consisted of 326 patients, of which 105 exhibit high microsatellite instability (MSI-H) and the remaining 221 are classified as either low (MSI-L) or stable (MSS). The COADREAD cohort included 261 samples, with 37 MSI-H and the remaining 224 classified as either MSI-L or MSS. In both tumor types, we trained models to distinguish MSI-H patients from MSI-L+MSS patients on 50 phenotype-stratified partitions of 75% training and 25% test folds. We then computed area under the ROC curve (AUROC) for each set of test fold predictions.

We compared AKLIMATE to Bayesian Multiple Kernel Learning (BMKL) because it performed well in several DREAM challenges, in particular winning the NCI-DREAM Drug Sensitivity Prediction Challenge [9]. Furthermore, its pathway-informed extension [7] shares several similarities with AKLIMATE—it is a multiple kernel learning method that operates on pathway-derived kernels. In particular, BMKL uses expression-based Gaussian kernels computed on features from the Pathway Interaction Database (PID) pathway collection [10]. We tested four versions of BMKL—sparse single-task BMKL(SBMKL), dense single-task BMKL(DBMKL), sparse multi-task BMKL (SBMTMKL), and dense multi-task BMKL (DBMTMKL). Sparse BMKL models are the focus of [7]—they use sparsity-inducing priors to train models with few non-zero kernel weights (we used the hyperparameters specified in [7]). We added DBMKL and DBMTMKL to the comparison because dense MKL models (almost all kernels receive non-zero weights) tend to produce higher predictive accuracy in many experimental settings (e.g. [11]). Their parameters were identical to the ones for the sparse models except for (ζ_κ, η_κ), which were set to (999, 1) in the dense models ((ζ_κ, η_κ) = (1, 999) in the sparse ones). Finally, the single-task models were trained separately on the UCEC and COADREAD cohorts, while the multitask versions learned MSI status on the two TCGA cohorts jointly, with each cohort representing a separate task. All train/test splits were matched across methods. All BMKL models used 196 PID gene sets; model parameters, kernel computations and data filtering steps matched the setup in [7].

Since AKLIMATE uses a much larger gene set compendium (S = 17, 273 feature sets, see S1 Text), we controlled for this prior information imbalance as a possible source of performance bias. For that purpose, we created a reduced version (AKLIMATE-reduced) that is restricted to the same set of 196 input PID sets as used by the original publication of BMKL. We refer to the full unrestricted model of AKLIMATE as simply AKLIMATE in this comparison.

AKLIMATE predicted MSI status in the UCEC cohort significantly better than AKLIMATE-reduced or any of the BMKL models (Fig 1A, mean AUROCs 0.885±0.051, 0.92±0.029, 0.903±0.039, 0.92±0.03 for SBMKL, DBMKL, SBMTMKL, and DBMTMKL respectively). In particular, AKLIMATE achieved mean AUROC of 0.962±0.019 compared to 0.938±0.031 for AKLIMATE-reduced (P < 4.1e − 08; paired Wilcoxon signed rank test), suggesting a predictive benefit to AKLIMATE’s larger collection of gene sets. A larger gene set collection is both more likely to contain sets derived specifically to describe the MSI process and more flexible in terms of the possible combinations of component gene sets. Indeed, the most informative feature set according to AKLIMATE was “MSI Colon Cancer” (16.9% relative contribution to model explanatory power)—a gene expression signature for MSI-H vs MSI-L+MSS in COADREAD cohorts [12] (Fig 1B). Futhermore, the next two most informative sets were “GO DNA Binding” (4.55% relative contribution) and “REACTOME Meiotic Recombination” (2.4% relative contribution), both of which are strongly relevant to DNA mismatch repair (MMR). Of note, the MutL Homolog 1 (MLH1) gene was the top-ranked single feature in AKLIMATE (26.7% relative contribution) and was present in the ten top ranked gene set kernels (Fig 1B). MLH1 is a key MMR gene involved in meiotic cross-over [13]—loss of MLH1 expression, usually through DNA methylation, is known to cause microsatellite instability. This demonstrates that AKLIMATE was able to pinpoint individual causal genes as it sifts through thousands of gene sets. In contrast, the meta-pathway constructed by AKLIMATE-reduced represents a poorer approximation to the underlying biological process, as evidenced by its lower AUROC. This is likely due to the fact that, out of the 196 PID pathways used, only “PID P53 Downstream” (35.2% relative contribution) contained MLH1 S2 Fig), limiting the influence of this key gene on the prediction task.

Fig 1

AKLIMATE performance on predicting MSI in UCEC TCGA.

(A) Performance of AKLIMATE and BMKL on classifying MSI-H vs MSI-L+MSS. AUROC computed for 50 75%/25% stratified train/test splits. P-values for paired Wilcoxon signed rank test. Methods compared: aklimate—AKLIMATE on full collection of feature sets; aklimate-reduced—AKLIMATE with 196 PID pathways; sbmkl—sparse single-task BMKL; dbmkl—dense single-task BMKL; sbmtmkl—sparse multi-task BMKL; dbmtmkl—dense multi-task BMKL. Multi-task BMKL models are trained to simultaneously predict MSI status on UCEC and COADREAD cohorts. (B) Top 10 predictive AKLIMATE feature sets and top 50 predictive features. Expression of top 50 features (left heatmap); Membership of most predictive features in most predictive feature sets (right heatmap). Features are organized by KNN clustering into 3 groups, followed by hierarchical clustering within each cluster. Feature set model weights scaled to sum up to 1 (barplot, top of right heatmap). Feature model weights scaled to sum up to 1 (barplot, right of right heatmap). Feature and feature set weights averaged across 50 train/test splits.

Interestingly, even though AKLIMATE-reduced used the same feature sets as BMKL, it achieved a statistically significant improvement over all BMKL varieties (Fig 1A). This included both the non-sparse and multi-task BMKL varieties, despite the fact that the latter drew benefit from an entire additional training set of COADREAD data. In this case, the difference in kernel representations may have contributed to the improved performance of AKLIMATE-reduced (see Discussion). For the COADREAD MSI classification task, we found that all methods performed equally well, with AKLIMATE-reduced having marginally lower accuracy compared to DBMKL & DBMTMKL and the full AKLIMATE model showing marginally higher mean AUROC compared to all BMKL varieties (S3 Fig). We suspect the strong MSI signal in the feature data makes this prediction task somewhat more facile, limiting our ability to discriminate among methods based on their predictive performance.

Breast cancer survival

For our second benchmark, we considered the task of predicting survival in the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort [14]. This problem is much more challenging than predicting MSI status, as demonstrated by the DREAM Breast Cancer Challenge [15] and elsewhere [16]. Another difference between this task and MSI inference is that the METABRIC cohort is annotated with curated clinical data. In fact, the clinical features are quite informative for survival prediction—pre-competition benchmarking by the DREAM Challenge organizers found that models that used exclusively clinical features significantly outperformed ones that used only genomic features, and performed only marginally worse than models in which clinical features were augmented by a subset of molecular features selected through prior domain-specific knowledge [17]. In addition, the best pre-competition clinical and molecular feature model had better accuracy than all but the top 5 models in the actual challenge [15].

A breast cancer model that foregoes the use of clinical data would clearly suffer from inferior performance as well as reduced relevance in real-world medical settings. To achieve clinical data integration, AKLIMATE introduces a special category of “global” features, which are added to the “local” features of each feature set prior to the construction of its corresponding Random Forest. Global features can therefore be interpreted as a uniform conditioning step applied to all AKLIMATE component Random Forests. In our METABRIC analysis, all clinical features are treated as global.

We compared AKLIMATE to two state-of-the art METABRIC survival predictors. The first one, which we refer to as the “Breast Cancer Challenge” (BCC) method, was the top-performer in the Sage Bionetworks–DREAM Breast Cancer Prognosis Challenge [15, 18]—an ensemble of Cox regression, gradient boosting regression, and K-nearest neighbors trained on different combinations of clinical variables and molecular-feature derived metagenes. The second one is a multiple kernel learning (MKL) method—Feature Selection MKL (FSMKL) [16]—that implements a pathway-informed extension of SimpleMKL [19] using linear and polynomial kernels created from clinical data and molecular features in pathways of the Kyoto Encyclopedia of Genes and Genomes (KEGG) [20]. In FSMKL, pathway features from different data types lead to the construction of separate kernels—in the case of METABRIC, each pathway produces distinct expression and copy number kernels (in contrast, AKLIMATE learns one kernel matrix from the combined pathway features across all data types). In addition, FSMKL treats each clinical feature as a singleton pathway that produces an individual kernel. To make our results directly comparable to FSMKL and BCC as presented in [16], we used a subset of the patient cohort (N = 639) and a reduced set of clinical variables to match the dataset used in that publication (see S1 Text).

We cast the problem as a classification task where molecular and clinical input features are used to predict whether a patient is alive or not at the 2000 day mark. Based on the 2000 day cutoff, there were 387 survivors and 252 non-survivors in the reduced cohort. Similar to the MSI analysis, we performed 50 stratified repeats of 80% train and 20% test partitions; AKLIMATE was trained on each training split and its accuracy computed on the respective test samples. To decrease computation time, AKLIMATE’s kernel construction step used 1000 trees instead of the default 2000 (see S1 Text for all other hyperparameter settings).

The full AKLIMATE model had higher mean accuracy than BCC (74.1±3.3% vs 73.3±0.2%) and was on par with FSMKL (74.1±3.3% vs 74.2±1.8%) with a slightly higher standard error across cross-validation folds compared with FSMKL (S4 Fig). The AKLIMATE accuracy was largely unchanged when it was restricted to use the same KEGG pathways as FSMKL (S5 Fig). All three algorithms selected clinical variables as important for prediction.

Of note, two of AKLIMATE’s main advantages were not fully utilized under these experimental settings. First, the benefit of incorporating prior knowledge is reduced because of the relatively low information content of genomic features. Second, the clinical variables may lack complex interactions, which may explain why modeling them with simpler linear techniques is just as effective. In fact, the most explanatory feature in the full AKLIMATE model (16% relative contribution) was the Nottingham Prognostic Index (NPI) [21]—a linear combination of tumor size, tumor grade, and number of lymph nodes involved. Even with these disadvantages, however, AKLIMATE achieved performance on par with two state-of-the-art methods that were specifically optimized for the METABRIC data.

In this task, clinical information proved to be the most influential data type for survival prediction. AKLIMATE models with clinical features alone were more accurate than AKLIMATE models with genomic features alone (p-val = 1.9e-08, paired Wilcoxon signed rank test, Fig 2A). This is even more striking considering the models used tens of thousands of genomic features versus only 15 clinical ones. However, the AKLIMATE models using both clinical and genomic features outperformed each single-component model (p-val = 1.4e-09 for full versus genomic, p-val = 6.7e-04 for full versus clinical, paired Wilcoxon signed rank test, Fig 2A), suggesting that the inclusion of genomic features contained complementary signals with respect to the clinical variables. The mean relative contributions of each data type in the full models were 61.6% for clinical, 29.6% for expression, and 8.8% for copy number.

Fig 2

AKLIMATE performance on predicting survival at 2000 days in the METABRIC cohort.

(A) Performance of AKLIMATE under different data type combinations. EXP+CNV—AKLIMATE with genomic features only; clinical—an RF model run with the clinical variables only; EXP+CNV+CLINICAL—AKLIMATE with genomic features as “local” variables and clinical features as “global” ones. FSMKL and BCC dashed lines show mean performances for the two models under 5-fold cross-validation as shown in [16]. (B) AKLIMATE results highlighting the top 10 predictive feature sets and top 50 predictive features. Figure organized as Fig 1. Clinical variables shown as column annotations; included only if among the 50 most informative features in the model. Clinical variables are ranked from top to bottom by relative predictive contribution. Survival status is a binary variable representing survival at 2000 days (labels) while days survived shows actual duration of survival. Samples sorted by days survived within the two classes. Feature and feature set weights averaged across 50 train/test splits.

Importantly, while the full AKLIMATE and FSMKL models achieved similar mean accuracy, the influence on prediction of individual features and pathways in each model were quite different. AKLIMATE heavily favored clinical variables, with 54.5% of the model’s relative explanatory power carried by just five features—NPI, tumor size, lymph node involvement, age and treatment (Fig 2B). FSMKL ranked tumor size as the most important clinical feature (second overall), with other clinical variables also generally considered relevant—e.g. NPI(9th), age (11th), histological type(14th), tumor group(34th) and PAM50(40th) [16]. Most clinical variables, however, were considered less informative than FSMKL’s top ranked KEGG pathway-based kernels. FSMKL’s highest weighted kernel was “Intestinal Immune Response for IgA production”, followed in importance by “Arachidonic acid metabolism”, “Systemic lupus erythematosus”, “Glycerophospholipid metabolism”, and “Homologous recombination”—all of these KEGG-based kernels scored higher than any clinical variable with the exception of tumor size [16].

AKLIMATE’s most informative feature sets were enriched for breast cancer progression and response to treatment, with 3 of the top 10 and 8 of the top 20 (S1 Table) related to these functional groups. For example, the most informative feature set (“BREAST CANCER ER-/PR- DN”) represents a signature that is correlated with reduced protein abundance of the estrogen (ER) and progesterone (PR) hormone receptors [22]. Similarly, the 9th most informative pathway (“BREAST CANCER HER2 ENDOCRINE THERAPY RESISTANCE”) [23] captures transcriptome changes associated with the development of resistance to targeted therapies. Both of these signatures serve as proxies for highly relevant information not available for the METABRIC cohort (protein activity for PR and therapy resistance). Furthermore, the latest American Joint Committee on Cancer breast cancer staging manual introduces tumor grade and ER/PR/HER2 receptor status among the key breast cancer biomarkers, which already include tumor size and lymph node engagement [24]. All of these appeared as highly informative AKLIMATE model features, either directly or via a proxy genomic signature.

Breast cancers are typically subtyped into one of several main pathological groups that heavily influence treatment selection (e.g. luminal / ER-positive, HER2-amplified, basal / triple-negative, normal-like). Appreciable differences in AKLIMATE’s METABRIC performance can be seen when the samples are divided up by their subtypes (based on PAM50 signatures [25], S6 Fig). The basals were associated with the least accurate (65.6% ± 8%, n = 131) and the luminal-As with the most accurate (80.8% ± 5.4%, n = 200) results, a trend that is expected given the known higher heterogeneity levels of less differentiated basal relative to luminal-A tumors.

shRNA knockdown viability

We tested AKLIMATE’s ability to integrate multiple data types and solve regression tasks to predict cell viability post shRNA knockdown. In this case, the prediction labels were continuous values representing cell line survival after the shRNA-mediated mRNA degradation of a particular gene. Gene profiles were computed with the Analytic Technique for Assessment of RNAi by Similarity (ATARiS) method [26] that creates consensus scores by combining viability phenotypes from multiple shRNAs targeting the same gene. We selected 37 such profiles for different genes across 216 cancer cell lines from the Cancer Cell Line Encyclopedia (CCLE) [27]. We chose these 37 tasks (out of 5711 available consensus profiles) because they had at least 10 cell lines showing strong (>2 standard deviations from mean) knockdown viability response and were in the top quartile by variance of all consensus profiles (see S1 Text). Based on the results of the DREAM9 gene essentiality prediction challenge [28], we expected this task to be more difficult than the other case studies.

We used expression and copy number measurements from CCLE [29] as predictive features. We augmented these two data types by adding discrete gene copy number calls made by GISTIC2 [30] and activities for 447 transcriptional and post-transcriptional regulators inferred by hierarchical VIPER [31, 32]. We focused on the 206 cell lines for which we had knockdown profiles, copy number and expression features.

We compared AKLIMATE’s performance to three of the five top performing methods in DREAM9 [28] as well as three baseline algorithms. We briefly describe the DREAM9 subchallenge 1 top performers next. Multiple Pathway Learning (MPL) [31] took 5th place (see https://www.synapse.org/#!Synapse:syn2384331/wiki/64760 for challenge results)—it used elastic-net regularized Multiple Kernel Learning with Gaussian kernels based on feature sets from the Molecular Signature Database (MSIGDB) [3]. MPL and Random Forest Ensemble (MPL-RF) took 2nd place—its prediction was based on averaging a Random Forest classifier with MPL. MPL and MPL-RF were our contributions to the DREAM9 challenge. Both methods were run with the same hyperparameters and pathway collections as in DREAM9 [28, 31]. Kernelized Gaussian Process Regression (GPR) took 3rd place—it used extensive filtering steps to reduce the input feature dimensionality, followed by principal component analysis, and finally Gaussian Process regression with covariance computed from the principal components [28] (see also https://www.synapse.org/#!Synapse:syn2664852/wiki/68499 for implementation and model description). We downloaded the code from the Synapse URL and ran it with the published DREAM9 hyperparameters.

To provide performance baselines, the DREAM9 winners were augmented by standalone Random Forest (RF), Generalized Linear Model (GLM) with lasso penalty (GLM-sparse), and GLM with L2 regularization (GLM-dense). RF was run with the ranger R package [33] with the following hyperparameters—sampling without replacement with 70% of the samples used for tree construction, minimum node size of 10, 1500 trees and 10% of the features randomly sampled for each node split. GLM-dense and GLM-sparse were run using the glmnet R package [34] with the response family set to “gaussian” and the strength of regularization λ optimized through cross-validation. The elastic net tradeoff between the lasso and ridge penalties α was set to α = 0.8 (GLM-sparse) and α = 0.001 (GLM-dense).

We did not include the nominal first place winner of DREAM9—an ensemble of four kernel ridge regression models with kernels trained through Kernel Canonical Correlation Analysis and Kernel Target Alignment—because we could not re-run the source code supplied with the challenge submission. We feel this omission is not material as the top 3 methods were declared joint co-winners—their results were shown to be statistically indistinguishable from each other but separable from the rest of the entries [28]. Furthermore, experiments in [31] suggest that this method underperforms MPL, MPL-RF and GPR when only high-quality shRNA knockdown profiles are considered.

To save computational time, we compared methods on a single stratified train/test split for each ATARiS profile (a different split for each profile) where 67% of the cell lines were used for training and 33% were withheld for testing. Each method was run with its recommended parameters and filtering steps—if no filtering steps were specified, the method used all available features. AKLIMATE’s prediction binarization quantile was set to its default value of q = 0.05 (see Methods, also S1 Text for all other settings).

AKLIMATE achieved average Root Mean Squared Error (RMSE) of 1.031 vs 1.048 for GPR, 1.056 for RF, 1.066 for GLM-dense, 1.07 for MPL, 1.071 for GLM-sparse and 1.081 for MPL-RF. While the improvements were modest, the mean RMSE difference was statistically significant in all but one case (Fig 3A). AKLIMATE was also the top performer when we considered the number of times an algorithm achieved the best RMSE on an individual prediction task (Fig 3B, AKLIMATE retained top spot under non-RMSE metrics as well—S7 Fig). AKLIMATE performed better than average across nearly all tasks (S8 Fig); its advantage was particularly pronounced in predicting the essentiality of key regulators—Beta-Catenin (CTNNB1), Forkead Box A1 (FOXA1), Mouse Double Minute 4 Homolog Regulator of P53 (MDM4), Phosphatidylinositol-4,5-Bisphosphate 3-Kinase Catalytic Subunit Alpha (PIK3CA)—or housekeeping genes—Proteasome 26S Subunit ATPase 2 (PSMC2), and the fifth ATPase subunit (PSMC5). As these gene classes are heavily studied and thus over-represented in our pathway compendium, AKLIMATE’s enhanced accuracy may be due to the relatively higher abundance of relevant prior knowledge.

Fig 3

Prediction of cell line viability after shRNA gene knockdowns.

(A) RMSEs of AKLIMATE and competing methods on 37 consensus viability profiles from the Achilles dataset. P-values for RMSE comparisons from paired Wilcoxon signed rank test. Methods: Random Forest (RF), Gaussian Process Regression (GPR), Multiple Pathway Learning (MPL), ensemble of MPL and Random Forest (MPL-RF), L2 regularized linear regression (GLM-dense), L1 regularized linear regression (GLM-sparse). (B) Number of times an algorithm produced the best RMSE on a prediction task. To prevent small relative RMSE differences from having a biasing effect on the win counts, for each task we consider all algorithms with RMSE within 1% of the min RMSE to be joint winners. For that reason, the total win counts add up to more than the number of regression tasks. (C) AKLIMATE top 10 informative feature sets and top 50 informative features for the task of predicting MDM4 shRNA knockdown viability. Organized as Fig 1. (D) RMSEs of KRAS AKLIMATE models with and without the use of mutation profiles for eight key regulators. Results shown for 10 matched stratified train/test splits where 80% of the cohort is used for training and 20% for testing. (E) AKLIMATE top 10 informative feature sets and top 50 informative features for the task of predicting KRAS shRNA knockdown viability when mutation features are used (feature and feature set weights averaged over 10 train/test splits). Organized as Fig 1.

For example, AKLIMATE’s ability to predict MDM4 knockdowns benefited from MDM4’s function as a p53 inhibitor—many of the most informative feature sets in the MDM4 model relate to p53’s regulome and its role in controlling apoptosis, hypoxia and DNA damage response (Fig 3C). The top features were also functionally linked to p53—Cyclin Dependent Kinase Inhibitor 1A (CDKN1A, 28.6% relative contribution) is a kinase that controls G1 cell cycle progession and is tightly regulated by p53; Mouse Double Minute 2 homolog P53 binding protein (MDM2, 8% relative contribution) participates in a regulatory feedback loop with p53; Zinc Finger Matrin-Type 3 (ZMAT3, 3.7% relative contribution of expression; 3.5% of copy number) interacts with p53 to control cell growth. In addition, the four p53 features from each data type were all present among the 50 most informative ones. Individually they carried little signal (relative contribution: 0.5% expression, 0.3% copy number, 0.3% inferred protein activity, 0.2% GISTIC score) but taken together they clearly implicated p53 as one of the top 10 most informative genes. The ability to capture such multi-omic interactions is one of AKLIMATE’s main strengths, made possible by the use of RF kernels that fully integrate all data types. The synergy between different p53-based features is impossible to observe in methods that assign individual kernels to each data type (e.g. BMKL, FSMKL and MPL). Encouragingly, AKLIMATE dominated MPL/MPL-RF on almost all tasks even though they share the same MKL solver (Figs 3A and S8 Fig).

FOXA1 knockdown offered another example of AKLIMATE showing superior predictive performance on a well-studied gene. FOXA1 dysregulation is an essential event in breast cancer progression and subtype characterization. Due to breast cancer’s prevalence and clinical importance, there is an extensive list of relevant signatures in our feature set compendium. Eight of the top 10 (and 12 of 14 overall) feature sets in the FOXA1 AKLIMATE model were directly related to breast cancer experiments under different conditions (S10 Fig). As expected, FOXA1 features were the most informative (relative contribution: 31.5% FOXA1 expression; 3.4% FOXA1 inferred activity), with the androgen receptor (AR) also among the top 5 most informative genes (relative contribution: 2.8% AR expression).

Out of the 37 shRNA prediction tasks we considered, shRNA knockdown of the Kirsten Rat Sarcoma Viral Oncogene Homolog (KRAS) was the most obvious example of a task involving a well-characterized gene that did not experience discernible accuracy improvement over competing methods (S8 Fig). Our hypothesis is that a true biological “driver” is absent from the set of molecular features presented to AKLIMATE. To test this, we added a mutation data type containing information for 8 key regulators (KRAS, NRAS, PIK3CA, BRAF, PTEN, APC, CTNNB1 and EGFR), 3 of which (KRAS, CTNNB1, PIK3CA) have knockdown profiles among the 37 shRNA prediction tasks. Even with the addition of these limited 8 features we observed dramatic improvement in KRAS shRNA prediction accuracy across all metrics (Fig 3D, mean RMSE 0.918±0.025 vs 1.02±0.024; mean Pearson 0.665±0.023 vs 0.529±0.025; mean Spearman 0.57±0.039 vs 0.453±0.034). Furthermore, the KRAS mutation feature was by far the most informative (23.5% relative contribution, Fig 3E), followed by KRAS copy number (6.9%), KRAS GISTIC (3.3%) and KRAS expression (3.1%). While KRAS expression, copy number and GISTIC features all appeared among the 50 most informative ones in the “no mutation” run, their combined relative contribution was only 3.33% (1.5% copy number, 1% expression, 0.8% GISTIC, S11 Fig). Thus, the addition of the KRAS mutation feature was not only key to improving the predictive performance of the model on its own, but it also helped “lift” the influence of KRAS features form other data types.

Our training features are measured at the gene level, but AKLIMATE has no constraints on the inclusion of more granular data. For example, it can evaluate the importance of mutations at particular amino acid positions or the type of resulting substitutions. In the case of KRAS, glycine replacement by either aspartic acid or valine in the 12th amino acid position appears to have the biggest negative effect on cell viability post-shRNA knockdown (Fig 3E). G12 is a well-known KRAS mutation hotspot [35]—AKLIMATE’s ability to prioritize relevant hotspots can be a key advantage in modeling drug response or recommending treatment strategies.

The addition of mutation data did not yield any predictive benefit in modeling PIK3CA or CTNNB1 post-knockdown viability (S12 and S13 Figs). CTNNB1 had only nine mutations in the cohort, all of them in different codons. PIK3CA had 30 mutations, with some hotspots—the lack of improvement in this case may be due to the change in protein sequence not being biologically relevant, the ability of other genomic features to fully capture the mutation signal, or the fact that the “no mutation” models were quite accurate to begin with.

Background

Before describing AKLIMATE in detail, we briefly summarize existing approaches to which AKLIMATE’s design and performance can be compared. There are three main challenges inherent to bioinformatics algorithms for supervised and unsupervised learning—prior knowledge integration, heterogeneous data interrogation, and ease of interpretation. New methods try to address some aspects of these three challenges [41–43]. Several common approaches emerge for supervised learning tasks. A popular way of increasing interpretability is to train models with regularization terms that constrain the number of included features—sparse models are considered easier to interpret than ones containing thousands of variables. Common regularization penalties are the lasso [44] and the elastic net [45]. More sophisticated regularization schemes can control the model behavior at the feature set level (e.g. the group lasso (GL) [46] and the overlap group lasso (OGL) [47]), allowing the incorporation of prior knowledge in the form of feature sets. However, both have drawbacks that make them less suitable for biological analysis where genes often participate in multiple processes. GL requires that if a feature’s weight is zero in one group, its coefficients in all other groups must necessarily be zero. OGL tends to assign positive coefficients to entire feature sets, making it less suitable in situations where only some of a set’s features are relevant.

Network-regularized methods [48, 49] represent a different approach where gene-gene interaction networks are used as regularization terms on the L₂-norm of feature weights. While they do use pathway-level information and are straightforward to interpret, they generally focus on an individual data type.

Multiple Kernel Learning (MKL) approaches [7, 16, 19, 50–52] can incorporate heterogeneous data by mapping each set of features through a kernel function and learning a linear combination of the kernel representations. Each kernel represents distinct sample-sample similarities providing flexible and powerful transformations to access either explicit or implicit feature combinations. To prevent overfitting, MKL methods generally include a regularization term—e.g., L₁ sparsity-inducing norm on the kernel weights [19] or the elastic net [51]. Prior knowledge can also be integrated by constructing individual kernels from a pathway’s member features within each data type [7, 16, 28, 31, 52]. Indeed, MKL methods with prior knowledge integration [9, 28] have won several Dialogue on Reverse-Engineering Assessment and Methods (DREAM) [53] challenges, including a predecessor of the approach described here [28, 31]. Nevertheless, MKL suffers drawbacks when the contributions of individual input features need to be evaluated. Except in trivial cases, it is generally impossible to assign importance to the original features once the method is trained in the kernel function feature space, thus limiting interpretability of solutions. In addition, feature heterogeneity necessitates the construction of separate kernels for each data type, limiting the ability of MKL to capture cross-data-type interactions.

All of these methods can be used as components in more complex ensemble learning models. Ensembles combine predictions from multiple algorithms into a more robust, and often more accurate, “wisdom of crowds” final prediction. The simplest and most common ensemble technique is averaging the predictions of component models. Averaging over uncorrelated models or models with complementary information can improve performance—it is one of the main reasons for the emergence of collaborative competitions such as the DREAM challenges [9, 54]. In fact, such ensemble methods often win DREAM challenges [18] or outperform competitors in genomic prediction tasks [55]. While they provide a boost in predictive accuracy, their interpretation is quite challenging—ensembles often combine different model types, making the computation of input feature importance impossible in the large majority of cases. A notable exception is the Random Forest (RF) [56]—an ensemble of decision trees that has been widely applied to bioinformatics problems.

A more general approach to combining the predictions of multiple models is model stacking [57]. In stacking, component (base) learner predictions are used as inputs to an overall (stacked) model which produces the final calls. To reduce overfitting and improve generalization, base learner predictions are generated in a cross-validation manner—a sample’s predicted label comes from a model trained on all folds except the one that includes the sample in question. Importantly, if the stacked model is a weighted combination of the base models (a Super Learner), it is asymptotically guaranteed to perform at least as well as the best base learner or any conical combination of the base learners [58, 59]. Stacked models exhibit identical pros and cons as standard ensemble models, with diminished interpretability traded off for increased accuracy. In some cases, however, interpretability can be tractable—e.g. [60] uses RF base learners with a stacked least square regression to compute a final weighted average of the base RF predictions. Although not examined in [60], we demonstrate that a similar setup can lead to an intuitive derivation of feature importance scores.

Overview

AKLIMATE is a stacked learning algorithm with RF base learners and an MKL meta-learner. Each base learner is trained using only the features from a pre-defined feature set representing a known biological concept or process—e.g. a biological pathway, a chromosomal region, a drug response or shRNA knockdown signature, or disease subtype biomarkers (Fig 4). Each feature in a feature set can be associated with multiple data modalities—i.e. a gene will have individual features corresponding to its copy number, mutation status, mRNA expression level, or protein abundance. Furthermore, although feature sets normally consist of genes, they are easily extendable to a more comprehensive membership. For example, they can be augmented with features for mutation hotspots, different splice forms, or relevant miRNA or methylation measurements. AKLIMATE’s MKL meta-learner trains on RF kernel matrices (see the RF Kernel section), each of which is derived from a corresponding base learner (Fig 5). An RF kernel captures a proximity measure between training samples based on the similarity of their predicted labels and decision paths in the trees of the RF. The MKL learning step finds the optimal weighted combination of the RF kernels. The MKL optimal solution can be interpreted as the meta-kernel associated with the most predictive meta-feature set derived from all interrogated feature sets.

Fig 4

Overview of AKLIMATE.

AKLIMATE takes as inputs multiple data types and a collection of feature sets. (A) Each feature set is overlapped with the available data types to identify all features that map to it (some feature sets might not extract features from all data types). (B) A feature-set-specific Random Forest model is trained on the multi-modal data identified in A. (C) RF models are ranked based on their predictive performance (metrics for mock RF models in A have matching color). (D) The top ranked RF models are converted to RF kernels (mock kernels shown for RF models with top 3 performance scores). (E) An elastic net MKL meta-learner finds the optimal combination of RF kernels and computes the final predictions. Elastic net hyperparameters are optimized via cross-validation.

Fig 5

RF kernel matrix construction.

The RF Kernel matrix is a geometric mean of the Hadamard product of three component similarity matrices. Each component captures a different aspect of a RF model. K₁ captures the similarity over RF tree predictions for sample labels (in the case of classification, probability of belonging to a class); two samples (A and B) show different predicted probabilities of belonging to the positive class (probability estimate in parentheses). K₂ represents the similarity over RF tree leaf indices to which samples are assigned; predicted leaf indices shown in parentheses for the two samples. K₃ reflects the proportion of times samples are assigned to the same RF tree leaf; e.g. the two samples end up together in one out of three example trees.

The key contributions of AKLIMATE are

The introduction of an integrative empirical kernel function that combines similarity in predicted labels with proximity in the space of RF trees (RF kernel).

The extension of kernel learning to a stacking framework. To our knowledge, AKLIMATE is the first stacked learning formulation to incorporate base kernels.

In the next sections we give a detailed description of each AKLIMATE component.

Random forest

An RF learner is an ensemble of decision trees each of which operates on a perturbed version of the training set. The perturbation is achieved by two randomization techniques—on one hand, each tree trains on a bootstrapped data set generated by drawing samples with replacement. On the other hand, at each node in a decision tree, only a subset of all features are considered as candidates for the next split. Tree randomization tends to reduce correlation among predictions of individual trees at the expense of increased variance of error estimates [61]. However, due to the averaging effect of ensembles of decorrelated models, it generally reduces the overall error variance of the ensemble RF [61]. That effect is partially negated by an increase in prediction bias, but broadly speaking the more decorrelated the trees the better the RF performance [61].

As each tree is constructed from a subsample of the training data, there is a tree-specific set of excluded, or out-of-bag (OOB) samples. The OOB predictions for the training set of a regression task are computed by walking each sample along the subset of trees in which that sample is OOB and averaging their predictions:

It has been shown that the OOB error rate provides a very good approximation to generalization error [56]. For that reason, AKLIMATE uses RF kernels based on OOB predictions as level-one data for the MKL meta-learner training (see RF Kernel and Stacked Learning below).

Kernel learning

Kernel learning is an approach that allows linear discriminant methods to be applied to problems with non-linear decision boundaries [62]. It relies on the “kernel trick”—a transformation of the input space to a different (potentially infinite dimensional) feature space where the train set samples are linearly separable. The utility of the “kernel trick” stems from the fact that the feature map between the input space and the feature space does not need to be explicitly stated—if the kernel function is positive definite (i.e. it has an associated Reproducing Kernel Hilbert Space), we only need to know the functional form of the feature space dot product in terms of the input space variables [62, 63]. A kernel function represents such a generalized dot product (Note: We define positive definiteness as in [62]—a kernel is positive definite if $$ {\sum }_{i=1}^n{\sum }_{j=1}^n{c}_i{c}_{j}K(x_i,x_j)\ge 0$$ for $$ \forall x_i,x_j\in \mathcal{X}$$ and $$ c_1,\dots ,c_n\in \mathbb{R}$$).

Commonly used kernel functions generally have an explicit closed form, such as a polynomial or a radial basis function. Kernel functions that encode more complex object relationships also exist, particularly if the objects can be defined in a recursive manner (ANOVA kernels, string kernels, graph kernels, see [64]). However, knowing the closed form kernel function is not a necessary condition—if the kernel positive definiteness can be verified and the kernel matrix of pairwise dot products (similarities) can be computed, we can still utilize the “kernel trick”. In fact, the Representer Theorem [62, 65] guarantees that the optimal solution to a large collection of optimization problems can be computed by kernel matrix evaluations on the finite-dimensional training data, ensuring the applicability of kernel-based alogirthms to real-world problems. More formally, for an arbitrary loss function L(f(x₁)…f(x_N)), the minimizer f* of the regularized risk function

AKLIMATE uses the dependency structure of a RF trained on a (possibly multi-modal) feature set to define an implicit empirical kernel function. It then computes an RF kernel matrix (see RF Kernel below) of pairwise training set similarities to use in a kernel learning algorithm. As different feature sets can contribute complementary information for a given prediction task, AKLIMATE utilizes a Multiple Kernel Learning approach to integrate the available RF kernels into an optimal predictor.

Multiple kernel learning

Kernel learning can lead to large improvements in accuracy, provided that an optimal kernel is selected. That choice can be difficult, however—the best kernel function for a multi-modal data set may be non-obvious, may entail tuning several hyperparameters, or may even lack a closed form expression. MKL addresses the problem of kernel selection by constructing a composite kernel that is a data-driven optimal combination of candidate kernels [50, 66–68]. Linear combinations of kernels of the form $$ \mathit{K}\left(\alpha \right)={\sum }_{i=1}^M{\alpha }_i{\mathit{K}}_i$$ with either a conical (∀i, α_i ≥ 0) or convex (∀i, α_i ≥ 0 and $$ {\sum }_{i=1}^M{\alpha }_i=1$$) sum constraint on the kernel weights are by far the most commonly used because of their many desirable properties (although algorithms that allow non-linear combinations do exist—see [67]). Some important advantages are:

Conical/convex combinations of positive definite kernels are positive definite.

The composite kernel is associated with a feature space that is the concatenation of all individual kernel feature spaces.

The Representer Theorem is readily extendable to the the conical/convex combination case—the optimal solution f* takes the form $$ f^{*}={\sum }_{m=1}^M{\sum }_{i=1}^N{k}_m(x,x_i){\alpha }_{m,i}$$ [51, 62].

MKL algorithms admit different forms of regularization depending on what norm of the kernel weights is chosen. AKLIMATE uses an elastic-net [45] regularizer on the norm of the individual kernel weights of the form $$ {\lambda }_1{\sum }_{i=1}^M\Vert {\alpha }_m{\Vert }_{K_m}+{\lambda }_2{\sum }_{i=1}^M{\Vert {\alpha }_m\Vert }_{K_m}^2$$, where λ₁, λ₂ ≥ 0, α_m = (α_m,1…α_m,N)^T, and $$ \Vert {\alpha }_m{\Vert }_{K_m}=\sqrt{{\alpha }_m{K}_m{\alpha }_m}$$ is the definition of a kernel norm as in [51] and [69]. The elastic net regularization provides the flexibility to find both sparse and dense solutions with appropriately tuned λ’s—i.e. λ₁ > > λ₂ gives few non-zero kernel weights (sparse) while λ₁ < λ₂ gives many non-zero kernel weights (dense). Also note that λ₁ = 0 and λ₂ > 0 allows for a uniform kernel weight solution.

The explicit form of AKLIMATE’s optimization problem is:

RF kernel

The most common approach to kernel matrix evaluation is to specify an explicit data dependency model for the kernel function—e.g. polynomial, Gaussian, ANOVA or graph kernels [64]. However, choosing the dependency structure a priori can lead to lack of robustness, particularly when it is not obvious what the right dependency structure is. AKLIMATE uses a RF to create an empirical approximation of the true dependency model and to evaluate the (implicit) kernel function associated with it (the RF kernel). RF kernels are robust to overfitting, generalize well to new data, and facilitate the integration of signals across data types.

Defining a kernel through a RF is not a new idea—in fact the concept was introduced at the same time as RFs [70]. In particular, for a random forest of trees with an equal number of leaves and uniform predictions in each leaf, a positive definite kernel can be defined based on the probability two samples share a leaf [70]:

AKLIMATE’s RF kernel extends the original definition (Eq (6)) by incorporating two additional RF-derived statistics. The intuition behind the first one is that two samples predicted to have the same label within a tree should be considered alike even if they fall into different leaves (Fig 5, K₁ kernel). The reasoning for the second one is that earlier node splits in a tree generally separate more distinct sample groups, while later splits tend to fine tune the decision boundary, highlighting more subtle differences. Thus, irrespective of the predicted labels, two samples that end up in different tree depths—one early- and one late-split leaf—should be considered less similar than samples landing in two late-split leaves (Fig 5, K₂ kernel). Incorporating these three patterns, AKLIMATE’s RF kernel matrix is computed using the following steps:

Calculate similarity over predictions across RF trees:


where p_i = (p_i,1, p_i,2, …, p_{i, M}) is a vector of predictions for data point i from the M trees in the RF. p_i always has continuous entries—either the actual predictions in a regression setting, or the probabilities of class membership for classification problems. The ‖p_i − p_j‖² distances are divided by the scaling constant σ = max_i,j‖p_i − p_j‖² so that they map to the [0, 1] range. Exponentiation of the negative distances converts them to similarities.

Compute similarity over leaf node indices across RF trees:


where t_i = (t_i,1, t_i,2, …, t_i,M) is a vector of the leaf indices of sample i across the RF. Tree nodes are indexed starting from the base node and then sequentially across each depth level of the tree. The scaling constant σ is set in the same manner as in K₁.

Calculate similarity over leaf node co-occurrence using the finite approximation of Eq (6) with variable number of leaves per tree:


where the exponential transformation is added to keep K₃ on a similar scale to K₁ and K₂.

Calculate the RF kernel as the geometric mean of the element-wise product of K₁,K₂ and K₃:

Importantly, since the components K₁, K₂, and K₃ are positive definite, so too is K. K₁ and K₂ are Gaussian kernels over $$ \mathcal{P}\times \mathcal{P}$$,$$ \mathcal{P}\subseteq {\mathcal{R}}^m$$ and $$ \mathcal{T}\times \mathcal{T}$$,$$ \mathcal{T}\subseteq {\mathcal{I}}^m$$ respectively, which ensures their positive definiteness [64]. K₃ is the exponent of a positive definite kernel, i.e. a positive definite kernel as well [64]. Finally, K is the Hadamard product of three positive definite kernels raised to a positive power—both of these operations preserve positive definiteness [64].

Stacked learning

Stacked learning is a generalization of ensemble learning in which the stacked model (meta-learner) uses the prediction output of its components (base learners) as training data for the computation of the final predictions [57, 74]. What makes stacked learning unique is that the base learner predictions (level-one data) are generated in a cross-validated manner that excludes each sample from the training set that produces its predicted label. More specifically, if our training set (level-zero data) is $$ X=\{X_i:i=1\dots N,X_i\in {\mathcal{R}}^p\}$$ with labels $$ Y=\{Y_i:i=1,\dots ,N,Y_i\in \mathcal{R}\}$$ and we have a collection of base learners Δ = {Δ₁, …, Δ_S} then the stacked generalization proceeds as follows [59, 75]:

Randomly split the level-zero data into V folds of roughly equal size—h₁, …, h_V (V-fold cross validation). Note that each h_v defines a set of indices that select a subset of the samples in X.

For each Δ_s base learner, train V models, collectively denoted as $$ {\widehat{\Delta }}_s=\{{\widehat{\Delta }}_{s,1}\dots {\widehat{\Delta }}_{s,V}\}$$. Each model $$ {\widehat{\Delta }}_{s,v}$$ uses $$ X\backslash X_{h_v}$$ to train on and generates predictions for $$ X_{h_v}$$.

Concatenate the V sets of predictions from $$ {\widehat{\Delta }}_s$$ into a vector z of length N. The N × S matrix Z of such vectors for all S base learners becomes the feature matrix for level-one training.

Train a meta-learner Ψ on Z with hyperparameter tuning if necessary, again using Y as the labels.

The base learners used in the final stacked model are created using all training samples, so one final (non-cross-validated) training round is needed. To this end, train each Δ_s base learner on the full level-zero training set X.

Form the stacked model using the base-learners trained on the full data set and the meta-learner to obtain ({Δ_s: s = 1, …, S}, Ψ). To predict on a new data point X_new, compute a 1 × S vector Z_new = {Δ_s(X_new): s = 1, …, S} and use Ψ(Z_new) as the stacked model’s prediction.

The predictive performance of the meta-learner is generally improved when the base learner predictions are maximally uncorrelated. This is achieved either by using different algorithms, or by varying the parameters of a particular modeling approach [58, 74]. AKLIMATE generates diversity through the use of feature subsets built around distinct biological concepts and processes.

Super learner

Stacked learning imposes no conditions on the choice of meta-learner Ψ. The disadvantage of such flexibility is the lack of theoretical results for the improved empirical performance of stacking. A Super Learner [58] is a type of stacked learner with restrictions on Ψ that give provable desirable properties. The main such constraint is that the optimal meta-learner Ψ* is the minimizer of a bounded loss function. A Super Learner for which {Δ₁, …, Δ_S} and Ψ have uniformly bounded loss functions exhibits the “oracle” property—Ψ* is asymptotically guaranteed to perform as well as the optimal base learner $$ {\Delta }_s^{*}$$ under the true data-generating distribution [76, 77]. Furthermore, if we constrain the choice of Ψ to $$ \Psi ={\sum }_{i=1}^S{\alpha }_i{\Delta }_i$$,∀α_i ≥ 0, Ψ* asymptotically converges to the performance of the optimal conical combination of {Δ₁, …, Δ_S} [58, 59]. This is the main reason why Ψ often takes the form of regularized linear or logistic regression.

For example, if the aim is to predict a continuous variable, one can set $$ \Psi ={\sum }_{i=1}^S{\alpha }_i{\widehat{\Delta }}_i$$ and solve the regression problem:

which can be regularized or subjected to a convex sum constraint on the α weights (e.g. ∀α_i ≥ 0, $$ {\sum }_{i=1}^S{\alpha }_i=1$$) [58, 59]. Similarly, the squared error loss of Eq (11) can be replaced with the logistic loss to solve a classification problem:

maintaining all theoretical Super Learner results.

AKLIMATE

To our knowledge, AKLIMATE is the first instance of a kernel-based stacked learner. The base learners {Δ₁, …, Δ_S} are RFs, each of which is used to produce an associated kernel. The meta-learner Ψ is an elastic-net regularized MKL that can be interpreted as the kernel learning counterpart of linear regression. Before describing the algorithm, we first discuss how the level-one RF kernels are constructed using OOB samples.

AKLIMATE selection of relevant RFs

AKLIMATE’s selection step for the best RFs Δ* in Algorithm 1 can filter the full collection of feature sets down to a subgroup of relevant sets two or more orders of magnitude smaller in size. This makes it possible to incorporate appreciably more feature sets than standard MKL algorithms which require kernel evaluation for all feature sets. When the number of such sets is in the thousands, MKL can be computationally very slow, even for data sets of small sample size. However, usually only a small proportion of the feature set collection is truly explanatory for a given prediction task. Thus, filtering out the non-relevant parts of the compendium does not impact MKL accuracy yet drastically improves computation time.

Algorithm 1 demonstrates the Δ* discovery process for a classification task. However, if Y is continuous (i.e. a regression problem), the predictions and labels are not directly comparable for equality. One approach is to take the squared error of the prediction-label differences and use that metric to re-rank RFs for each sample. In our experience, this leads to the selection of suboptimal RFs due to overfitting. Instead, AKLIMATE uses a more robust scheme that shows better results in practice—the vector of predictions for each sample across all RFs $$ {\left\{{\widehat{Y}}_n^s\right\}}_{s=1}^S$$ is binarized into matching and non-matching predictions and then the standard classification case selection rule is applied. The binarization is done as follows:

where q is a user-specified quantile of the empirical distribution of absolute prediction errors $$ {\left\{|{\widehat{Y}}_n^s-Y_n|\right\}}_{s=1}^S$$ (default q = 0.05). This setup prioritizes RFs that perform near-optimally on individual data points and optimally when the training set is considered as a whole.