Identifying interactions in omics data for clinical biomarker discovery using symbolic regression

Abstract Motivation The identification of predictive biomarker signatures from omics and multi-omics data for clinical applications is an active area of research. Recent developments in assay technologies and machine learning (ML) methods have led to significant improvements in predictive performance. However, most high-performing ML methods suffer from complex architectures and lack interpretability. Results We present the application of a novel symbolic-regression-based algorithm, the QLattice, on a selection of clinical omics datasets. This approach generates parsimonious high-performing models that can both predict disease outcomes and reveal putative disease mechanisms, demonstrating the importance of selecting maximally relevant and minimally redundant features in omics-based machine-learning applications. The simplicity and high-predictive power of these biomarker signatures make them attractive tools for high-stakes applications in areas such as primary care, clinical decision-making and patient stratification. Availability and implementation The QLattice is available as part of a python package (feyn), which is available at the Python Package Index (https://pypi.org/project/feyn/) and can be installed via pip. The documentation provides guides, tutorials and the API reference (https://docs.abzu.ai/). All code and data used to generate the models and plots discussed in this work can be found in https://github.com/abzu-ai/QLattice-clinical-omics. Supplementary information Supplementary material is available at Bioinformatics online.


Background
The rapid increase in biological data obtained through highthroughput technologies offers new opportunities to unravel the networks of molecular interactions that underlie health and disease (Perkel, 2021). An important contribution to this is made by genomics, transcriptomics, proteomics, lipidomics and metabolomics studies, which generate thousands of measurements per sample and offer the unique opportunity to uncover molecular signatures associated with a particular condition or phenotype. These signatures have the potential to act as biomarkers, i.e. a biological characteristic used in the evaluation of normal, abnormal or pathogenic conditions. Biomarker profiles have been found to be particularly useful for medical decision making, where use cases such as surrogate endpoints, exposure, diagnosis and disease management have been identified (Ghosh and Poisson, 2009). Although the large amount of omics data contains extensive information, it is not always trivial to extract actionable insights from it. Challenges include the high dimensionality of datasets where the number of variables far exceeds the number of samples, unbalanced measured outcomes (target variables), heterogeneous molecular profiles with multiple subtypes of patients and diseases, and instrumental and experimental biases (Libbrecht and Noble, 2015;Podgó rski, 2021;Whalen et al., 2021).
Classical statistical modelling has long been the gold standard for the analysis of genomics and transcriptomics data analysis. As a result, a significant amount of post-processing is required to condense information into meaningful results, e.g. through manual searching, enrichment or pathway analysis. An inherent challenge in the wide data matrices typical of omics is the existence of dependencies between features. This phenomenon is called 'multicollinearity' or 'concurvity' when linear and non-linear dependencies are involved, respectively (Buja et al., 1989). The increasing availability of affordable computing power and high-throughput omics data has led to the increasing use of machine learning (ML) in the life sciences and pharmaceutical industries. In addition, ML methods have been used for biomarker discovery based on omics data, where they are beginning to outperform state-of-the-art assays (Mann et al., 2021). Feature selection methods such as minimum redundancy maximum relevance have also shown great utility in identifying parsimonious sets of features that act as simple, predictive and robust signatures (Radovic, 2017).
Due to the inherent noise of biological data and the 'curse of dimensionality' (Altman and Krzywinski, 2018;Hastie et al., 2001) (more features than observations), it is a non-trivial task to perform traditional ML without misleading or overfitting the model during training, such that it is unable to robustly predict outcomes on unseen samples (Domingos, 2012). In addition, most state-of-the-art ML models are difficult to interpret and are therefore often considered complex 'black boxes' (Lundberg and Lee, 2017). Applying black-box ML models such as random forests and neural networks to omics data has proven effective in identifying predictive biomarkers (Chen et al., 2020), but the underlying relationships between features remain hidden, and especially for decisions where the stakes are high, it has been argued that interpretable methods should be used wherever possible (Rudin, 2019).

Symbolic regression and parsimonious models
Recently, the QLattice, a new ML method based on symbolic regression (SR), has shown promising results in terms of performance and interpretability (https://doi.org/10.48550/arXiv.2103.15147). The goal of any implementation of SR is to model a relationship between one or more independent variables X and a dependent variable y by finding a suitable combination of mathematical operators and parameters. Even when considering only expressions with finite length, the search space is usually too large for any kind of brute force strategy, and thus, alternative methods are required. All SR algorithms can be thought of as methods of searching this combinatorial space effectively. SR is an active field of research and there are multiple examples of recent implementations (Burlacu et al., 2020;Udrescu et al., 2020) (https://gplearn.readthedocs.io/en/stable/). SR is particularly suitable for scenarios where the number of features in the model should be small and their interpretation and interactions are of primary interest. Furthermore, it seeks to solve problems where the mathematical form of the data generating process cannot be assumed, or approximated, a priori. This is in contrast to the typical regression problem where parameters are fitted to a presupposed model, like linear models or polynomials. Thanks to its unconstrained nature, SR can usually attain higher performances while keeping the number of explicit parameters as low as possible.
It is well known that most functions can be approximated by using an arbitrarily large number of coefficients and functions belonging to a complete set (e.g. Fourier series, Chebyshev polynomials, etc.). Analogously, one can theoretically build a model that explains y in terms of X with arbitrarily low train error, even if the approximated mathematical model is ostensibly different from the data-generating process. This does not necessarily pose a problem to types of research where the primary objective is to produce a working model that fits well the data, but vital information may be lost along with interpretability as model complexity grows. The most well-known example of this is deep neural networks, where e.g. modern language models contain billions of parameters, inevitably trading off interpretability for performance.
In contrast, the aspiration of SR is that domain knowledge can be applied and extracted more efficiently by seeking simpler mathematical models to preserve explainability from a human perspective. In principle, this increases the likelihood of discovering driving mechanisms in data, and inclines SR methods towards maximum information gathering, which is vital in (e.g. life) sciences where both performance and interpretability is important. In practice, SR methods achieve this by using parsimonious models that explain the data with a minimal number of parameters. Additionally, one can use complexity measures such as Bayesian information criterion (BIC) and Akaike information criterion (AIC) to ensure that the resulting models generalise well from train to test set.
Here, we applied the QLattice to four different omics problems to identify biomarker signatures that predict clinical outcomes while also revealing new interactions in the data. We demonstrate how highly complex problems can quickly be condensed into a set of simple models that can be reasoned and used as hypotheses for potential mechanisms underlying the problem at hand.

The QLattice
The QLattice is a SR engine that aims to solve the optimisation problem of finding the functions that best fit the data. It applies an evolutionary algorithm framework to find the combinations of inputs, operators and parameters that minimise the fitting error in a supervised learning problem (see Koza (1992) for seminal work on genetic programming, and a practical guide in Poli et al. (2008)).
The QLattice algorithm works as follows: first it generates an initial sample of functions, fits them with gradient descent, and evaluates them for fitness. This initial sample is formed using a set of estimated priors assigned to each input based on its mutual information with the output. Then, the best performing functions are used to create a new generation of functions consisting of three groups: (i) the best performers from the previous generation, (ii) mutated versions of the best performers from the previous generation and (iii) a completely new set of sampled functions. Instead of sampling mutations and new functions from a uniform distribution of inputs and operators, the QLattice draws from a probability distribution that is learned thanks to a mapping between the functional space and a lattice. Thus, with each generation the QLattice improves the probability distribution estimation. As this iterative process continues, the QLattice expands the search for the best fitting functions. The result of a training run is a list of functions sorted by a userdefined quality metric. These functions serve as hypotheses that each serve as their own solution to the problem. A more extensive description of the methodology can be found in Larsen (2021).
The QLattice can be used in both regression and classification tasks for supervised learning problems. In the case of classification, the algorithm is designed to work with binary problems, although it can be easily extended to multi-class targets using a one-versus-rest approach (see Bishop (2006) for detailed description of the method). All the QLattice models discussed in this manuscript are trained to perform binary classification tasks. The target variables are encoded as 0 or 1, and the output of the models is to be interpreted as a probability. In order to keep the outputs between 0 and 1, all the mathematical expressions are wrapped with the logistic regression function 1=ð1 þ exp ðÀf ðXÞÞÞ, expressed throughout the text as logreg.
The Feyn Python library (https://doi.org/10.48550/ arXiv.2104.05417) is the interface between the user and the QLattice, and it is used to train and analyse new models. Its highlevel train function returns a list of ten models sorted by a criterion of choice (see documentation at https://docs.abzu.ai). The default sorting option is the BIC, which amounts to the training loss plus a complexity penalty, and allows selection of the most generalisable models without compromising training speed (Hastie et al., 2001).

Cross-validation
Overfitting and spurious correlations are major concerns when applying ML to the wide datasets typical of many areas of computational biology such as genetics, transcriptomics and proteomics (i.e. when the number of features is much larger than the number of observations). For these kinds of datasets, simple models with complexity penalties tend to offer competitive performances (Hastie et al., 2001). This is the case of the models selected by the QLattice when the BIC criterion is enabled.
The BIC criterion used for model selection, however, does not provide an unbiased estimate of the test performance. Therefore, we use a standard k-fold cross-validation scheme to estimate the performance of the QLattice and determine what one can expect from the models selected by it. We use a scheme with five folds: four folds as a train set and one as a test set. In each of the five training loops, we reset the QLattice and call the train function to avoid "data leakage" in the feature selection. Individual models' performances are estimated using single train/test splits.
Finally, in the Supplementary material, we present a comprehensive benchmark of the QLattice along with other ML algorithms in combination with different feature selection techniques on all four datasets discussed in the manuscript.

Selection of models for further analysis
In ML, the emphasis is usually put on test-set performance. In most instances, model selection is done with the sole goal of finding the models that will generalise best on new data. BIC is a good example of such a model selection tool, and the QLattice uses it to explore and find the models with strongest signal-both in the train and test sets. In the cases where the user is only interested in prediction performance, one should select the model with the lowest BIC score. This is the selection criterion we followed in the benchmarking section.
However, interpretable algorithms like the QLattice have more goals than performance. They are also used to generate hypotheses about the features involved in a process, and their specific relations. When using BIC as a criterion, all the models returned by the QLattice can be expected to highlight robust patterns in the data. Although the performance of the models might differ, one should consider all of them valid. The list of models returned by the QLattice might highlight a combination of patterns: different mechanisms, the same one represented by multicollinear features, already known biomarkers, or completely new candidates. The evaluation of the user (human in the loop) is then necessary to extract the relevant learnings from the models and put them in the context of the question at hand.
For the sake of clarity, we only discuss one model in each of the first three cases. We selected the models according to performance and interpretability: when models had very similar performances, we chose the simplest models first. In the insulin response case we used evidence from previous studies to choose a model, where the gene features could be easily interpreted.
Note, that the biomarker candidates selected by the models on the different cases analysed in this article should be further investigated-both in relation to the disease mechanism and with regards to confounding factors (e.g. cohort dependencies).

Data preparation
Proteomics: Alzheimer's disease. The data were taken from Bader et al. (2020) and consist of 1166 protein expression of the cerebrospinal fluid of 137 subjects, collected in three sample groups (we address the possible confounding factors in Section 3). We used the QLattice to predict whether a patient would develop Alzheimer (dependent variable ¼ 1) or not (dependent variable ¼ 0).
Gene expression: relevant genes for insulin response in obese and never-obese women. The data were retrieved from Mileti et al. (2021). The dataset consists of gene expression from a total of 23 never obese and 23 obese women sequenced before and 2 years after bariatric surgery (post-obese) using RNA sequencing (CAGE) (Mileti et al., 2021). The only pre-processing done was to normalise the data from raw counts to (tag-per-million normalisation, the gold standard for CAGE data; Mileti et al., 2021). The data are balanced for MValue, body mass index and age across sample groups. We used to QLattice to model the response to insulin based on geneexpression measurements and predicted whether an individual is in a fasting (target variable ¼ 0) or hyperinsulinemic (target variable ¼ 1) state.
Epigenomics: hepatocellular carcinoma. The data were processed to contain only the 1712 most important features, filtered for variance. The curated dataset contained 1712 CpG island (CGI) features with a binary target of 55 cancer-free (target variable ¼ 0) and 36 cancer (target variable ¼ 1) individuals coming from a single sample group (plasma samples) Wen (2015). The CGI features cover the methylated alleles per million mapped reads.
Multi-omics: breast cancer. The dataset was obtained from Ciriello et al. (2015) and contains multi-omics data from 705 breast tumour samples of different patients. We use the QLattice to predict outcomes; a survival outcome is encoded with target variable ¼ 0 (611 patients) and a fatal outcome with target variable ¼ 1 (94 patients). The data was extracted from The Cancer Genome Atlas through the R-package 'curatedTCGAData' (Ramos et al., 2020) and included four data types: somatic mutations, copy number variations, gene expressions and protein expressions. The raw data were pre-processed with a variance threshold limiting each type of input to the highest variance features. The data were stratified for lobular and ductal subtypes in each train/ test split and the models were assessed for potential confounding influences from factors such as age, stage and treatment regimen.

Results
In the following cases, we showcase different aspects of the QLattice using four different omics data types: • Interpretability: In the proteomics case, we show how the QLattice finds high-performing models that can be easily interpreted. • Feature combinations: In the gene-expression case, we demonstrate how the QLattice finds biomarker signatures that together explain the data better than any single feature on its own. • Multicollinearity: In the epigenomics case, we show how the QLattice deals with multicollinearity typical of omics data by choosing the combination of features that best explains the target while minimising complexity of the model. • Multi-omics and non-linear interactions: In the multiomics case, we highlight how the QLattice can find non-linear interactions within and across omics data types that help to stratify patient populations.

Proteomics: Alzheimer's disease
Background Despite many decades of research, neurodegenerative diseases remain a major threat to human health and are a substantial cause of mortality. Alzheimer's disease (AD) is the most common type of dementia, and currently no therapeutics can halt or significantly slow its fatal progression (van der Schaar et al., 2021). Furthermore, short of an autopsy, there is no definitive way to diagnose AD, and it is in general impossible to predict who will develop the disease.
Here, we demonstrate how the QLattice can be used to discover protein biomarkers for AD working with the data from Bader et al. (2020). We will use this example as an introduction to the QLattice functionality and capabilities.

Model analysis
After splitting the dataset into 80% train and 20% test partitions, we ran the QLattice on the train partition to obtain 10 best unique models from the QLattice (Table 1). Each model points to a relation that serves as a data-derived hypothesis. Thus, all 10 models potentially hold insights into the mechanisms involved in AD.
We chose the model with the lowest BIC-score for thorough analysis. This model uses MAPT, age at CSF collection and LILRA2 as inputs combined with additions to predict the probability of AD for a given patient. As a ML model, it can be analysed in terms of test prediction metrics (Fig. 1) to verify that the relations found are not spurious (see Walsh (2021) for a review on the matter).
We ran the cross validation scheme outlined in Section 2. The estimated test performance of the QLattice top models was AUC ¼ 0.94 (mean of the five folds, with a standard deviation of 0.05). We note that the predictive power might be overestimated due to the presence of confounders in the data.

Model interpretation
The known AD biomarker MAPT (tau protein) was consistently found in the highest scoring QLattice models, while the additional features varied between models. Figure 2 shows how MAPT contributes prominently to the signal of the chosen model. The plot shows the signal flow in the model, and the colour represents the strength of the association of each node to the clinical outcome. The association measure used is mutual information (Cover and Thomas, 2006). Thus, the features age at CSF collection and LILRA2 are both secondary to MAPT but both improve the model as made clear by the rising mutual information numbers displayed on top of the nodes. MAPT on its own has a mutual information score of 0.37, but this number rises to 0.56 when applying the additional features and the right mathematical operators-in this case additions.
The partial dependence plot (Fig. 3) shows that at fixed LILRA2, higher levels of MAPT leads to positive AD prediction. When the MAPT level reaches around 25 000 the model starts to predict AD-positive. In addition, the effect of age is displayed in the plot. Unsurprisingly, at a higher age comparably lower MAPT levels trigger the model to predict AD-positive (when the predicted probability rises above 0.5), as displayed by the coloured curves.
It should be mentioned that the protein levels of LILRA2 only present a significant difference between AD and non-AD patients for one of the cohorts depicted in Bader et al. (2020) (P-value of 0.008 with a Student's t-test). Since we fix LILRA2 in Figure 3, the model response to MAPT is cohort-independent. Moreover, the three age values depicted in Figure 3 are all between 62 and 72 years old. This ensures that the observed model response is free from any bias arising from the younger control group present in one of the cohorts (Bader et al., 2020). Including suspected confounding factors in the model is a standard practice to statistically control for confounders when using linear and logistic regression (see Angrist and Pischke (2008) for an extensive review on the topic).
The QLattice models provide data-derived hypotheses that can quickly provide an overview of possible explanations to a given question. For instance, the first model in Table 1 may be translated into the following hypothesis: 'MAPT is a main driver of AD since it is positively correlated with AD status' and 'MAPT interacts both with the age of the patient and with the protein LILRA2'. Thus, the mathematical simplicity of the QLattice models allows direct translation into hypotheses that can be readily understood and tested. This marks a significant departure from black-box ML models, where the inner working of the models is usually more opaque.
3.2 Gene expression: relevant genes for insulin response in obese and never obese women Background Obesity is a major public health problem, and obese people are at higher risk of heart disease, stroke and type 2 diabetes. Obesity is considered a medical condition caused by eating more calories than  2. Model signal path for AD. A prominent signal contribution from MAPT was found in all 10 models (green). The signal is expressed in terms of mutual information and displayed above the nodes in the model (see Cover and Thomas (2006)) (A color version of this figure appears in the online version of this article.) Fig. 3. Partial dependence plot for the AD model: marginal effect of MAPT on AD-risk necessary, but it can also be caused by a decreased response to insulin.
To shed light on this, a recent publication (Mileti et al., 2021) focused on white adipose tissue, which is one of the main insulinresponsive tissues. In this study, obese subjects underwent gastric bypass surgery and lost weight. Weight loss can support the subsequent restoration of the insulin response. In Mileti et al. (2021), insulin sensitivity was determined using the hyperinsulinaemic euglycemic clamp, while the insulin response was measured using cap analysis of gene expression (CAGE) from 23 obese women before and 2 years after bariatric surgery. To control for the effects of surgery, 23 never obese women were also included.
The experiment was designed to understand the effects of insulin on the expression of different genes. In traditional differential gene expression (DGE) analysis, the individual genes with the strongest and most consistent changes between conditions are highlighted. Here, we propose a complementary approach to DGE analysis that uses the QLattice to identify sets of genes and their interactions that best separate two groups of samples.
Specifically, we modelled the response to insulin based on geneexpression measurements and predicted whether an individual is in a fasting or hyperinsulinemic state. As well as being a predictive algorithm, the QLattice looks for different interactions between genes that describe the insulin response in two classes of individuals.

Model analysis
We inspect the 10 models returned by the QLattice in Table 2 after we ran it on the training set (80-20% split). We select the second model for further analysis because it contains PDK4, an established insulin target (Mileti et al., 2021); C2CD2L, a positive regulator of insulin secretion during glucose stimulus; and PHF23 a negative regulator of autophagy. To our best knowledge, defects in autophagy homeostasis are also implicated in metabolic disorders such as obesity and insulin resistance as discussed in Zhang et al. (2018). The high performance of this model is summarised in Figure 4 for both the training and test sets. In addition, the QLattice identified other genes known to be insulin targets or found in the article such as C19orf80 and LDLR (Mileti et al., 2021).
Qlattice as complementary approach to differential gene expression Differential gene expression analysis (DGE) is generally used to detect quantitative changes in expression levels between experimental groups based on normalized read-count data. There are several methods for differential expression analysis based on negative binomial distributions (Love et al., 2014;Robinson et al., 2010) or based on a negative binomial model (Bayesian approaches) (Hardcastle, 2021;Leng and Kendziorski, 2021;Smyth, 2004). Differential expression tools can perform pairwise comparisons or multiple comparisons.
Alternatively, DGE can be used to identify candidate biomarkers, as it provides a robust method for selecting genes that offer the greatest biological insight into the processes influenced by the condition(s) under investigation. However, this robustness can sometimes translate into rigidity. Signatures expressed in linear combinations, interactions or through non-linear relationships may be overlooked when using DGE. SR-based ML models offer a complementary view on the data and highlight predictive signatures. The advantage of this approach is that even simple feature combinations can lead to a high predictive performance. As we can see in the model decision boundaries of Figure 5, a linear combination of the features PDK4, PHF23 and C2CD2L can characterise the insulin response for almost all individuals in the sample. The strength of the signal is found as well in the test set (see Fig. 4).
Although PDK4 and PHF23 are reported as significant in the DGE analysis (according to FDR), they do not appear at the top of the list ordered by log-fold change (the one used in Mileti et al. (2021)). This apparent discrepancy between the DGE and the QLattice choice can be explained by the fact that the DGE only considers the univariate distributions. From the density plots in Figure 6, we can indeed see a considerable overlap between the two classes when we look at the univariate distributions of PDK4 and PHF23, which is smaller for the distributions of the linear combination of genes. The effect can also be seen in the mutual information between the variables or their combinations, and the target variable.
In summary, we find that the QLattice can be used as a complementary method to DGE, as it is very good at finding feature combinations that carry strong signals, and as it efficiently explores the feature space without requiring an exhaustive exploration of all features. There have been efforts in this direction using mutual information and partial information decomposition (Chan et al., 2017). Consequently, the QLattice can suggest specific operations for the proposed combinations and help to better understand biologically relevant interactions that were previously hidden.
Epigenomics: hepatocellular carcinoma. Background. Primary liver cancer is a major health burden worldwide and develops in response to chronic inflammation of the liver. This can be caused by a number of insults such as viral infections as well as both alcoholic steatohepatitis and non-alcoholic steatohepatitis. The most common form of liver cancer is hepatocellular carcinoma (HCC), which accounts for 90% of liver cancers and is the third leading cause of cancer mortality worldwide (Llovet et al., 2021;Yang and Roberts, 2010).
In this HCC diagnosis example, we explore how the QLattice performs on highly multicollinear data. The dataset was taken from a study by Wen et al. (2015) and contains data generated by methylated CpG tandems amplification and sequencing, a method that can detect thousands of hypermethylated CpG islands (CGIs) simultaneously in circulating cell-free DNA (ccfDNA). The aim is to explain liver cancer occurrence using methylation biomarkers as features. After pre-processing (see Section 2) the curated dataset contained nearly 2000 features. As demonstrated below, the QLattice gave highly predictive models using only a few key interactions.
Model analysis. As in the previous case, we split the dataset into train and test partitions (80-20%) and ran the QLattice with default settings on the training set. We inspected the ten models returned in Table 3 balancing simplicity and performance. The 10 models all perform equally well, and we therefore chose the model with the least features for further examination (n. 5). The models is shown in Figure 7 and its performance metrics are summarised in Figure 8.
Using 5-fold cross-validation, the QLattice top models yielded an average performance of AUC ¼ 0.93 (mean of the five folds, with a standard deviation of 0.01).
Model interpretation. As can be seen in Figure 9, the primary separator of the two features in the selected model is chr17_59473060_59483266. Individuals who do not have cancer have stable, low levels of methylated alleles, while individuals with cancer generally have higher, more variable levels of this trait. In addition, we find that some cancer individuals have low levels of chr17_59473060_59483266. Furthermore, from the 2d plot of partial dependence in Figure 9 we can also see that low values of both chr17_59473060_59483266 and chr3_9987895_9989619 can be used to identify cancer individuals. This dynamic is captured well in the 2d partial dependence plot of Figure 9. This is an easily understood model, two genes interacting, generating a top performing model. The model's AUC and confusion matrices on the train and test set are shown in Figure 10.
The models that were generated (Table 3) perform equally well. Aside from chr17_59473060_59483266 all models contain different secondary features and thus there could be molecular substitutes among the other features. To show whether there is multicollinearity, an overview of other correlated features is given in the correlation heatmap, Supplementary Figure S1. The figure shows whether the selected model feature belong to a group of highly correlated features. If this is the case, we can most likely replace this one feature with another from the same group and achieve similar model performance. In this case, the QLattice achieves high performance by selecting one feature from each main variance group in the dataset.
Instead of using dimensionality reduction such as PCA to group features with similar variance into single features, the QLattice selects representatives from each variance group. The representative that performs best in combination with the other features in the training dataset is selected.

Multi-omics: breast cancer
Background Breast cancer is the most common cancer in women, worldwide. There are two main types of breast cancer, ductal and lobular The signal is expressed in terms of mutual information and displayed above the nodes in the model (Cover and Thomas, 2006) (Katz and Morrow, 2013;van Seijen et al., 2019). In this case study, we explore a multi-omics dataset and identify potential regulatory interactions across omics-types [copy numbers (cn), somatic mutations (mu), gene expression (rs) and protein expression (pp)] that could explain and predict survival outcomes of breast-cancer patients. We benchmark the QLattice models with a random forest and show that in addition to revealing interactions the QLattice performs as well as complex 'black-box' models. The dataset was obtained from Ciriello et al. (2015) and contains multiomics data from 705 breast tumour samples. Two-feature model analysis. Upon running the QLattice on different partitions of the data, one can expect different models being selected. These models bring similar albeit complementary insights, as they are able to see different sub-samples of the data. In this case, we obtained diverse models by keeping the lowest BIC-scoring ones from each partition of our cross-validation scheme.
To maximise interpretability, we started by exploring simple models that allowed for a maximum of two features. The mean test AUC for the best models of all folds was 0.635, with a standard deviation of 0.070. Equations (1) Two of the expressions correspond to a bivariate normal distribution, while the others have a multiplicative interaction as shown in equations (1). All the chosen features in the models above are measurements of gene expression.
From the Pearson correlation heatmap in Figure 11, we observe that all five models contain a gene-expression feature from the group with highest pairwise Pearson correlation: rs_PCK1, rs_MRAP, rs_LOC283392, rs_APOB and rs_ACVR1C; their correlation values range from 0.774 to 0.835. Then these features are each combined in a non-linear interaction with the remaining gene expression variables.
Pairwise correlation gives a measure of the similarity between the input features. In addition, one can calculate the correlation between input features and the output variable, as shown in Table 4. The latter gives a measure of the relevancy of the input features relative to the output. Note in Table 4 that rs_PCK1, rs_MRAP, rs_LOC283392, rs_APOB and rs_ACVR1C are the features with highest relevance in this group. Therefore, akin to the HCC case, the models yielded by the QLattice combine a gene expression variable with high relevance with another gene expression with low similarity score.
The models' decisions boundaries are depicted in Figure 12. The bivariate Gaussian function and the product between the gene  expression features identify a 'hotspot', i.e. there is a particular range for these gene expressions that indicate whether a breastcancer patient died or survived. Strikingly, these patients were predominantly suffering from ductal breast cancers. There seems to be a putative molecular interaction that is an important biomarker for ductal breast-cancer survival. The model's AUC and confusion matrices on the train and test sets are shown in Figure 13.
The genes in the models depicted in equations (1) were found to have no relation to the age of the patients. The Pearson correlation coefficient between gene expression and age was computed and the highest absolute Pearson correlation coefficient value between gene expression and age was 0.153 (P-value of 4 Â 10 À5 ).
Comparison with multi-omics models Allowing models with higher complexity-more features and operations-can potentially unlock better performing models that mix different omics. To this end, we ran the same 5-fold cross-validation scheme allowing a maximum of five features. This should allow for any signal beyond the gene expression 'hotspot' to be captured by the QLattice. The resulting average test AUC score on the best models is 0.671 with standard deviation of 0.040. This average result is certainly larger than test AUC of the two feature models, although both scores could be considered statistically compatible. The average train AUC of the models on Eqs. 2 is 0.743 (r ¼ 0:020), which is significantly higher than the average train AUC of the two-features models, 0.683 (r ¼ 0:043). Since their mean test AUC scores are on par, the discrepancy in the training set implies that the more complex models depicted above tend to overfit when compared to the simpler gene expression models presented before. When it comes to the functional form of the models on Eqs. 2, it is interesting to note that they all possess a non-linear interaction between gene expression features (prefix rs). For most, this interaction is multiplicative, while for the model from Fold 3, the nonlinear boundary is set by the tanh function of rs_APOB.  Note: Values computed using SciPy (Virtanen et al., 2020). Comparison with random forest In order to get a better sense of the performance of the QLattice, we compare it with a widely used 'black-box' model: the random forest. We use the implementation by scikit-learn, and tune its hyperparameters and estimate its performance using a 'nested' cross-validation scheme (Cawley and Talbot, 2010). The best parameters lay around n estimators ¼ 50 and max depth ¼ 4 for the different folds, and the average performance is an AUC of 0.604 with standard deviation of 0.106, on par with the QLattice. This is a very remarkable result, considering that the QLattice is only using two features while the random forest uses potentially all of them. A benchmark with other algorithms and feature selection techniques can be found in the Supplementary material. Taking into consideration how the multi-omics models in Eqs 2 tend to overfit and the random forest result in comparison to the QLattice, we can conclude that the models in Eq. (1) reveal the core patterns in the data. In summary, the interaction of two gene expression variables allows for the identification of a 'hotspot' where the probability of a poor outcome of the disease is high. One of the genes in the model belongs to a group of genes with pairwise Pearson correlation above 0.7, while the other is taken from the remaining pool of variables. A possible next step in the study of this data is to pinpoint the combinations of gene expression variables that best predict vital.status.

Conclusion
Given the large amounts of data being generated and a need for more efficient treatment regimens, predictive analytics in the clinic is gaining traction. A range of methods exist that can predict a certain outcome based on omics data; however, there is a scarcity of interpretable alternatives. Here, we showed that we can identify simple yet highly predictive and explainable biomarker signatures by combining sophisticated feature selection with a powerful model search algorithm. Due to the small number of features, the models are robust and can be readily interpreted. This makes them a valuable starting point for researchers and clinicians who are looking to find new and biomarker signatures while learning about the underlying interactions in the data.