PaperlayerView exact revision on arXiv

Source-linked AI summary

Iterative Random Forests to detect predictive and stable high-order interactions

Sumanta Basu, Karl Kumbier, James B. Brown, Bin Yu

arXiv:1706.08457v4stat.MLq-bio.GN

TL;DR

High-order genomic interactions are difficult to identify because assays measure individual elements within larger molecular systems. The paper introduces iRF, which combines iteratively feature-weighted Random Forests with generalized Random Intersection Trees and bootstrap stability scoring. Across simulations and genomic case studies, iRF recovered stable high-order rules while maintaining comparable predictive performance, including known and novel regulatory interactions.

  • Problem

    Genomic assays measure individual molecular elements, but understanding gene regulation requires identifying their stable high-order interactions.

  • Method

    iRF sequentially trains feature-weighted Random Forests, extracts decision-path interactions with generalized Random Intersection Trees, and scores their bootstrap stability.

  • Results

    iRF generally matched Random Forest predictive performance while recovering generating rules up to order 8 and identifying known and novel enhancer interactions.

  • Takeaways & Limitations

    iRF provides a computationally efficient framework for generating hypotheses about high-order mechanisms in genome biology.

  • Takeaways & Limitations

    Recovered interactions are candidate hypotheses requiring domain knowledge and experimental evaluation to determine whether they represent domain-relevant or physical interactions.

Abstract

from arXiv · show

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 (RF), 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 same order of computational cost as 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, novel 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 re-discovered a central role of H3K36me3 in chromatin-mediated splicing regulation, and identified novel 5th and 6th order interactions, indicative of multi-valent nucleosomes with specific roles in splicing regulation. By decoupling the order of interactions from the computational cost of identification, iRF opens new avenues of inquiry into the molecular mechanisms underlying genome biology.

1 Introduction

Genomic assays measure individual molecular interactions, while gene regulation depends on larger combinatorial systems. iRF addresses the computational challenge of identifying stable high-order interactions without explicitly enumerating them.

  • Genomic assays often measure single protein targets even though regulatory elements interact as components of larger molecular machines.
  • Existing interaction methods focus on pairwise relationships or require computationally infeasible enumeration of higher-order combinations.
  • iRF uses sequentially feature-weighted Random Forests to stabilize decision paths and generalized Random Intersection Trees to identify prevalent high-order combinations.
  • The method seeks stable, interpretable interactions while retaining Random Forest predictive accuracy and computational efficiency.

2 Our method: iterative Random Forests

iRF combines iterative feature re-weighting, Random Forest decision-path transformations, generalized Random Intersection Trees, and bootstrap stability scoring to discover class-specific high-order interactions.

  • Generalized RIT: Generalized RIT applies randomized intersections to binary feature-index sets derived from Random Forest decision paths and leaf labels.
  • Feature-weighted RF: Feature-weighted forests sample split variables with probabilities proportional to non-negative feature weights, rather than uniformly.
  • RIT: RIT searches for interactions prevalent in one class by intersecting active-feature sets from randomly selected observations and comparing class prevalences.
  • Purpose: The supervised-learning framing links recovered high-order interactions to accurate prediction in feature-weighted forests.
  • Iteratively re-weighted RF: iRF iteratively re-weights Random Forest feature sampling using feature importance from the preceding iteration.
  • Bagged stability scores: Bootstrap stability scores quantify the proportion of bootstrap samples in which each recovered interaction appears as an RIT output.

3 Simulation experiments

Biologically inspired simulations show that iRF preserves predictive performance while recovering high-order generating rules more reliably than standard Random Forests and at lower interaction-search cost.

  • Simulations used Boolean AND, OR, and XOR rules across datasets ranging from 100–5000 observations and 50–2500 features.
  • The workflow combines iterative forests, decision-path interaction extraction, and bootstrap stability assessment.
  • iRF predictive performance with K > 1 was generally comparable with Random Forests using K = 1.
  • iRF recovered generating rules up to order 8 as the most stable interaction in settings where Random Forests rarely recovered interactions above order 2.
  • Iterative re-weighting made early-important features more likely to appear among the first splits later, enabling conditionally important active variables to be selected.

4 Case study I: enhancer elements in Drosophila

In early Drosophila embryos, iRF predicted enhancer status and recovered stable pairwise and higher-order transcription-factor interactions, including an AND-like Zelda–Giant–Krüppel relationship.

  • The study uses early Drosophila embryonic enhancers as a simplified system for examining combinatorial transcription-factor binding and patterned gene expression.
  • iRF achieved an AUC-PR of 0.5 on held-out enhancer-status test data with K = 5.The corresponding MCC was 0.43 with PPV 0.71 after thresholding predicted probabilities to maximize training MCC.
  • 80% of stable pairwise transcription-factor findings with stability score > 0.5 had previously been reported as physical interactions.Highly stable examples included Gt–Kr at 1.0, Gt–Hb at 0.93, and Hb–Kr at 0.73.
  • iRF identified stable higher-order interactions involving Zelda, including Zld–Gt–Twi at 1.0 and Zld–Gt–Kr at 0.7.
  • Response surfaces supported an order-3 AND rule in which sufficiently high binding of Zld, Gt, and Kr marked enhancer activation in a subset of sequences.When Kr occupancy was above its median, active-enhancer proportions reached 60% when Zld and Gt were sufficiently bound; below-median Kr, they remained uniformly low at <10%.

5 Case study II: alternative splicing in a human-derived cell line

iRF identifies stable high-order combinations of chromatin marks and RNA Pol II states associated with alternative splicing in human-derived cells. These interactions include an order-6 pattern and a fully stable order-5 subset linked to high splicing rates.

  • Alternative splicing produces multiple distinct mRNAs from the same gene, diversifying the proteome through regulated spliceosome activity.
  • The extent to which chromatin state and DNA-binding factors interact en masse to modulate co-transcriptional splicing remains unknown.
  • 0.51 AUC-PR was achieved on held-out test data for K = 2.The corresponding held-out-test MCC was 0.47, with PPV 0.72 after thresholding probabilities using training data.
  • An order-6 interaction surrounding H3K36me3 and S2 phospho-Pol II had stability score 0.5, while two order-5 subsets had stability scores 1.0.The six-feature enrichment formed an AND-type rule associated with highly spliced exons.
  • The recovered interactions involved H3K36me3, other chromatin marks, and post-translationally modified RNA Pol II states.The impact of serine 2 phosphorylation of Pol II appeared highly dependent on local chromatin state.

6 Discussion

The discussion presents iRF as a hypothesis-generating tool for nonlinear genomic interactions, while emphasizing experimental validation and methodological boundaries. It also identifies follow-up questions about multivalent histone-mark dependencies and the precise form of discovered rules.

  • iRF identified known and novel interactions in Drosophila enhancer activation and proposed new high-order interactions in human-derived splicing regulation.
  • Experimental validation of complex biological interactions remains necessary and challenging despite emerging genome and epigenome engineering tools.
  • The reported interactions raise questions about proteins mediating dependencies and the role of Phospho-S2 Pol II in splicing regulation.The paper suggests Proteomics on ChIP samples and Co-ChIP as possible approaches for resolving these factors and marks.
  • iRF handles data heterogeneity only implicitly, and detectable interaction order depends directly on tree depth, which is on the order of log2(n).
  • iRF does not distinguish between additive and non-additive interaction forms.The authors are exploring rule-structure tests to clarify rule-response relationships.

S1 Algorithms

The paper presents Random Intersection Trees and iterative Random Forests as algorithmic components for recovering feature combinations from decision-tree structures. The procedures use class-labeled feature sets, tree sampling, bootstrap data, and iterative weighted forest fitting.

  • Random Intersection Trees: Random Intersection Trees take class-labeled feature sets as input, with each feature set paired with a binary label.
  • Random Intersection Trees: The RIT tree has depth D and nchild children at each node across levels 0 through D − 1.
  • Random Intersection Trees: Each RIT node samples an observation uniformly from the specified class, and the procedure initializes S1 from the first sampled feature set.
  • Iterative Random Forests: Iterative Random Forests begin with data D, class C, iteration and bootstrap parameters, and uniform feature weights w(1) = (1/p, . . . , 1/p).
  • Iterative Random Forests: At each iteration, the algorithm fits a feature-weighted RF on D and then generates bootstrap samples for further forest fitting.
  • Iterative Random Forests: For each bootstrap sample, the fitted forest records feature-set and label pairs for observations falling in tree leaf nodes.

S2 Remarks on iRF

iRF combines iterative feature re-weighting, Random Forests, and generalized Random Intersection Trees to recover stable high-order interactions while retaining predictive performance. Re-weighting improves recovery of higher-order rules and supports stability assessment across bootstrap samples, but recovered interactions remain hypothesis-generating until experimentally evaluated.

  • S2.1 Iterative re-weighting: Feature re-weighting acts as soft dimension reduction, encouraging Random Forests to select stable features without removing marginally unimportant features involved in high-order interactions.Iterative re-weighting is also described as probabilistic regularization of the base Random Forest learner.
  • S2.1 Iterative re-weighting: Order-8 data-generating rules become highly stable with feature re-weighting, whereas iRF rarely recovers interactions above order 2 without it under comparable simulations.The comparison is between k > 1 and k = 1 settings.
  • S2.1 Iterative re-weighting: 9 order-3 enhancer interactions exceed stability score 0.5 with k = 5, compared with none without iterative re-weighting.The fourth iteration also recovers additional lower-stability interactions of orders 3, 4, and 5.
  • S2.1 Iterative re-weighting: Iteratively re-weighted and unweighted Random Forests achieve similar predictive accuracy on held-out test data.This supports using re-weighting to alter feature selection and interaction recovery without a reported loss in held-out prediction accuracy.
  • S2.2 Bootstrap sampling: Inner bootstrap sampling robustifies tree-derived binary representations, while outer bootstrap sampling assesses uncertainty in recovered interactions.The final iteration fits weighted Random Forests on different bootstrap samples to estimate interaction stability.
  • S2.6 Sensitivity to tuning parameters: Prediction accuracy is highly consistent across tuning choices except at k = 1 in the splicing data, while recovered-interaction stability is generally correlated across RIT settings.The examined prediction metric is AUC-PR; interaction stability was evaluated across RIT parameters.
  • S2.9 Interaction evaluation through prediction: Candidate high-order interactions are intended for hypothesis generation, and determining which are domain-relevant requires subject-matter knowledge and experimental follow-up.The paper notes that non-additivity tests may help prioritize interactions for biological examination.

S3 Data processing

The case studies transform enhancer activity and chromatin-mediated splicing measurements into supervised classification datasets. The processing choices define labels, feature summaries, data splits, and class-imbalance handling for iRF analysis.

  • Enhancer data: 7809 genomic sequences were evaluated for enhancer activity using a stable-integration transgenic reporter assay.Each tested sequence was 100–3000 nt and inserted into a targeted genomic site before reporter detection.
  • Enhancer data: Approximately 10% of enhancer sequences were labeled active, and the 3912 training and 3897 test observations were randomly divided.The sequences were selected using prior biological knowledge rather than sampled randomly from the genome.
  • Splicing data: The splicing dataset combines exon PSI values with ChIP-seq measurements of chromatin marks and transcription-factor binding in K562 cells.The data were used to identify interactions associated with chromatin-mediated splicing.
  • Splicing data: 270 features were produced by taking the maximum normalized fold-enrichment for each ChIP assay over each exon’s genomic region.Only internal exons with at least 100 RPKM were retained; PSI above 70% defined class 1 and PSI below 30% defined the other class.
  • Splicing data: Approximately 90% of splicing observations belonged to class 1, so AUC-PR was reported for class 0 using balanced training and test sets of 11911 and 11912 observations.The analysis used iRF with B = 30 and K = 2.

S4 Evaluating Drosophila enhancer interactions

The Drosophila enhancer analysis compares iRF-recovered transcription-factor interactions with previously reported physical interactions. Stable pairwise findings show substantial agreement with prior evidence, while comprehensive validation of high-order interactions remains limited.

  • S4 Evaluating Drosophila enhancer interactions: 16 (80%) of the 20 pairwise TF interactions with stability score > 0.5 had been previously reported in physical or regulatory interaction forms.The prior evidence included one TF regulating the other or joint binding associated with increased expression of target genes.
  • S4 Evaluating Drosophila enhancer interactions: Only a small number of high-order interactions have been studied previously, preventing comprehensive comparison of the high-order interactions recovered by iRF.The paper specifically notes that existing examples are insufficient for evaluating the full set of recovered high-order interactions.
  • S4 Evaluating Drosophila enhancer interactions: Table S1 lists previously identified pairwise TF interactions recovered by iRF with stability score > 0.5.The table is used to compare stable iRF findings with prior reports of physical interactions.

S5 Simulation experiments

Simulation studies show that iterative re-weighting improves recovery and stability of high-order interactions while generally preserving predictive performance, but performance degrades under higher noise, strong feature correlation, and competing interaction rules.

  • Boolean-rule simulations: RF never recovers the true order-4 interaction in the Boolean simulations, whereas iRF identifies it as the most stable order-4 interaction given enough training observations.The improvement is accompanied by more low-stability false positives, allowing true positives to remain distinguishable by stability ranking.
  • Overall performance: iRF generally maintains predictive performance comparable to RF while improving recovery of high-order interactions through iterative re-weighting.The authors attribute improved recovery to stabilized decision paths and adaptive regularization.
  • Noise robustness: Increasing iteration improves predictive accuracy and interaction recovery at 15% and 20% noise, but higher noise substantially reduces both absolute performance and improvement rates.Weight separation between active and inactive variables develops more slowly at higher noise levels, indicating that additional iterations may be needed in low signal-to-noise settings.
  • Mixture models: When one rule dominates, iRF fails to recover the other, while at π = 0.75 it recovers subsets of both XOR and AND interactions.The limitation arises because global feature-weight updates do not distinguish features important for different observations or feature-space regions.
  • Correlated features: At ρ = 0.75, iRF recovers the full order-8 interaction only at k = 5 and simultaneously recovers many more false positives, especially under block covariance.Distinguishing highly correlated features is difficult for feature-selection methods generally; grouping related variables is suggested when prior knowledge is available.
  • Enhancer-data simulations: Iterations k > 1 recover true high-order interactions more frequently at fixed sample sizes while maintaining comparable predictive accuracy across the enhancer-data noise settings.True positives are consistently more stable, although later iterations also recover many low-stability false high-order interactions.

S6 Computational cost of detecting high-order interaction

iRF scales more favorably than Rulefit3 for high-order interaction detection: runtime grows linearly with feature count for iRF but exponentially for Rulefit3, while iRF can search arbitrarily high-order interactions.

  • Runtime scaling: Rulefit3 runtime increases exponentially with the number of features, whereas iRF runtime increases linearly in the enhancer-data comparison.Rulefit3 searches interactions up to order 3, while iRF leverages deep RF decision trees to search arbitrarily high-order interactions.
  • Runtime scaling: The computational difference follows from Rulefit3’s exponentially growing interaction search space rather than an implementation optimization.The benchmark used balanced enhancer-data samples with 7809 observations and up to 80 features.
  • Resource profile: iRF completed the splicing-dataset profile in 26 minutes 59 seconds with a reported memory consumption of 499910 Mb.The profile used n = 11911, p = 270, B = 30, K = 3, and 24 CPU cores.

S7 List of datasets

The supplementary materials provide scripts, processed datasets, and figures documenting iRF analyses, simulations, tuning, stability, and interaction validation.

  • Scripts: Scripts cover iRF analyses for enhancer, splicing, Boolean simulations, enhancer simulations, and runtime comparisons with Rulefit3.The accompanying Rulefit3 package is configured for Linux systems.
  • Datasets: The enhancer dataset contains 7,809 genomic regions and 80 assay features, with 3,912 training and 3,897 testing observations.The response vector marks active elements.
  • Datasets: The splicing dataset contains 23,823 exons and 270 assay features, divided into 11,911 training and 11,912 testing observations.The response vector marks highly spliced exons, and assay metadata include identifiers and collapsed replicate groups.
  • Supplementary analyses: Supplementary figures evaluate prediction accuracy, interaction recovery, feature weights, decision-path entry depth, and sensitivity to RF and RIT parameters.The evaluations include enhancer, splicing, Boolean-rule, noise, and mixture-model settings.
Loading 1706.08457v4…