End-to-End Workflow: From FASTQ to Regulatory Elements

Overview

This vignette walks through the complete PRIME toolkit pipeline, from raw CAGE sequencing reads to a genome-wide map of predicted regulatory elements. The pipeline involves three tools:

FASTQ files
    │
    ▼
┌──────────┐
│ PRIMEprep │  Shell: QC → trim → rRNA filter → STAR map → G-correction → BigWig
└────┬─────┘
     │ BigWig files (plus/minus strand)
     ▼
┌──────────┐
│  PRIME   │  R package (extends CAGEfightR)
│          │  ● quantify CTSSs (via CAGEfightR)
│          │  ● pool replicates, subsample, normalize
│          │  ● call tag clusters & divergent loci
│          │  ● estimate noise
│          │  ● complexity & saturation analysis
└────┬─────┘
     │ CTSS RangedSummarizedExperiment objects
     ▼
┌────────────┐
│ PRIMEmodel │  R + Python (LightGBM)
│            │  ● genome-wide prediction of regulatory elements
│            │  ● focal prediction on defined regions
│            │  ● post-processing & BED output
└────────────┘

Each step is described in detail in the dedicated vignettes:

vignette("preprocessing") — PRIMEprep
vignette("ctss-processing") — CTSS quantification and QC
vignette("divergent-loci") — divergent loci calling
vignette("normalization-batches") — normalization and subsampling
vignette("noise-estimation") — genomic background noise
vignette("prediction") — regulatory element prediction

This vignette provides a concise end-to-end example.

Step 1: Preprocess CAGE data with PRIMEprep

Run the PRIMEprep shell pipeline on each FASTQ file:

cd PRIMEprep
./PRIMEprep.sh \
    -f /data/sample1.fastq.gz \
    -g /ref/genome.fa \
    -b /ref/STAR_index \
    -d /ref/rRNAdust_db \
    -t 16 \
    -o /results/sample1

The key output is bw_files/ containing strand-specific BigWig files:

/results/sample1/bw_files/sample1.plus.bw
/results/sample1/bw_files/sample1.minus.bw

See vignette("preprocessing") for full details on PRIMEprep.

Step 2: Quantify CTSSs with CAGEfightR and PRIME

library(CAGEfightR)
library(PRIME)
library(SummarizedExperiment)
library(GenomeInfoDb)
library(ggplot2)

Bundled example data path

dir_design <- system.file("extdata", "design_matrix_first10pct.tsv",
                          package = "PRIME")
dir_bw <- system.file("extdata", "cage_bw", package = "PRIME")

Load the design matrix using :

design <- read.table(dir_design, header = TRUE, sep = "\t")
rownames(design) <- design$Name

# Required columns: Name, BigWigPlus, BigWigMinus
# Optional columns: CellLine, Type, Replicate, Batch, ...

Build the CTSS object

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) <- design$Name

ctss <- CAGEfightR::quantifyCTSSs(plusStrand = bw_plus,
                                  minusStrand = bw_minus,
                                  design      = design)
ctss <- CAGEfightR::calcTotalTags(ctss)
ctss <- CAGEfightR::calcTPM(ctss)
ctss <- CAGEfightR::calcPooled(ctss)

# Restrict to standard chromosomes
ctss <- GenomeInfoDb::keepStandardChromosomes(ctss, pruning.mode = "coarse")

Pool replicates and subsample

# Pool biological replicates
ctss_pooled <- PRIME::poolReplicates(
    ctss,
    replicates = paste(colData(ctss)$CellLine,
                       colData(ctss)$Type,
                       sep = "_")
)
ctss_pooled <- CAGEfightR::calcTPM(ctss_pooled)

# Subsample to the minimum library depth for fair comparison
min_depth <- min(colData(ctss_pooled)$totalTags)
ctss_pooled_sub  <- PRIME::subsampleTarget(ctss_pooled, target = min_depth)
ctss_pooled_sub  <- CAGEfightR::calcTPM(ctss_pooled_sub)

# Remove singletons
ctss_clean <- PRIME::rmSingletons(ctss_pooled_sub)

Step 3: Assess complexity and noise

See vignette("ctss-processing") and vignette("noise-estimation") for full details.

Assess complexity on the replicated data (before pooling)

complexity <- PRIME::calcCTSSComplexity(ctss,
                                        step           = 2e6,
                                        minCTSSsupport = 1,
                                        CTSSunexpressed = 0)

# Set up the dataframe for plotting
df_complexity <- data.frame(
  totalTags   = complexity$totalTags,
  numberCTSSs = complexity$numberCTSSs,
  sample      = complexity$sample
)

# Remove the missing points
df_complexity <- na.omit(df_complexity)

# Plot
ggplot(df_complexity, aes(x = totalTags,
                          y = numberCTSSs,
                          color = sample,
                          group = sample)) +
  geom_line(linewidth = 0.8, alpha = 0.8) +
  geom_point(size = 1.2) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Total Tags", y = "Number of CTSSs", color = NULL,
       title = "Library complexity saturation curves") +
  theme_bw()

Step 4: Call divergent loci

Call divergent loci

CTSSs_by_chr <- split(ctss_clean, GenomicRanges::seqnames(ctss_clean))

DLs <- lapply(names(CTSSs_by_chr), function(chr) {
    PRIME::divergentLociTCsSummit(
        ctss         = CTSSs_by_chr[[chr]],
        callingAssay = "counts.noSingletons"
    )
})

# Combine results across chromosomes
DLs <- do.call(c, DLs)

Step 5: Predict regulatory elements with PRIMEmodel

Typically, you can run PRIMEmodel::predict() on replicated CTSS data “without any prior filtering.” However, this depends entirely on the specific downstream analysis steps you want the PRIMEmodel to predict.

As an example continuing from the previous steps, the ctss_pooled_sub object will be passed through the model to predict the regions of interest between the subsampled pooled replicates of K562 whole-cell CAGE with the singleton removed (ctss_clean) from step 2-3.

Then, those prediction results will be used by PRIMEmodel::predictFocal() to call each individual replicate within the predicted regions.

This is a downsized version of the analysis performed with pooled samples on the FANTOM5 dataset. Further details can be found in the manuscript.

To ensure the process remains computationally efficient, only chromosomes 16 and 17 will be selected as an example.

pool_target_samples <- c("K562_C")
ctss_clean_slt <- ctss_clean[, pool_target_samples]
ctss_clean_slt <- GenomeInfoDb::keepSeqlevels(ctss_clean_slt,
                                              c("chr16", "chr17"),
                                              pruning.mode = "coarse")

The full setup for the reticulate Python environment is described in vignette("prediction"). The full PRIMEmodel prediction can be run on this example data when that environment is available:

library(reticulate)
library(PRIMEmodel)
library(GenomicRanges)

region_result <- PRIMEmodel::predict(
    ctss_clean_slt,
    python_path     = reticulate::py_config()$python,
    score_threshold = 0.75,
    score_diff      = 0.1,
    num_cores       = NULL,
    keep_tmp        = FALSE
)

For a lightweight, self-contained example, load the precomputed prediction regions bundled with PRIMEmodel. This also keeps the focal prediction example independent of the optional genome-wide prediction chunk above.

region_result <- readRDS(system.file(
    "extdata",
    "predicted_regions_gr.rds",
    package = "PRIMEmodel"
))

Run focal prediction on the predicted region for each K562 whole-cell CAGE:

target_samples <- c("K562_C1", "K562_C2")
ctss_slt <- ctss[, target_samples]
ctss_slt <- GenomeInfoDb::keepSeqlevels(ctss_slt,
                                        c("chr16", "chr17"),
                                        pruning.mode = "coarse")

focal_result <- PRIMEmodel::predictFocal(
    ctss_slt,
    region_result,
    python_path = reticulate::py_config()$python
  )

Full reproducibility

The complete analysis code for the associated publication (Einarsson, Navamajiti, et al. 2026) is available at: https://github.com/anderssonlab/nucCAGE_PRIME_paper

That repository contains step-by-step R Markdown notebooks for all analyses described in the paper.