bio-gene-regulatory-networks-grn-inference

What this skill does

## 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

| Family | Tool | Citation | Mechanism | Structural bias |
|--------|------|----------|-----------|-----------------|
| Mutual information | ARACNe-AP | Lachmann 2016 *Bioinformatics* | MI + data-processing-inequality pruning of indirect edges | good on feed-forward loops; deletes the direct leg of true FFLs |
| Tree ensemble (RF) | GENIE3 | Huynh-Thu 2010 *PLoS ONE* | per-target random-forest variable importance | good on cascades; trades away FFLs |
| Tree ensemble (GBM) | GRNBoost2 | Moerman 2019 *Bioinformatics* | per-target gradient boosting; fast/scalable | as GENIE3; stochastic without a seed |
| Info-theoretic | CLR / MRNET | Faith 2007; Meyer 2007 | MI z-scored against per-gene background / mRMR | suppress hub artifacts |
| TF activity | VIPER / msVIPER | Alvarez 2016 *Nat Genet* | regulon-enrichment (aREA) on a signature | needs a regulon + Mode of Regulation |

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

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Bulk RNA-seq, want a TF -> target network | GENIE3 or ARACNe-AP | tree-ensemble or MI edge inference |
| Find master regulators of a phenotype | ARACNe regulon -> msVIPER | activity inference is robust to edge errors |
| Score per-sample TF activity for stratification | `viper()` per-sample matrix | turns expression into an activity readout |
| Want robustness / no single best method | ensemble multiple inferences | DREAM5 wisdom-of-crowds |
| Single-cell data with motif resources | -> scenic-regulons | motif pruning adds directness SCENIC-style |
| Just co-expression modules (no direction) | -> coexpression-networks | WGCNA modules, no TF privileging |
| Compare TF activity between conditions | msViper on a 2-group signature | differential activity, not differential edges |

## 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.

```r
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.

```bash
# 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).

```r
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) * s