/arboreto-grn-inference

name: "arboreto-grn-inference"
description: "GRN inference from expression via GRNBoost2 (gradient boosting) or GENIE3 (Random Forest). Load matrix, filter by TFs, infer TF-target-importance links, save network. Dask-parallelized to single-cell scale. Core SCENIC component."
license: "BSD-3-Clause"

Arboreto GRN Inference

Overview

Arboreto infers gene regulatory networks (GRNs) from gene expression data using parallelized tree-based regression. For each target gene, it trains a regression model with all other genes (or a specified TF list) as features and emits TF-target-importance triplets. It provides two interchangeable algorithms -- GRNBoost2 (gradient boosting, fast) and GENIE3 (Random Forest, classic) -- sharing identical input/output formats. Computation is Dask-parallelized, scaling from laptop cores to HPC clusters.

When to Use

• Inferring transcription factor-to-target gene regulatory relationships from bulk RNA-seq expression data
• Building gene regulatory networks from single-cell RNA-seq count matrices (cells as rows, genes as columns)
• Generating the adjacency matrix (Step 1) of the pySCENIC regulatory analysis pipeline
• Comparing regulatory network structure across experimental conditions (e.g., control vs treatment)
• Producing consensus regulatory networks by running inference across multiple random seeds
• Validating GRN results by comparing GRNBoost2 and GENIE3 outputs on the same dataset
• For downstream regulon identification and activity scoring, use arboreto output with pySCENIC
• For single-cell preprocessing (QC, normalization, clustering) before GRN inference, use **scanpy-scrna-seq**

Prerequisites

• **Python packages**: `arboreto`, `pandas`, `numpy`, `dask`, `distributed`, `scikit-learn`, `scipy`
• **Data requirements**: Gene expression matrix (genes as columns, observations as rows) in TSV/CSV; optionally a TF list file (one gene name per line)
• **Environment**: Python 3.8+; optional `networkx`, `matplotlib` for visualization

pip install arboreto distributed networkx matplotlib

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: algorithm
    kind: required
    source: user
    ask: "Infer the network with gradient boosting (faster) or random forests (the original GENIE3 formulation)?"
    default: "GRNBoost2"

  - id: D2
    param: expressionMatrix
    kind: required
    source: upstream
    ask: "Which processed expression matrix, and which cells or samples within it, should the network be inferred from?"
    default: null

  - id: D3
    param: regulatorList
    kind: required
    source: literature
    ask: "Restrict candidate regulators to a curated transcription-factor list, or let every gene be a candidate regulator?"
    default: "a species-matched TF list"

  - id: D4
    param: edgeThreshold
    kind: required
    source: user
    ask: "The output ranks every regulator-target pair - how many of the top edges should be kept as the network?"
    default: null

  - id: D5
    param: seed
    kind: never_ask
    source: data
    reason: "Fixes the draw for reproducibility; does not change what the data supports"
    default: 0

  - id: D6
    param: daskClient
    kind: never_ask
    source: data
    reason: "Execution backend; affects runtime only"
    default: "local"

D3 is the decision that separates a usable network from a hairball. Leaving regulators unrestricted lets any correlated gene appear as a regulator, and the result still returns a fully populated ranked table. D4 has no safe default: the algorithm scores all pairs, so where the list is cut is the user's call about the network they intend to interpret.

Quick Start

Complete GRN inference in a single block. The `if __name__ == '__main__':` guard is required because Dask spawns worker processes via multiprocessing.

import pandas as pd
from arboreto.algo import grnboost2
from arboreto.utils import load_tf_names

if __name__ == '__main__':
    # Load expression matrix (observations x genes)
    expression_matrix = pd.read_csv('expression_data.tsv', sep='\t')
    tf_names = load_tf_names('tf_list.txt')  # optional TF filter

    # Infer GRN (uses all local cores by default)
    network = grnboost2(expression_data=expression_matrix,
                        tf_names=tf_names, seed=777)

    # Filter top links and save
    top_network = network[network['importance'] > 1.0]
    top_network.to_csv('grn_output.tsv', sep='\t', index=False, header=False)
    print(f"Inferred {len(network)} links, kept {len(top_network)} above threshold")
    # Example: Inferred 185432 links, kept 12876 above threshold

Workflow

Step 1: Load Expression Data

Arboreto accepts a pandas DataFrame (recommended) or NumPy array. Rows are observations (cells or samples), columns are genes. Gene names must be column headers.

import pandas as pd

# From TSV (genes as columns, observations as rows)
expression_matrix = pd.read_csv('expression_data.tsv', sep='\t')
print(f"Shape: {expression_matrix.shape}  "
      f"({expression_matrix.shape[0]} observations x {expression_matrix.shape[1]} genes)")
# Example: Shape: (5000, 18654)  (5000 observations x 18654 genes)

# From AnnData (e.g., after scanpy preprocessing)
import anndata as ad
adata = ad.read_h5ad('preprocessed.h5ad')
expression_matrix = pd.DataFrame(
    adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
    columns=adata.var_names.tolist()
)
print(f"Converted AnnData: {expression_matrix.shape}")
# Example: Converted AnnData: (5000, 18654)

Step 2: Load Transcription Factor List

Providing a TF list restricts regulators to known transcription factors, reducing computation time and improving biological relevance. If omitted, all genes are treated as potential regulators.

from arboreto.utils import load_tf_names

# From file (one TF name per line)
tf_names = load_tf_names('human_tfs.txt')