--- title: "From Signal Tracks to Replicate-Supported Epigenomic Domain Phenotypes" author: "Ying ZHANG" date: "`r Sys.Date()`" package: epiPortrait output: BiocStyle::html_document: toc: true toc_depth: 2 fig_width: 7 fig_height: 5 vignette: > %\VignetteIndexEntry{From Signal Tracks to Replicate-Supported Epigenomic Domain Phenotypes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) library(SummarizedExperiment) library(GenomicRanges) # rtracklayer BigWig I/O is unsupported on Windows; the core workflow uses the # shipped tiny BigWigs on Unix/macOS and falls back to example_se on Windows. use_bw <- .Platform$OS.type != "windows" ``` # Introduction Epigenomic domains are often summarized by signal enrichment, yet signal magnitude is only one aspect of domain organization. A focal, high-amplitude H3K27ac region and a broad, moderately enriched region may carry similar integrated signal while having different signal geometries. Likewise, broad H3K4me3 or H3K27me3 domains can change in measured breadth without a proportional change in signal. A single ranking statistic can obscure these patterns. Their biological meaning requires mark- and experiment-specific validation. `epiPortrait` represents these measurements separately. It distinguishes **integrated signal magnitude** from **replicate-specific native domain breadth**, retains within-domain signal dispersion as a complementary structural descriptor, and aggregates evidence only after each biological replicate has been evaluated independently. The resulting domain portrait is **not** intended to replace peak calling or differential-enrichment statistics. It provides a downstream phenotype layer for asking whether a domain is predominantly intensity-extreme, breadth-extreme, jointly extreme, or undergoing a transition between these states. ## Design principles The design of `epiPortrait` rests on five principles: | # | Principle | |:--|:----------| | 1 | **Shared coordinates are for comparison, not Breadth truth.** The consensus candidate-domain universe is a comparison frame; it never defines a sample's breadth. | | 2 | **Breadth comes from replicate-specific native geometry.** Each replicate's own peak/domain calls carry the breadth evidence. | | 3 | **Intensity and Breadth are classified separately.** They are not merged into one rank, but they may remain correlated because integrated signal depends partly on interval length. | | 4 | **Replicate evidence is retained before group aggregation.** Per-replicate calls are aggregated with explicit support rules; `Uncertain` is an abstention when evidence is insufficient. | | 5 | **Biological interpretation depends on the histone mark.** The calculations are shared across marks, but the labels do not imply the same biology for active and repressive marks. | ## A domain can change in different ways Across conditions, a domain's architecture state can remodel in distinct directions: | Ref | Target | Interpretation | |:----|:-------|:---------------| | Typical | Intensity-Super | acquisition of a relative intensity-extreme state | | Typical | Breadth-Super | acquisition of a relative broad state | | Typical | Dual-Super | acquisition of both relative states | | Intensity-Super | Dual-Super | additional broad-state call | | Breadth-Super | Dual-Super | additional intensity-state call | | Breadth-Super | Typical | loss of the relative broad-state call | 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 | Tool / family | Primary question | Where epiPortrait differs | |:--------------|:-----------------|:--------------------------| | MACS2 / SICER / epic2 | Where are enriched peaks/domains? | epiPortrait starts after domain calling | | DiffBind / csaw | Where does enrichment differ statistically? | epiPortrait describes domain phenotype and remodeling mode | | ROSE | Which stitched enhancers are signal-extreme? | epiPortrait separates Intensity from native Breadth and adds replicate support | | ChIPseeker | Where is a peak relative to genes/features? | epiPortrait uses annotation after quantitative phenotyping | | deepTools | How do signal tracks look across regions/samples? | epiPortrait converts track measurements into domain-level phenotypes | `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. # Installation ```{r install, eval = FALSE} if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c( "GenomicRanges", "SummarizedExperiment", "GenomicFeatures", "org.Hs.eg.db", "org.Mm.eg.db" )) if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes") remotes::install_github("zying106/epiPortrait", build_vignettes = TRUE) ``` # Quick Start with Example Data `epiPortrait` ships a simulated `SummarizedExperiment` (`example_se`) with 500 domains across 6 samples (3 Control, 3 Treatment). Its assays contain Intensity, SignalDispersion, NativeMaxPeakWidth, NativeOccupiedWidth, and NativePeakCount. The canonical Breadth call is derived from replicate-specific native peak widths rather than from a single domain-level assay. ```{r quickstart} library(epiPortrait) data(example_se) example_se ``` ```{r quickstart2} # Normalization is off by default: input BigWigs should already be # quantitatively comparable (CPM/RPGC/spike-in). Post-hoc TotalSignal # rescaling would assume conserved total signal and may remove real global # biological shifts. se <- normalize_portrait(example_se, method = "None") # Main analysis for a multi-condition study: per_group mode, default elbow # inflection and default log convention. se <- call_super_domains(se, feature = "Intensity", mode = "per_group", group_var = "Condition", verbose = FALSE) table(rowData(se)$Intensity_Call__Control) ``` # Module 1: Building the Portrait Matrix The rest of this vignette runs on a tiny synthetic BigWig dataset shipped in `inst/extdata` (4 samples: 2 Control, 2 Treatment; 3 domains on chr21), which exercises the real BigWig I/O and portrait-matrix construction path. Subsequent inflection-based calling is demonstrated on the 500-domain `example_se`, because three domains do not provide an informative ranked distribution for cutoff estimation. ## 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 ```{r consensus} 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 ``` ## 1.3 Optionally stitch into macro-domains ```{r stitch} 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.4 Build the portrait matrix ```{r build} 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 ``` 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. ### Custom features You can extend the portrait with arbitrary per-domain statistics. Example pattern for real BigWig inputs: ```r 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)) } ) ) ``` # Module 2: Normalization ```{r normalize} se <- normalize_portrait(example_se, method = "None") ``` | Method | Use case | |:-------|:---------| | `None` (default) | Recommended — input BigWigs already quantitatively comparable (CPM/RPGC/spike-in) | | `TotalSignal` | Rescale to mean domain-set total signal; **assumes conserved total signal and may remove genuine global biological shifts** | | `Quantile` | Force identical distributions across samples (limma); not recommended for between-condition comparison | | `TMM` | Count-based scaling (`force_TMM = TRUE` for continuous signal); advanced | Row-wise Z-score is intentionally **not** an analysis normalization method in epiPortrait: it would destroy the cross-domain magnitude ranking that Super calling relies on. Plotting / PCA layers perform display-only scaling internally. # Module 3: Super-Domain Calling `call_super_domains()` ranks domains by a feature and separates super-domains from typical domains via a hockey-stick inflection point. ```{r calling} se <- call_super_domains(example_se, feature = "Intensity", method = "elbow", verbose = FALSE) table(rowData(se)$Intensity_Domain_Type) ``` ## 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](https://doi.org/10.1016/j.cell.2013.03.035)). 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. ```{r hockey-methods} 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) ``` **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](https://doi.org/10.1016/j.cell.2013.03.035)), 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. ```{r per-group} 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) ``` ## 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: ```{r breadth-calling} se_b <- call_super_domains(example_se, feature = "Breadth", mode = "per_group", group_var = "Condition", verbose = FALSE) table(rowData(se_b)$Breadth_Call__Control) ``` 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): ```{r combine} 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) ``` # Module 4: Condition Transitions `compare_superdomains()` compares super-domain status between conditions using `mode = "per_group"` calls (which are replicate-aware within each group): ```{r transitions} se_t <- call_super_domains(example_se, feature = "Intensity", mode = "per_group", group_var = "Condition", verbose = FALSE) se_t <- compare_superdomains(se_t, group_var = "Condition", ref_group = "Control", target_group = "Treatment", feature = "Intensity") table(rowData(se_t)$Intensity_Relative_Transition) ``` Transition classes depend on `cutoff_scope`: - **`"relative"`** (default): each group is called with its own cutoff. The output classes `Relative_Prominence_Up` / `Relative_Prominence_Down` express a *relative rank-state* change and **must not** be interpreted as absolute signal gain/loss — a domain whose signal is unchanged can still switch class if the global background of the treatment shifts the group cutoff. - **`"reference"`** / **`"pooled"`**: a common cutoff is applied to both groups, producing common-scale state `Gain` / `Loss` labels. These labels still do not establish an absolute biochemical change and require quantitatively comparable BigWigs. ```{r transitions-reference} se_t2 <- compare_superdomains(se_t, group_var = "Condition", ref_group = "Control", target_group = "Treatment", feature = "Intensity", cutoff_scope = "reference") table(rowData(se_t2)$Intensity_Transition__reference) ``` `Uncertain` is assigned when a per-group call was not reliable (no-call). # Module 5: Annotation & Visualization ## 5.1 Domain-aware annotation ```{r annotate} 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) # long-format links: metadata(se)$domain_gene_links candidates <- get_domain_genes(se) head(candidates) ``` ## 5.2 Hockey-stick plot ```{r hockey, fig.width = 8, fig.height = 5} 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: ```{r portrait, fig.width = 7, fig.height = 5} plot_domain_feature_profile(example_se, peak_id = rownames(example_se)[1], group_var = "Condition") ``` ## 5.4 PCA and correlation QC ```{r qc-plots, fig.width = 8, fig.height = 5} 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: ```{r academic-viz, fig.width = 8, fig.height = 6} 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: ```{r save-analysis} 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")]) prov <- get_call_provenance(se2, "Intensity") prov[c("method", "log_transform_used", "cutoff", "call_status")] 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: ```{r differential} se_diff <- analyze_differential_domains( se_viz, feature = "Intensity", ref_group = "Control", target_group = "Treatment") table(rowData(se_diff)$Intensity_DiffStatus) # Gain / Loss / NS ``` Paired, blocked, and multi-factor designs are supported through `design` and a valid design coefficient or contrast: ```{r differential2, eval=FALSE} 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: ```{r differential-viz} 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: ```{r bedpe-score} 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 get_domain_genes(se_bp, rank_by = "bedpe_score") # dominant-loop-first prioritization ``` - `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: ```{r per-sample-transitions} 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 head(rowData(se_ps)$Intensity_SampleTransition__Control_1_vs_Treatment_1) ``` Labels match `compare_superdomain_classes()` semantics: `Persistent_`, `_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. # Module 6: Mark-Specific Interpretation The same calculations can be applied across marks, but the biological meaning of the resulting labels changes with the mark. ## 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](https://doi.org/10.1016/j.cell.2014.06.027)). 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: | Pattern | Design | What epiPortrait reports | |:--------|:-------|:-------------------------| | Disease-state remodeling | Normal → Cancer | class transitions and width remodeling | | Drug / perturbation response | Vehicle → Drug, WT → KO | per-condition domain states and transitions | | Therapy resistance | Sensitive → Resistant | acquisition / loss of super states | | Differentiation / lineage | State A → State B | state reorganization along the trajectory | | Cohort recurrence | Independent patients / cell-line models | cross-sample recurrence of a domain phenotype | 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: ```{r cross-mark, eval=FALSE} # se_h3k27ac carries Intensity_Call__; se_h3k4me3 carries Breadth_Call__ 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 # _Call__ 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 | Claim | Requirement / caveat | |:------|:----------------------| | Intensity (between-condition magnitude) | BigWigs must be quantitatively comparable (CPM/RPGC/spike-in) | | Breadth | Depends on the upstream native peak/domain caller and boundary quality | | Dual | A conjunction of the two axes, not a new independent measurement | | SignalDispersion | Describes spatial signal organization, not physical chromatin conformation | | Annotation | Nearest / overlapping genes are candidates, not causal targets | | Cohort support | Independent patients are recurrence evidence, not conventional experimental replicates | | Transition | Default `"relative"` scope is a rank-state change, not absolute biochemical gain/loss | ## 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. ### 6.7.1 Mark-aware preset workflow (e.g. H3K27me3) `get_mark_preset()` is a configuration helper: it returns mark-appropriate defaults (for broad repressive marks, `stitch_distance = 0`, i.e. no extra stitching of the already caller-defined domains, and repressive display terminology) but the core pipeline is unchanged. A `stitch_distance` of `0` is a sentinel meaning "do not re-stitch": only call `stitch_epi_peaks()` when the preset returns a positive distance. ```r 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") ``` The internal classes remain `Intensity-Super` / `Breadth-Super` / `Dual-Super`; the `mark` argument only changes surface labels. ## 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). ```{r bootstrap} se_boot <- call_super_domains(example_se, feature = "Intensity", n_bootstrap = 25, seed = 1) metadata(se_boot)$superdomain_calls$Intensity$cutoff_stability_interval ``` # Troubleshooting | Problem | Likely cause | Fix | |:--------|:-------------|:----| | `Feature 'X' not found` | Feature name typo | Use `assayNames(se)` / `colnames(rowData(se))` | | `Breadth requires native peaks` | `feature = "Breadth"` without `peak_path` in the sample sheet | Re-run `build_portrait_matrix()` with per-sample native peak files | | Negative-signal error | BigWig contains negative values | Prefer a compatible non-negative track. Use `negative_policy = "clip_zero"` only when negative values are known technical artefacts, and report that choice. | | `Constant feature values` | All domains identical for a feature (e.g., `IntervalWidth` in a single-sample set) | Use a dynamic feature like `Intensity` | | Super count unexpectedly high | log-transformed scale + tangent | Try `log_transform = FALSE` | | No reliable inflection (`no_call`) | Curve lacks a clear break | Inspect the hockey-stick plot. If a fixed top fraction is scientifically justified, set `quantile_cutoff` and report it; this option does not assess cutoff quality. The default data-driven path retains `no_call` and never switches automatically. | # Session Info ```{r session} sessionInfo() ```