Bio Gene Regulatory Networks Grn Inference

Version Compatibility

Reference examples tested with: VIPER 1.36+ (Bioconductor), GENIE3 1.24+ (Bioconductor), ARACNe-AP (Java, build from source), arboreto 0.1.6+ (Python GRNBoost2).

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

GENIE3 expects the expression matrix as genes-in-rows, samples-in-columns (the transpose of the WGCNA convention); a transposed matrix silently produces a meaningless network.

GRN Inference and TF Activity

"Infer a gene regulatory network from my bulk expression data and find the master regulators" -> Reverse-engineer TF -> target edges from an expression matrix, assemble them into regulons, then score the protein activity of each TF from a gene-expression signature.

• R: GENIE3() (tree-ensemble) or ARACNe-AP (mutual information) for edges
• R: viper::aracne2regulon() -> msviper() / viper() for TF activity

The Single Most Important Modern Insight -- Activity, Not Edges: A Regulon Is a Multiplexed Reporter

A GRN inferred from observational expression is, by default, an undirected statistical association graph: correlation and mutual information are symmetric and cannot distinguish TF -> target from target -> TF or from a shared upstream driver. Tree-ensemble methods (GENIE3/GRNBoost2) appear directed only because they restrict predictors to a TF list -- that direction is an input assumption, not an inference. Individual inferred edges are therefore unreliable, and benchmarks confirm it (below).

The paradigm that survives this (the Califano-lab lineage: ARACNe -> VIPER) is to stop trusting individual edges and instead read TF protein activity from the regulon as a whole. VIPER treats a regulon (a TF and its inferred targets, each carrying a Mode-of-Regulation sign) as a multiplexed reporter assay: even if many edges are wrong, the coordinated up/down shift of the targets in a signature is a robust estimate of the regulator's activity (Alvarez 2016 Nat Genet 48:838). This is why VIPER can identify an active master regulator whose own mRNA is unchanged -- because the protein is regulated post-transcriptionally. Master-regulator analysis (MARINa/VIPER) is a fundamentally different computation from "the most-connected node," and edge-level precision matters less than activity inference. The Mode-of-Regulation signs are load-bearing: without them VIPER collapses to a plain enrichment test.

Method Taxonomy

DREAM5 (Marbach 2012 Nat Methods 9:796): no single method dominates; families make complementary errors, so an ensemble ("wisdom of crowds") is the most robust. And methods that excel on synthetic data collapse on real eukaryotic data (yeast was near-random) because TF and target mRNA decorrelate -- synthetic AUPRC does not transfer.

Decision Tree by Scenario

Edge Inference with GENIE3 (R)

Goal: Reverse-engineer a ranked TF -> target network from a bulk expression matrix.

Approach: Fit a per-target random forest predicting each gene from candidate regulators (TFs); the regulator's variable importance is the edge weight. Restrict predictors to a TF list to orient edges.

library(GENIE3)

# GENIE3 convention: genes in ROWS, samples in COLUMNS (transpose of WGCNA).
expr <- as.matrix(read.csv('normalized_counts.csv', row.names = 1))
regulators <- readLines('tf_list.txt')                  # candidate TFs only

set.seed(42)                                            # tree ensembles are stochastic
weight_matrix <- GENIE3(expr, regulators = regulators, treeMethod = 'RF',
                        K = 'sqrt', nTrees = 1000, nCores = 8)
link_list <- getLinkList(weight_matrix)                 # ranked edge list (NOT thresholded)
head(link_list)

Edge Inference with ARACNe-AP (Java CLI)

Goal: Build a mutual-information network with indirect edges pruned.

Approach: ARACNe-AP is a two-phase Java pipeline: compute the MI threshold, run many bootstrap reconstructions, then consolidate them (with data-processing-inequality pruning) into a final network. Running a single bootstrap or skipping consolidation is the classic misuse.

# Phase 1: MI threshold at a chosen p-value (needs the expression matrix + TF list).
java -Xmx32G -jar aracne.jar -e expr.txt -o out/ --tfs tf_list.txt \
    --pvalue 1E-8 --seed 1 --calculateThreshold

# Phase 2: many bootstraps (vary --seed) -- 100 is conventional.
for s in $(seq 1 100); do
    java -Xmx32G -jar aracne.jar -e expr.txt -o out/ --tfs tf_list.txt \
        --pvalue 1E-8 --seed $s
done

# Phase 3: consolidate bootstraps into the final network (DPI + a Poisson edge-significance
# test with Bonferroni correction across bootstraps).
java -Xmx32G -jar aracne.jar -o out/ --consolidate

TF Activity with VIPER / msVIPER (R)

Goal: Infer transcription-factor protein activity from a regulon and an expression signature, and rank master regulators.

Approach: Convert an ARACNe network into a regulon object (assigning each target a Mode-of-Regulation sign and likelihood), build a null model by sample permutation, then run msVIPER on a two-group signature (master regulators) or VIPER per sample (activity matrix).

library(viper)

# Build the regulon from the ARACNe network + matched expression (assigns Mode of Regulation).
# ARACNe-AP network.txt has a header + 4 columns (Regulator, Target, MI, p-value); viper's
# '3col' reader wants Regulator/Target/MI with no header, so strip them first (shell):
#   tail -n +2 out/network.txt | cut -f1-3 > net_3col.txt
# ('adj' is the legacy ARACNE adjacency-matrix format, not ARACNe-AP.)
regulon <- aracne2regulon('net_3col.txt', eset, format = '3col')

# msVIPER: master regulators of a two-group contrast.
signature <- rowTtest(eset, pheno = 'group', group1 = 'tumor', group2 = 'normal')
sig_z <- (qnorm(signature$p.value / 2, lower.tail = FALSE) * sign(signature$statistic))[, 1]
nullmodel <- ttestNull(eset, pheno = 'group', group1 = 'tumor', group2 = 'normal', per = 1000)
mra <- msviper(sig_z, regulon, nullmodel)
summary(mra)                                            # top master regulators by NES

# VIPER: a per-sample TF-activity matrix for clustering/stratification.
activity <- viper(eset, regulon, method = 'scale')

For single cells or tissues lacking a matched network, metaVIPER integrates multiple interactomes; DIGGIT then intersects master regulators with genetic alterations to nominate causal drivers.