Design principles

The design of epiPortrait rests on five principles:

A domain can change in different ways

Across conditions, a domain’s architecture state can remodel in distinct directions:

The default cutoff_scope = "relative" transition expresses a relative rank-state change, not absolute biochemical gain/loss (see Module 4).

Conceptual model

epiPortrait distinguishes the shared coordinate frame (the consensus candidate-domain universe, one interval per domain across all samples) from the sample-specific native peak geometry (each replicate’s own peak/domain calls). The Breadth-Super call is derived from native peak geometry, not from the width of the shared consensus interval.

Relationship to existing methods

epiPortrait complements these tools as a replicate-aware quantitative domain-phenotyping layer (Intensity × native Breadth decomposition, plus within-domain dispersion), rather than a peak caller or count-based differential method. The primary use case is histone-mark ChIP-seq or CUT&Tag; continuous-track workflows also apply to ATAC-seq-derived domains.

Because the input is normalized continuous signal (BigWig) rather than raw alignments, epiPortrait does not perform read alignment, duplicate removal, or count-based differential binding. It is a downstream domain architecture analysis tool.

1.1 The candidate-domain universe is required

Two different “peak files” are involved — do not confuse them.

1. Candidate-domain universe (consensus_peaks, one shared set for all samples) defines which regions you analyze. It comes from your peak caller’s bed/narrowPeak output (MACS2, SICER, epic2, …), the same upstream input ROSE consumes. Every analysis needs this — including single-sample Intensity-only.
   • Single sample: consensus_peaks is that sample’s own peak file.
   • Multiple samples: it is the consensus/union via get_consensus_peaks().
   • Exception: you may import an externally defined domain set (fixed tiling, database interval table, or published domains) as a GRanges instead.
2. Per-sample native peaks (peak_path, one file per sample) are each sample’s own peak calls, used to measure native enrichment width per sample. Only the Breadth axis needs this; Intensity-only analyses use the candidate universe alone and omit peak_path.

Bottom line: peak files are always required (to build the candidate universe) — the only difference is how many and what for. Intensity-only uses them to define the domains; Breadth additionally uses each sample’s own peaks for width measurement. There is no analysis mode that runs from BigWig alone.

1.2 Consensus peaks from your peak-caller output

if (use_bw) {
  peak_list <- list(
    rep1 = rtracklayer::import(system.file("extdata", "peaks.bed", package = "epiPortrait")),
    rep2 = rtracklayer::import(system.file("extdata", "peaks.bed", package = "epiPortrait"))
  )
  consensus_peaks <- get_consensus_peaks(peak_list, min_reps = 1)
} else {
  consensus_peaks <- rowRanges(example_se)
}
consensus_peaks

## GRanges object with 3 ranges and 1 metadata column:
##       seqnames            ranges strand | Support_Reps
##          <Rle>         <IRanges>  <Rle> |    <numeric>
##   [1]    chr21 43999000-44002000      * |            2
##   [2]    chr21 44019000-44023000      * |            2
##   [3]    chr21 44029000-44044000      * |            2
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

1.3 Optionally stitch into macro-domains

if (use_bw) {
  # For super-enhancer analysis, 12.5 kb is the ROSE-standard stitching distance
  macro_domains <- stitch_epi_peaks(consensus_peaks, stitch_distance = 12500)
} else {
  macro_domains <- consensus_peaks
}
length(macro_domains)

## [1] 2

1.4 Build the portrait matrix

if (use_bw) {
  extdata <- system.file("extdata", package = "epiPortrait")
  samples <- data.frame(
    SampleID  = c("C1", "C2", "T1", "T2"),
    Condition = c("Control", "Control", "Treatment", "Treatment"),
    Replicate = c(1, 2, 1, 2),
    bw_path   = file.path(extdata, c("C1.bw", "C2.bw", "T1.bw", "T2.bw")),
    peak_path = file.path(extdata, c("C1_peaks.bed", "C2_peaks.bed",
                                     "T1_peaks.bed", "T2_peaks.bed"))
  )
  se <- build_portrait_matrix(
    sample_sheet    = samples,
    consensus_peaks = macro_domains,
    workers         = 1
  )
} else {
  se <- example_se
}
se

## class: RangedSummarizedExperiment 
## dim: 2 4 
## metadata(6): signal_contract domain_provenance ...
##   native_peak_available provenance
## assays(5): Intensity SignalDispersion NativeMaxPeakWidth
##   NativeOccupiedWidth NativePeakCount
## rownames(2): epiDomain_000001 epiDomain_000002
## rowData names(2): Constituent_Peaks IntervalWidth
## colnames(4): C1 C2 T1 T2
## colData names(6): SampleID Condition ... peak_path NegativeFraction

The output SummarizedExperiment contains the core assays Intensity and SignalDispersion, plus — when per-sample native peak files are provided via peak_path — the sample-specific breadth geometry (NativeMaxPeakWidth, NativeOccupiedWidth, NativePeakCount). The static genomic interval length is stored in rowData(se)$IntervalWidth (not an assay — it is identical across samples).

peak_path is optional. Without it, native breadth geometry and canonical Breadth-Super calling are unavailable; Intensity-Super remains available, and SignalDispersion can still be explored as a secondary architecture descriptor. When peak_path is supplied, the imported native peaks are retained in metadata(se)$native_peaks so that peak-level breadth calls can be inspected.

For your own data, replace the sample sheet paths with your real normalized BigWigs (and native peak files), and consensus_peaks with your consensus peak set.

se <- build_portrait_matrix(
  sample_sheet    = samples,
  consensus_peaks = macro_domains,
  custom_features = list(
    Entropy = function(x) {
      p <- (x + 1) / sum(x + 1)
      -sum(p * log2(p))
    }
  )
)

3.1 Inflection methods

Two heuristics are available via the method argument:

• "elbow" (default): the point with the maximum perpendicular distance to the first-last ranked line — classic knee detection.
• "tangent": a ROSE-inspired tangent-optimization inflection (Whyte et al., 2013). It slides a line whose slope is (max - min)/n and minimizes the number of points below it. Note that this is a geometric variant of ROSE’s calculate_cutoff() step, not a bit-exact reproduction; it is only an inflection step, not the full ROSE pipeline.

se_elbow <- call_super_domains(example_se, feature = "Intensity",
                               method = "elbow", verbose = FALSE)
se_tangent <- call_super_domains(example_se, feature = "Intensity",
                                 method = "tangent", verbose = FALSE)
table(rowData(se_tangent)$Intensity_Domain_Type)

## 
## Intensity_Super_Element       Intensity_Typical 
##                      50                     450

Explicit top-fraction cutoff (quantile_cutoff). Both axes offer an explicit top-fraction alternative that bypasses inflection detection: for Intensity it is the given quantile of the ranked domain values; for Breadth it is the same quantile of each replicate’s genome-wide eligible native PeakWidth distribution, so the top 100*(1-quantile_cutoff)% widest native peaks are Broad (e.g. 0.95 = top 5% — a common broad-domain practice). This option does not assess cutoff quality. It never automatically falls back to, or from, the data-driven inflection path, so the no-call semantics of the default path are unchanged.

3.2 log-transform is optional

By default Intensity and SignalDispersion use the package’s automatic log10(x + 1) ranking transform (log_transform = NULL, auto per feature). Breadth uses the native peak-width calling pipeline described below, and the static IntervalWidth is not auto-log-transformed. Set log_transform = FALSE to rank on the raw scale, or log_transform = TRUE to force it.

Optional H3K27ac benchmark against ROSE. If you plan a head-to-head comparison with the original ROSE method (Whyte et al., 2013), use method = "tangent" with log_transform = FALSE (raw signal, ROSE-style scale). Because the tangent implementation is a geometric variant rather than a bit-exact port, validate any ROSE-reproduction numbers against the reference calculate_cutoff() before drawing quantitative conclusions. This is an optional benchmark setting, not the package’s primary workflow. The log-transformed default is retained as the epiPortrait convention because it stabilises the skewed intensity distribution for the elbow heuristic and downstream analysis.

All Super labels are relative to the candidate-domain universe and selected cutoff. They should not be interpreted as functional super-enhancers without independent evidence. Intensity is the area under the signal track, so it can increase with signal amplitude, interval length, or both.

3.3 Modes

• mode = "global_consensus" (default; alias "consensus"): calls each sample independently, computes replicate support per condition group, and requires the support rule to be met in EVERY group (persistent consensus).
• mode = "per_group" (recommended for multi-condition studies): each replicate is called independently using its own ranked feature distribution; replicate-level calls are then aggregated within each condition using the selected support rule. One call column per group. The per-domain x per-replicate call matrix is available via get_replicate_calls(). Group-mean ranks are stored for visualization only and never determine the group call.
• mode = "per_sample": calls each sample independently.

se_pg <- call_super_domains(example_se, feature = "Intensity",
                            mode = "per_group", group_var = "Condition",
                            verbose = FALSE)
grep("Intensity_Call__", colnames(SummarizedExperiment::rowData(se_pg)), value = TRUE)

## [1] "Intensity_Call__Control"   "Intensity_Call__Treatment"

3.4 Cutting on breadth: peak-level Breadth-Super

Breadth is the peak-level breadth feature. It does not rank a domain x sample assay: instead each replicate’s genome-wide eligible native peak width distribution is cut by an elbow/inflection, broad peaks are mapped to the shared domains by unique assignment, and the evidence is aggregated across replicates:

se_b <- call_super_domains(example_se, feature = "Breadth",
                           mode = "per_group", group_var = "Condition",
                           verbose = FALSE)
table(rowData(se_b)$Breadth_Call__Control)

## 
## Breadth_Super_Element       Breadth_Typical 
##                    54                   440

Key parameters for feature = "Breadth":

• min_peak_overlap_fraction (default 0.5): a native peak is assigned to at most one shared domain (the one with maximum overlap) only if at least this fraction of its extent lies inside it; peaks below the threshold are Unmapped, and tied maxima are Ambiguous — preventing a broad peak from double-counting into neighbouring domains.
• valid_chroms: optional vector of allowed chromosomes for the genome-wide eligible native peak set.

The peak-level call (PeakWidth / PeakBroadCall / PeakRank) is stored in metadata(se)$breadth_peak_calls, and the unique-assignment provenance in metadata(se)$breadth_peak_mapping.

IntervalWidth (alias Width) remains available for static region-level broadness ranking only — it is constant across samples and never defines a sample-specific breadth.

3.5 Combined taxonomy

Combine Intensity and Breadth calls into the four-class taxonomy (plus Uncertain for domains where a call was not reliable):

se <- call_super_domains(example_se, feature = "Intensity",
                         mode = "per_group", group_var = "Condition",
                         verbose = FALSE)
se <- call_super_domains(se, feature = "Breadth",
                         mode = "per_group", group_var = "Condition", verbose = FALSE)
se <- combine_superdomain_calls(se, group_var = "Condition")
table(rowData(se)$Combined_Class__Control)

## 
##   Breadth-Super      Dual-Super Intensity-Super         Typical       Uncertain 
##              39              15              34             406               6

5.1 Domain-aware annotation

se <- annotate_epi_domains(se, genome = "hg38")
# rowData: primary_genomic_context, nearest_tss_gene_id,
#          nearest_tss_gene_symbol, nearest_tss_distance_bp,
#          promoter_overlap_gene_count, gene_body_overlap_gene_count,
#          fully_contained_gene_count
head(SummarizedExperiment::rowData(se)$primary_genomic_context)

## [1] "Intergenic" "Intergenic" "Intergenic" "Intergenic" "Intergenic"
## [6] "Intergenic"

# long-format links: metadata(se)$domain_gene_links
candidates <- get_domain_genes(se)
head(candidates)

##          domain_id   gene_id  gene_symbol
## 1 epiDomain_000053 102724200 LOC102724200
## 2 epiDomain_000057 105379487    LINC01670
## 3 epiDomain_000076 105379499 LOC105379499
## 4 epiDomain_000098 105372730    LINC01666
## 5 epiDomain_000115 105372733 LOC105372733
## 6 epiDomain_000160 112268283 LOC112268283
##                                   relation_types bedpe_supported
## 1 gene_body_overlap;nearest_tss;promoter_overlap           FALSE
## 2 gene_body_overlap;nearest_tss;promoter_overlap           FALSE
## 3                   nearest_tss;promoter_overlap           FALSE
## 4 gene_body_overlap;nearest_tss;promoter_overlap           FALSE
## 5                   nearest_tss;promoter_overlap           FALSE
## 6                   nearest_tss;promoter_overlap           FALSE
##   bedpe_support_count bedpe_contact_score nearest_tss_distance_bp
## 1                   0                   0                    2695
## 2                   0                   0                    2603
## 3                   0                   0                    2096
## 4                   0                   0                     946
## 5                   0                   0                     476
## 6                   0                   0                    2464
##   expression_value expression_status expression_rank candidate_priority
## 1               NA              <NA>              NA                  1
## 2               NA              <NA>              NA                  2
## 3               NA              <NA>              NA                  3
## 4               NA              <NA>              NA                  4
## 5               NA              <NA>              NA                  5
## 6               NA              <NA>              NA                  6

5.2 Hockey-stick plot

se <- call_super_domains(example_se, feature = "Intensity",
                         method = "tangent", log_transform = FALSE, verbose = FALSE)
plot_hockey_stick(se, feature = "Intensity", label_genes = FALSE)

5.3 Single-domain feature profile

Quantitative per-domain view across conditions: each feature is faceted with biological replicates as points and group means as open points; a dashed line connects the group means so the across-condition direction of change is read directly. When the object carries combined architecture classes (combine_superdomain_calls()), the domain’s class in each condition is annotated under the axis, linking the raw values to its architectural identity:

plot_domain_feature_profile(example_se, peak_id = rownames(example_se)[1],
                            group_var = "Condition")

5.4 PCA and correlation QC

plot_portrait_pca(example_se, feature = "Intensity", group_var = "Condition")

plot_portrait_correlation(example_se, feature = "Intensity")

5.5 Result visualizations

The plotting functions summarize Intensity versus native breadth, replicate support, and condition-level state transitions:

se_viz <- call_super_domains(example_se, feature = "Intensity",
                             mode = "per_group", group_var = "Condition",
                             verbose = FALSE)
se_viz <- call_super_domains(se_viz, feature = "Breadth",
                             mode = "per_group", group_var = "Condition",
                             verbose = FALSE)
se_viz <- combine_superdomain_calls(se_viz, group_var = "Condition")
se_viz <- compare_superdomain_classes(se_viz,
                                      ref_group = "Control",
                                      target_group = "Treatment")

plot_domain_landscape(se_viz, group = "Control")            # no domain labels by default

plot_class_transition(se_viz, ref_group = "Control", target_group = "Treatment")

plot_replicate_support(se_viz, feature = "Intensity", group = "Control")   # support >= 1/3

plot_domain_class_composition(se_viz)                       # default: super-class focus

plot_domain_class_composition(se_viz, focus = "all")        # full composition incl. Typical

# label the representative domain of each super class (opt-in)
plot_domain_landscape(se_viz, group = "Control", label_n = 5)

# per-group hockey with a descriptive median replicate cutoff
plot_hockey_stick(se_viz, feature = "Intensity", group = "Control",
                  label_genes = FALSE)

The landscape places each domain in the Intensity by native-breadth plane and colours it by architecture class (pass label_n > 0 to name the most representative domain of each super class). The transition matrix expresses a relative architecture-state change (per-group cutoffs) — it must not be read as absolute gain/loss. Because non-super domains usually dominate the total, the class-composition plot defaults to a focused view of the super classes (with their share of all domains annotated above each bar); pass focus = "all" for the full composition.

5.6 Saving and reusing an analysis

Save the SummarizedExperiment to preserve the assays, calls, and provenance. Flat tables and plots can then be regenerated without repeating BigWig I/O:

f <- tempfile(fileext = ".rds")
saveRDS(se_viz, f)
se2 <- readRDS(f)

# flat tables for Excel / Python / reviewers
outdir <- tempfile("epiPortrait_results_")
dir.create(outdir)
export_epiportrait_results(se2, outdir = outdir)
# a few representative columns (the full table is written to
# epiPortrait_results/domain_results.tsv)
head(get_domain_results(se2, group_var = "Condition")[, c(
  "Domain_ID", "chr", "start", "end", "IntervalWidth",
  "Intensity_Call__Control", "Breadth_Call__Control",
  "Combined_Class__Control", "Combined_Class__Treatment")])

##                         Domain_ID   chr   start     end IntervalWidth
## epiDomain_000001 epiDomain_000001 chr21 1000000 1000999          1000
## epiDomain_000002 epiDomain_000002 chr21 1080000 1080999          1000
## epiDomain_000003 epiDomain_000003 chr21 1160000 1179999         20000
## epiDomain_000004 epiDomain_000004 chr21 1240000 1240999          1000
## epiDomain_000005 epiDomain_000005 chr21 1320000 1320999          1000
## epiDomain_000006 epiDomain_000006 chr21 1400000 1400999          1000
##                  Intensity_Call__Control Breadth_Call__Control
## epiDomain_000001       Intensity_Typical       Breadth_Typical
## epiDomain_000002       Intensity_Typical       Breadth_Typical
## epiDomain_000003       Intensity_Typical       Breadth_Typical
## epiDomain_000004       Intensity_Typical       Breadth_Typical
## epiDomain_000005       Intensity_Typical       Breadth_Typical
## epiDomain_000006       Intensity_Typical       Breadth_Typical
##                  Combined_Class__Control Combined_Class__Treatment
## epiDomain_000001                 Typical                   Typical
## epiDomain_000002                 Typical                   Typical
## epiDomain_000003                 Typical                   Typical
## epiDomain_000004                 Typical                   Typical
## epiDomain_000005                 Typical                   Typical
## epiDomain_000006                 Typical                   Typical

prov <- get_call_provenance(se2, "Intensity")
prov[c("method", "log_transform_used", "cutoff", "call_status")]

## $method
## [1] "elbow"
## 
## $log_transform_used
## [1] TRUE
## 
## $cutoff
## [1] NA
## 
## $call_status
## [1] NA

unlink(f)
unlink(outdir, recursive = TRUE)

5.7 Differential domain analysis on continuous signal

Where compare_superdomains() reports relative state transitions between discrete phenotypes, analyze_differential_domains() adds a statistical P-value for a continuous signal change between condition groups, computed directly on the integrated BigWig intensity with limma:

se_diff <- analyze_differential_domains(
  se_viz, feature = "Intensity",
  ref_group = "Control", target_group = "Treatment")
table(rowData(se_diff)$Intensity_DiffStatus)   # Gain / Loss / NS

## 
##  NS 
## 500

Paired, blocked, and multi-factor designs are supported through design and a valid design coefficient or contrast:

se_diff <- analyze_differential_domains(
  se_diff, feature = "Intensity",
  # Patient and Condition must be columns of colData(se_diff)
  design = ~ Patient + Condition,
  contrast = "ConditionTreatment")

The volcano plot can label the top-ranked domains:

plot_differential_volcano(se_diff, feature = "Intensity", label_n = 5)

Scope caveat. The test uses log-transformed integrated continuous signal, not raw fragment counts. The simple two-group design requires at least two samples per group to estimate residual variance. For count-based inference, obtain per-domain counts and use DiffBind, csaw, edgeR, or DESeq2.

5.8 Score-aware BEDPE/3D annotation

BEDPE contact evidence can be attached two ways. The default is contact-only: any record linking a domain anchor to a gene promoter counts toward that domain-gene pair. In many HiChIP pipelines (e.g. FitHiChIP TMM output), however, the loop set is shared across conditions and only the per-record score (contact frequency) varies with the perturbation. Collapsing scores to binary contact flags then discards the quantitative dynamics. Passing bedpe_score_col keeps the strength. This small self-contained example builds a synthetic SummarizedExperiment over two tiny hg38-like genes, with two BEDPE records (scores 12.3 and 8.9) bridging domain 1 to the gene-2 promoter:

gids  <- c("gene1", "gene2")
gmod  <- GRanges(rep("chr1", 4),
                 IRanges(c(1000, 1000, 20000, 20000),
                         c(3000, 3000, 26000, 26000)),
                 strand = rep(c("+", "-"), each = 2))
mcols(gmod)$type   <- rep(c("gene", "transcript"), 2)
mcols(gmod)$ID     <- c(rbind(gids, paste0(gids, "_tx1")))
mcols(gmod)$Parent <- c(rbind(rep(NA, 2), gids))
txdb <- txdbmaker::makeTxDbFromGRanges(gmod)

se_small <- SummarizedExperiment(
  assays  = list(Intensity = matrix(1, 3, 2)),
  rowRanges = GRanges("chr1", IRanges(c(1500, 21000, 35000), c(5000, 30000, 90000))))
rownames(se_small) <- paste0("epiDomain_", 1:3)
colnames(se_small) <- c("C1", "T1")

bedpe_df <- data.frame(
  chrom1 = "chr1", start1 = 1499,  end1 = 2500,   # inside epiDomain_1
  chrom2 = "chr1", start2 = 21999, end2 = 23000,  # gene2 promoter anchor
  score  = c(12.3, 8.9),                           # two records -> summed, not maxed
  stringsAsFactors = FALSE)
se_bp <- annotate_epi_domains(se_small, genome = txdb,
                              bedpe = bedpe_df, bedpe_score_col = "score")
rowData(se_bp)$bedpe_contact_score                  # per-domain summed score

## [1] 21.2 21.2  0.0

get_domain_genes(se_bp, rank_by = "bedpe_score")    # dominant-loop-first prioritization

##     domain_id gene_id gene_symbol         relation_types bedpe_supported
## 1 epiDomain_1   gene2        <NA> bedpe_promoter_contact            TRUE
## 2 epiDomain_2   gene1        <NA> bedpe_promoter_contact            TRUE
##   bedpe_support_count bedpe_contact_score nearest_tss_distance_bp
## 1                   2                21.2                      NA
## 2                   2                21.2                      NA
##   expression_value expression_status expression_rank candidate_priority
## 1               NA              <NA>              NA                  1
## 2               NA              <NA>              NA                  2

• bedpe_score_col accepts a numeric column index (>= 7) or a column name.
• Scores are summed over unique supporting records per domain-gene pair (by record id, not score value — TMM scores are often quantized and repeat).
• Semantics: no BEDPE evidence → 0; evidence without a usable score (or a run without bedpe_score_col) → NA, so “no 3D evidence” and “unscored 3D evidence” stay distinguishable.
• get_domain_genes(rank_by = "bedpe_score") ranks candidate genes by descending contact strength, with evidence tier and expression as tie-breakers. This ranking does not establish a causal target gene.

5.9 Per-sample time-course transitions

compare_superdomain_classes() requires per-group call columns. In time-course designs with one sample per condition (e.g. 0 min / 30 min / 4 h / 16 h / 72 h), call in per_sample mode and use transition_matrix_per_sample() to tabulate pairwise transitions. With the simulated example_se the “samples” are Control_1..3 / Treatment_1..3; a single reference sample demonstrates the ref-to-target form:

se_ps <- call_super_domains(example_se, feature = "Intensity",
                            mode = "per_sample", verbose = FALSE)
# ref = one sample (time point); all others are targets
#   a real time course would pass e.g. ref = "0m", timepoints = c("0m","30m","4h","16h","72h")
se_ps <- transition_matrix_per_sample(se_ps, feature = "Intensity",
                                      ref = "Control_1",
                                      timepoints = colnames(example_se),
                                      verbose = FALSE)
names(S4Vectors::metadata(se_ps)$transitions)   # one entry per (ref, target) pair

## [1] "Control_1_vs_Control_2"   "Control_1_vs_Control_3"  
## [3] "Control_1_vs_Treatment_1" "Control_1_vs_Treatment_2"
## [5] "Control_1_vs_Treatment_3"

head(rowData(se_ps)$Intensity_SampleTransition__Control_1_vs_Treatment_1)

## [1] "Persistent_Intensity_Typical" "Persistent_Intensity_Typical"
## [3] "Persistent_Intensity_Typical" "Persistent_Intensity_Typical"
## [5] "Persistent_Intensity_Typical" "Persistent_Intensity_Typical"

Labels match compare_superdomain_classes() semantics: Persistent_<state>, <from>_to_<to>, and Uncertain when either time point has a no-call / unreliable call. Omit ref to compute all ordered time-point pairs. With one sample per time point these transitions are descriptive; they do not estimate biological variability.

6.1 H3K4me3: promoter breadth and transcriptional state

Broad H3K4me3 domains have been associated with cell-identity genes and transcriptional consistency (Benayoun et al., 2014). With epiPortrait, one can test whether a broad promoter is also intensity-extreme and whether a breadth change accompanies an intensity change:

• Breadth captures broad promoter domains;
• Intensity captures high-magnitude promoters;
• Dual-Super identifies promoters that are both broad and high-magnitude.

epiPortrait does not claim to discover broad H3K4me3; it decomposes this known breadth biology into separable quantitative remodeling modes.

6.2 H3K27ac: enhancer-domain magnitude versus regulatory territory

A high-ranking enhancer domain can have concentrated signal, a broad native enriched span, or both. Intensity-Super, Breadth-Super, and Dual-Super express these enhancer domain phenotypes. For an optional head-to-head comparison with the ROSE super-enhancer pipeline, use method = "tangent" with log_transform = FALSE (a geometric variant, not a bit-exact ROSE reproduction).

6.3 H3K27me3 / H3K9me3: broad repressive domains

For broad repressive marks, “Super” is an internal class label; the mark-appropriate surface terminology is Intensity-Extreme, Extended-Domain, and Dual-Extreme. Use the upstream broad-domain caller (e.g. epic2 / SICER) as the native peak source, and keep the preset’s stitch_distance = 0 sentinel (no re-stitching of the caller-defined domains).

6.4 Application patterns

Organize your analysis around the biological question rather than the mark:

For cohort-style datasets, independent patients provide population-level recurrence evidence, not repeated measurements of one biological sample. Report the support denominator and choose the threshold to match the recurrence question. With one sample per condition, use mode = "per_sample" rather than making a group-level support claim.

6.5 Cross-mark integration across domain universes

A SummarizedExperiment object represents one mark. When an analysis combines H3K27ac Intensity-Super enhancer domains with H3K4me3 Breadth-Super promoter domains, the objects need not share a domain universe. Stitching and promoter filtering may differ by mark. integrate_cross_mark() projects calls from one object onto overlapping domains in the other:

# se_h3k27ac carries Intensity_Call__<tp>; se_h3k4me3 carries Breadth_Call__<tp>
se_h3k4me3 <- integrate_cross_mark(
  se_h3k4me3, se_h3k27ac,                        # left: promoter domains annotated
  call_fmt_a = "Breadth_Call__%s",               #    by the enhancer-mark calls
  call_fmt_b = "Intensity_Call__%s",
  mark_a = "Breadth", mark_b = "Intensity",
  timepoints = c("0m", "30m", "4h", "16h", "72h"),
  aggregate_ov = "any_super")                    # Super if any overlapping b-domain is Super
rowData(se_h3k4me3)$Intensity_Call__4h           # <mark_b>_Call__<tp> added
rowData(se_h3k4me3)$Intensity_NOverlaps__4h      # # overlapping H3K27ac domains (0 = outside)

• aggregate_ov: "any_super" (any overlapping B domain is super) or "majority" (strict > half of the non-NA overlapping calls are super).
• min_overlap_bp: require at least this many overlapping base pairs.
• Domains with no overlap in B stay NA for the overlaid call, flagging the two-universe “boundary” regions (e.g. promoter-proximal enhancers that the enhancer-mark pipeline’s ROSE-style promoter filter deliberately excludes).
• Promoter-proximal regions absent from an enhancer-mark universe can therefore remain represented in the promoter-mark object.

6.6 Interpretation guardrails

6.7 Input compatibility and parameter guidance

epiPortrait expects normalized BigWig files that share the same genome build, bin size, and normalization strategy (e.g., RPGC-, CPM-, or spike-in-normalized tracks from the same pipeline). Mixing BigWigs from different sources is risky for cross-sample comparisons; within-sample ranking and super-domain calling remain interpretable. Positions with no coverage in a BigWig are treated as zero signal during Intensity integration, so domains with coverage gaps are downward-biased; prefer tracks with complete, consistent coverage across the analyzed region set and check the ZeroFraction / NegativeFraction QC columns from check_signal_compatibility().

CPM and RPGC normalize library scale but do not correct a true global change in histone-mark occupancy. Use spike-in controls when absolute global changes are part of the biological claim, or restrict interpretation to relative within-dataset ranks. Native breadth also depends on sequencing depth, peak caller, and caller parameters; keep these choices constant across samples.

preset <- get_mark_preset("H3K27me3")

# upstream broad-domain caller output imported as GRanges (e.g. epic2/SICER)
domains <- rtracklayer::import("H3K27me3_domains.bed")
# stitch_distance = 0 is a sentinel: keep the upstream caller's geometry
if (preset$stitch_distance > 0) {
  domains <- stitch_epi_peaks(domains, stitch_distance = preset$stitch_distance)
}

se <- build_portrait_matrix(sample_sheet, consensus_peaks = domains)
se <- call_super_domains(se, feature = "Intensity", mode = "per_group")
se <- call_super_domains(se, feature = "Breadth", mode = "per_group")
se <- combine_superdomain_calls(se, group_var = "Condition")

# continuous width remodeling (descriptive, not a significance test)
se <- compute_width_transition(se, ref_group = "Control", target_group = "Treatment")

# repressive display terminology (Extended-Domain / Dual-Extreme)
plot_domain_landscape(se, group = "Control", mark = "H3K27me3")

6.8 Parameter guidance

• Start with method = "elbow" (default) for the main analysis; use method = "tangent" with log_transform = FALSE only for an optional H3K27ac / ROSE-style benchmark (a geometric variant, not a bit-exact ROSE reproduction).
• Use n_bootstrap (e.g., 200) for a cutoff stability interval (a resampling-stability diagnostic, not a classical CI) and success rate. For Intensity the ranked set is the candidate domains; for Breadth it is each replicate’s eligible native PeakWidth distribution (a per-replicate interval is stored in get_call_provenance(se, "Breadth")$replicates[[s]]). Bootstrap resampling is performed on the raw scale with the same transform scale as the main analysis. Its information content shrinks with the size of the ranked set; at genome-wide scale, check per-replicate cutoff consistency (and, for Breadth, peak-set perturbation) instead.
• For multi-condition studies, use mode = "per_group" for the main analysis, then compare_superdomains(). Use mode = "global_consensus" (alias "consensus") only when the goal is to identify domains reproducibly super in every condition; consensus support is computed per condition group (condition-aware).
• Choose support_rule deliberately: "majority" (default, needs > n/2 replicates) is stricter than "fraction" 0.5 for n=2 (needs 2/2 vs 1/2).

se_boot <- call_super_domains(example_se, feature = "Intensity",
                              n_bootstrap = 25, seed = 1)
metadata(se_boot)$superdomain_calls$Intensity$cutoff_stability_interval

## NULL