Predicting Regulatory Elements with PRIMEmodel

Overview

PRIMEmodel uses a LightGBM model trained on CAGE data to predict the locations of active regulatory elements (enhancers and promoters) genome-wide or in user-defined focal regions.

PRIMEmodel takes a CTSS RangedSummarizedExperiment as input and outputs the GRangesList with predicted regulatory element coordinates and scores.

Two prediction modes are available:

Genome-wide prediction (PRIMEmodel::predict()) — scans the whole genome using a sliding-window approach over tag clusters.
Focal prediction (PRIMEmodel::predictFocal()) — restricts prediction to user-defined regions of interest (e.g. FANTOM5 enhancers or called divergent loci, pooled prediction from PRIMEmodel).

Setup

Setup a PRIME-Compatible Virtual Environment

txt_file <- file.path(find.package("PRIMEmodel"), "envfile", "environment.txt")
txt_content <- readLines(txt_file)

library(reticulate)

# Install Python locally
reticulate::install_python(version = "3.9.22")

# Create a new virtual environment for running PRIMEmodel
virtualenv_create("PRIMEmodel", python = "3.9.22")

# Activate the environment and install dependencies
use_virtualenv("PRIMEmodel")
virtualenv_install(envname = "PRIMEmodel", packages = txt_content)

# Optional (if not already covered by environment.txt)
virtualenv_install(envname = "PRIMEmodel", packages = "matplotlib")

library(PRIMEmodel)
library(GenomicRanges)

# Activate virtual environment
use_virtualenv("PRIMEmodel", required = TRUE)

PRIMEmodel requires a Python environment (>= 3.9) with LightGBM. See the installation guide for details. Full function documentation is available on the PRIMEmodel pkgdown site.

Obtaining CTSS input

# design_matrix: data frame with columns Name, BigWigPlus, BigWigMinus
# e.g. from PRIME
dir_design <- system.file("extdata", "design_matrix_first10pct.tsv", package = "PRIME")
design <- read.table(dir_design, header = TRUE, sep = "\t")
rownames(design) <- design$Name

# dir_bw: directory containing BigWig files from PRIMEprep
# e.g. from PRIME
dir_bw <- system.file("extdata", "cage_bw", package = "PRIME")

PRIMEmodel accepts CTSS data produced by either of two approaches: ### Option A: PRIMEmodel convenience function

ctss <- PRIMEmodel::plc_get_ctss_from_bw(dir_bw, design)
ctss <- GenomeInfoDb::keepStandardChromosomes(ctss, pruning.mode = "coarse")

Option B: Standard CAGEfightR / PRIME workflow.

library(CAGEfightR)

bw_plus  <- rtracklayer::BigWigFileList(file.path(dir_bw, design$BigWigPlus))
bw_minus <- rtracklayer::BigWigFileList(file.path(dir_bw, design$BigWigMinus))
names(bw_plus) <- names(bw_minus) <- rownames(design)

ctss <- CAGEfightR::quantifyCTSSs(plusStrand = bw_plus,
                                  minusStrand = bw_minus,
                                  design      = design)
ctss <- CAGEfightR::calcTotalTags(ctss)
ctss <- GenomeInfoDb::keepStandardChromosomes(ctss, pruning.mode = "coarse")

Option B is useful when you have already performed CTSS QC steps with PRIME (e.g., replicate pooling, subsampling, singleton removal) before running PRIMEmodel. See vignette("ctss-processing") for the full CTSS workflow.

Due to the computationally heavy nature of these processing pipelines, we highly recommend running this function on a server. By default, it utilizes half of the available CPU cores for parallelization. If parallelization is unavailable, the function will automatically fall back to sequential processing, which will significantly increase the execution time.

If you are running this example on a personal or local computer, we highly recommend using the pre-prepared .rds object, which contains only data for chromosomes 16 and 17 to test the function time-efficiently.

ctss <- readRDS(system.file("extdata", "ctss_rse_chr16to17.rds", package = "PRIMEmodel"))

Genome-wide prediction

PRIMEmodel::predict() runs the complete 6-step prediction pipeline:

Extract CTSS — retrieve CTSS counts from the RSE object
Identify tag clusters (TCs) — call unidirectional TCs genome-wide
Slide through TCs — apply a sliding window over each TC
Normalized profiles — normalize CTSS profiles within windows
Predict probabilities — apply the LightGBM model
Post-processing — merge windows, apply score threshold, write BED

result <- PRIMEmodel::predict(
    ctss,
    python_path     = reticulate::py_config()$python,
    score_threshold = 0.75,
    score_diff      = 0.1,
    num_cores       = NULL,   # NULL uses half of available cores (max 25 cores)
    keep_tmp        = FALSE
)

The result is a GRangesList object with predicted regulatory element coordinates and prediction scores.

Running via the bash script

For very large datasets or for reproducibility on HPC clusters, the bash script interface is recommended:

# Genome-wide prediction
cd PRIMEmodel/genomewide_prediction
./PRIMEmodel.sh --config bash_config_predict.sh --predict

# Run individual steps
./PRIMEmodel.sh --config bash_config_predict.sh -1 -2 -3 -4 -5 -6

The individual step flags (-1 through -6) correspond to the six pipeline steps described above, allowing re-runs from any intermediate checkpoint.

Focal prediction

PRIMEmodel::predictFocal() restricts prediction to a set of pre-defined genomic regions (e.g., called divergent loci, candidate enhancers, pooled prediction from PRIMEmodel)

The regions object is a GRanges object containing the regions of interest. These should either have a uniform width of 401 bp or include a ‘thick’ column to allow internal resizing/extension within the function.

# regions: a GRanges object with regions of interest
regions <- readRDS(system.file("extdata", "predicted_regions_gr.rds", package = "PRIMEmodel"))

result_focal <- PRIMEmodel::predictFocal(
    ctss,
    regions,
    python_path = reticulate::py_config()$python
)

Focal prediction via bash

cd PRIMEmodel/genomewide_prediction
./PRIMEmodel.sh --config bash_config_predictFocal.sh --predictFocal

Choosing between genome-wide and focal prediction

Mode	Use case
Genome-wide	Discovery of all active regulatory elements in the genome
Focal	Scoring/classifying a specific set of candidate regions

Focal prediction is faster and is appropriate when you already have a set of candidate regions (e.g., from divergent loci calling or an external database).

Output

Both predict() and predictFocal() return a GRanges object and write BED files to the specified output directory. Each predicted element has:

Genomic coordinates (chr, start, end, strand)
A prediction score (0–1; higher = more confident regulatory element)
Element type annotation (enhancer / promoter)