Iterative random forests to discover predictive and stable high-order interactions

Abstract

Genomics has revolutionized biology, enabling the interrogation of whole transcriptomes, genome-wide binding sites for proteins, and many other molecular processes. However, individual genomic assays measure elements that interact in vivo as components of larger molecular machines. Understanding how these high-order interactions drive gene expression presents a substantial statistical challenge. Building on random forests (RFs) and random intersection trees (RITs) and through extensive, biologically inspired simulations, we developed the iterative random forest algorithm (iRF). iRF trains a feature-weighted ensemble of decision trees to detect stable, high-order interactions with the same order of computational cost as the RF. We demonstrate the utility of iRF for high-order interaction discovery in two prediction problems: enhancer activity in the early Drosophila embryo and alternative splicing of primary transcripts in human-derived cell lines. In Drosophila, among the 20 pairwise transcription factor interactions iRF identifies as stable (returned in more than half of bootstrap replicates), 80% have been previously reported as physical interactions. Moreover, third-order interactions, e.g., between Zelda (Zld), Giant (Gt), and Twist (Twi), suggest high-order relationships that are candidates for follow-up experiments. In human-derived cells, iRF rediscovered a central role of H3K36me3 in chromatin-mediated splicing regulation and identified interesting fifth- and sixth-order interactions, indicative of multivalent nucleosomes with specific roles in splicing regulation. By decoupling the order of interactions from the computational cost of identification, iRF opens additional avenues of inquiry into the molecular mechanisms underlying genome biology.

Keywords: genomics; high-order interaction; interpretable machine learning; random forests; stability.

Fig. 1.

iRF workflow. Iteratively reweighted RFs (blue boxes) are trained on full data $\mathcal{D}$ and pass Gini importance as weights to the next iteration. In iteration $K$ (red box), feature-weighted RFs are grown using $w^{(K)}$ on $B$ bootstrap samples of the full data ${\mathcal{D}}_{(1)},\dots ,{\mathcal{D}}_{(B)}$. Decision paths and predicted leaf node labels are passed to the RIT (green box), which computes prevalent interactions in the RF ensemble. Recovered interactions are scored for stability across (outer-layer) bootstrap samples.

Fig. 2.

(A) Accuracy of iRF (AUC-PR) in predicting active elements from TF binding and histone modification data. (B) The 20 most stable interactions recovered by iRF after five iterations. Interactions that are a strict subset of another interaction with stability score $\ge 0.5$ have been removed for cleaner visualization. iRF recovers known interactions among Gt, Kr, and Hb and interacting roles of master regulator Zld. (C) Surface maps demonstrating the proportion of active enhancers by quantiles of Zld, Gt, and Kr binding (held-out test data). On the subset of data where Kr binding is lower than the median Kr level, the proportion of active enhancers does not change with Gt and Zld. On the subset of data with Kr binding above the median level, the structure of the response surface reflects an order-3 AND interaction: Increased levels of Zld, Gt, and Kr binding are indicative of enhancer status for a subset of observations. (D) Quantiles of Zld, Gt, and Kr binding grouped by enhancer status (balanced sample of held-out test data). The block of active elements highlighted in red represents the subset of observations for which the AND interaction is active.

Fig. 3.

(A) Accuracy of iRF (AUC-PR) in classifying included exons from excluded exons in held-out test data. iRF shows $7\%$ increase in AUC-PR over RF. (B) An order-6 interaction recovered by iRF (stability score 0.5) displayed on a superheat map which juxtaposes two separately clustered heat maps of exons with high and low splicing rates. Coenrichment of all six plotted features reflects an AND-type rule indicative of high splicing rates for the exons highlighted in red (held-out test data). The subset of Pol II, S2 phospho-Pol II, H3K36me3, H3K79me2, and H4K20me1 was recovered as an order-5 interaction in all bootstrap samples (stability score $1.0$). (C) The 20 most stable interactions recovered in the second iteration of iRF. Interactions that are a strict subset of another interaction with stability score $\ge 0.5$ have been removed for cleaner visualization.