Purpose and signature

get_data() reads an HDF5-backed contact map, extracts a genomic interval, rebins contacts to the requested resolution, and returns a complete long-format bin-pair grid with optional genomic annotations.

get_data(
    file_path, chr, start, end, resolution,
    genome_package = NULL,
    acc_wig = NULL, chain_file = NULL,
    te_granges = NULL
)

Arguments

file_path

Path to a .cool file or supported .h5 schema. Binary .hic is not HDF5 and is rejected; convert it to .cool or extract counts with strawr::straw(). For .mcool, select/export the required resolution as .cool.

chr

Chromosome name exactly as represented by the input and annotations, for example "chr4".

start

One-based start coordinate of the requested genomic interval.

end

End coordinate of the requested interval. It must be greater than or equal to start.

resolution

Positive target bin width in base pairs. Native bins are merged or sliced to construct the requested grid.

genome_package

Optional installed BSgenome package name used to calculate GC content. With NULL, GC is returned as NA.

acc_wig

Optional accessibility bedGraph/wig path. Accessibility is imported only when chain_file is also supplied.

chain_file

Optional liftOver chain used to map the accessibility track into the contact-map genome.

te_granges

Optional TE annotation path readable by rtracklayer::import() or an existing GRanges object.

Internal processing

The function checks required packages, rejects unsupported .hic input, dispatches between supported Cooler and interval/matrix HDF5 schemas, filters both contact ends to the requested chromosome and interval, mirrors an upper-triangle Cooler representation, rebins counts, creates every bin pair, and fills unobserved pairs with zero counts. When the optional resources are available it then calculates combined-bin GC, mean accessibility, and summed TE overlaps.

The returned data frame contains start, end.i., start.j., end, chrom, GC, ACC, TES, and interactions. Coordinate columns describe both ends of each pair; interactions contains non-negative observed counts.

Example

cool_file <- system.file(
    "extdata", "BG3_WT_merged_hic_matrix_chr4_100Kb.cool",
    package = "HiCPotts"
)

imported_contacts <- get_data(
    file_path = cool_file,
    chr = "chr4",
    start = 1L,
    end = 500000L,
    resolution = 100000L,
    genome_package = NULL
)

head(imported_contacts)
#>    start end.i. start.j.   end chrom GC ACC TES interactions
#> 1      1  1e+05        1 1e+05  chr4 NA  NA  NA        19023
#> 2 100001  2e+05        1 1e+05  chr4 NA  NA  NA         9927
#> 3 200001  3e+05        1 1e+05  chr4 NA  NA  NA         2259
#> 4 300001  4e+05        1 1e+05  chr4 NA  NA  NA         1665
#> 5 400001  5e+05        1 1e+05  chr4 NA  NA  NA         1086
#> 6      1  1e+05   100001 2e+05  chr4 NA  NA  NA         9927

Purpose and signature

process_data() validates a long-format interaction table and reshapes counts and covariates into one or more square matrices accepted by run_chain_betas().

process_data(
    data, N, scale_max = NA_real_,
    standardization_y = FALSE,
    pad_with_zero = FALSE,
    mirror = FALSE
)

Arguments

data

Data frame containing start, end, interactions, GC, TES, and ACC. start.j. is used for the second-bin start when present. Missing values are rejected.

N

Positive integer number of rows and columns in each output matrix. Full blocks contain N^2 rows.

scale_max

Upper scaling value used only when standardization_y = TRUE.

standardization_y

If FALSE, preserve raw counts. If TRUE, min-max scale and round them. Raw non-negative integer counts are required for the package’s Poisson/NB likelihoods, so FALSE is recommended for inference.

pad_with_zero

If TRUE, append zeros to an incomplete final block. This creates unobserved zero cells and should be used only when that scientific assumption is justified.

mirror

If TRUE, interpret input as one complete triangle and reconstruct the opposite triangle. Use this for triangular symmetric Hi-C input, not for intentionally ordered/asymmetric data.

Internal processing and return value

The function checks columns, missingness, count validity and block dimensions. It optionally mirrors triangular input or pads an incomplete block, calculates distance as abs(start.j. - start) when start.j. exists (otherwise the legacy abs(end - start) definition), and reshapes values in the package’s column-major order.

It returns:

• x_vars: a named list with distance, GC, TES, and ACC; each entry is a list of N by N matrices;
• y: a list of matching N by N count matrices.

Examples

# A small reproducible symmetric dataset with simulated annotations.
set.seed(4921)
bins <- seq.int(1L, by = 100000L, length.out = 5L)
hic <- expand.grid(start = bins, start.j. = bins)
hic$end.i. <- hic$start + 99999L
hic$end <- hic$start.j. + 99999L
hic$chrom <- "chr4"
symmetric_values <- function(values) {
    mat <- matrix(values, 5L, 5L)
    mat[lower.tri(mat)] <- t(mat)[lower.tri(mat)]
    as.vector(mat)
}
hic$GC <- symmetric_values(runif(25L, 0.3, 0.7))
hic$ACC <- symmetric_values(runif(25L))
hic$TES <- symmetric_values(rpois(25L, 2))
hic$interactions <- symmetric_values(rpois(25L, 8))

processed <- process_data(
    hic,
    N = length(unique(hic$start)),
    standardization_y = FALSE,
    mirror = FALSE
)
#> Genomic distance: using abs(start.j. - start)

# One triangle of the same matrix gives the same processed counts.
triangular <- hic[hic$start <= hic$start.j., ]
processed_triangular <- process_data(
    triangular,
    N = 5L,
    standardization_y = FALSE,
    mirror = TRUE
)
#> Genomic distance: using abs(start.j. - start)

y_matrix <- processed$y[[1L]]
x_matrices <- lapply(processed$x_vars, function(x) x[[1L]])
stopifnot(identical(y_matrix, processed_triangular$y[[1L]]))

Purpose and signature

run_chain_betas() is the only exported fitting entry point. It validates the data contract, prepares priors and a shared ABC tolerance, constructs independently seeded starting allocations, runs the native sampler, optionally relabels the components, and returns either raw chains or the robust diagnostic workflow.

run_chain_betas(
    N, x_vars, y, dist = "ZIP", gamma_start = 0.3,
    iterations = 20000L, burnin = NULL, n_chains = 1L,
    seeds = NULL, robust = FALSE, initialization = "auto",
    theta_start = NULL, size_start = NULL,
    use_data_priors = TRUE, user_fixed_priors = NULL,
    gamma_prior_shape = c(1, 1), epsilon = NULL,
    abc_potts_sweeps = 0L, abc_sim_reps = 4L,
    distance_metric = "manhattan", relabel = TRUE,
    mc_cores = 1L, verbose = FALSE, progress_interval = 50L,
    slope_sd_standardized = 0.5,
    abc_epsilon_quantile = 0.10,
    gamma_update_interval = 5L,
    z_probability_burnin = NULL,
    mcse_stop = TRUE, mcse_min_iterations = 10000L,
    mcse_check_interval = 500L,
    mcse_relative_threshold = 0.05,
    diagnostic_control = list(), gamma_prior = NULL
)

Data, model and chain arguments

N

Dimension shared by every square response and covariate matrix.

x_vars

Named list containing distance, GC, TES, and ACC. For one dataset entries may be matrices; for several datasets each entry is a list matching y.

y

One non-negative integer count matrix, or a list of matrices fitted independently.

dist

Count family: "Poisson", "NB", "ZIP", or "ZINB".

gamma_start

Initial Potts spatial-coupling value strictly between zero and one. It is not the gamma prior.

iterations

Positive maximum number of production MCMC updates per chain. It is exact when mcse_stop = FALSE.

burnin

Initial production iterations excluded by robust diagnostics. NULL uses half of iterations, capped at 5,000.

n_chains

Number of independent chains per dataset. Four or more are recommended for convergence assessment.

seeds

Integer seed per chain. If omitted, seq_len(n_chains) is used; if n_chains is omitted, seed-vector length determines it.

robust

TRUE returns a mode-screened hicpotts_robust_fit with diagnostics. FALSE returns the raw fitted chains.

initialization

"auto", or one method per chain from "likelihood_informed", "count_quantile", "distance_adjusted", "noise_anchored_random", and "random".

theta_start

Starting zero-inflation probability for ZIP and ZINB. Omit for Poisson and NB.

size_start

Positive length-three NB2 dispersion starting vector for NB and ZINB. Omit for Poisson and ZIP.

Prior and ABC arguments

use_data_priors

With robust = TRUE, run a separate multi-chain empirical-Bayes pilot, pool one prior, and freeze it for production. With FALSE, use user_fixed_priors or robust automatically constructed scaled priors.

user_fixed_priors

Named component1/component2/component3 list used when use_data_priors = FALSE. Each component supplies meany, meanx1:meanx4, sdy, and sdx1:sdx4. No empirical-Bayes pilot is needed.

gamma_prior_shape

Positive length-two vector containing the Beta prior shapes for gamma. c(1, 1) is uniform.

epsilon

Positive ABC kernel bandwidth. NULL performs one deterministic prior-predictive calibration shared by every chain of the dataset.

abc_potts_sweeps

Gibbs sweeps per auxiliary Potts field. Zero selects max(100, 100 * N^2 / 400). Reducing it speeds fitting but can under-equilibrate fields and bias gamma upward.

abc_sim_reps

Auxiliary fields averaged per gamma proposal. Cost is approximately linear in this value.

abc_epsilon_quantile

Quantile strictly between zero and one used during automatic ABC tolerance calibration.

gamma_update_interval

MCMC iterations between gamma proposals. Larger values are cheaper but provide fewer gamma transitions.

slope_sd_standardized

Standardized slope prior SD used only when the robust fixed-prior workflow constructs priors automatically. Reported coefficients remain on the original log1p covariate scale.

Execution, probability and diagnostic arguments

distance_metric

Compatibility argument retained for older code. It is ignored.

relabel

Apply the biological component-2/3 orientation before returning the fit.

mc_cores

Number of chains that may execute concurrently. Windows uses PSOCK workers; Unix-like systems use forked workers.

verbose

Print sampler progress.

progress_interval

Iterations between progress messages when verbose = TRUE.

z_probability_burnin

Iteration after which sampled latent states contribute to membership probabilities. NULL uses burnin.

mcse_stop

Permit early stopping when every monitored parameter satisfies the relative MCSE rule. Set FALSE for an exact iteration count.

mcse_min_iterations

Earliest iteration at which MCSE stopping can occur.

mcse_check_interval

Iterations between MCSE evaluations.

mcse_relative_threshold

Required maximum batch-means MCSE divided by posterior SD.

diagnostic_control

Named overrides for robust reliability and mode-screen thresholds; listed below. These affect reporting and selection, not the posterior target.

gamma_prior

Deprecated alias for gamma_start. Do not supply both.

Recognised diagnostic_control entries are:

minimum_component_cells

Default: 0. Minimum cells required in each component.

minimum_ess

Default: 200. Minimum bulk and tail effective sample size.

maximum_rhat

Default: 1.01. Maximum split-R-hat.

gamma_boundary_tolerance

Default: 0.01. Distance from zero or one treated as a boundary.

maximum_gamma_boundary_fraction

Default: 0.95. Largest permitted retained gamma fraction near a boundary.

minimum_gamma_unique

Default: 20. Minimum distinct retained gamma values per chain.

mode_agreement_threshold

Default: 0.8. Allocation agreement used when screening replicated modes.

minimum_mode_chains

Default: 2. Minimum chains required to retain a replicated mode.

Internal processing and return value

The function first standardizes one-dataset and multiple-dataset layouts. It validates counts, covariates and starting values, chooses the requested initialization for each seed, and calibrates one shared ABC tolerance when epsilon = NULL.

With robust empirical-Bayes priors, a pilot of at most 5,000 iterations per chain estimates one shared regression-prior object; pilot draws never enter production summaries. With fixed user priors there is no prior-estimation pilot, but the same deterministic ABC calibration is still shared across all chains. Production chains therefore target the same posterior in either mode.

For robust = FALSE, the return value is a list of raw fits. A raw fit includes component coefficient chains, gamma, family-dependent theta/size draws, final and checkpoint allocations, membership frequencies, batch summaries, sampler settings and performance information. For robust = TRUE, the hicpotts_robust_fit additionally contains fits, all_fits, mode_selection, diagnostics, covariate_diagnostics, priors, empirical_bayes_pilot, and settings.

Example

fit <- run_chain_betas(
    N = nrow(y_matrix),
    x_vars = x_matrices,
    y = y_matrix,
    dist = "ZINB",
    gamma_start = 0.3,
    gamma_prior_shape = c(1, 1),
    theta_start = 0.05,
    size_start = rep(11, 3L),
    iterations = 60L,
    burnin = 20L,
    n_chains = 2L,
    seeds = c(101L, 202L),
    robust = TRUE,
    initialization = "auto",
    use_data_priors = FALSE,
    epsilon = 0.2,
    abc_potts_sweeps = 10L,
    abc_sim_reps = 1L,
    mcse_stop = FALSE,
    mc_cores = 1L,
    verbose = FALSE
)
#> Warning in (function (N, x_vars, y, dist = "ZIP", gamma_prior = 0.3, iterations
#> = 20000L, : Fewer than 15,000 iterations requested; rely on ESS and R-hat
#> before interpreting estimates.
#> Warning in (function (N, x_vars, y, dist = "ZIP", gamma_prior = 0.3, iterations
#> = 20000L, : Fewer than four independent chains requested; four are recommended
#> for general use.

The fixed seed, short chains and explicit ABC settings keep this example quick to rebuild. For a production analysis, use the longer template below and assess the diagnostic output for your dataset. Setting use_data_priors = TRUE enables the shared empirical-Bayes pilot.

Purpose, argument and return value

relabel_hicpotts(x)

x is one raw HiCPotts fit or a list of raw fits. The function preserves component 1 as low-mean noise and makes one global component-2/3 orientation decision for each complete retained chain. Component 3 is the elevated component whose standardized covariate slopes most closely resemble component 1 and whose intercept is above component 1; component 2 is the remaining unrestricted true-signal component.

Internally it scores both possible component-2/3 branches, decides whether a global swap is needed, and applies the same permutation to coefficient chains, dispersion, allocations, checkpoints, membership probabilities, batch summaries and stored priors. The return value has the same outer structure as x and records the permutation, relationship probability and relabelling rule. run_chain_betas(relabel = TRUE) normally performs this step automatically.

oriented_chains <- relabel_hicpotts(fit$all_fits)
oriented_chains[[1L]]$relabel_rule
#> [1] "Component 1 retained as low-mean noise. Components 2/3 left in the internally selected orientation; component 3 is tested as the elevated component most consistent with component 1's standardised slopes and component 2 remains unrestricted. Component-3 relationship probability = 0.8065; no hard probability threshold is imposed by the package."
oriented_chains[[1L]]$noise_relationship_probability
#> [1] 0.8064516

Purpose and signature

This is the official classification function. It pools sampled post-burn-in latent-state frequencies; it does not refit, evaluate posterior-mean densities, or add new smoothing.

classify_hicpotts(
    fit, data = NULL,
    component_names = hicpotts_component_definition()$label,
    use = c("selected", "all"), relabel = TRUE,
    min_draws = 100L,
    reflect = c("auto", "always", "never")
)

Arguments

fit

Raw fit, list of chains, robust fit, or block-aware fit.

data

Optional original data frame in the same column-major cell order. When omitted, lattice indices are returned.

component_names

Three distinct output labels. Changing them changes wording only, not component identities.

use

For robust fits, pool mode-screened "selected" chains or "all" retained chains. Selected results are conditional on that replicated mode.

relabel

Apply biological chain orientation and cross-chain alignment before pooling.

min_draws

Minimum pooled post-burn-in allocation draws. Too few draws cause an error rather than a final-state fallback.

reflect

"auto" averages mirrored probabilities only when symmetric input is confirmed; "always" forces it and "never" disables it.

Internal processing and return value

The function selects the requested robust chains, applies biological relabelling, aligns allocations, sums the stored membership frequencies, and optionally averages mirrored cells. It returns the original columns or lattice indices plus prob1, prob2, prob3, the MAP component and label, winning probability, probability margin, normalized entropy, and classification. Attributes record chains, draws, scope and reflection behavior.

classification <- classify_hicpotts(
    fit,
    data = hic,
    use = "selected",
    min_draws = 5L,
    reflect = "auto"
)
#> Warning in classify_hicpotts(fit, data = hic, use = "selected", min_draws = 5L,
#> : Review the recorded diagnostic criteria in fit$diagnostics$reliability_flags.
#> Criteria below their configured thresholds: replicated_allocation_mode,
#> gamma_movement, effective_sample_size, split_rhat, independent_chains.

table(classification$classification)
#> 
#>        noise       signal false signal 
#>            8           15            2
head(classification[c("prob1", "prob2", "prob3", "probability_margin")])
#>    prob1   prob2   prob3 probability_margin
#> 1 0.0875 0.50000 0.41250            0.08750
#> 2 0.3000 0.36250 0.33750            0.02500
#> 3 0.3000 0.53125 0.16875            0.23125
#> 4 0.4625 0.35000 0.18750            0.11250
#> 5 0.3875 0.25000 0.36250            0.02500
#> 6 0.3000 0.36250 0.33750            0.02500

Purpose and signature

This function is a secondary parameter-based sensitivity analysis. The official allocation remains classify_hicpotts(). The integrated method evaluates membership probabilities over retained posterior draws; the plug-in method evaluates one set of posterior-mean parameters.

compute_HMRFHiC_probabilities(
    data = NULL, chain_betas = NULL, iterations = NULL,
    dist = "ZINB", max_interactions = NA_integer_,
    consistent_dist = FALSE, relabel = FALSE,
    N = NULL, potts_iterations = 5L,
    method = c("integrated", "plugin"), n_draws = 200L,
    component_definition = hicpotts_component_definition()
)

Arguments

data

Required long-format data frame with start, end, interactions, GC, TES, and ACC; may also be block data accepted by combine_hicpotts_blocks().

chain_betas

Raw fit, chain list, robust fit, or block-aware fit.

iterations

Compatibility argument that must be a scalar of at least two. Retained indices are derived from each chain’s actual stored length, so early-stopped chains are safe.

dist

Fitted family: Poisson, NB, ZIP, or ZINB. It must match the fit.

max_interactions

Optional integer cap applied to counts before density evaluation; NA applies no cap.

consistent_dist

FALSE matches the fitted model: zero inflation applies only to component 1. TRUE applies the zero-inflated family to all components and is retained only for backwards comparison.

relabel

Orient chains before extracting parameters.

N

Optional lattice dimension. Supply it with N^2 == nrow(data) to include the mean-field Potts neighbour term. Without it, spatial-free probabilities are returned with a warning.

potts_iterations

Mean-field spatial smoothing sweeps when N is valid.

method

"integrated" averages normalized probabilities across draws; "plugin" evaluates posterior-mean parameters and understates uncertainty.

n_draws

Posterior draws used by method = "integrated"; runtime is approximately linear in this value.

component_definition

Three-row data frame with component, label, and definition. Custom values rename output only.

Internal processing and return value

After validating and optionally relabelling the fits, the function transforms the four covariates with log1p, selects retained draws, evaluates the component-specific count densities, optionally iterates the mean-field Potts term over four-neighbour lattice adjacency, and normalizes each row. It returns the supplied data plus prob1, prob2, and prob3.

parameter_probabilities <- compute_HMRFHiC_probabilities(
    data = hic,
    chain_betas = fit,
    iterations = 60L,
    dist = "ZINB",
    relabel = TRUE,
    N = nrow(y_matrix),
    potts_iterations = 5L,
    method = "integrated",
    n_draws = 10L
)
#> Genomic distance: using abs(start.j. - start)

Signature and arguments

diagnose_hicpotts_fit(
    fit, burnin = NULL, ci_level = 0.95,
    prob_result = NULL, minimum_component_cells = 0L,
    minimum_ess = 200, maximum_rhat = 1.01,
    minimum_chains = 4L,
    gamma_boundary_tolerance = 0.01,
    maximum_gamma_boundary_fraction = 0.95,
    minimum_gamma_unique = 20L,
    covariate_diagnostics = NULL, relabel = TRUE
)

fit

Raw fit, independent-chain list, robust fit, block-aware fit, or compatible diagnostic input.

burnin

Initial draws discarded; NULL uses the fit setting when available or half the chain.

ci_level

Posterior interval level.

prob_result

Optional probability/classification data frame used to calculate classification entropy and occupancy summaries.

minimum_component_cells

Minimum internal cells per component. Zero records occupancy without rejecting an empty component.

minimum_ess

Required bulk and tail ESS for parameter resolution.

maximum_rhat

Largest permitted split-R-hat.

minimum_chains

Recommended minimum independent chain count.

gamma_boundary_tolerance

Distance from zero or one treated as gamma’s boundary region.

maximum_gamma_boundary_fraction

Maximum retained gamma fraction permitted in the boundary region for every chain.

minimum_gamma_unique

Minimum distinct retained gamma values in every chain.

covariate_diagnostics

Optional precomputed conditioning report. Robust fits reuse their stored report; absence means no demonstrated problem, not proven good conditioning.

relabel

Orient components before pooling chains.

Internal processing and return value

The function extracts comparable retained draws, checks shared posterior targets and coefficient scales, computes posterior estimates, intervals, sign probabilities, bulk/tail ESS and split-R-hat, and evaluates gamma movement, coefficient movement, occupancy, chain count, relabelling and conditioning.

It returns parameters, reliability_flags, component_occupancy, classification_entropy, mcse_status, gamma_diagnostics, beta_mixing_status, noise_relationship_probabilities, the fits used, and warnings. Interpret each row of reliability_flags; there is deliberately no single overall pass/fail word.

diagnostics <- diagnose_hicpotts_fit(
    fit,
    prob_result = classification,
    minimum_ess = 200,
    maximum_rhat = 1.01,
    minimum_chains = 4L
)

diagnostics$reliability_flags[c("criterion", "passed", "threshold")]
#>                  criterion passed
#> 1           mcse_precision   TRUE
#> 2            beta_movement   TRUE
#> 3  common_posterior_target   TRUE
#> 4       gamma_not_boundary   TRUE
#> 5           gamma_movement  FALSE
#> 6      component_occupancy   TRUE
#> 7    effective_sample_size  FALSE
#> 8               split_rhat  FALSE
#> 9       independent_chains  FALSE
#> 10              relabelled   TRUE
#> 11       coefficient_scale   TRUE
#> 12  covariate_conditioning   TRUE
#>                                                        threshold
#> 1      batch-means MCSE/SD threshold before requested iterations
#> 2  >0 retained ordinary beta acceptances per component and chain
#> 3    identical frozen regression prior in every production chain
#> 4                                    <0.95 within 0.01 of 0 or 1
#> 5                         >=20 distinct retained draws per chain
#> 6                                                      >=0 cells
#> 7                                                          >=200
#> 8                                                         <=1.01
#> 9                                                            >=4
#> 10                                                      required
#> 11                                     original manuscript scale
#> 12                              correlation/condition thresholds