Learning epistatic polygenic phenotypes with Boolean interactions

Ethics statement

This research was covered by the UK Biobank ethics agreement. UK Biobank has approval from the North West Multi-centre Research Ethics Committee (MREC). Complete details of all UK Biobank procedures are available at www.ukbiobank.ac.uk. All participants provided written informed consent.

Phenotype and genotype data from the UK Biobank

As positive control case study for epiTree, we analyzed epistatic interactions in UK Biobank data associated with a previously-studied, genetically-driven phenotype: red hair [15]. We considered a total of 337,535 unrelated white British individuals and their self-reported hair color (Data field 1747; 15,326 with Red hair: positive cases and 322,209 individuals with Blonde, Light Brown, Dark Brown, Black, or Other hair color: controls). We searched for interactions in a balanced, random sample of cases and controls (15,000 red hair individuals and 15,000 controls). To assess the generalizability of our results to unobserved data, we performed a random sample split into training and test sets with 26K training samples and 4K test samples. We note that, in general, one might prefer non-randomly to randomly sampled test sets—e.g., a data from an external experiment. Such an external data set was non available to us in this setting.

We modeled hair color using ∼ 15, 000, 000 common variants imputed to the Haplotype Reference Consortium (HRC) and UK10K reference panels from ∼800, 000 directly genotyped variants, which were obtained by UK Biobank using one of two similar arrays [28]. Details regarding the ascertainment and quality control of these genotypes have been previously described [15, 28]. In brief, genotyped variants were subject to outlier-based filtration on effects due to batch, plate, sex, array, as well as discordance across control replicates. Samples with excess heterozygosity or missingness were excluded from the data release. Imputed variants were further subject to filtration based on Hardy-Weinburg equilibrium, missingness in white British individuals, minor allele frequency (> 10⁻⁴), and imputation quality. Respective details on phenotype and genotype data for the MS case study can be found in Section S1.8. in S1 File.

Two-step interaction screening with biologically inspired dimension reduction

We developed a two-step procedure for screening interactions to improve computational efficiency and help stabilize of results. First, we aggregated SNP-level information by gene to search for gene-level interactions associated with the target phenotype. This initial step (i) substantially reduced dimensionality of the data and (ii) detected interactions at the level of a biologically interpretable unit. Second, variants surrounding putatively interacting genes were selected to search for SNP-level interactions.

iRF non-linear model selection

Given a set of features (either estimated gene expression or individual genetic variants), we applied iRF [31, 32] as a non-linear model selection step to extract candidate interactions from the training data. The iRF algorithm fits a series of feature-weighted random forests (RF) [33] to iteratively stabilize feature selection along the decision paths of trees in the forest. After the final iteration, iRF identifies stable candidate interactions by searching for features that frequently co-occur on decision tree paths using random intersection trees (RIT) [34]. In addition to the interactions explicitly identified by iRF, we further expanded our set of candidate interactions for PCS epistasis inference (see below) by taking all inter-chromosome pairwise interactions among the top 50 iRF genes (with respect to Gini importance).

A key benefit of iRF is that the computational complexity of searching for interactions does not grow exponentially with the size of an interaction. This allowed us to identify high-order (i.e. beyond pairwise) Boolean interactions in a computationally tractable manner. Moreover, accuracy of iRF predictions on the hold-out test data provides a measure of generalizability for our non-linear model selection method, and correspondingly, the identified interactions. Finally, iRF detects known interactions in both simulation and biological case studies [31, 32, 35, 36] and analytically tractable versions of iRF yield provably consistent interaction discovery [37].

Predictability, Computability, Stability (PCS) inference for epistatic interactions based on Boolean models

The significance of a candidate interaction for a binary phenotype is traditionally evaluated through logistic regression (see Section S1.1.2 in S1 File for a recap). However, the model imposed by logistic regression—including a response rescaling, transformations of gene expression data, and the polynomial interaction term—does not arise from a known biological observation or mechanism. Indeed, additive or multiplicative forms may not capture the stochastic nature of living systems or threshold effects commonly observed in molecular and cellular phenomena [38–41]. A mathematically inaccurate representation of the biological mechanism of action between multiple genetic loci may result in unstable and irreproducible results that provide a poor fit of the data.

To address these issues, we developed a PCS epistasis p-value to evaluate. The PCS framework, recently proposed in [20], is based on the three core principles of data science: Predictability, Computability and Stability (PCS). It unifies and expands ideas from machine learning and statistics by using predictivity as a universal reality (or model) check, appropriate computational strategies for every step of an analysis process including simulating realistic reference/null distributions, and stability analysis to assess reproducibility at every stage of an analysis—from problem formulation/data collection, to modeling including inference, and to post-hoc model interpretations. PCS epistasis p-values are based on Boolean models of non-epistasis (null) and epistasis (alternative), learned adaptively from training data, and evaluated on hold-out test data. This approach has the advantage of being able to learn flexible interaction forms from the training data, rather than assuming multiplicative interactions that lack biological justification [8, 9, 11, 42].

Red hair

To assess how well genotype/phenotype relationships learned by iRF (e.g., predicting ‘red hair’ phenotype from the gene and variant features, respectively) generalized to new samples, we evaluated prediction accuracy on the hold-out test set. Fig 3 reports ROC curves for the prediction accuracy of iRF and competitors on the gene level (left) and variant level (right) (all fit using training data): penalized logistic regression with L1 penalty term and random forests (RF). With an area under the ROC curve (AUROC) of 0.93 (95% bootstrap confidence interval was [0.922, 0.936]) on hold-out test data, iRF demonstrates high prediction accuracy and outperforms its closest competitor (penalized logistic regression AUROC = 0.9, 95% bootstrap confidence interval of [0.894, 0.913]; ranger AUROC = 0.88, 95% bootstrap confidence interval of [0.871, 0.892]).

Before conducting inference on interactions recovered by iRF, we applied two filtering steps. First, we performed n.bootstrap = 50 bootstrap replicates of the iRF search and removed interactions appearing in fewer than 50% of these replicates. Second, we removed interactions among genes/variants that were in high linkage disequilibrium (LD) as it has been well-documented that LD can lead to spurious estimates of epistasis [19, 44, 45]. We note that many of the gene-level interactions we detected were among genes in the vicinity of MC1R (S23 Fig in S1 File; LD R² of at least 40% for all removed interactions); our LD filtering was equivalent to restricting results to interactions among genes/variants on different chromosomes. This resulted in 10 order-2, 6 order-3, and 2 order-4 interactions (gene-level analysis, Fig 5, see also Fig 4) as well as 25 order-2 and 3 order-3 interactions (SNP-level analysis, Fig 6), representing a massive reduction in the ∼ 10¹⁶ possible order-4 interactions among ∼ 10⁴ genes. See also S23 Fig in S1 File, which shows the number of discovered interactions per interaction-order within the individual screening steps. Interactions primarily centered around genes in the MC1R and ASIP regions, which corroborates previously reported epistatic interactions surrounding red hair, see [15]. We stress that while other works a priori restrict their search to interactions involving MC1R, our approach does not require such pre-selection.

We evaluated the significance of the individual candidate interactions using PCS epistasis inference (Figs 5 and 6; comparisons with logistic regression reported in S15, S16 Figs in S1 File). Interactions with significant PCS p-value (< 0.05 with Bonferroni correction) correspond to those for which the epistasis model results in a stable decrease in prediction error relative to the no epistasis model. For several of these interactions, CART tree-based models show a considerable improvement in prediction error relative to logistic regression (e.g. ASIP, TUBB3; ASIP, VPS9D1; S15, S16 Figs in S1 File). In contrast, prediction accuracy of the epistasis model shows little improvement over the no epistasis model for interactions that are significant with respect to the logistic regression p-value. This finding suggests that the thresholding, Boolean form of interactions used by CART is a more accurate model of interaction behavior between genes comapred with the polynomial interaction model used by logistic regression. We hypothesize that this is a result of thresholding behavior that has been previously observed in biomolecular interactions [38–41].

PCS detects different types of epistasis

PCS epistasis and logistic regression p-values can differ substantially. This reflects the fact that each captures a different form of interaction, and thus, different genotype/phenotype relationships. The PCS p-values, based on CART models, operate directly on the penetrance scale and capture thresholding-type behavior. Logistic regression considers epistasis on a logit scale and is restricted to polynomial interaction terms. A key advantage of having CART as the building block for PCS p-values in our epiTree test is that interaction behavior is interpretable and easily visualized via the respective decision trees (Fig 2 for an example). We highlight two examples from the red hair case study below.

First, we considered interactions where PCS epistasis p-values were significant and logistic regression p-values were not. For example, the interaction between ASIP and TUBB3 shows the greatest difference between PCS and classical p-values (see S15 Fig in S1 File and Fig 5). Fig 7 reports response surfaces, indicating penetrance (P(red hair)) as a function of interacting features (top/bottom rows: logistic regression/CART, left/right columns: non-epistasis/epistasis models). For comparison, we report the smoothed distribution of test data as ground truth (Fig 7e). Here, the polynomial interaction term assumed by logistic regression fails to capture the non-monotic relationship between TUBB3 and P(red hair). On the other hand, the CART model more accurately captures this relationship, which is confirmed by the improvement in prediction accuracy relative to logistic regression (S15 Fig in S1 File). In other words, the PCS p-value detects a non-additive, epistatic relation between the two genes that cannot be captured by the logistic regression model. Such non-monotonic behaviour, as it is observed in Fig 7, is common among PrediXcan estimated gene expression features. This is due to the fact that estimated expression of a single gene corresponds to a weighted linear combination of discrete variants. When a small number of variants are strongly associated with the response and have large PrediXcan weights for a given gene, the response surface exhibits marked transitions that correspond to different values of the highly weighted variants. This is exactly the case for TUBB3, with 25% of its weight mass on the variant rs8048449, which shows strong marginal association with red hair. A similar example for the interaction ASIP—DBNDD1 is shown in Fig 8.

Second, we considered interactions where PCS epistasis p-values were not significant and logistic regression p-values were. For example, the interaction between DEF8 and ASIP has highly significant logistic regression p-value (S1 File). The response surfaces indicate that these models differ most in the region where DEF8 is small and ASIP is large. While the logistic regression model estimates a high probability of red hair in this region, there is limited data to support how responses behave here. Indeed, all four models (logistic regression non-epistasis, logistic regression epistasis, CART non-epistasis, CART epistasis) achieve nearly identical prediction accuracy (with average cross entropies of 0.624, 0.620, 0.622, 0.630, see S15 Fig in S1 File and Fig 5). In other words, there is no evidence with respect to prediction accuracy that a particular model should be preferred, despite the significance of the classical logistic regression p-value. S12 Fig in S1 File shows another example for the interaction ZNF276-RPL36P4 where the PCS p-value is not significant but the logistic regression p-value is. Here, the data clearly shows a non-linear relationship in one of the variables, which cannot be captured by the linear components of the null model for logistic regression. CART can capture such a relation well and outperforms logistic regression with a simple additive relationship. That is, the logistic regression significant p-values appear to be driven by a main effect that is not captured by the linear form of main effect in logistic regression.

We provide a similar discussion for the remaining three pairwise interactions for which either the PCS p-value or the logistic regression p-value was significant (with Bonferroni correction) in Section S1.2. in S1 File.

Biological findings

PCS epistasis p-values detect diverse interaction forms in simulation studies

To evaluate the performance of epistasis inference methods, we conducted a series of data-based simulations comparing (i) epiTree PCS epistasis p-values (ii) logistic regression p-values (iii) permutation random forest (pRF) p-values [55]. We generated responses using four different models, each corresponding to a different form of interaction, applied to the gene expression data from the red hair case study. For each response generating model, we generated a set of candidate interactions using iRF and computed p-values for each inference method. For three of the four response generating models, epiTree both consistently recovered true interaction(s) (along with subsets of interacting features) and ranked these interactions as most significant across all candidates. In contrast, pRF consistently recovered true interactions in 2/4 response generating models. However, pRF tended to rank any candidate interaction that included the true interaction as highly significant (even if the candidate included inactive features). Finally, logistic regression recovered the true interaction in 2/4 response generating models, though at a lower significance level and with greater variability. For a detailed description of our simulation study and results, see Section S1.7 in S1 File.

Summary of findings

Here we introduce the epiTree pipeline, building on a novel tree-based method and PCS framework, to detect epistasis. epiTree is grounded in a biologically relevant dimensional-reduction schema derived from tissue-specific gene expression to derive variant-level interactions, and additionally, may detect interactions of between more than two loci. A particular advantage of our methodology is that it allow for flexible forms of Boolean interactions and additive components. For many interactions, we find that this flexibility provides better fit to the data compared with multiplicative terms used in logistic regression. Moreover, by modeling responses directly on the penetrance scale, our approach avoids the need for transforming responses as, e.g., performed in logistic regression.

Without prior knowledge of the phenotype of red-hair coloration, our approach detects known interactions between MC1R (and nearby genes related to pigmentation) and ASIP, as well as novel interactions between CDK10 and SIAH2. Our findings may suggest a previously unrecognized relationship between SIAH2 and the red hair phenotype. Interestingly, SIAH2 is known to be associated with breast cancer risk and cellular response to hypoxia -a phenotype that has been associated with melanoma and pigmentation in general [56]. The red hair phenotype, which was considered in this manuscript, mainly served as an illustrative example for the epiTree pipeline. We selected this phenotype because of its well established epistasis that provides a positive control for interaction discovery. We stress that for further validation regarding new epistasis findings for the red hair phenotype independent data will be needed. We also consider another case study for the more complex phenotype of multiple sclerosis (MS). While epiTree detects interactions among genes related to processes dysregulated in MS, we are not aware of previously reported findings on gene-level interactions among unlinked genes associated with MS to benchmark against.

Limitations and directions for future work

Our application of iRF is based upon biologically-relevant dimensional reduction in the form of estimated tissue-specific expression values derived from GTeX via PrediXcan [29]. While this approach detects both known and previously undescribed epistasis, it depends on the accuracy of gene expression estimates derived from bulk RNA sequencing on heterogeneous cell types. We make the assumption that tissue-specific data is sufficiently representative of the biology from which the phenotype of interest is derived. This assumption of representativeness is not likely to hold for traits deriving from rare cell types or from tissues and cell-types, which are not sampled in GTeX in the first place. Additionally, we note that the number of individuals in GTeX is not currently representative of all population strata, which serves to limit the application of our approach for detection of epistasis to human phenotypes by excluding diverse genetic populations. Finally, our use of PrediXcan for dimensionality reduction models epistatic relationships at the level of transcription, which may miss variants that are not associated with changes in gene expression levels. Nevertheless, our epiTree test is a stand-alone test for epistasis of any given interaction relative to a phenotype that arises from other knowledge rather than the PrediXcan.

In principle, the epiTree pipeline can detect interactions of arbitrary order. However, the maximal detectable interaction order depends on tree depth—for the iRF screening step an interaction corresponds to subsets of features on individual decision paths. In general, the depths of a decision tree learned from data scales logarithmically with the number of samples. This implies that, in general, we cannot detect interactions of order higher than $\mathcal{O}\left(\text{ln}\right(n\left)\right)$. From a practical perspective, this limitation is fairly minor as even order-3 interactions remain unexplored in many biological settings. A corollary of this limitation is that the number of samples falling into an individual leaf node splitting on d features is $\mathcal{O}(n/2^d)$ (assuming balanced splits), where n is the sample size. Thus for an interaction represented by a single leaf node, n/2^d needs to be large enough to detect deviation from additivity. Thus, detection of higher-order interactions becomes exponentially more difficult as the size of the interaction increases. This problem applies not just to the epiTree pipeline, but to interaction discovery in general, where the size of all possible candidate interactions of order d grows exponentially. However, while traditional methods that screen through all possible interactions are typically restricted to pairwise interaction due to computational reasons, the epiTree pipeline can detect interactions of higher order if the interaction signal is strong enough. More mathematical details on how the order of an interaction enters the complexity of the interaction discovery problem for Boolean interaction recovery with a theoretically tractable version of iRF can be found in [37].

Another potential limitation of the iRF interaction screening is that it does not explicitly take the genomic location of individual features into account. Even on the gene expression level, features that are close-by on the genome are typically highly correlated. Such features might appear exchangeably in interactions, as observed for the red hair phenotype with MC1R, see Fig 5. More generally, it is well known that for random forest based algorithms, correlated features can lead to masking effects, where the effect of a feature might be hidden by another dominating, correlated feature [62]. Finally, as discussed earlier, genetic features that are close-by on the genome might show artificial epistasis that is caused by hidden linkage with unobserved causal variants. In fact, for the red hair phenotype iRF reported various interactions among genes in the vicitinity of MC1R (recall Results section) that we filtered out subsequently. One idea to overcome these limitations, would be to analyze by haplotypes of merged close-by genes into a single hyper-feature. However, two correlated genetic features do not necessarily show the same epistatic behaviour. Thus, there is a general tradeoff. Aggregating features in this way may increase the ability to detect epistasis by concentrating the forest on a single feature. On the other hand, aggregation may lead to loss of specificity that decreases the chance of detecting epistasis. For the red hair analysis, all interactions were either among genes between different chromosomes, in which case there is no correlation or linkage between them, or among genes that were all in the vicinity of MC1R and contained variants that were highly linked to each other (recall Results Section). Therefore, we followed the most simple approach for the red hair analysis and filtered out interactions on the same chromosome. In general, we note that more advanced approaches are possible by taking drop-off of linkage explicitly into account.

Although, the CART tree-based model that we propose for the PCS p-value (see Materials and method and Supplementary Information) can describe significantly more flexible additive and non-additive relationships for the epistasis and non-epistasis model than a single linear or multiplicative term, we also emphasize that not every functional relationship might be well described by a relatively simple CART model. In particular, when a relationship is truly linear or multiplicative, clearly, working with a linear model will be more powerful than a CART approximation. The response surface plots, as shown in e.g., in Fig 7 and discussed in the Results section, are a powerful tool to investigate whether CART or linear models might be more appropriate for the given data at hand. Moreover, we stress that PCS p-values are, in principle, not restricted to CART models, but could be combined with any learning algorithm that appears appropriate for the data. As with any machine learning algorithms there is typically a tradeoff between interpretability and prediction accuracy. We find the CART tree-based models provide a good balance between both, being sufficiently flexible and very straight-forward to interpret, while capturing the thresholding behavior of interactions among many biomolecules.

The epiTree pipeline involves several screening steps in order to result in a set of candidate interactions, such as estimating gene expression features from the SNP data, applying iRF, and fitting individual tree-based null and alternative models. We stress, however, that all these screening steps are based on the training data only, while the final evaluation of the PCS p-values is based on the hold-out test data. With this sample-split approach, which is well accepted as a post-selection tool, we obtain valid inference that avoids overfitting from the screening procedure.

Here, we obtained our training and test data sets from a random sample split. While such an approach results in valid inference, we recognize that using an independent or external data set as test data will be preferred. At the same time, we also recognize that in many cases such an independent data set will not be available (also for our analysis an independent data set for the red hair phenotype was not available to us). Therefore, our PCS p-values are designed in such a way that they explore the full data distribution via bootstrap sampling of the test-data, which makes them more stable to the somewhat arbitrary choice of the random sample split.

The sample splitting paradigm that we follow in our approach holds the advantage of separating the discovery from the inference stage, which prevents overfitting and thus, makes findings more reliable, in general. We note, however, that compared to classical p-values, the need of a separate training set, will generally imply that there is less data available to evaluate the p-values, thus, a decrease in power. However, we note that due to the iRF interaction selection step, in general our approach requires the testing of orders of magnitude fewer hypotheses than any brute-force p-value approach. Thus, multiplicity correction of p-values will be much less of a problem in our pipeline, i.e., less stringent significance thresholds for p-values apply, which can compensate for this effect.