diff --git a/DESCRIPTION b/DESCRIPTION
index 95aa9de..6f74d86 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -4,13 +4,13 @@ Title: Single Cell Oriented Reconstruction of PANDA Individually Optimized Netwo
 Version: 1.3.3
 Authors@R: c(person("Daniel","Osorio", email = "daniecos@uio.no", role = c("aut","cre"), comment = c(ORCID = "0000-0003-4424-8422")),
   person("Marieke L.","Kuijjer", email = "marieke.kuijjer@ncmbm.uio.no", role = c("aut"), comment = c(ORCID = "0000-0001-6280-3130")))
-Description: Constructs cell-type–specific gene regulatory networks from single-cell RNA-sequencing data. The method implements the SCORPION algorithm, which first aggregates individual cells into super-cells and then applies PANDA (Passing Attributes between Networks for Data Assimilation) to infer transcription factor–target regulatory relationships. It also provides statistical methods for differential edge analysis.
+Description: Constructs cell-type-specific gene regulatory networks from single-cell RNA-sequencing data. The method implements the SCORPION algorithm, which first aggregates individual cells into super-cells and then applies PANDA (Passing Attributes between Networks for Data Assimilation) to infer transcription factor-target regulatory relationships. It also provides statistical methods for differential edge analysis.
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
 Depends: R (>= 3.5.0)
 Imports: cli, methods, irlba, igraph, RANN, Matrix, pbapply, dplyr, furrr, future
-Suggests: RhpcBLASctl, testthat, mori
+Suggests: RhpcBLASctl, testthat, mori, circlize, biomaRt, fgsea
 URL: https://github.com/kuijjerlab/SCORPION
 BugReports: https://github.com/kuijjerlab/SCORPION/issues
 RoxygenNote: 7.3.3
diff --git a/NAMESPACE b/NAMESPACE
index 141030c..07d209d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
+export(circosEdges)
+export(enrichEdges)
+export(maEdges)
 export(regressEdges)
 export(runSCORPION)
 export(scorpion)
@@ -11,6 +14,7 @@ importFrom(RANN,nn2)
 importFrom(cli,cli_abort)
 importFrom(cli,cli_alert_info)
 importFrom(cli,cli_alert_success)
+importFrom(cli,cli_alert_warning)
 importFrom(cli,cli_h1)
 importFrom(cli,cli_progress_along)
 importFrom(dplyr,"%>%")
@@ -23,6 +27,9 @@ importFrom(furrr,future_map)
 importFrom(future,multisession)
 importFrom(future,plan)
 importFrom(future,sequential)
+importFrom(grDevices,adjustcolor)
+importFrom(grDevices,hcl.colors)
+importFrom(graphics,legend)
 importFrom(igraph,cluster_louvain)
 importFrom(igraph,cluster_walktrap)
 importFrom(igraph,contract)
diff --git a/R/circosEdges.R b/R/circosEdges.R
new file mode 100644
index 0000000..fa73687
--- /dev/null
+++ b/R/circosEdges.R
@@ -0,0 +1,670 @@
+#' @title Circos plot of differential network edges
+#' @description Draws a circular (Circos) plot of transcription factor to target
+#'   links from a \code{\link{testEdges}} two-sample result. Genes are placed on
+#'   their genomic coordinates, links are coloured continuously by significance,
+#'   flagged as known or novel against an optional a priori network, and genes
+#'   belonging to supplied gene sets can be labelled around the circle.
+#' @author Daniel Osorio <daniecos@uio.no>
+#' @param edgesDF A data.frame produced by \code{\link{testEdges}} (two-sample or
+#'   paired). Must contain the columns \code{tf}, \code{target},
+#'   \code{log2FoldChange} and \code{pAdj} (\code{pValue} is used as a fallback
+#'   when \code{pAdj} is absent).
+#' @param species Ensembl dataset name passed to \pkg{biomaRt} when
+#'   \code{geneCoords} is \code{NULL}, e.g. \code{"hsapiens_gene_ensembl"} for
+#'   human or \code{"mmusculus_gene_ensembl"} for mouse. See
+#'   \code{biomaRt::listDatasets()} for the full multi-species list.
+#' @param geneCoords Optional data.frame supplying gene coordinates from any
+#'   source (overrides the \pkg{biomaRt} download). Must have columns
+#'   \code{gene}, \code{chr}, \code{start} and \code{end}.
+#' @param priorNet Optional a priori TF-target network whose first two columns
+#'   are the TF and target. Links present here are labelled \code{"known"}, all
+#'   others \code{"novel"}. Accepts a data.frame or a matrix.
+#' @param geneSets Optional gene-set annotation used to label genes around the
+#'   circle: either a path to a GMT file or a named list of character vectors.
+#' @param pAdjThreshold Numeric significance cutoff applied to \code{pAdj}
+#'   (or \code{pValue} when \code{pAdj} is missing). Default \code{0.05}.
+#' @param log2FCThreshold Numeric minimum absolute \code{log2FoldChange} required
+#'   to draw a link. Default \code{0}.
+#' @param maxEdges Integer cap on the number of links drawn; when exceeded, the
+#'   most significant edges are kept. Default \code{500}.
+#' @param colorBy Name of the \code{edgesDF} column mapped to the link colour
+#'   ramp. Default \code{"log2FoldChange"}, giving a continuous diverging colour
+#'   scale centred at zero. When the default is used but \code{log2FoldChange} is
+#'   absent (e.g. single-sample \code{\link{testEdges}} output), it falls back to
+#'   \code{meanEdge}.
+#' @param linkColors Length-3 vector of colours for the low, mid and high ends of
+#'   \code{colorBy}. Default blue-white-red; a diverging ramp is used when
+#'   \code{colorBy} has negative values, otherwise a sequential low-to-high ramp.
+#' @param lwdRange Length-2 numeric giving the minimum and maximum link line
+#'   width; each link's thickness is scaled linearly within this range by its
+#'   \code{-log10} adjusted p-value. Default \code{c(0.5, 4)}.
+#' @param nmaxTF,nmaxTarget Integers giving how many TFs and targets to label,
+#'   selected by the largest absolute out-degree and in-degree respectively.
+#'   Use \code{NULL} or \code{Inf} to label all. Defaults \code{20}.
+#' @param knownColor,novelColor Border colours distinguishing known from novel
+#'   links. Defaults grey and orange.
+#' @param geneSetColors Optional named vector mapping gene-set names to colours.
+#'   When \code{NULL}, colours are generated automatically.
+#' @param chromosomes Optional character vector restricting and ordering the
+#'   chromosomes shown. When \code{NULL}, all chromosomes present are used.
+#' @param mainChromosomesOnly Logical; when \code{TRUE} (the default) and
+#'   \code{chromosomes} is \code{NULL}, only the main chromosomes (numbered,
+#'   plus X, Y and MT) are kept and unplaced scaffolds/contigs are dropped.
+#' @param ensemblMirror \pkg{biomaRt} mirror to query: one of \code{"www"},
+#'   \code{"useast"} or \code{"asia"}. Default \code{"www"}.
+#' @param transparency Numeric link transparency in \code{[0, 1]} (0 is opaque).
+#'   Default \code{0.5}.
+#' @param hRatio Numeric in \code{[0, 1]} controlling how far link ribbons bend
+#'   toward the circle centre; smaller values give flatter, less tangled links.
+#'   Default \code{0.6}.
+#' @param fontFamily Font family used for all plot text, e.g. \code{"sans"}
+#'   (Helvetica/Arial, the default) for a publication look.
+#' @param legend Logical; whether to draw legends for effect size, novelty and
+#'   gene sets. Default \code{TRUE}.
+#' @return Invisibly, a list with \code{edges} (the plotted links annotated with
+#'   coordinates and novelty) and \code{coords} (the gene coordinate table with
+#'   \code{outDegree}, \code{inDegree} and total \code{degree} columns, each the
+#'   sum of \code{log2FoldChange} over a gene's outgoing / incoming links). The
+#'   function is called for the side effect of drawing the plot.
+#' @seealso \code{\link{testEdges}}, \code{\link{runSCORPION}}
+#' @details Requires the \pkg{circlize} package, and \pkg{biomaRt} when gene
+#'   coordinates are downloaded automatically (\code{geneCoords = NULL}). Genes
+#'   without coordinates, and links whose TF or target lacks coordinates, are
+#'   dropped with a message.
+#' @examples
+#' \dontrun{
+#' data(scorpionTest)
+#' nets <- runSCORPION(
+#'   gexMatrix = scorpionTest$gex,
+#'   tfMotifs = scorpionTest$tf,
+#'   ppiNet = scorpionTest$ppi,
+#'   cellsMetadata = scorpionTest$metadata,
+#'   groupBy = c("donor", "region")
+#' )
+#' res <- testEdges(
+#'   networksDF = nets,
+#'   testType = "two.sample",
+#'   group1 = grep("--T$", colnames(nets), value = TRUE),
+#'   group2 = grep("--N$", colnames(nets), value = TRUE)
+#' )
+#'
+#' # Human coordinates auto-downloaded from Ensembl, known/novel vs a prior net
+#' circosEdges(
+#'   edgesDF = res,
+#'   species = "hsapiens_gene_ensembl",
+#'   priorNet = scorpionTest$tf,
+#'   geneSets = "hallmark.gmt"
+#' )
+#' }
+#' @export
+#' @importFrom cli cli_abort cli_alert_info cli_alert_success cli_alert_warning
+#' @importFrom grDevices adjustcolor hcl.colors
+#' @importFrom graphics legend
+circosEdges <- function(edgesDF,
+                        species = "hsapiens_gene_ensembl",
+                        geneCoords = NULL,
+                        priorNet = NULL,
+                        geneSets = NULL,
+                        pAdjThreshold = 0.05,
+                        log2FCThreshold = 0,
+                        maxEdges = 500L,
+                        colorBy = "log2FoldChange",
+                        linkColors = c("#2166AC", "#F7F7F7", "#B2182B"),
+                        lwdRange = c(0.5, 4),
+                        nmaxTF = 20L,
+                        nmaxTarget = 20L,
+                        knownColor = "grey60",
+                        novelColor = "#D95F02",
+                        geneSetColors = NULL,
+                        chromosomes = NULL,
+                        mainChromosomesOnly = TRUE,
+                        ensemblMirror = "www",
+                        transparency = 0.5,
+                        hRatio = 0.6,
+                        fontFamily = "sans",
+                        legend = TRUE) {
+
+  if (!requireNamespace("circlize", quietly = TRUE)) {
+    cli::cli_abort(c(
+      "Package {.pkg circlize} is required to draw Circos plots.",
+      "i" = "Install it with {.code install.packages('circlize')}."
+    ))
+  }
+
+  # Validate input columns
+  if (!is.data.frame(edgesDF)) {
+    cli::cli_abort("{.arg edgesDF} must be a data.frame from {.fn testEdges}.")
+  }
+  # Single-sample testEdges output has no log2FoldChange; fall back to meanEdge
+  if (missing(colorBy) && !("log2FoldChange" %in% colnames(edgesDF)) &&
+      "meanEdge" %in% colnames(edgesDF)) {
+    colorBy <- "meanEdge"
+    cli::cli_alert_info("No {.field log2FoldChange} column; colouring by {.field meanEdge}.")
+  }
+  required_cols <- c("tf", "target")
+  if (!identical(colorBy, "log10Padj")) required_cols <- c(required_cols, colorBy)
+  missing_cols <- setdiff(required_cols, colnames(edgesDF))
+  if (length(missing_cols)) {
+    cli::cli_abort("{.arg edgesDF} is missing required column{?s}: {.val {missing_cols}}")
+  }
+  if (length(linkColors) != 3) {
+    cli::cli_abort("{.arg linkColors} must be a vector of exactly 3 colours.")
+  }
+  if (transparency < 0 || transparency > 1) {
+    cli::cli_abort("{.arg transparency} must be between 0 and 1.")
+  }
+
+  # Significance column (pAdj preferred, pValue fallback)
+  sig_col <- if ("pAdj" %in% colnames(edgesDF)) "pAdj" else if ("pValue" %in% colnames(edgesDF)) "pValue" else NA_character_
+
+  edges <- edgesDF
+  edges$tf <- as.character(edges$tf)
+  edges$target <- as.character(edges$target)
+
+  # Filter by significance and effect size
+  if (!is.na(sig_col)) {
+    keep <- !is.na(edges[[sig_col]]) & edges[[sig_col]] < pAdjThreshold
+    edges <- edges[keep, , drop = FALSE]
+  }
+  # Keep every significant edge for degree computation (before fold-change
+  # filtering and the maxEdges cap, which only affect which links are drawn)
+  edges_sig <- edges
+  if (log2FCThreshold > 0 && "log2FoldChange" %in% colnames(edges)) {
+    keep <- !is.na(edges$log2FoldChange) & abs(edges$log2FoldChange) >= log2FCThreshold
+    edges <- edges[keep, , drop = FALSE]
+  }
+  if (identical(colorBy, "log10Padj")) {
+    if (is.na(sig_col)) {
+      cli::cli_abort("{.arg colorBy = 'log10Padj'} requires a {.field pAdj} or {.field pValue} column.")
+    }
+    edges$log10Padj <- -log10(pmax(edges[[sig_col]], .Machine$double.xmin))
+  }
+  edges <- edges[!is.na(edges[[colorBy]]), , drop = FALSE]
+  if (nrow(edges) == 0) {
+    cli::cli_abort("No edges pass the significance / fold-change thresholds.")
+  }
+
+  # Cap to the most significant edges
+  if (!is.na(sig_col) && nrow(edges) > maxEdges) {
+    edges <- edges[order(edges[[sig_col]]), , drop = FALSE][seq_len(maxEdges), , drop = FALSE]
+    cli::cli_alert_info("Showing top {maxEdges} edges by {sig_col}.")
+  }
+
+  genes <- unique(c(edges$tf, edges$target))
+
+  # Resolve gene coordinates
+  if (is.null(geneCoords)) {
+    coords <- .fetchGeneCoords(genes, species, ensemblMirror)
+  } else {
+    coords <- .validateGeneCoords(geneCoords)
+  }
+  coords <- coords[coords$gene %in% genes, , drop = FALSE]
+  coords <- coords[!duplicated(coords$gene), , drop = FALSE]
+  if (nrow(coords) == 0) {
+    cli::cli_abort("None of the edge genes could be mapped to coordinates.")
+  }
+
+  # Restrict to main chromosomes (drop unplaced scaffolds/contigs)
+  if (isTRUE(mainChromosomesOnly) && is.null(chromosomes)) {
+    bare <- sub("^chr", "", coords$chr, ignore.case = TRUE)
+    is_main <- grepl("^[0-9]+$", bare) | toupper(bare) %in% c("X", "Y", "MT", "M")
+    dropped <- unique(coords$chr[!is_main])
+    coords <- coords[is_main, , drop = FALSE]
+    if (length(dropped)) {
+      cli::cli_alert_info("Excluding {length(dropped)} non-main sequence{?s} (scaffolds/contigs).")
+    }
+  }
+
+  # Optionally restrict chromosomes
+  if (!is.null(chromosomes)) {
+    coords <- coords[coords$chr %in% as.character(chromosomes), , drop = FALSE]
+  }
+  if (nrow(coords) == 0) {
+    cli::cli_abort("No genes remain after chromosome filtering.")
+  }
+  coord_idx <- match(genes, coords$gene)
+  names(coord_idx) <- genes
+
+  # Keep only links whose endpoints both have coordinates
+  tf_i <- match(edges$tf, coords$gene)
+  tg_i <- match(edges$target, coords$gene)
+  mappable <- !is.na(tf_i) & !is.na(tg_i)
+  if (!all(mappable)) {
+    cli::cli_alert_warning("Dropping {sum(!mappable)} link{?s} with unmapped TF or target.")
+    edges <- edges[mappable, , drop = FALSE]
+    tf_i <- tf_i[mappable]
+    tg_i <- tg_i[mappable]
+  }
+  if (nrow(edges) == 0) {
+    cli::cli_abort("No links remain after coordinate mapping.")
+  }
+
+  # Flag known vs novel against the a priori network
+  novelty <- rep("novel", nrow(edges))
+  if (!is.null(priorNet)) {
+    prior <- as.data.frame(priorNet, stringsAsFactors = FALSE)
+    if (ncol(prior) < 2) {
+      cli::cli_abort("{.arg priorNet} must have at least two columns (TF, target).")
+    }
+    prior_key <- paste(as.character(prior[[1]]), as.character(prior[[2]]), sep = "\r")
+    edge_key <- paste(edges$tf, edges$target, sep = "\r")
+    novelty[edge_key %in% prior_key] <- "known"
+  }
+  edges$novelty <- novelty
+
+  # Per-gene degree: sum of log2FoldChange over outgoing (TF) and incoming
+  # (target) links, using every edge under pAdjThreshold (not just drawn links)
+  deg_col <- if ("log2FoldChange" %in% colnames(edges_sig)) "log2FoldChange" else colorBy
+  out_deg <- tapply(edges_sig[[deg_col]], edges_sig$tf, sum, na.rm = TRUE)
+  in_deg <- tapply(edges_sig[[deg_col]], edges_sig$target, sum, na.rm = TRUE)
+  coords$outDegree <- as.numeric(out_deg[coords$gene])
+  coords$outDegree[is.na(coords$outDegree)] <- 0
+  coords$inDegree <- as.numeric(in_deg[coords$gene])
+  coords$inDegree[is.na(coords$inDegree)] <- 0
+  coords$degree <- coords$outDegree + coords$inDegree
+
+  # Build chromosome layout (BED-like ranges spanning each chromosome)
+  chr_order <- if (!is.null(chromosomes)) as.character(chromosomes) else .orderChr(coords$chr)
+  chr_order <- chr_order[chr_order %in% coords$chr]
+  chr_layout <- do.call(rbind, lapply(chr_order, function(ch) {
+    s <- coords[coords$chr == ch, , drop = FALSE]
+    rng <- range(c(s$start, s$end), na.rm = TRUE)
+    data.frame(chr = ch, start = rng[1], end = rng[2], stringsAsFactors = FALSE)
+  }))
+
+  # Link fill colours: diverging (signed values) or sequential (all non-negative)
+  vals <- edges[[colorBy]]
+  diverging <- any(vals < 0, na.rm = TRUE)
+  if (diverging) {
+    maxabs <- max(abs(vals), na.rm = TRUE)
+    if (!is.finite(maxabs) || maxabs == 0) maxabs <- 1
+    value_range <- c(-maxabs, maxabs)
+    col_fun <- circlize::colorRamp2(c(-maxabs, 0, maxabs), linkColors)
+  } else {
+    rng <- range(vals, na.rm = TRUE)
+    if (!all(is.finite(rng)) || diff(rng) == 0) rng <- c(0, max(1, rng[2]))
+    value_range <- rng
+    col_fun <- circlize::colorRamp2(c(rng[1], mean(rng), rng[2]), linkColors)
+  }
+  fill_cols <- grDevices::adjustcolor(col_fun(vals), alpha.f = 1 - transparency)
+  border_cols <- ifelse(edges$novelty == "known", knownColor, novelColor)
+
+  # Link thickness scaled linearly by -log10 adjusted p-value
+  if (!is.na(sig_col)) {
+    sig_vals <- -log10(pmax(edges[[sig_col]], .Machine$double.xmin))
+    srng <- range(sig_vals, na.rm = TRUE)
+    if (!all(is.finite(srng)) || diff(srng) == 0) {
+      link_lwd <- rep(mean(lwdRange), nrow(edges))
+    } else {
+      link_lwd <- lwdRange[1] + (sig_vals - srng[1]) / (srng[2] - srng[1]) * (lwdRange[2] - lwdRange[1])
+    }
+  } else {
+    link_lwd <- rep(lwdRange[1], nrow(edges))
+  }
+
+  region1 <- data.frame(
+    chr = coords$chr[tf_i], start = coords$start[tf_i], end = coords$start[tf_i],
+    stringsAsFactors = FALSE
+  )
+  region2 <- data.frame(
+    chr = coords$chr[tg_i], start = coords$start[tg_i], end = coords$start[tg_i],
+    stringsAsFactors = FALSE
+  )
+
+  # Gene-set colour assignment (optional)
+  set_colors <- NULL
+  gene_set_list <- NULL
+  if (!is.null(geneSets)) {
+    sets <- if (is.character(geneSets) && length(geneSets) == 1) .parseGMT(geneSets) else geneSets
+    if (!is.list(sets) || is.null(names(sets))) {
+      cli::cli_abort("{.arg geneSets} must be a GMT file path or a named list of gene vectors.")
+    }
+    gene_set_list <- sets
+    set_names <- names(sets)
+    if (is.null(geneSetColors)) {
+      set_colors <- stats::setNames(grDevices::hcl.colors(length(set_names), "Dark3"), set_names)
+    } else {
+      set_colors <- geneSetColors
+    }
+  }
+
+  # Select the top TFs and targets to label by |out-degree| / |in-degree|,
+  # styling TFs in bold (font 2) and targets in italics (font 3).
+  tf_genes <- unique(edges$tf)
+  target_genes <- unique(edges$target)
+  pick_top <- function(cands, score_col, n) {
+    if (!length(cands)) return(character(0))
+    if (is.null(n) || !is.finite(n)) return(cands)
+    sc <- abs(coords[[score_col]][match(cands, coords$gene)])
+    cands[order(-sc)][seq_len(min(as.integer(n), length(cands)))]
+  }
+  label_genes <- union(pick_top(tf_genes, "outDegree", nmaxTF),
+                       pick_top(target_genes, "inDegree", nmaxTarget))
+  lab_idx <- which(coords$gene %in% label_genes)
+  label_font <- ifelse(coords$gene[lab_idx] %in% tf_genes, 2L, 3L)
+  # Colour TF labels by out-degree sign: blue if negative, red if positive
+  # (fully saturated endpoints, no fade). Targets are black.
+  # Cap the ring colour scale at the 90th percentile of non-zero degrees so
+  # typical genes show saturated colour instead of washing out near the max.
+  deg_vals <- c(coords$outDegree, coords$inDegree)
+  nz <- abs(deg_vals[deg_vals != 0 & is.finite(deg_vals)])
+  deg_absmax <- if (length(nz)) as.numeric(stats::quantile(nz, 0.9, names = FALSE)) else 0
+  if (!is.finite(deg_absmax) || deg_absmax == 0) {
+    deg_absmax <- max(abs(deg_vals), na.rm = TRUE)
+  }
+  if (!is.finite(deg_absmax) || deg_absmax == 0) deg_absmax <- 1
+  deg_col_fun <- circlize::colorRamp2(c(-deg_absmax, 0, deg_absmax), linkColors)
+  is_tf_lab <- coords$gene[lab_idx] %in% tf_genes
+  label_col <- rep("black", length(lab_idx))
+  label_col[is_tf_lab] <- ifelse(
+    coords$outDegree[lab_idx][is_tf_lab] < 0, linkColors[1], linkColors[3]
+  )
+  label_bed <- data.frame(
+    chr = coords$chr[lab_idx], start = coords$start[lab_idx], end = coords$start[lab_idx],
+    gene = coords$gene[lab_idx], font = label_font, col = label_col,
+    stringsAsFactors = FALSE
+  )
+
+  # Draw
+  op <- graphics::par(family = fontFamily, xpd = NA, mar = c(0, 0, 0, 0))
+  on.exit(graphics::par(op), add = TRUE)
+  circlize::circos.clear()
+  on.exit(circlize::circos.clear(), add = TRUE)
+  n_sectors <- nrow(chr_layout)
+  gaps <- rep(1.5, n_sectors)
+  # Seam gap widens with the longest track label so names always fit
+  seam_labels <- c("in-degree", "out-degree", names(gene_set_list))
+  gaps[n_sectors] <- min(90, max(15, max(nchar(seam_labels)) * 2.1))
+  circlize::circos.par(
+    cell.padding = c(0, 0, 0, 0), points.overflow.warning = FALSE,
+    gap.after = gaps, start.degree = 90, track.margin = c(0.006, 0.004),
+    canvas.xlim = c(-1.2, 1.2), canvas.ylim = c(-1.2, 1.2)
+  )
+  circlize::circos.genomicInitialize(chr_layout, plotType = NULL)
+
+  # Outermost ring: gene / TF labels
+  if (!is.null(label_bed) && nrow(label_bed) > 0) {
+    circlize::circos.genomicLabels(
+      label_bed[, c("chr", "start", "end", "gene")],
+      labels.column = 4,
+      col = label_bed$col,
+      line_col = label_bed$col,
+      font = label_bed$font,
+      side = "outside",
+      cex = 0.5
+    )
+  }
+
+  # Shared tile geometry for the degree rings (deg_col_fun defined above)
+  span_max <- max(chr_layout$end - chr_layout$start, na.rm = TRUE)
+  tile_w <- span_max / 300
+
+  ring_bed <- data.frame(
+    chr = coords$chr,
+    start = pmax(0, coords$start - tile_w), end = coords$start + tile_w,
+    stringsAsFactors = FALSE
+  )
+  seam_label <- function(name, col = "grey30") {
+    if (circlize::CELL_META$sector.numeric.index == 1) {
+      circlize::circos.text(
+        circlize::CELL_META$cell.xlim[1], 0.5, labels = name,
+        facing = "downward", adj = c(1, 0.5), cex = 0.5, col = col
+      )
+    }
+  }
+
+  # Thin heatmap ring of a per-gene value on the diverging degree scale
+  circlize::circos.par(track.margin = c(0.006, 0.004))
+  deg_ring <- function(v, name) {
+    circlize::circos.genomicTrack(
+      cbind(ring_bed, value = v), ylim = c(0, 1), track.height = 0.03,
+      bg.border = "grey85", bg.lwd = 0.5,
+      panel.fun = function(region, value, ...) {
+        # Draw weakest first so the strongest tiles sit on top when they overlap
+        ord <- order(abs(value[[1]]), decreasing = FALSE, na.last = TRUE)
+        circlize::circos.genomicRect(
+          region[ord, , drop = FALSE], value[ord, , drop = FALSE],
+          ytop = 1, ybottom = 0, border = NA,
+          col = deg_col_fun(value[[1]][ord])
+        )
+        seam_label(name)
+      }
+    )
+  }
+
+  # In-degree ring: always the outermost data track, just inside the labels
+  deg_ring(coords$inDegree, "in-degree")
+
+  # One thin presence ring per gene set (member genes coloured by set),
+  # drawn just inside the in-degree ring
+  if (!is.null(gene_set_list) && length(gene_set_list)) {
+    for (sn in names(gene_set_list)) {
+      scol <- if (!is.null(set_colors)) unname(set_colors[[sn]]) else "grey40"
+      member <- as.integer(coords$gene %in% gene_set_list[[sn]])
+      local({
+        col_s <- scol
+        nm <- sn
+        circlize::circos.genomicTrack(
+          cbind(ring_bed, value = member), ylim = c(0, 1), track.height = 0.022,
+          bg.border = "grey85", bg.lwd = 0.5,
+          panel.fun = function(region, value, ...) {
+            keep <- value[[1]] == 1
+            if (any(keep)) {
+              circlize::circos.genomicRect(
+                region[keep, , drop = FALSE], value[keep, , drop = FALSE],
+                ytop = 1, ybottom = 0, col = col_s, border = NA
+              )
+            }
+            seam_label(nm)
+          }
+        )
+      })
+    }
+  }
+
+  # Chromosome name band, then the out-degree ring innermost
+  circlize::circos.track(
+    ylim = c(0, 1), track.height = 0.045, bg.col = "grey92", bg.border = NA,
+    panel.fun = function(x, y) {
+      circlize::circos.text(
+        circlize::CELL_META$xcenter, circlize::CELL_META$ycenter,
+        labels = circlize::CELL_META$sector.index,
+        cex = 0.7, font = 2, facing = "inside", niceFacing = TRUE
+      )
+    }
+  )
+  deg_ring(coords$outDegree, "out-degree")
+
+  # Draw strongest links last so the signal sits on top of faint edges
+  draw_order <- order(abs(vals), decreasing = FALSE, na.last = TRUE)
+  circlize::circos.genomicLink(
+    region1[draw_order, , drop = FALSE], region2[draw_order, , drop = FALSE],
+    col = fill_cols[draw_order], border = border_cols[draw_order],
+    lwd = link_lwd[draw_order], h.ratio = hRatio
+  )
+
+  if (isTRUE(legend)) {
+    present_sets <- names(gene_set_list)
+    color_label <- if (identical(colorBy, "log10Padj")) "-log10(padj)" else colorBy
+    .drawCircosLegends(
+      col_fun = col_fun, valueRange = value_range, colorBy = color_label,
+      hasPrior = !is.null(priorNet), knownColor = knownColor, novelColor = novelColor,
+      set_colors = if (length(present_sets)) set_colors[present_sets] else NULL,
+      lwdRange = lwdRange, widthShown = !is.na(sig_col),
+      degCol_fun = deg_col_fun, degRange = c(-deg_absmax, deg_absmax),
+      degColBy = deg_col
+    )
+  }
+
+  cli::cli_alert_success("Drew {nrow(edges)} link{?s} across {nrow(chr_layout)} chromosome{?s}.")
+  invisible(list(edges = edges, coords = coords))
+}
+
+#' Parse a GMT gene-set file into a named list.
+#' @keywords internal
+.parseGMT <- function(path) {
+  if (!file.exists(path)) {
+    cli::cli_abort("GMT file not found: {.path {path}}")
+  }
+  lines <- readLines(path, warn = FALSE)
+  lines <- lines[nzchar(trimws(lines))]
+  sets <- lapply(lines, function(l) {
+    parts <- strsplit(l, "\t", fixed = TRUE)[[1]]
+    genes <- if (length(parts) > 2) parts[-(1:2)] else character(0)
+    genes[nzchar(genes)]
+  })
+  names(sets) <- vapply(lines, function(l) strsplit(l, "\t", fixed = TRUE)[[1]][1], character(1))
+  sets
+}
+
+#' Validate a user-supplied gene coordinate table.
+#' @keywords internal
+.validateGeneCoords <- function(geneCoords) {
+  if (!is.data.frame(geneCoords)) {
+    cli::cli_abort("{.arg geneCoords} must be a data.frame.")
+  }
+  needed <- c("gene", "chr", "start", "end")
+  missing_cols <- setdiff(needed, colnames(geneCoords))
+  if (length(missing_cols)) {
+    cli::cli_abort("{.arg geneCoords} is missing column{?s}: {.val {missing_cols}}")
+  }
+  data.frame(
+    gene = as.character(geneCoords$gene),
+    chr = as.character(geneCoords$chr),
+    start = as.numeric(geneCoords$start),
+    end = as.numeric(geneCoords$end),
+    stringsAsFactors = FALSE
+  )
+}
+
+#' Download gene coordinates from Ensembl via biomaRt.
+#' @keywords internal
+.fetchGeneCoords <- function(genes, species, mirror) {
+  if (!requireNamespace("biomaRt", quietly = TRUE)) {
+    cli::cli_abort(c(
+      "Package {.pkg biomaRt} is required to download gene coordinates.",
+      "i" = "Install it with {.code BiocManager::install('biomaRt')}, or pass {.arg geneCoords}."
+    ))
+  }
+  cli::cli_alert_info("Querying Ensembl ({species}) for {length(genes)} genes via biomaRt.")
+  mart <- biomaRt::useEnsembl(biomart = "genes", dataset = species, mirror = mirror)
+  attrs <- c("external_gene_name", "chromosome_name", "start_position", "end_position")
+  bm <- biomaRt::getBM(
+    attributes = attrs, filters = "external_gene_name",
+    values = genes, mart = mart
+  )
+  bm <- bm[!is.na(bm$start_position) & nzchar(bm$chromosome_name), , drop = FALSE]
+  data.frame(
+    gene = as.character(bm$external_gene_name),
+    chr = as.character(bm$chromosome_name),
+    start = as.numeric(bm$start_position),
+    end = as.numeric(bm$end_position),
+    stringsAsFactors = FALSE
+  )
+}
+
+#' Natural-sort chromosome names (1..N, then X, Y, MT, then the rest).
+#' @keywords internal
+.orderChr <- function(chr) {
+  u <- unique(as.character(chr))
+  bare <- sub("^chr", "", u, ignore.case = TRUE)
+  num <- suppressWarnings(as.numeric(bare))
+  special <- match(toupper(bare), c("X", "Y", "MT", "M"))
+  # numeric first (by value), then X/Y/MT, then everything else alphabetically
+  key1 <- ifelse(!is.na(num), 1L, ifelse(!is.na(special), 2L, 3L))
+  key2 <- ifelse(!is.na(num), num, ifelse(!is.na(special), special, 0))
+  u[order(key1, key2, u)]
+}
+
+#' Draw effect-size, novelty and gene-set legends on a Circos plot.
+#' @keywords internal
+.drawCircosLegends <- function(col_fun, valueRange, colorBy, hasPrior,
+                               knownColor, novelColor, set_colors,
+                               lwdRange = NULL, widthShown = FALSE,
+                               degCol_fun = NULL, degRange = NULL,
+                               degColBy = "log2FoldChange") {
+  # Render known measures as tidy plotmath; fall back to the raw string
+  nice_label <- function(txt) {
+    if (!is.character(txt)) return(txt)
+    switch(txt,
+      "log2FoldChange" = expression(log[2] ~ "fold change"),
+      "-log10(padj)"   = expression(-log[10] * "(" * p[adj] * ")"),
+      txt
+    )
+  }
+
+  y_top <- -1.05     # common top edge shared by every legend title
+  gap   <- 0.08      # horizontal gap between blocks
+  barw  <- 0.45      # colour-bar length
+  bw    <- 0.03      # colour-bar thickness (thin)
+  # Measure text so colour bars line up with the legend swatch rows
+  th     <- graphics::strheight("Ag", cex = 0.68, font = 2)  # title height
+  kh     <- graphics::strheight("Ag", cex = 0.60)            # key-row height
+  bar_cy <- y_top - th - 0.5 * kh - 0.015                    # first-key centre
+
+  # A continuous colour-bar block, self-centred within its measured width
+  bar_block <- function(cfun, rng, title) {
+    tw <- graphics::strwidth(title, cex = 0.68, font = 2)
+    w  <- max(barw, tw)
+    list(w = w, draw = function(xl) {
+      cx  <- xl + w / 2
+      bxl <- cx - barw / 2
+      bxr <- cx + barw / 2
+      ytt <- bar_cy + bw / 2
+      ybb <- bar_cy - bw / 2
+      n <- 128L
+      xs <- seq(bxl, bxr, length.out = n + 1L)
+      gvals <- seq(rng[1], rng[2], length.out = n)
+      graphics::rect(xs[-(n + 1L)], ybb, xs[-1L], ytt,
+                     col = cfun(gvals), border = NA, xpd = NA)
+      graphics::rect(bxl, ybb, bxr, ytt, border = "grey40", lwd = 0.6, xpd = NA)
+      breaks <- pretty(rng, n = 4)
+      breaks <- breaks[breaks >= rng[1] & breaks <= rng[2]]
+      tick_x <- bxl + (breaks - rng[1]) / diff(rng) * (bxr - bxl)
+      graphics::segments(tick_x, ybb, tick_x, ybb - 0.018, lwd = 0.6, xpd = NA)
+      graphics::text(tick_x, ybb - 0.04, labels = format(breaks, digits = 2),
+                     adj = c(0.5, 1), cex = 0.5, xpd = NA)
+      graphics::text(cx, y_top, labels = title, adj = c(0.5, 1),
+                     cex = 0.68, font = 2, xpd = NA)
+    })
+  }
+
+  # A discrete graphics::legend block: measure its width, then draw at an anchor
+  legend_block <- function(...) {
+    m <- graphics::legend(x = 0, y = 0, ..., plot = FALSE)
+    list(w = m$rect$w, draw = function(xl)
+      graphics::legend(x = xl, y = y_top, ..., xjust = 0, yjust = 1, xpd = NA))
+  }
+
+  # Assemble only the blocks that apply to this plot. The gene-set ring seam
+  # labels already name each set, so no separate gene-set legend is drawn.
+  blocks <- list(bar_block(col_fun, valueRange, nice_label(colorBy)))
+  if (!is.null(degCol_fun) && !is.null(degRange)) {
+    deg_title <- if (identical(degColBy, "log2FoldChange")) {
+      expression("degree" ~ (Sigma ~ log[2] * "FC"))
+    } else {
+      bquote("degree" ~ (Sigma ~ .(degColBy)))
+    }
+    blocks <- c(blocks, list(bar_block(degCol_fun, degRange, deg_title)))
+  }
+  if (isTRUE(hasPrior)) {
+    blocks <- c(blocks, list(legend_block(
+      title = "Link", legend = c("known", "novel"),
+      col = c(knownColor, novelColor), lwd = 3, bty = "n", cex = 0.6,
+      y.intersp = 0.5)))
+  }
+  if (isTRUE(widthShown) && !is.null(lwdRange)) {
+    blocks <- c(blocks, list(legend_block(
+      title = expression(-log[10] * "(" * p[adj] * ")"),
+      legend = c("low", "high"), lwd = lwdRange, col = "grey30",
+      bty = "n", cex = 0.6, y.intersp = 0.5)))
+  }
+
+  # Pack the blocks left-to-right and centre the whole row under the circle
+  widths <- vapply(blocks, function(b) b$w, numeric(1))
+  x <- -(sum(widths) + gap * (length(blocks) - 1)) / 2
+  for (b in blocks) {
+    b$draw(x)
+    x <- x + b$w + gap
+  }
+}
diff --git a/R/data-scorpionTest.R b/R/data-scorpionTest.R
index 63e96d7..4fa6812 100644
--- a/R/data-scorpionTest.R
+++ b/R/data-scorpionTest.R
@@ -1,13 +1,27 @@
 #' @name scorpionTest
 #' @docType data
-#' @title Example single-cell gene expression, motif, and ppi data
-#' @description This data is a list containing three objects. The motif \code{data.frame} describes a set of pairwise connections where a specific known sequence motif of a transcription factor was found upstream of the corresponding gene. The expression \code{dgCMatrix} is a set of 230 gene expression levels measured across 80 cells. Finally, the ppi \code{data.frame} describes a set of known pairwise protein-protein interactions.
+#' @title Example single-cell colorectal cancer data for SCORPION
+#' @description
+#' A list bundling the inputs required to build and compare gene regulatory
+#' networks with SCORPION, derived from a colorectal cancer single-cell
+#' RNA-sequencing experiment. It contains a gene expression matrix, a
+#' transcription factor motif prior, a protein-protein interaction prior, and
+#' cell-level metadata.
 #' @usage data(scorpionTest)
-#' @format A list containing three datasets.
+#' @format A named list with four elements:
 #' \describe{
-#' \item{\code{gex}}{A subsetted version of 10X Genomics' 3k PBMC dataset provided by the \code{Seurat} package.}
-#' \item{\code{tf}}{Subset of the transcription-factor and target gene list provided by the \code{dorothea} package for Homo sapiens.}
-#' \item{\code{ppi}}{The known protein-protein interactions and the combined score downloaded from the STRING database}
+#'   \item{\code{gex}}{A \code{dgCMatrix} gene expression matrix with 300 genes
+#'     (rows) and 1,954 cells (columns).}
+#'   \item{\code{tf}}{A \code{data.frame} of transcription factor-target motif
+#'     pairs from DoRothEA with columns \code{source_genesymbol},
+#'     \code{target_genesymbol} and \code{weight} (371,738 rows).}
+#'   \item{\code{ppi}}{A \code{data.frame} of protein-protein interactions with
+#'     columns \code{source_genesymbol}, \code{target_genesymbol} and
+#'     \code{weight} (4,076 rows).}
+#'   \item{\code{metadata}}{A \code{data.frame} of cell-level annotations with
+#'     columns \code{cell_id}, \code{donor}, \code{region} and \code{cell_type}
+#'     (1,954 rows). Region codes are \code{T} (tumor), \code{B} (border) and
+#'     \code{N} (normal).}
 #' }
 #' @examples
 #' # Loading example data
@@ -15,23 +29,4 @@
 #'
 #' # The structure of the data
 #' str(scorpionTest)
-#'
-#' # List of 3
-#' # $ gex:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
-#' # .. ..@ i       : int [1:4456] 1 5 8 11 22 30 33 34 36 38 ...
-#' # .. ..@ p       : int [1:81] 0 47 99 149 205 258 306 342 387 423 ...
-#' # .. ..@ Dim     : int [1:2] 230 80
-#' # .. ..@ Dimnames:List of 2
-#' # .. .. ..$ : chr [1:230] "MS4A1" "CD79B" "CD79A" "HLA-DRA" ...
-#' # .. .. ..$ : chr [1:80] "ATGCCAGAACGACT" "CATGGCCTGTGCAT" "GAACCTGATGAACC" "TGACTGGATTCTCA" ...
-#' # .. ..@ x       : num [1:4456] 1 1 3 1 1 4 1 5 1 1 ...
-#' # .. ..@ factors : list()
-#' # $ tf :'data.frame':	4485 obs. of  3 variables:
-#' #   ..$ tf    : chr [1:4485] "ADNP" "ADNP" "ADNP" "AEBP2" ...
-#' # ..$ target: chr [1:4485] "PRF1" "TMEM40" "TNFRSF1B" "CFP" ...
-#' # ..$ mor   : num [1:4485] 1 1 1 1 1 1 1 1 1 1 ...
-#' # $ ppi:'data.frame':	12754 obs. of  3 variables:
-#' #   ..$ X.node1       : chr [1:12754] "ADNP" "ADNP" "ADNP" "AEBP2" ...
-#' # ..$ node2         : chr [1:12754] "ZBTB14" "NFIA" "CDC5L" "YY1" ...
-#' # ..$ combined_score: num [1:12754] 0.769 0.64 0.581 0.597 0.54 0.753 0.659 0.548 0.59 0.654 ...
 NULL
diff --git a/R/enrichEdges.R b/R/enrichEdges.R
new file mode 100644
index 0000000..3dec935
--- /dev/null
+++ b/R/enrichEdges.R
@@ -0,0 +1,325 @@
+#' @title Gene set enrichment analysis of TF-target edges
+#' @description
+#' Performs gene set enrichment analysis separately for each transcription
+#' factor (TF) using the edge-level values supplied in \code{numericValue}.
+#' Enrichment is performed with the multilevel implementation of \pkg{fgsea}.
+#' Calculations for individual TFs are performed in parallel.
+#' @author Daniel Osorio <daniecos@uio.no>
+#' @param edgesDF A data.frame of TF-target edges, typically produced by
+#'   \code{\link{testEdges}}. Must contain a \code{tf} column, a \code{target}
+#'   column, and the numeric column named by \code{numericValue}.
+#' @param geneSets A named list of gene sets. The names of the list elements are
+#'   used as gene set identifiers.
+#' @param numericValue Character string naming the column in \code{edgesDF} used
+#'   as the ranking statistic for enrichment analysis.
+#' @param nCores Integer specifying the number of parallel workers to use.
+#'   Default 3.
+#' @param seed Integer specifying the random seed used by the parallel
+#'   enrichment calculations. Default 1.
+#' @return A data.frame of enrichment results with one row per TF-gene set pair:
+#'   \itemize{
+#'     \item{tf: Transcription factor}
+#'     \item{geneSet: Gene set identifier}
+#'     \item{pValue: Raw enrichment p-value}
+#'     \item{pAdj: Benjamini-Hochberg adjusted p-value}
+#'     \item{log2Err: Expected log2 error of the p-value estimate}
+#'     \item{ES: Enrichment score}
+#'     \item{NES: Normalized enrichment score}
+#'     \item{geneSetSize: Number of genes from the set found among the targets}
+#'   }
+#' @details
+#' For each TF, the values in \code{numericValue} are used as ranked statistics
+#' for its target genes. Edges with missing targets or missing or non-finite
+#' values in the selected numeric column are excluded before enrichment
+#' analysis. If a target occurs more than once for a TF, only the observation
+#' with the largest absolute value of the selected ranking statistic is
+#' retained.
+#'
+#' Gene set enrichment is performed using \code{fgsea::fgseaMultilevel}; the
+#' \pkg{fgsea} package (Bioconductor) is required. The \code{leadingEdge} column
+#' returned by \code{fgseaMultilevel} is not included in the output. P-values
+#' are adjusted across all TF-gene set enrichment tests using the
+#' Benjamini-Hochberg procedure.
+#' @seealso \code{\link{testEdges}}, \code{\link{maEdges}}
+#' @examples
+#' \dontrun{
+#' data(scorpionTest)
+#' nets <- runSCORPION(
+#'   gexMatrix = scorpionTest$gex,
+#'   tfMotifs = scorpionTest$tf,
+#'   ppiNet = scorpionTest$ppi,
+#'   cellsMetadata = scorpionTest$metadata,
+#'   groupBy = c("donor", "region")
+#' )
+#' res <- testEdges(
+#'   networksDF = nets,
+#'   testType = "two.sample",
+#'   group1 = grep("--T$", colnames(nets), value = TRUE),
+#'   group2 = grep("--N$", colnames(nets), value = TRUE)
+#' )
+#'
+#' geneSets <- list(SetA = c("ACKR1", "ACTA2"), SetB = c("ACTG2", "ADAMDEC1"))
+#' enr <- enrichEdges(
+#'   edgesDF = res,
+#'   geneSets = geneSets,
+#'   numericValue = "log2FoldChange"
+#' )
+#' }
+#' @export
+#' @importFrom stats p.adjust
+enrichEdges <- function(edgesDF,
+                        geneSets,
+                        numericValue,
+                        nCores = 3,
+                        seed = 1) {
+  if (!requireNamespace("fgsea", quietly = TRUE)) {
+    stop(
+      "Package 'fgsea' (Bioconductor) is required for enrichEdges(). ",
+      "Install it with BiocManager::install('fgsea').",
+      call. = FALSE
+    )
+  }
+  # ============================================================
+  # Check input
+  # ============================================================
+
+  stopifnot(
+    is.data.frame(edgesDF),
+    is.list(geneSets),
+    length(geneSets) > 0L,
+    is.character(numericValue),
+    length(numericValue) == 1L,
+    is.numeric(nCores),
+    length(nCores) == 1L,
+    is.finite(nCores),
+    nCores >= 1,
+    nCores == as.integer(nCores),
+    is.numeric(seed),
+    length(seed) == 1L,
+    is.finite(seed),
+    seed == as.integer(seed)
+  )
+
+  required_columns <- c("tf", "target", numericValue)
+
+  missing_columns <- setdiff(required_columns, colnames(edgesDF))
+
+  if (length(missing_columns)) {
+    stop(
+      "Missing required column(s): ",
+      paste(missing_columns, collapse = ", "),
+      call. = FALSE
+    )
+  }
+
+  if (!is.numeric(edgesDF[[numericValue]])) {
+    stop("`numericValue` must specify a numeric column.", call. = FALSE)
+  }
+
+  if (is.null(names(geneSets)) ||
+      any(!nzchar(names(geneSets)))) {
+    stop("`geneSets` must be a named list.", call. = FALSE)
+  }
+
+  # ============================================================
+  # TFs
+  # ============================================================
+
+  tf_list <- unique(edgesDF$tf[!is.na(edgesDF$tf)])
+
+  # ============================================================
+  # Enrichment for one TF
+  # ============================================================
+
+  enrich_tf_edges <- function(selected_tf) {
+    tf_edges <- edgesDF[
+      !is.na(edgesDF$tf) & edgesDF$tf == selected_tf,
+      ,
+      drop = FALSE
+    ]
+
+    # ----------------------------------------------------------
+    # Keep only edges with valid ranking values and targets.
+    # ----------------------------------------------------------
+
+    valid <- (
+      !is.na(tf_edges$target) &
+        is.finite(tf_edges[[numericValue]])
+    )
+
+    tf_edges <- tf_edges[valid, , drop = FALSE]
+
+    if (!nrow(tf_edges)) {
+      return(NULL)
+    }
+
+    # ----------------------------------------------------------
+    # Construct ranked statistics.
+    #
+    # If a target occurs more than once for a TF, retain the
+    # observation with the largest absolute statistic.
+    # ----------------------------------------------------------
+
+    tf_edges <- tf_edges[
+      order(
+        abs(tf_edges[[numericValue]]),
+        decreasing = TRUE
+      ),
+      ,
+      drop = FALSE
+    ]
+
+    tf_edges <- tf_edges[
+      !duplicated(tf_edges$target),
+      ,
+      drop = FALSE
+    ]
+
+    weights <- tf_edges[[numericValue]]
+
+    names(weights) <- tf_edges$target
+
+    # ----------------------------------------------------------
+    # fgsea
+    # ----------------------------------------------------------
+
+    enrichment <- suppressWarnings(
+      fgsea::fgseaMultilevel(
+        pathways = geneSets,
+        stats = weights
+      )
+    )
+
+    if (!nrow(enrichment)) {
+      return(NULL)
+    }
+
+    enrichment <- as.data.frame(enrichment)
+
+    # ----------------------------------------------------------
+    # Keep the relevant output columns.
+    # fgsea returns padj, which is replaced below by the global
+    # adjustment across all TF-gene set tests.
+    # ----------------------------------------------------------
+
+    enrichment <- enrichment[
+      setdiff(colnames(enrichment), "leadingEdge"),
+      drop = FALSE
+    ]
+
+    enrichment <- enrichment[
+      ,
+      c(
+        "pathway",
+        "pval",
+        "log2err",
+        "ES",
+        "NES",
+        "size"
+      ),
+      drop = FALSE
+    ]
+
+    colnames(enrichment) <- c(
+      "geneSet",
+      "pValue",
+      "log2Err",
+      "ES",
+      "NES",
+      "geneSetSize"
+    )
+
+    enrichment$tf <- selected_tf
+
+    enrichment <- enrichment[
+      ,
+      c(
+        "tf",
+        "geneSet",
+        "pValue",
+        "log2Err",
+        "ES",
+        "NES",
+        "geneSetSize"
+      ),
+      drop = FALSE
+    ]
+
+    enrichment
+  }
+
+  # ============================================================
+  # Parallel enrichment
+  # ============================================================
+
+  old_plan <- future::plan()
+
+  on.exit(
+    future::plan(old_plan),
+    add = TRUE
+  )
+
+  future::plan(
+    future::multisession,
+    workers = nCores
+  )
+
+  enrichment <- furrr::future_map_dfr(
+    tf_list,
+    enrich_tf_edges,
+    .progress = TRUE,
+    .options = furrr::furrr_options(seed = seed)
+  )
+
+  # ============================================================
+  # Handle case where no TF has valid information
+  # ============================================================
+
+  if (!nrow(enrichment)) {
+    return(
+      data.frame(
+        tf = character(),
+        geneSet = character(),
+        pValue = numeric(),
+        pAdj = numeric(),
+        log2Err = numeric(),
+        ES = numeric(),
+        NES = numeric(),
+        geneSetSize = integer(),
+        stringsAsFactors = FALSE
+      )
+    )
+  }
+
+  # ============================================================
+  # BH correction across all TF-gene set tests
+  # ============================================================
+
+  enrichment$pAdj <- p.adjust(
+    enrichment$pValue,
+    method = "BH"
+  )
+
+  # ============================================================
+  # Final column order
+  # ============================================================
+
+  enrichment <- enrichment[
+    ,
+    c(
+      "tf",
+      "geneSet",
+      "pValue",
+      "pAdj",
+      "log2Err",
+      "ES",
+      "NES",
+      "geneSetSize"
+    ),
+    drop = FALSE
+  ]
+
+  rownames(enrichment) <- NULL
+
+  enrichment
+}
diff --git a/R/maEdges.R b/R/maEdges.R
new file mode 100644
index 0000000..c466bdf
--- /dev/null
+++ b/R/maEdges.R
@@ -0,0 +1,687 @@
+#' @title Meta-analysis of TF-target edges across studies
+#' @description
+#' Performs a meta-analysis of TF-target edges across multiple studies using
+#' either a fixed-effect or DerSimonian-Laird random-effects model. Missing or
+#' non-finite effect sizes and p-values are excluded from the corresponding
+#' study. A TF-target pair is only counted as contributing to a study when both
+#' its effect size and p-value are valid.
+#' @author Daniel Osorio <daniecos@uio.no>
+#' @param edgesList A list of data.frames, one per study, typically produced by
+#'   \code{\link{testEdges}}. Each data.frame must contain the columns
+#'   \code{tf}, \code{target}, \code{log2FoldChange} and \code{pValue}.
+#' @param method Meta-analysis model. Either \code{"random"} (DerSimonian-Laird
+#'   random-effects) or \code{"fixed"} (inverse-variance fixed-effect). Default
+#'   \code{"random"}.
+#' @param minStudies Minimum number of studies with valid numeric information
+#'   required for a TF-target pair to be included. Default 2.
+#' @return A data.frame of meta-analysis results with one row per TF-target pair:
+#'   \itemize{
+#'     \item{tf: Transcription factor}
+#'     \item{target: Target gene}
+#'     \item{k: Number of contributing studies}
+#'     \item{log2FoldChange: Meta-analytic effect size}
+#'     \item{se: Standard error of the effect size}
+#'     \item{ciLow, ciHigh: Lower and upper bounds of the 95\% confidence interval}
+#'     \item{zStatistic: Z statistic}
+#'     \item{pValue: Raw p-value}
+#'     \item{Q: Cochran's Q heterogeneity statistic}
+#'     \item{iSquared: I-squared heterogeneity (percentage)}
+#'     \item{tauSquared: Between-study variance estimate}
+#'     \item{pAdj: Benjamini-Hochberg adjusted p-value}
+#'   }
+#' @details
+#' Standard errors are estimated from the reported log2 fold change and
+#' two-sided p-value. Fixed-effect estimates use inverse-variance weighting.
+#' Random-effects estimates use DerSimonian-Laird estimates of between-study
+#' variance.
+#'
+#' Missing or non-finite values in \code{log2FoldChange} or \code{pValue} do not
+#' contribute to the meta-analysis and do not count toward \code{minStudies}.
+#' TF-target pairs with fewer than \code{minStudies} valid studies are excluded.
+#' Benjamini-Hochberg correction is performed once across all TF-target
+#' meta-analysis tests.
+#' @seealso \code{\link{testEdges}}, \code{\link{enrichEdges}}
+#' @examples
+#' \dontrun{
+#' data(scorpionTest)
+#' nets <- runSCORPION(
+#'   gexMatrix = scorpionTest$gex,
+#'   tfMotifs = scorpionTest$tf,
+#'   ppiNet = scorpionTest$ppi,
+#'   cellsMetadata = scorpionTest$metadata,
+#'   groupBy = c("donor", "region")
+#' )
+#'
+#' # Two independent comparisons treated as separate studies
+#' studyA <- testEdges(nets, "two.sample",
+#'   group1 = "P31--T", group2 = "P31--N")
+#' studyB <- testEdges(nets, "two.sample",
+#'   group1 = "P32--T", group2 = "P32--N")
+#'
+#' meta <- maEdges(list(studyA, studyB), method = "random")
+#' }
+#' @export
+#' @importFrom stats p.adjust qnorm pnorm
+maEdges <- function(edgesList,
+                    method = c("random", "fixed"),
+                    minStudies = 2L) {
+  method <- match.arg(method)
+
+  n_studies <- length(edgesList)
+
+  stopifnot(is.list(edgesList), n_studies >= 2L)
+
+  # ============================================================
+  # Build the UNION of TFs across all studies.
+  # ============================================================
+
+  tfs <- unique(unlist(lapply(edgesList, function(x) {
+    unique(x$tf[!is.na(x$tf)])
+  }), use.names = FALSE))
+
+  n_tf <- length(tfs)
+
+  cli_alert_info("TFs to process: {format(n_tf, big.mark = ',')}")
+
+  # ============================================================
+  # Output list
+  # ============================================================
+
+  result_list <- vector("list", n_tf)
+
+  # ============================================================
+  # Progress
+  # ============================================================
+
+  pb <- cli_progress_bar("TF-level meta-analysis",
+                         total = n_tf,
+                         clear = FALSE)
+
+  # ============================================================
+  # Process one TF at a time
+  # ============================================================
+
+  for (tf_i in seq_along(tfs)) {
+    tf <- tfs[tf_i]
+
+    # --------------------------------------------------------
+    # Extract this TF from each study.
+    #
+    # Numeric NAs are retained here so that the target universe
+    # remains complete. They are filtered before contributing to
+    # study counts or meta-analysis estimates.
+    # --------------------------------------------------------
+
+    study_data <- vector("list", n_studies)
+
+    for (s in seq_len(n_studies)) {
+      x <- edgesList[[s]]
+
+      idx_tf <-
+        !is.na(x$tf) &
+        x$tf == tf
+
+      if (!any(idx_tf)) {
+        next
+      }
+
+      study_data[[s]] <- data.frame(
+        target =
+          x$target[idx_tf],
+        effect =
+          x$log2FoldChange[idx_tf],
+        p =
+          x$pValue[idx_tf],
+        stringsAsFactors = FALSE
+      )
+    }
+
+    # --------------------------------------------------------
+    # Get target universe for this TF.
+    # --------------------------------------------------------
+
+    target_union <- unique(unlist(lapply(study_data, function(x) {
+      if (is.null(x)) {
+        NULL
+      } else {
+        x$target
+      }
+    }), use.names = FALSE))
+
+    target_union <- target_union[!is.na(target_union)]
+
+    if (!length(target_union)) {
+      cli_progress_update(id = pb, status = sprintf("TF %d/%d", tf_i, n_tf))
+
+      next
+    }
+
+    n_targets <- length(target_union)
+
+    # ========================================================
+    # PASS 1 ACCUMULATORS
+    # ========================================================
+
+    k <- integer(n_targets)
+
+    sum_w <- numeric(n_targets)
+    sum_wy <- numeric(n_targets)
+    sum_wy2 <- numeric(n_targets)
+    sum_w2 <- numeric(n_targets)
+
+    # ========================================================
+    # PASS 1
+    # ========================================================
+
+    for (s in seq_len(n_studies)) {
+      x <- study_data[[s]]
+
+      if (is.null(x)) {
+        next
+      }
+
+      idx <- match(x$target, target_union)
+
+      ok <- !is.na(idx)
+
+      if (!any(ok)) {
+        next
+      }
+
+      idx <- idx[ok]
+
+      effect <- x$effect[ok]
+      p <- x$p[ok]
+
+      # ----------------------------------------------------
+      # A study contributes only if BOTH numeric quantities
+      # are valid.
+      # ----------------------------------------------------
+
+      valid_input <-
+        is.finite(effect) &
+        is.finite(p) &
+        p > 0 &
+        p <= 1
+
+      if (!any(valid_input)) {
+        next
+      }
+
+      idx <- idx[valid_input]
+      effect <- effect[valid_input]
+      p <- p[valid_input]
+
+      # ----------------------------------------------------
+      # Convert p-value to z
+      # ----------------------------------------------------
+
+      p <- pmin(pmax(p, .Machine$double.xmin), 1 - .Machine$double.eps)
+
+      z <- qnorm(1 - p / 2)
+
+      # ----------------------------------------------------
+      # SE from p-value + log2FC
+      # ----------------------------------------------------
+
+      se <-
+        abs(effect) / z
+
+      valid <-
+        is.finite(effect) &
+        is.finite(z) &
+        z > 0 &
+        is.finite(se) &
+        se > 0
+
+      if (!any(valid)) {
+        next
+      }
+
+      idx <- idx[valid]
+      effect <- effect[valid]
+      se <- se[valid]
+
+      w <- 1 / se^2
+
+      valid_w <-
+        is.finite(w) &
+        w > 0
+
+      if (!any(valid_w)) {
+        next
+      }
+
+      idx <- idx[valid_w]
+      effect <- effect[valid_w]
+      w <- w[valid_w]
+
+      # ----------------------------------------------------
+      # Accumulate
+      # ----------------------------------------------------
+
+      k[idx] <-
+        k[idx] + 1L
+
+      sum_w[idx] <-
+        sum_w[idx] + w
+
+      sum_wy[idx] <-
+        sum_wy[idx] +
+        w * effect
+
+      sum_wy2[idx] <-
+        sum_wy2[idx] +
+        w * effect^2
+
+      sum_w2[idx] <-
+        sum_w2[idx] +
+        w^2
+    }
+
+    # ========================================================
+    # Keep targets with enough VALID studies
+    # ========================================================
+
+    valid_meta <-
+      k >= minStudies &
+      is.finite(sum_w) &
+      sum_w > 0
+
+    if (!any(valid_meta)) {
+      cli_progress_update(id = pb, status = sprintf("TF %d/%d", tf_i, n_tf))
+
+      next
+    }
+
+    # --------------------------------------------------------
+    # Restrict all vectors to meta-analyzed targets
+    # --------------------------------------------------------
+
+    target_meta <-
+      target_union[valid_meta]
+
+    k <-
+      k[valid_meta]
+
+    sum_w <-
+      sum_w[valid_meta]
+
+    sum_wy <-
+      sum_wy[valid_meta]
+
+    sum_wy2 <-
+      sum_wy2[valid_meta]
+
+    sum_w2 <-
+      sum_w2[valid_meta]
+
+    n_meta <-
+      length(target_meta)
+
+    # ========================================================
+    # Fixed effect
+    # ========================================================
+
+    fixed_log2FC <-
+      sum_wy / sum_w
+
+    fixed_SE <-
+      sqrt(1 / sum_w)
+
+    # ========================================================
+    # Cochran Q
+    # ========================================================
+
+    Q <-
+      pmax(0, sum_wy2 -
+             sum_wy^2 / sum_w)
+
+    df <-
+      k - 1L
+
+    # ========================================================
+    # DerSimonian-Laird tau²
+    # ========================================================
+
+    C <-
+      sum_w -
+      sum_w2 / sum_w
+
+    tau2 <-
+      numeric(n_meta)
+
+    tau_ok <-
+      k > 1L &
+      is.finite(C) &
+      C > 0
+
+    tau2[tau_ok] <-
+      pmax(0, (Q[tau_ok] - df[tau_ok]) /
+             C[tau_ok])
+
+    # ========================================================
+    # PASS 2: random effects
+    # ========================================================
+
+    if (method == "random") {
+      sum_wr <-
+        numeric(n_meta)
+
+      sum_wry <-
+        numeric(n_meta)
+
+      for (s in seq_len(n_studies)) {
+        x <- study_data[[s]]
+
+        if (is.null(x)) {
+          next
+        }
+
+        idx <- match(x$target, target_meta)
+
+        ok <- !is.na(idx)
+
+        if (!any(ok)) {
+          next
+        }
+
+        idx <- idx[ok]
+
+        effect <- x$effect[ok]
+        p <- x$p[ok]
+
+        # ----------------------------------------------------
+        # Same validity criteria as PASS 1.
+        # ----------------------------------------------------
+
+        valid_input <-
+          is.finite(effect) &
+          is.finite(p) &
+          p > 0 &
+          p <= 1
+
+        if (!any(valid_input)) {
+          next
+        }
+
+        idx <- idx[valid_input]
+        effect <- effect[valid_input]
+        p <- p[valid_input]
+
+        p <- pmin(pmax(p, .Machine$double.xmin), 1 - .Machine$double.eps)
+
+        z <- qnorm(1 - p / 2)
+
+        se <-
+          abs(effect) / z
+
+        valid <-
+          is.finite(effect) &
+          is.finite(z) &
+          z > 0 &
+          is.finite(se) &
+          se > 0
+
+        if (!any(valid)) {
+          next
+        }
+
+        idx <- idx[valid]
+        effect <- effect[valid]
+        se <- se[valid]
+
+        wr <-
+          1 /
+          (se^2 +
+             tau2[idx])
+
+        valid_w <-
+          is.finite(wr) &
+          wr > 0
+
+        if (!any(valid_w)) {
+          next
+        }
+
+        idx <- idx[valid_w]
+        effect <- effect[valid_w]
+        wr <- wr[valid_w]
+
+        sum_wr[idx] <-
+          sum_wr[idx] + wr
+
+        sum_wry[idx] <-
+          sum_wry[idx] +
+          wr * effect
+      }
+
+      # ------------------------------------------------------
+      # Random-effects estimate
+      # ------------------------------------------------------
+
+      meta_log2FC <-
+        numeric(n_meta)
+
+      meta_SE <-
+        numeric(n_meta)
+
+      random_ok <-
+        is.finite(sum_wr) &
+        sum_wr > 0
+
+      meta_log2FC[random_ok] <-
+        sum_wry[random_ok] /
+        sum_wr[random_ok]
+
+      meta_SE[random_ok] <-
+        sqrt(1 / sum_wr[random_ok])
+
+    } else {
+      meta_log2FC <-
+        fixed_log2FC
+
+      meta_SE <-
+        fixed_SE
+    }
+
+    # ========================================================
+    # Meta-analysis statistics
+    # ========================================================
+
+    valid_effect <-
+      is.finite(meta_log2FC) &
+      is.finite(meta_SE) &
+      meta_SE > 0
+
+    if (!any(valid_effect)) {
+      cli_progress_update(id = pb, status = sprintf("TF %d/%d", tf_i, n_tf))
+
+      next
+    }
+
+    meta_z <-
+      numeric(n_meta)
+
+    meta_z[valid_effect] <-
+      meta_log2FC[valid_effect] /
+      meta_SE[valid_effect]
+
+    pValue <-
+      numeric(n_meta)
+
+    pValue[valid_effect] <-
+      2 *
+      pnorm(-abs(meta_z[valid_effect]))
+
+    # --------------------------------------------------------
+    # Confidence intervals
+    # --------------------------------------------------------
+
+    CI_low <-
+      numeric(n_meta)
+
+    CI_high <-
+      numeric(n_meta)
+
+    CI_low[valid_effect] <-
+      meta_log2FC[valid_effect] -
+      qnorm(.975) *
+      meta_SE[valid_effect]
+
+    CI_high[valid_effect] <-
+      meta_log2FC[valid_effect] +
+      qnorm(.975) *
+      meta_SE[valid_effect]
+
+    # --------------------------------------------------------
+    # I²
+    # --------------------------------------------------------
+
+    I2 <-
+      numeric(n_meta)
+
+    i2_ok <-
+      is.finite(Q) &
+      Q > 0 &
+      df > 0
+
+    I2[i2_ok] <-
+      pmax(0, (Q[i2_ok] - df[i2_ok]) /
+             Q[i2_ok]) * 100
+
+    # ========================================================
+    # Output
+    # ========================================================
+
+    out <- data.frame(
+      tf =
+        rep(tf, sum(valid_effect)),
+      target =
+        target_meta[valid_effect],
+      k =
+        k[valid_effect],
+      log2FoldChange =
+        meta_log2FC[valid_effect],
+      se =
+        meta_SE[valid_effect],
+      ciLow =
+        CI_low[valid_effect],
+      ciHigh =
+        CI_high[valid_effect],
+      zStatistic =
+        meta_z[valid_effect],
+      pValue =
+        pValue[valid_effect],
+      Q =
+        Q[valid_effect],
+      iSquared =
+        I2[valid_effect],
+      tauSquared =
+        tau2[valid_effect],
+      stringsAsFactors = FALSE
+    )
+
+    # --------------------------------------------------------
+    # Defensive final filter
+    # --------------------------------------------------------
+
+    out <- out[is.finite(out$log2FoldChange) &
+                 is.finite(out$se) &
+                 is.finite(out$pValue), , drop = FALSE]
+
+    if (nrow(out) > 0L) {
+      result_list[[tf_i]] <- out
+    }
+
+    # --------------------------------------------------------
+    # Cleanup
+    # --------------------------------------------------------
+
+    rm(
+      study_data,
+      target_union,
+      target_meta,
+      k,
+      sum_w,
+      sum_wy,
+      sum_wy2,
+      sum_w2,
+      fixed_log2FC,
+      fixed_SE,
+      Q,
+      df,
+      C,
+      tau2,
+      meta_log2FC,
+      meta_SE,
+      meta_z,
+      pValue,
+      CI_low,
+      CI_high,
+      I2,
+      out
+    )
+
+    gc(FALSE)
+
+    cli_progress_update(id = pb, status = sprintf("TF %d/%d", tf_i, n_tf))
+  }
+
+  cli_progress_done(id = pb)
+
+  # ============================================================
+  # Combine
+  # ============================================================
+
+  cli_alert_info("Combining results")
+
+  valid_results <- result_list[vapply(result_list, function(x) {
+    !is.null(x) &&
+      nrow(x) > 0L
+  }, logical(1))]
+
+  if (!length(valid_results)) {
+    return(
+      data.frame(
+        tf = character(),
+        target = character(),
+        k = integer(),
+        log2FoldChange = numeric(),
+        se = numeric(),
+        ciLow = numeric(),
+        ciHigh = numeric(),
+        zStatistic = numeric(),
+        pValue = numeric(),
+        Q = numeric(),
+        iSquared = numeric(),
+        tauSquared = numeric(),
+        pAdj = numeric(),
+        stringsAsFactors = FALSE
+      )
+    )
+  }
+
+  result <- do.call(rbind, valid_results)
+
+  rownames(result) <- NULL
+
+  # ============================================================
+  # BH correction
+  # ============================================================
+
+  result$pAdj <-
+    p.adjust(result$pValue, method = "BH")
+
+  # ============================================================
+  # Sort safely
+  # ============================================================
+
+  result <- result[order(result$tf, result$target, na.last = TRUE), , drop = FALSE]
+
+  rownames(result) <- NULL
+
+  gc()
+
+  return(result)
+}
diff --git a/R/runSCORPION.R b/R/runSCORPION.R
index a495769..835f536 100644
--- a/R/runSCORPION.R
+++ b/R/runSCORPION.R
@@ -1,5 +1,6 @@
 #' @title Run SCORPION across cell groups and return combined networks
 #' @description Builds per-group regulatory networks by running \code{\link{scorpion}} on subsets of cells defined by \code{cellsMetadata} and combining the resulting networks into a wide-format data frame where each column corresponds to a network.
+#' @author Daniel Osorio <daniecos@uio.no>
 #' @param gexMatrix An expression dataset with genes in the rows and barcodes (cells) in the columns.
 #' @param tfMotifs A motif dataset, a data.frame or a matrix containing 3 columns. Each row describes a motif associated with a transcription factor (column 1) a gene (column 2) and a score (column 3).
 #' @param ppiNet A Protein-Protein-Interaction dataset, a data.frame or matrix containing 3 columns. Each row describes a protein-protein interaction between transcription factor 1 (column 1), transcription factor 2 (column 2) and a score (column 3).
@@ -24,6 +25,7 @@
 #' @param scaleByPresent Boolean to indicate scaling of correlations by percentage of positive samples. Default FALSE.
 #' @param filterExpr Boolean to indicate whether or not to remove genes with 0 expression across all cells. Default FALSE.
 #' @return A data.frame in wide format where rows represent TF-target pairs (union across all networks) and columns represent network identifiers. Cell values are edge weights from the corresponding network.
+#' @seealso \code{\link{scorpion}}, \code{\link{testEdges}}, \code{\link{regressEdges}}
 #' @details
 #' This function is a wrapper around \code{\link{scorpion}} that groups cells according to metadata columns, filters out groups with insufficient cells, runs network inference on each remaining group independently, and finally combines all resulting networks into a single wide-format data frame.
 #' @examples
diff --git a/R/testEdges.R b/R/testEdges.R
index 22f9338..955106b 100644
--- a/R/testEdges.R
+++ b/R/testEdges.R
@@ -4,6 +4,7 @@
 #' @description Performs statistical testing of network edges from runSCORPION output.
 #' Supports single-sample tests (testing if edges differ from zero) and two-sample
 #' tests (comparing edges between two groups).
+#' @author Daniel Osorio <daniecos@uio.no>
 #' @param networksDF A data.frame output from \code{\link{runSCORPION}} containing
 #'   TF-target pairs as rows and network identifiers as columns.
 #' @param testType Character specifying the test type. Options are:
@@ -57,6 +58,7 @@
 #'     \item{pAdj: Adjusted p-value}
 #'     \item{For two-sample tests: meanGroup1, meanGroup2, cohensD, log2FoldChange (Group1 - Group2)}
 #'   }
+#' @seealso \code{\link{runSCORPION}}, \code{\link{regressEdges}}, \code{\link{maEdges}}, \code{\link{circosEdges}}
 #' @details
 #' For single-sample tests, the function tests whether the mean edge weight across
 #' replicates significantly differs from zero using a one-sample t-test.
@@ -644,6 +646,7 @@ testEdgesPaired <- function(networksDF, group1, group2, alternative, minLog2FC,
 #' @description Performs linear regression on network edges from runSCORPION output
 #' to identify edges that show significant trends across ordered conditions (e.g.,
 #' disease progression: Normal -> Border -> Tumor).
+#' @author Daniel Osorio <daniecos@uio.no>
 #' @param networksDF A data.frame output from \code{\link{runSCORPION}} containing
 #'   TF-target pairs as rows and network identifiers as columns.
 #' @param orderedGroups A named list where each element is a character vector of
@@ -670,6 +673,7 @@ testEdgesPaired <- function(networksDF, group1, group2, alternative, minLog2FC,
 #'     \item{meanEdge: Overall mean edge weight across all conditions}
 #'     \item{One column per condition showing mean edge weight in that condition}
 #'   }
+#' @seealso \code{\link{runSCORPION}}, \code{\link{testEdges}}
 #' @details
 #' This function performs simple linear regression for each edge, modeling edge weight
 #' as a function of an ordered categorical variable (coded as 0, 1, 2, ... for each
diff --git a/README.md b/README.md
index bfcb7d0..b321e1e 100644
--- a/README.md
+++ b/README.md
@@ -415,6 +415,250 @@ decreasing <- results_reg[results_reg$pAdj < 0.05 & results_reg$slope < 0, ]
 
 ---
 
+### maEdges
+
+Combines differential edge results from several independent studies into a
+single meta-analytic estimate per TF-target pair, using either a fixed-effect
+or DerSimonian-Laird random-effects model. Standard errors are recovered from
+each study's `log2FoldChange` and `pValue`.
+
+**Usage:**
+
+```r
+results <- maEdges(
+  edgesList,
+  method = c("random", "fixed"),
+  minStudies = 2
+)
+```
+
+**Parameters:**
+
+| Parameter | Description | Default |
+|-----------|-------------|---------|
+| `edgesList` | List of `testEdges()` result data frames, one per study | Required |
+| `method` | Meta-analysis model: `random` or `fixed` | `random` |
+| `minStudies` | Minimum studies with valid values required per TF-target pair | 2 |
+
+**Return value:**
+
+A data frame containing:
+
+| Column | Description |
+|--------|-------------|
+| `tf`, `target` | TF-target pair identifiers |
+| `k` | Number of contributing studies |
+| `log2FoldChange` | Meta-analytic effect size |
+| `se` | Standard error of the effect size |
+| `ciLow`, `ciHigh` | 95% confidence interval bounds |
+| `zStatistic` | Z statistic |
+| `pValue` | Raw p-value |
+| `Q` | Cochran's Q heterogeneity statistic |
+| `iSquared` | I-squared heterogeneity (%) |
+| `tauSquared` | Between-study variance estimate |
+| `pAdj` | Adjusted p-value (Benjamini-Hochberg) |
+
+**Example:**
+
+```r
+# Treat each patient's Tumor vs Normal comparison as a study
+studyA <- testEdges(nets, "two.sample", group1 = "P31--T", group2 = "P31--N")
+studyB <- testEdges(nets, "two.sample", group1 = "P32--T", group2 = "P32--N")
+
+meta <- maEdges(list(studyA, studyB), method = "random")
+```
+
+---
+
+## Functional Enrichment
+
+### enrichEdges
+
+Runs gene set enrichment analysis separately for each transcription factor,
+ranking its target genes by a chosen edge-level statistic (e.g.
+`log2FoldChange`) and testing enrichment with the multilevel `fgsea` algorithm.
+TFs are processed in parallel.
+
+Requires the `fgsea` package (Bioconductor):
+
+```r
+BiocManager::install("fgsea")
+```
+
+**Usage:**
+
+```r
+results <- enrichEdges(
+  edgesDF,
+  geneSets,
+  numericValue,
+  nCores = 3,
+  seed = 1
+)
+```
+
+**Parameters:**
+
+| Parameter | Description | Default |
+|-----------|-------------|---------|
+| `edgesDF` | Output from `testEdges()` with `tf`, `target` and `numericValue` columns | Required |
+| `geneSets` | Named list of gene sets | Required |
+| `numericValue` | Name of the `edgesDF` column used as the ranking statistic | Required |
+| `nCores` | Number of parallel workers | 3 |
+| `seed` | Random seed for the parallel enrichment calculations | 1 |
+
+**Return value:**
+
+A data frame containing:
+
+| Column | Description |
+|--------|-------------|
+| `tf` | Transcription factor |
+| `geneSet` | Gene set identifier |
+| `pValue` | Raw enrichment p-value |
+| `pAdj` | Adjusted p-value (Benjamini-Hochberg) |
+| `log2Err` | Expected log2 error of the p-value estimate |
+| `ES` | Enrichment score |
+| `NES` | Normalized enrichment score |
+| `geneSetSize` | Number of set genes found among the targets |
+
+**Example:**
+
+```r
+res <- testEdges(
+  networksDF = nets,
+  testType = "two.sample",
+  group1 = grep("--T$", colnames(nets), value = TRUE),
+  group2 = grep("--N$", colnames(nets), value = TRUE)
+)
+
+geneSets <- list(
+  SetA = c("ACKR1", "ACTA2"),
+  SetB = c("ACTG2", "ADAMDEC1")
+)
+
+enr <- enrichEdges(
+  edgesDF = res,
+  geneSets = geneSets,
+  numericValue = "log2FoldChange"
+)
+```
+
+---
+
+## Visualization
+
+### circosEdges
+
+Draws a circular (Circos) plot of differential TF–target links from a
+`testEdges()` two-sample result. Genes are placed on their genomic coordinates;
+link ribbons are colored continuously by `log2FoldChange` (diverging scale
+centered at zero) and their thickness is proportional to `-log10(pAdj)`. A
+track of per-gene degree (the summed `log2FoldChange` of each gene's outgoing /
+incoming links) is drawn as needles. Links are flagged known vs novel against
+an optional a priori network, and gene labels are styled by role — **TFs in
+bold, targets in *italics*** — showing only the top TFs/targets by absolute
+degree to avoid overlap. By default only the main chromosomes are shown.
+
+Requires the `circlize` package, and `biomaRt` when gene coordinates are
+downloaded automatically:
+
+```r
+install.packages("circlize")
+BiocManager::install("biomaRt")
+```
+
+**Usage:**
+
+```r
+circosEdges(
+  edgesDF,
+  species = "hsapiens_gene_ensembl",
+  geneCoords = NULL,
+  priorNet = NULL,
+  geneSets = NULL,
+  pAdjThreshold = 0.05,
+  log2FCThreshold = 0,
+  maxEdges = 500,
+  colorBy = "log2FoldChange",
+  linkColors = c("#2166AC", "#F7F7F7", "#B2182B"),
+  lwdRange = c(0.5, 4),
+  nmaxTF = 20,
+  nmaxTarget = 20,
+  knownColor = "grey60",
+  novelColor = "#D95F02",
+  geneSetColors = NULL,
+  chromosomes = NULL,
+  mainChromosomesOnly = TRUE,
+  ensemblMirror = "www",
+  transparency = 0.5,
+  legend = TRUE
+)
+```
+
+**Parameters:**
+
+| Parameter | Description | Default |
+|-----------|-------------|---------|
+| `edgesDF` | Output from `testEdges()` (two-sample or paired); needs `tf`, `target`, `log2FoldChange`, `pAdj` | Required |
+| `species` | Ensembl dataset for `biomaRt` when `geneCoords` is `NULL` (e.g. `"mmusculus_gene_ensembl"`) | `"hsapiens_gene_ensembl"` |
+| `geneCoords` | Optional data frame (`gene`, `chr`, `start`, `end`) overriding the download for any source | `NULL` |
+| `priorNet` | Optional a priori TF–target network (first two columns); listed links are labeled *known*, others *novel* | `NULL` |
+| `geneSets` | Optional GMT file path or named list of gene vectors to label around the circle | `NULL` |
+| `pAdjThreshold` | Significance cutoff on `pAdj` (falls back to `pValue`) | `0.05` |
+| `log2FCThreshold` | Minimum absolute `log2FoldChange` to draw a link | `0` |
+| `maxEdges` | Cap on links drawn; keeps the most significant | `500` |
+| `colorBy` | `edgesDF` column mapped to the continuous link color ramp | `"log2FoldChange"` |
+| `linkColors` | Length-3 low/mid/high color ramp (diverging for signed values) | blue–white–red |
+| `lwdRange` | Min/max link line width; thickness scales with `-log10(pAdj)` | `c(0.5, 4)` |
+| `nmaxTF`, `nmaxTarget` | Number of TFs/targets to label, by largest absolute out-/in-degree (`NULL`/`Inf` = all) | `20`, `20` |
+| `knownColor`, `novelColor` | Border colors for known vs novel links | grey, orange |
+| `geneSetColors` | Optional named colors for gene sets | auto |
+| `chromosomes` | Optional subset/ordering of chromosomes shown | all present |
+| `mainChromosomesOnly` | Drop unplaced scaffolds/contigs, keeping numbered chromosomes plus X, Y, MT | `TRUE` |
+| `ensemblMirror` | `biomaRt` mirror: `"www"`, `"useast"`, `"asia"` | `"www"` |
+| `transparency` | Link transparency in `[0, 1]` (0 = opaque) | `0.5` |
+| `legend` | Draw legends for effect size, novelty and gene sets | `TRUE` |
+
+**Return value:**
+
+Invisibly, a list with `edges` (plotted links annotated with coordinates and
+`novelty`) and `coords` (the gene coordinate table with `outDegree`, `inDegree`
+and total `degree` columns). Called for the side effect of drawing the plot.
+
+**Example:**
+
+```r
+# Two-sample comparison: Tumor vs Normal
+res <- testEdges(
+  networksDF = nets,
+  testType = "two.sample",
+  group1 = grep("--T$", colnames(nets), value = TRUE),
+  group2 = grep("--N$", colnames(nets), value = TRUE)
+)
+
+# Human coordinates auto-downloaded from Ensembl,
+# links flagged known/novel against a prior network,
+# and hallmark gene sets labeled around the circle
+circosEdges(
+  edgesDF = res,
+  species = "hsapiens_gene_ensembl",
+  priorNet = scorpionTest$tf,
+  geneSets = "hallmark.gmt"
+)
+
+# Offline / other species: supply your own coordinates
+coords <- data.frame(
+  gene = c("TF1", "GENE1"),
+  chr = c("1", "2"),
+  start = c(1000, 5000),
+  end = c(2000, 6000)
+)
+circosEdges(res, geneCoords = coords)
+```
+
+---
+
 ## Citation
 
 If you use SCORPION in your research, please cite:
diff --git a/data/scorpionTest.RData b/data/scorpionTest.RData
index 36041cf..baef709 100644
Binary files a/data/scorpionTest.RData and b/data/scorpionTest.RData differ
diff --git a/man/circosEdges.Rd b/man/circosEdges.Rd
new file mode 100644
index 0000000..4ffc6f3
--- /dev/null
+++ b/man/circosEdges.Rd
@@ -0,0 +1,162 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circosEdges.R
+\name{circosEdges}
+\alias{circosEdges}
+\title{Circos plot of differential network edges}
+\usage{
+circosEdges(
+  edgesDF,
+  species = "hsapiens_gene_ensembl",
+  geneCoords = NULL,
+  priorNet = NULL,
+  geneSets = NULL,
+  pAdjThreshold = 0.05,
+  log2FCThreshold = 0,
+  maxEdges = 500L,
+  colorBy = "log2FoldChange",
+  linkColors = c("#2166AC", "#F7F7F7", "#B2182B"),
+  lwdRange = c(0.5, 4),
+  nmaxTF = 20L,
+  nmaxTarget = 20L,
+  knownColor = "grey60",
+  novelColor = "#D95F02",
+  geneSetColors = NULL,
+  chromosomes = NULL,
+  mainChromosomesOnly = TRUE,
+  ensemblMirror = "www",
+  transparency = 0.5,
+  hRatio = 0.6,
+  fontFamily = "sans",
+  legend = TRUE
+)
+}
+\arguments{
+\item{edgesDF}{A data.frame produced by \code{\link{testEdges}} (two-sample or
+paired). Must contain the columns \code{tf}, \code{target},
+\code{log2FoldChange} and \code{pAdj} (\code{pValue} is used as a fallback
+when \code{pAdj} is absent).}
+
+\item{species}{Ensembl dataset name passed to \pkg{biomaRt} when
+\code{geneCoords} is \code{NULL}, e.g. \code{"hsapiens_gene_ensembl"} for
+human or \code{"mmusculus_gene_ensembl"} for mouse. See
+\code{biomaRt::listDatasets()} for the full multi-species list.}
+
+\item{geneCoords}{Optional data.frame supplying gene coordinates from any
+source (overrides the \pkg{biomaRt} download). Must have columns
+\code{gene}, \code{chr}, \code{start} and \code{end}.}
+
+\item{priorNet}{Optional a priori TF-target network whose first two columns
+are the TF and target. Links present here are labelled \code{"known"}, all
+others \code{"novel"}. Accepts a data.frame or a matrix.}
+
+\item{geneSets}{Optional gene-set annotation used to label genes around the
+circle: either a path to a GMT file or a named list of character vectors.}
+
+\item{pAdjThreshold}{Numeric significance cutoff applied to \code{pAdj}
+(or \code{pValue} when \code{pAdj} is missing). Default \code{0.05}.}
+
+\item{log2FCThreshold}{Numeric minimum absolute \code{log2FoldChange} required
+to draw a link. Default \code{0}.}
+
+\item{maxEdges}{Integer cap on the number of links drawn; when exceeded, the
+most significant edges are kept. Default \code{500}.}
+
+\item{colorBy}{Name of the \code{edgesDF} column mapped to the link colour
+ramp. Default \code{"log2FoldChange"}, giving a continuous diverging colour
+scale centred at zero. When the default is used but \code{log2FoldChange} is
+absent (e.g. single-sample \code{\link{testEdges}} output), it falls back to
+\code{meanEdge}.}
+
+\item{linkColors}{Length-3 vector of colours for the low, mid and high ends of
+\code{colorBy}. Default blue-white-red; a diverging ramp is used when
+\code{colorBy} has negative values, otherwise a sequential low-to-high ramp.}
+
+\item{lwdRange}{Length-2 numeric giving the minimum and maximum link line
+width; each link's thickness is scaled linearly within this range by its
+\code{-log10} adjusted p-value. Default \code{c(0.5, 4)}.}
+
+\item{nmaxTF, nmaxTarget}{Integers giving how many TFs and targets to label,
+selected by the largest absolute out-degree and in-degree respectively.
+Use \code{NULL} or \code{Inf} to label all. Defaults \code{20}.}
+
+\item{knownColor, novelColor}{Border colours distinguishing known from novel
+links. Defaults grey and orange.}
+
+\item{geneSetColors}{Optional named vector mapping gene-set names to colours.
+When \code{NULL}, colours are generated automatically.}
+
+\item{chromosomes}{Optional character vector restricting and ordering the
+chromosomes shown. When \code{NULL}, all chromosomes present are used.}
+
+\item{mainChromosomesOnly}{Logical; when \code{TRUE} (the default) and
+\code{chromosomes} is \code{NULL}, only the main chromosomes (numbered,
+plus X, Y and MT) are kept and unplaced scaffolds/contigs are dropped.}
+
+\item{ensemblMirror}{\pkg{biomaRt} mirror to query: one of \code{"www"},
+\code{"useast"} or \code{"asia"}. Default \code{"www"}.}
+
+\item{transparency}{Numeric link transparency in \code{[0, 1]} (0 is opaque).
+Default \code{0.5}.}
+
+\item{hRatio}{Numeric in \code{[0, 1]} controlling how far link ribbons bend
+toward the circle centre; smaller values give flatter, less tangled links.
+Default \code{0.6}.}
+
+\item{fontFamily}{Font family used for all plot text, e.g. \code{"sans"}
+(Helvetica/Arial, the default) for a publication look.}
+
+\item{legend}{Logical; whether to draw legends for effect size, novelty and
+gene sets. Default \code{TRUE}.}
+}
+\value{
+Invisibly, a list with \code{edges} (the plotted links annotated with
+  coordinates and novelty) and \code{coords} (the gene coordinate table with
+  \code{outDegree}, \code{inDegree} and total \code{degree} columns, each the
+  sum of \code{log2FoldChange} over a gene's outgoing / incoming links). The
+  function is called for the side effect of drawing the plot.
+}
+\description{
+Draws a circular (Circos) plot of transcription factor to target
+  links from a \code{\link{testEdges}} two-sample result. Genes are placed on
+  their genomic coordinates, links are coloured continuously by significance,
+  flagged as known or novel against an optional a priori network, and genes
+  belonging to supplied gene sets can be labelled around the circle.
+}
+\details{
+Requires the \pkg{circlize} package, and \pkg{biomaRt} when gene
+  coordinates are downloaded automatically (\code{geneCoords = NULL}). Genes
+  without coordinates, and links whose TF or target lacks coordinates, are
+  dropped with a message.
+}
+\examples{
+\dontrun{
+data(scorpionTest)
+nets <- runSCORPION(
+  gexMatrix = scorpionTest$gex,
+  tfMotifs = scorpionTest$tf,
+  ppiNet = scorpionTest$ppi,
+  cellsMetadata = scorpionTest$metadata,
+  groupBy = c("donor", "region")
+)
+res <- testEdges(
+  networksDF = nets,
+  testType = "two.sample",
+  group1 = grep("--T$", colnames(nets), value = TRUE),
+  group2 = grep("--N$", colnames(nets), value = TRUE)
+)
+
+# Human coordinates auto-downloaded from Ensembl, known/novel vs a prior net
+circosEdges(
+  edgesDF = res,
+  species = "hsapiens_gene_ensembl",
+  priorNet = scorpionTest$tf,
+  geneSets = "hallmark.gmt"
+)
+}
+}
+\seealso{
+\code{\link{testEdges}}, \code{\link{runSCORPION}}
+}
+\author{
+Daniel Osorio <daniecos@uio.no>
+}
diff --git a/man/dot-drawCircosLegends.Rd b/man/dot-drawCircosLegends.Rd
new file mode 100644
index 0000000..1771e75
--- /dev/null
+++ b/man/dot-drawCircosLegends.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circosEdges.R
+\name{.drawCircosLegends}
+\alias{.drawCircosLegends}
+\title{Draw effect-size, novelty and gene-set legends on a Circos plot.}
+\usage{
+.drawCircosLegends(
+  col_fun,
+  valueRange,
+  colorBy,
+  hasPrior,
+  knownColor,
+  novelColor,
+  set_colors,
+  lwdRange = NULL,
+  widthShown = FALSE,
+  degCol_fun = NULL,
+  degRange = NULL,
+  degColBy = "log2FoldChange"
+)
+}
+\description{
+Draw effect-size, novelty and gene-set legends on a Circos plot.
+}
+\keyword{internal}
diff --git a/man/dot-fetchGeneCoords.Rd b/man/dot-fetchGeneCoords.Rd
new file mode 100644
index 0000000..fe6d36a
--- /dev/null
+++ b/man/dot-fetchGeneCoords.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circosEdges.R
+\name{.fetchGeneCoords}
+\alias{.fetchGeneCoords}
+\title{Download gene coordinates from Ensembl via biomaRt.}
+\usage{
+.fetchGeneCoords(genes, species, mirror)
+}
+\description{
+Download gene coordinates from Ensembl via biomaRt.
+}
+\keyword{internal}
diff --git a/man/dot-orderChr.Rd b/man/dot-orderChr.Rd
new file mode 100644
index 0000000..48ebfe3
--- /dev/null
+++ b/man/dot-orderChr.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circosEdges.R
+\name{.orderChr}
+\alias{.orderChr}
+\title{Natural-sort chromosome names (1..N, then X, Y, MT, then the rest).}
+\usage{
+.orderChr(chr)
+}
+\description{
+Natural-sort chromosome names (1..N, then X, Y, MT, then the rest).
+}
+\keyword{internal}
diff --git a/man/dot-parseGMT.Rd b/man/dot-parseGMT.Rd
new file mode 100644
index 0000000..42c0cb9
--- /dev/null
+++ b/man/dot-parseGMT.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circosEdges.R
+\name{.parseGMT}
+\alias{.parseGMT}
+\title{Parse a GMT gene-set file into a named list.}
+\usage{
+.parseGMT(path)
+}
+\description{
+Parse a GMT gene-set file into a named list.
+}
+\keyword{internal}
diff --git a/man/dot-validateGeneCoords.Rd b/man/dot-validateGeneCoords.Rd
new file mode 100644
index 0000000..cd4ef36
--- /dev/null
+++ b/man/dot-validateGeneCoords.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circosEdges.R
+\name{.validateGeneCoords}
+\alias{.validateGeneCoords}
+\title{Validate a user-supplied gene coordinate table.}
+\usage{
+.validateGeneCoords(geneCoords)
+}
+\description{
+Validate a user-supplied gene coordinate table.
+}
+\keyword{internal}
diff --git a/man/enrichEdges.Rd b/man/enrichEdges.Rd
new file mode 100644
index 0000000..3c86b29
--- /dev/null
+++ b/man/enrichEdges.Rd
@@ -0,0 +1,89 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/enrichEdges.R
+\name{enrichEdges}
+\alias{enrichEdges}
+\title{Gene set enrichment analysis of TF-target edges}
+\usage{
+enrichEdges(edgesDF, geneSets, numericValue, nCores = 3, seed = 1)
+}
+\arguments{
+\item{edgesDF}{A data.frame of TF-target edges, typically produced by
+\code{\link{testEdges}}. Must contain a \code{tf} column, a \code{target}
+column, and the numeric column named by \code{numericValue}.}
+
+\item{geneSets}{A named list of gene sets. The names of the list elements are
+used as gene set identifiers.}
+
+\item{numericValue}{Character string naming the column in \code{edgesDF} used
+as the ranking statistic for enrichment analysis.}
+
+\item{nCores}{Integer specifying the number of parallel workers to use.
+Default 3.}
+
+\item{seed}{Integer specifying the random seed used by the parallel
+enrichment calculations. Default 1.}
+}
+\value{
+A data.frame of enrichment results with one row per TF-gene set pair:
+  \itemize{
+    \item{tf: Transcription factor}
+    \item{geneSet: Gene set identifier}
+    \item{pValue: Raw enrichment p-value}
+    \item{pAdj: Benjamini-Hochberg adjusted p-value}
+    \item{log2Err: Expected log2 error of the p-value estimate}
+    \item{ES: Enrichment score}
+    \item{NES: Normalized enrichment score}
+    \item{geneSetSize: Number of genes from the set found among the targets}
+  }
+}
+\description{
+Performs gene set enrichment analysis separately for each transcription
+factor (TF) using the edge-level values supplied in \code{numericValue}.
+Enrichment is performed with the multilevel implementation of \pkg{fgsea}.
+Calculations for individual TFs are performed in parallel.
+}
+\details{
+For each TF, the values in \code{numericValue} are used as ranked statistics
+for its target genes. Edges with missing targets or missing or non-finite
+values in the selected numeric column are excluded before enrichment
+analysis. If a target occurs more than once for a TF, only the observation
+with the largest absolute value of the selected ranking statistic is
+retained.
+
+Gene set enrichment is performed using \code{fgsea::fgseaMultilevel}; the
+\pkg{fgsea} package (Bioconductor) is required. The \code{leadingEdge} column
+returned by \code{fgseaMultilevel} is not included in the output. P-values
+are adjusted across all TF-gene set enrichment tests using the
+Benjamini-Hochberg procedure.
+}
+\examples{
+\dontrun{
+data(scorpionTest)
+nets <- runSCORPION(
+  gexMatrix = scorpionTest$gex,
+  tfMotifs = scorpionTest$tf,
+  ppiNet = scorpionTest$ppi,
+  cellsMetadata = scorpionTest$metadata,
+  groupBy = c("donor", "region")
+)
+res <- testEdges(
+  networksDF = nets,
+  testType = "two.sample",
+  group1 = grep("--T$", colnames(nets), value = TRUE),
+  group2 = grep("--N$", colnames(nets), value = TRUE)
+)
+
+geneSets <- list(SetA = c("ACKR1", "ACTA2"), SetB = c("ACTG2", "ADAMDEC1"))
+enr <- enrichEdges(
+  edgesDF = res,
+  geneSets = geneSets,
+  numericValue = "log2FoldChange"
+)
+}
+}
+\seealso{
+\code{\link{testEdges}}, \code{\link{maEdges}}
+}
+\author{
+Daniel Osorio <daniecos@uio.no>
+}
diff --git a/man/maEdges.Rd b/man/maEdges.Rd
new file mode 100644
index 0000000..0e63935
--- /dev/null
+++ b/man/maEdges.Rd
@@ -0,0 +1,82 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/maEdges.R
+\name{maEdges}
+\alias{maEdges}
+\title{Meta-analysis of TF-target edges across studies}
+\usage{
+maEdges(edgesList, method = c("random", "fixed"), minStudies = 2L)
+}
+\arguments{
+\item{edgesList}{A list of data.frames, one per study, typically produced by
+\code{\link{testEdges}}. Each data.frame must contain the columns
+\code{tf}, \code{target}, \code{log2FoldChange} and \code{pValue}.}
+
+\item{method}{Meta-analysis model. Either \code{"random"} (DerSimonian-Laird
+random-effects) or \code{"fixed"} (inverse-variance fixed-effect). Default
+\code{"random"}.}
+
+\item{minStudies}{Minimum number of studies with valid numeric information
+required for a TF-target pair to be included. Default 2.}
+}
+\value{
+A data.frame of meta-analysis results with one row per TF-target pair:
+  \itemize{
+    \item{tf: Transcription factor}
+    \item{target: Target gene}
+    \item{k: Number of contributing studies}
+    \item{log2FoldChange: Meta-analytic effect size}
+    \item{se: Standard error of the effect size}
+    \item{ciLow, ciHigh: Lower and upper bounds of the 95\% confidence interval}
+    \item{zStatistic: Z statistic}
+    \item{pValue: Raw p-value}
+    \item{Q: Cochran's Q heterogeneity statistic}
+    \item{iSquared: I-squared heterogeneity (percentage)}
+    \item{tauSquared: Between-study variance estimate}
+    \item{pAdj: Benjamini-Hochberg adjusted p-value}
+  }
+}
+\description{
+Performs a meta-analysis of TF-target edges across multiple studies using
+either a fixed-effect or DerSimonian-Laird random-effects model. Missing or
+non-finite effect sizes and p-values are excluded from the corresponding
+study. A TF-target pair is only counted as contributing to a study when both
+its effect size and p-value are valid.
+}
+\details{
+Standard errors are estimated from the reported log2 fold change and
+two-sided p-value. Fixed-effect estimates use inverse-variance weighting.
+Random-effects estimates use DerSimonian-Laird estimates of between-study
+variance.
+
+Missing or non-finite values in \code{log2FoldChange} or \code{pValue} do not
+contribute to the meta-analysis and do not count toward \code{minStudies}.
+TF-target pairs with fewer than \code{minStudies} valid studies are excluded.
+Benjamini-Hochberg correction is performed once across all TF-target
+meta-analysis tests.
+}
+\examples{
+\dontrun{
+data(scorpionTest)
+nets <- runSCORPION(
+  gexMatrix = scorpionTest$gex,
+  tfMotifs = scorpionTest$tf,
+  ppiNet = scorpionTest$ppi,
+  cellsMetadata = scorpionTest$metadata,
+  groupBy = c("donor", "region")
+)
+
+# Two independent comparisons treated as separate studies
+studyA <- testEdges(nets, "two.sample",
+  group1 = "P31--T", group2 = "P31--N")
+studyB <- testEdges(nets, "two.sample",
+  group1 = "P32--T", group2 = "P32--N")
+
+meta <- maEdges(list(studyA, studyB), method = "random")
+}
+}
+\seealso{
+\code{\link{testEdges}}, \code{\link{enrichEdges}}
+}
+\author{
+Daniel Osorio <daniecos@uio.no>
+}
diff --git a/man/regressEdges.Rd b/man/regressEdges.Rd
index 7d1f81e..a6c400e 100644
--- a/man/regressEdges.Rd
+++ b/man/regressEdges.Rd
@@ -104,3 +104,9 @@ strong_trends <- results_regression[results_regression$pAdj < 0.05 &
                                      abs(results_regression$meanEdge) > 0.1, ]
 }
 }
+\seealso{
+\code{\link{runSCORPION}}, \code{\link{testEdges}}
+}
+\author{
+Daniel Osorio <daniecos@uio.no>
+}
diff --git a/man/runSCORPION.Rd b/man/runSCORPION.Rd
index 919f168..7830f14 100644
--- a/man/runSCORPION.Rd
+++ b/man/runSCORPION.Rd
@@ -247,3 +247,9 @@ nets_batch_corrected <- runSCORPION(
 # 6 AEBP2_EED_EZH2_RBBP4_SUZ12  ACKR1 -0.3274770 -0.34475499 -0.12449908
 }
 }
+\seealso{
+\code{\link{scorpion}}, \code{\link{testEdges}}, \code{\link{regressEdges}}
+}
+\author{
+Daniel Osorio <daniecos@uio.no>
+}
diff --git a/man/scorpionTest.Rd b/man/scorpionTest.Rd
index 5ddbb44..b3bff89 100644
--- a/man/scorpionTest.Rd
+++ b/man/scorpionTest.Rd
@@ -3,20 +3,33 @@
 \docType{data}
 \name{scorpionTest}
 \alias{scorpionTest}
-\title{Example single-cell gene expression, motif, and ppi data}
+\title{Example single-cell colorectal cancer data for SCORPION}
 \format{
-A list containing three datasets.
+A named list with four elements:
 \describe{
-\item{\code{gex}}{A subsetted version of 10X Genomics' 3k PBMC dataset provided by the \code{Seurat} package.}
-\item{\code{tf}}{Subset of the transcription-factor and target gene list provided by the \code{dorothea} package for Homo sapiens.}
-\item{\code{ppi}}{The known protein-protein interactions and the combined score downloaded from the STRING database}
+  \item{\code{gex}}{A \code{dgCMatrix} gene expression matrix with 300 genes
+    (rows) and 1,954 cells (columns).}
+  \item{\code{tf}}{A \code{data.frame} of transcription factor-target motif
+    pairs from DoRothEA with columns \code{source_genesymbol},
+    \code{target_genesymbol} and \code{weight} (371,738 rows).}
+  \item{\code{ppi}}{A \code{data.frame} of protein-protein interactions with
+    columns \code{source_genesymbol}, \code{target_genesymbol} and
+    \code{weight} (4,076 rows).}
+  \item{\code{metadata}}{A \code{data.frame} of cell-level annotations with
+    columns \code{cell_id}, \code{donor}, \code{region} and \code{cell_type}
+    (1,954 rows). Region codes are \code{T} (tumor), \code{B} (border) and
+    \code{N} (normal).}
 }
 }
 \usage{
 data(scorpionTest)
 }
 \description{
-This data is a list containing three objects. The motif \code{data.frame} describes a set of pairwise connections where a specific known sequence motif of a transcription factor was found upstream of the corresponding gene. The expression \code{dgCMatrix} is a set of 230 gene expression levels measured across 80 cells. Finally, the ppi \code{data.frame} describes a set of known pairwise protein-protein interactions.
+A list bundling the inputs required to build and compare gene regulatory
+networks with SCORPION, derived from a colorectal cancer single-cell
+RNA-sequencing experiment. It contains a gene expression matrix, a
+transcription factor motif prior, a protein-protein interaction prior, and
+cell-level metadata.
 }
 \examples{
 # Loading example data
@@ -24,23 +37,4 @@ data(scorpionTest)
 
 # The structure of the data
 str(scorpionTest)
-
-# List of 3
-# $ gex:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
-# .. ..@ i       : int [1:4456] 1 5 8 11 22 30 33 34 36 38 ...
-# .. ..@ p       : int [1:81] 0 47 99 149 205 258 306 342 387 423 ...
-# .. ..@ Dim     : int [1:2] 230 80
-# .. ..@ Dimnames:List of 2
-# .. .. ..$ : chr [1:230] "MS4A1" "CD79B" "CD79A" "HLA-DRA" ...
-# .. .. ..$ : chr [1:80] "ATGCCAGAACGACT" "CATGGCCTGTGCAT" "GAACCTGATGAACC" "TGACTGGATTCTCA" ...
-# .. ..@ x       : num [1:4456] 1 1 3 1 1 4 1 5 1 1 ...
-# .. ..@ factors : list()
-# $ tf :'data.frame':	4485 obs. of  3 variables:
-#   ..$ tf    : chr [1:4485] "ADNP" "ADNP" "ADNP" "AEBP2" ...
-# ..$ target: chr [1:4485] "PRF1" "TMEM40" "TNFRSF1B" "CFP" ...
-# ..$ mor   : num [1:4485] 1 1 1 1 1 1 1 1 1 1 ...
-# $ ppi:'data.frame':	12754 obs. of  3 variables:
-#   ..$ X.node1       : chr [1:12754] "ADNP" "ADNP" "ADNP" "AEBP2" ...
-# ..$ node2         : chr [1:12754] "ZBTB14" "NFIA" "CDC5L" "YY1" ...
-# ..$ combined_score: num [1:12754] 0.769 0.64 0.581 0.597 0.54 0.753 0.659 0.548 0.59 0.654 ...
 }
diff --git a/man/testEdges.Rd b/man/testEdges.Rd
index b9ee5e3..c49ee3d 100644
--- a/man/testEdges.Rd
+++ b/man/testEdges.Rd
@@ -176,3 +176,9 @@ results_paired <- testEdges(
 )
 }
 }
+\seealso{
+\code{\link{runSCORPION}}, \code{\link{regressEdges}}, \code{\link{maEdges}}, \code{\link{circosEdges}}
+}
+\author{
+Daniel Osorio <daniecos@uio.no>
+}
diff --git a/tests/testthat/test_circos_edges.R b/tests/testthat/test_circos_edges.R
new file mode 100644
index 0000000..02465b7
--- /dev/null
+++ b/tests/testthat/test_circos_edges.R
@@ -0,0 +1,99 @@
+test_that("circosEdges draws with user-supplied coordinates (offline)", {
+  skip_on_cran()
+  skip_if_not_installed("circlize")
+
+  set.seed(1)
+  n <- 30
+  edges <- data.frame(
+    tf = paste0("TF", sample(1:5, n, replace = TRUE)),
+    target = paste0("G", sample(1:20, n, replace = TRUE)),
+    log2FoldChange = rnorm(n),
+    tStatistic = rnorm(n),
+    pValue = runif(n, 0, 0.01),
+    pAdj = runif(n, 0, 0.04),
+    stringsAsFactors = FALSE
+  )
+  edges <- edges[edges$tf != edges$target, ]
+
+  all_genes <- unique(c(edges$tf, edges$target))
+  coords <- data.frame(
+    gene = all_genes,
+    chr = sample(c("1", "2", "X"), length(all_genes), replace = TRUE),
+    start = sample(1:1e6, length(all_genes)),
+    end = sample(1:1e6, length(all_genes)) + 1000,
+    stringsAsFactors = FALSE
+  )
+
+  prior <- data.frame(tf = edges$tf[1:5], target = edges$target[1:5])
+  gene_sets <- list(SetA = all_genes[1:5], SetB = all_genes[6:10])
+
+  tmp <- tempfile(fileext = ".pdf")
+  grDevices::pdf(tmp)
+  on.exit({
+    grDevices::dev.off()
+    unlink(tmp)
+  }, add = TRUE)
+
+  out <- circosEdges(
+    edgesDF = edges,
+    geneCoords = coords,
+    priorNet = prior,
+    geneSets = gene_sets,
+    pAdjThreshold = 0.05,
+    nmaxTF = 3,
+    nmaxTarget = 3
+  )
+
+  expect_type(out, "list")
+  expect_true(all(c("edges", "coords") %in% names(out)))
+  expect_true(all(c("outDegree", "inDegree", "degree") %in% names(out$coords)))
+  expect_equal(out$coords$degree, out$coords$outDegree + out$coords$inDegree)
+  expect_true(all(out$edges$novelty %in% c("known", "novel")))
+  expect_true(any(out$edges$novelty == "known"))
+  expect_true(all(out$edges[["pAdj"]] < 0.05))
+})
+
+test_that("circosEdges parses a GMT file", {
+  skip_on_cran()
+  gmt <- tempfile(fileext = ".gmt")
+  writeLines(c(
+    "SET1\tdescription\tGENEA\tGENEB\tGENEC",
+    "SET2\thttp://example\tGENEC\tGENED"
+  ), gmt)
+  on.exit(unlink(gmt), add = TRUE)
+
+  parseGMT <- getFromNamespace(".parseGMT", "SCORPION")
+  sets <- parseGMT(gmt)
+  expect_named(sets, c("SET1", "SET2"))
+  expect_equal(sets$SET1, c("GENEA", "GENEB", "GENEC"))
+  expect_equal(sets$SET2, c("GENEC", "GENED"))
+})
+
+test_that("circosEdges errors on missing columns and empty selections", {
+  skip_on_cran()
+  skip_if_not_installed("circlize")
+
+  bad <- data.frame(tf = "A", target = "B")
+  expect_error(circosEdges(bad, colorBy = "log2FoldChange"), "log2FoldChange")
+
+  edges <- data.frame(
+    tf = "TF1", target = "G1", log2FoldChange = 1,
+    pValue = 0.5, pAdj = 0.9, stringsAsFactors = FALSE
+  )
+  coords <- data.frame(
+    gene = c("TF1", "G1"), chr = c("1", "1"),
+    start = c(1, 100), end = c(10, 110), stringsAsFactors = FALSE
+  )
+  expect_error(
+    circosEdges(edges, geneCoords = coords, pAdjThreshold = 0.05),
+    "thresholds"
+  )
+})
+
+test_that(".orderChr sorts chromosomes naturally", {
+  orderChr <- getFromNamespace(".orderChr", "SCORPION")
+  expect_equal(
+    orderChr(c("2", "X", "1", "MT", "10", "Y")),
+    c("1", "2", "10", "X", "Y", "MT")
+  )
+})