# Bayesian Analysis of Hi-C Interactions with HiCPotts

# Introduction

High-throughput chromosome conformation capture (Hi-C) technologies have revolutionized our understanding of 3D genome organization by mapping interactions between genomic loci. However, Hi-C data are inherently noisy and affected by experimental biases such as GC content, transposable elements, and DNA accessibility, which complicate the identification of biologically significant interactions.

The HiCPotts package provides a comprehensive framework for Bayesian analysis of Hi-C interaction data using a Hidden Markov Random Field (HMRF) model. Hi-C is a high-throughput sequencing technique that captures chromatin interactions across the genome, revealing spatial organization. This package models these interactions using a mixture of distributions (Poisson, Negative Binomial, Zero-Inflated Poisson, or Zero-Inflated Negative Binomial) while accounting for covariates, genomic distance, GC content, accessibility, and transposable element (TE) counts. The HMRF framework incorporates spatial dependencies via a Potts model, and the package employs Markov Chain Monte Carlo (MCMC) methods for parameter estimation.

An HMRF models spatial dependencies in a lattice (e.g., a matrix of Hi-C interactions) by assigning each site to one of several latent states (mixture components). The Potts model governs spatial interactions, encouraging neighboring sites to share the same state, controlled by an interaction parameter (\(\gamma\)).

The package uses MCMC to estimate model parameters (regression coefficients (\(\beta\)), zero-inflation parameter (\(\theta\)), dispersion parameter for Binomial distributions, and (\(\gamma\))) by sampling from their posterior distributions. Covariates like genomic distance and GC content adjust for biases in interaction counts. Each Hi-C interaction is modeled as belonging to one of three mixture components following a specified distribution.

The package performance was optimized through C++ implementations using Rcpp and RcppArmadillo. The package also supports parallel processing and flexible distribution choices, making it suitable for large-scale genomic analyses.

# Motivation

Most existing computational methods fail to adequately model the spatial dependencies and overdispersion in Hi-C contact matrices, limiting their ability to distinguish true signals from other components such as noise. The HiCPotts package addresses these challenges by providing a novel Bayesian framework to detect enriched interactions while accounting for experimental biases in Hi-C data. Its integration into Bioconductor is motivated by its robust statistical approach, computational efficiency, and ability to provide insights into bias sources, making it a valuable tool for researchers studying chromatin architecture in diverse biological contexts.

HiCPotts package also extends the existing knowledge of classifying interacting loci into two components to three components (true signal, false signal, and noise) and its focus on bias correction (DNA accessibility, transposable elements) enhance its utility for integrative genomic studies.

# Comparison with Existing Packages

There are several Bioconductor packages that addresses Hi-C data analysis, each with distinct functionalities and scopes. We highlight the uniqueness of the HiCPotts package with other existing packages:

1. diffHic: The diffHic package focuses on detecting differential interactions between biological conditions using the edgeR framework for statistical modeling. diffHic provides methods for read pair alignment, binning, filtering, and normalization of biases (e.g., trended or CNV-driven). While diffHic excels at differential analysis, it does not explicitly model spatial dependencies in Hi-C data or account for overdispersion as HiCPotts does. HiCPotts is better suited for identifying enriched interactions within a Hi-C experiment and exploring bias sources, whereas diffHic is ideal for comparative studies across conditions.

2. HiCcompare and multiHiCcompare: HiCcompare offers joint normalization and difference detection for multiple Hi-C datasets, operating on sparse chromatin interaction matrices. multiHiCcompare extends this to handle multiple groups and replicates using cyclic loess normalization and a general linear model (GLM) based on edgeR. Both packages emphasize comparative analysis and normalization but do not focus on detecting enriched interactions within a dataset or modeling spatial dependencies. HiCPotts’s Bayesian approach and bias correction make it complementary, as it prioritizes significant interaction detection and bias insight over differential analysis.

3. HiCDCPlus: HiCDCPlus enables significant interaction calling and differential analysis for Hi-C and HiChIP data using a negative binomial generalized linear model. It includes tools for topologically associating domain (TAD) and A/B compartment calling, integrating with visualization tools. Like HiCPotts, it calls significant interactions, but HiCPotts’s HMRF-based model and ABC approach provide superior handling of spatial dependencies and computational tractability for chromosome-wide analysis. HiCDCPlus requires GC content information (computable internally), while HiCPotts additionally corrects for transposable elements and DNA accessibility biases.

4. scHiCcompare: Designed for single-cell Hi-C data, scHiCcompare supports imputation, normalization, and differential interaction analysis across single-cell datasets. Its focus on single-cell data makes it distinct from HiCPotts, which targets bulk Hi-C data. HiCPotts’s ability to model overdispersion and spatial dependencies is not replicated in scHiCcompare, which prioritizes single-cell-specific challenges.

Other Tools: Packages like HiCdat and HiCExperiment provide preprocessing, visualization, or data manipulation for Hi-C data but lack the statistical rigor of HiCPotts for interaction detection. HiCdat offers a graphical interface for preprocessing and integrative analysis with other omics data, while HiCExperiment provides data structures for 3C-related experiments. Neither focuses on enriched interaction detection or bias correction like HiCPotts.

# Features of HiCPotts

HiCPotts allows researchers to identify significant intra-chromosomal interactions in Hi-C data while correcting for experimental biases. Its key features include:

1. Use of Zero-Inflated distributions to handle overdispersion.

2. An HMRF-based Bayesian framework with the Potts model for spatial dependency.

3. ABC for computationally efficient handling of the Potts model’s normalizing constant.

4. Bias correction for GC content, transposable elements, and DNA accessibility.

5. Unlike diffHic, HiCcompare, and multiHiCcompare, which focus on differential analysis, or HiCDCPlus, which balances interaction calling and differential analysis, HiCPotts prioritizes enriched interaction detection. Its bias correction and spatial modeling make it a powerful complement to existing tools, enhancing Bioconductor’s suite for 3D genome analysis. HiCPotts can be used for initial interaction detection, followed by diffHic or HiCcompare for differential studies, creating a comprehensive Hi-C analysis workflow.

# Installation

HiCPotts depends on several CRAN and Bioconductor packages.

Install them as follows:

Install CRAN dependencies

    install.packages(c("Rcpp", "RcppArmadillo", "parallel"))

    # Install Bioconductor dependencies
    if (!requireNamespace("BiocManager", quietly = TRUE)) {
      install.packages("BiocManager")
    }
    BiocManager::install(c("rhdf5", "strawr", "rtracklayer", "GenomicRanges", "BSgenome", "Biostrings"))

    # Loading the package:
    # install.packages(HiCPotts)
    library(HiCPotts)

# Workflow Overview

The HiCPotts workflow involves four main steps:

1. Data Loading: Use `get_data()` to read Hi-C contact matrices (`.hic`, `.cool`, `.h5`) and annotate bins with GC content, accessibility, and TE counts.

2. Data Processing: Use `process_data()` to convert the data into (N N) matrices of interactions and covariates.

3. MCMC Simulation: Use `run_chain_betas()` to run MCMC chains, estimating parameters and latent state assignments.

4. Probability Computation: Use `compute_HMRFHiC_probabilities()` to calculate posterior probabilities of component assignments.

We’ll demonstrate this workflow using synthetic Hi-C data for simplicity.

## Step 1: Loading Hi-C Data

The `get_data()` function reads Hi-C contact matrices and annotates bins with covariates. For this example, we simulate a small dataset instead of using real Hi-C files, which require specific file formats and genome annotations.

Simulate a 10x10 Hi-C dataset

In practice, the `get_data()` function is used to load real Hi-C data from .cool, .mcool, or .hic files, similar to the HiCExperiment package. Both convert these files into structured outputs, but they differ in their results. HiCExperiment generates a HiCExperiment object that includes a contact matrix, genomic regions, metadata (e.g., resolution, chromosome), and pairwise interactions. In contrast, `get_data()` produces a data frame from the same file types, incorporating optional calculations (if available) to identify sources of experimental biases, such as GC content or DNA accessibility. Also the function support loading profiles as bigwig or bedgraph and will be imported as GRanges.

## Step 2: Processing Data

The `process_data()` function converts the data frame generated through the `prepare_data()` function into a list of (N \(\times\) N) matrices for interactions (`y`) and covariates (`x_vars`), optionally scaling interaction counts. Also, if a HiCExperiment object is already available and the sources of biases are known and organized in a separate data frame, these can be combined into a single data frame using base R functions and the `prepare_data()` function to convert into a list for analysis.

This produces matrices for distance, GC, TEs, ACC, and interaction counts, ready for MCMC.

## Step 3: Running MCMC Simulations

The `run_chain_betas()` function runs MCMC simulations to estimate parameters ((\(\beta\)), (\(\gamma\)), (\(\theta\) (for Zero-Inflated distributions)), size(for Negative binomial distributions)) and latent state assignments ((z)). In the example below, we use the Zero-Inflated Negative Binomial (ZINB) distribution to model the simulated sparse Hi-C data.

The output includes chains for regression parameters (`chains`), the Potts interaction parameter (`gamma`), zero-inflation parameter (`theta`), and dispersion parameters (`size`).

## Step 4: Computing Posterior Probabilities

The `compute_HMRFHiC_probabilities()` function calculates posterior probabilities for each interaction belonging to one of three mixture components, using the MCMC chains.

# Worked Examples: Hi-C and Micro-C Data

The synthetic example above is useful for illustrating the API but is not biologically meaningful. This section walks through two fully realistic analyses — one for Hi-C and one for Micro-C — using the exact same four-step workflow. The Hi-C data files are shipped with the package under `inst/extdata`; the Micro-C paths are illustrative and should be replaced with your own files.

## Example 1: Hi-C Data

We analyse a 100 kb-resolution Hi-C contact map of D. melanogaster chr4 in BG3 cells, integrating DNase-I accessibility and transposable-element annotations for bias correction.

### Step 1 — Locate input files

    # Hi-C contact matrix (.cool, 100 kb, BG3 cells, chr4)
    hic_file   <- system.file("extdata",
                              "BG3_WT_merged_hic_matrix_chr4_100Kb.cool",
                              package = "HiCPotts")

    # DNase-I accessibility bedGraph (dm3 genome coordinates)
    dnase_file <- system.file("extdata",
                              "DNaseI_BG3_gr_chr4.bedGraph",
                              package = "HiCPotts")

    # LiftOver chain — maps dm3 accessibility up to dm6 assembly
    chain_file <- system.file("extdata",
                              "dm3ToDm6_chr4_only.chain",
                              package = "HiCPotts")

    # Transposable-element annotations (GTF, dm6 chr4)
    te_file    <- system.file("extdata",
                              "dm6_TEs_chr4.gtf",
                              package = "HiCPotts")

    stopifnot(file.exists(hic_file),
              file.exists(dnase_file),
              file.exists(chain_file),
              file.exists(te_file))

### Step 2 — Load Hi-C data and compute covariates

    # get_data() does four things in one call:
    #   1. Reads the .cool file and rebins to the requested resolution
    #   2. Calculates GC content from the BSgenome package
    #   3. Lifts over the DNase-I signal (dm3 -> dm6) and aggregates per bin
    #   4. Counts transposable elements overlapping each bin pair
    hic_data <- get_data(
      file_path      = hic_file,
      chr            = "chr4",
      start          = 1,
      end            = 400000,
      resolution     = 200000,
      genome_package = "BSgenome.Dmelanogaster.UCSC.dm6",
      acc_wig        = dnase_file,
      chain_file     = chain_file,
      te_granges     = te_file
    )

    head(hic_data)
    ##    start end.i. start.j.   end chrom        GC      ACC TES interactions
    ## 1      1  2e+05        1 2e+05  chr4 0.3358900 67.78359 276        72158
    ## 2 200001  4e+05        1 2e+05  chr4 0.3385825 66.73802 287        29660
    ## 3      1  2e+05   200001 4e+05  chr4 0.3385825 66.73802 287        29660
    ## 4 200001  4e+05   200001 4e+05  chr4 0.3412750 65.69244 298        89911

### Step 3 — Process data into matrices

    # (end - start) / resolution, rounded up for the edge bin.
    N_hic <- 14

    processed_hic <- process_data(
      hic_data,
      N                 = N_hic,
      scale_max         = 500,
      standardization_y = TRUE
    )

    x_vars_hic <- processed_hic[["x_vars"]]
    y_hic      <- processed_hic[["y"]]

### Step 4 — Run MCMC with a realistic sampler configuration

    # ZINB is recommended for Hi-C because of structural zeros and
    # over-dispersion typical of ligation-based count data.
    set.seed(2025)

    hic_results <- run_chain_betas(
      N                = N_hic,
      gamma_prior      = 0.3,
      iterations       = 5000,
      x_vars           = x_vars_hic,
      y                = y_hic,
      theta_start      = 0.5,
      size_start       = c(1, 1, 1),
      use_data_priors  = TRUE,
      distance_metric  = "manhattan",
      dist             = "ZINB",
      mc_cores         = min(4L, parallel::detectCores() - 1L)
    )

    # run_chain_betas() returns a list with one element per input dataset.
    # Unpack the first (and here only) chain set:
    hic_chains <- hic_results[[1]][["chains"]]
    hic_gamma  <- hic_results[[1]][["gamma"]]
    hic_theta  <- hic_results[[1]][["theta"]]
    hic_size   <- hic_results[[1]][["size"]]

### Step 5 — Compute posterior component probabilities

    # Pass the full results list, not just the chains.
    hic_probs <- compute_HMRFHiC_probabilities(
      data        = hic_data,
      chain_betas = hic_results,
      iterations  = 5000,
      dist        = "ZINB"
    )

    # hic_probs is a data.frame with prob1, prob2, prob3 columns that sum to 1
    # per row. prob1 is the noise/baseline component; prob2 and prob3 are the
    # elevated-signal regimes.
    head(hic_probs[, c("start", "end", "interactions", "prob1", "prob2", "prob3")])

    # Bin-pairs confidently classified as elevated signal (prob2 > 0.9). This value can be increase or reduce based on user's justification.
    genuine_signal <- hic_probs[hic_probs$prob2 > 0.9,
                                c("start", "end", "interactions", "prob2")]
    head(genuine_signal)

## Example 2: Micro-C Data

Micro-C uses micrococcal nuclease (MNase) instead of a restriction enzyme, yielding nucleosome-resolution contact maps (~150 bp). The HiCPotts workflow is identical to Hi-C; practical differences:

- Input file: typically a higher-resolution `.cool` or `.mcool`.
- Resolution: 1 kb – 5 kb rather than 10 – 100 kb.
- Structural zeros: even more abundant at fine resolution — ZINB strongly recommended.
- Iterations: larger lattices need longer chains to converge.

### Step 1 — Locate input files

    # Replace these with the actual paths on your machine.
    microc_file <- "path/to/MicroC_BG3_chr4_1kb.cool"
    dnase_file  <- "path/to/DNaseI_BG3_gr_chr4.bedGraph"
    chain_file  <- "path/to/dm3ToDm6_chr4_only.chain"
    te_file     <- "path/to/dm6_TEs_chr4.gtf"

### Step 2 — Load Micro-C data

    # Aggregate from the native 1 kb to 5 kb to keep the lattice manageable.
    microc_data <- get_data(
      file_path      = microc_file,
      chr            = "chr4",
      start          = 1,
      end            = 1350000,
      resolution     = 5000,
      genome_package = "BSgenome.Dmelanogaster.UCSC.dm6",
      acc_wig        = dnase_file,
      chain_file     = chain_file,
      te_granges     = te_file
    )

    head(microc_data)

### Step 3 — Process into matrices

    # 1,350,000 / 5,000 = 270 bins -> a 270 x 270 lattice.
    N_microc <- 270

    processed_microc <- process_data(
      microc_data,
      N                 = N_microc,
      scale_max         = 500,
      standardization_y = TRUE
    )

    x_vars_microc <- processed_microc[["x_vars"]]
    y_microc      <- processed_microc[["y"]]

### Step 4 — Run MCMC

    set.seed(2025)

    microc_results <- run_chain_betas(
      N                = N_microc,
      gamma_prior      = 0.3,
      iterations       = 10000,
      x_vars           = x_vars_microc,
      y                = y_microc,
      theta_start      = 0.5,
      size_start       = c(1, 1, 1),
      use_data_priors  = TRUE,
      distance_metric  = "manhattan",
      dist             = "ZINB",
      mc_cores         = min(4L, parallel::detectCores() - 1L)
    )

### Step 5 — Posterior probabilities and classification

    microc_probs <- compute_HMRFHiC_probabilities(
      data        = microc_data,
      chain_betas = microc_results,
      iterations  = 10000,
      dist        = "ZINB"
    )

    # Hard component assignment: the column with the highest posterior probability.
    probs_mat    <- as.matrix(microc_probs[, c("prob1", "prob2", "prob3")])
    component    <- max.col(probs_mat, ties.method = "first")
    table(component)

## Practical Notes on Hi-C vs. Micro-C

Aspect | Hi-C | Micro-C
--- | --- | ---
Typical resolution | 10 – 100 kb | 1 – 5 kb
Lattice size `N` | 10 – 1000 | 100 – 10000
Recommended `dist` | ZINB or NB | ZINB
Recommended `iterations` | 2000 – 5000 | 5000 – 20000
`mc_cores` | 2 – 4 | 4 – 8

The same `get_data() → process_data() → run_chain_betas() → compute_HMRFHiC_probabilities()` pipeline applies to both; only the argument values and compute budget change.

# Advanced Usage

## Custom Priors

To use user-specified priors instead of data-driven priors:

## Parallel Processing

For multiple Hi-C matrices, use mc_cores to parallelize:

## Other Distributions

The package supports Poisson, NB, ZIP, and ZINB distributions. For example, to use Poisson:

# Technical Notes

Performance: The following functions, `pz_123()`, `run_metropolis_MCMC_betas()`, and `Neighbours_combined()`, are implemented in C++ using Rcpp and RcppArmadillo for efficiency, crucial for large Hi-C matrices.

Utility Functions: The following functions, `proposaldensity_combined()`, `likelihood_gamma()`, `gamma_prior_value()`, `posterior_combined()`, `prior_combined()`, `likelihood_combined()`, `proposalfunction()`, and `size_prior()`, support the MCMC process and are typically not called directly by users.

# Conclusion

The HiCPotts package offers a powerful tool for Bayesian analysis of Hi-C data, integrating spatial dependencies, information from sources of bias associated with Hi-C data, and flexible mixture models. This page covered the core workflow, but the package’s functions can be customized for specific research needs, such as different genomic regions or distribution assumptions.

Users trying to bring in extra covariates; beyond genomic distance, GC content, ACC-score and TEs, can do so by classifying each new factor according to how it relates to those four (for example, as a distance-like term, a sequence-composition term, or an interaction‐score term) and then inputting it into the same model framework. Also, the package currently supports only intrachromosomal analyses. In the future, we plan to extend to interchromosomal contacts.

For further details, consult the package documentation (`?HiCPotts`) or contact the package maintainers. We hope HiCPotts facilitates your genomic research!

    ## R version 4.6.0 (2026-04-24)
    ## Platform: x86_64-pc-linux-gnu
    ## Running under: Ubuntu 24.04.4 LTS
    ##
    ## Matrix products: default
    ## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
    ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
    ##
    ## locale:
    ##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
    ##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
    ##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
    ##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
    ##
    ## time zone: Etc/UTC
    ## tzcode source: system (glibc)
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] HiCPotts_1.2.0   BiocStyle_2.40.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] sass_0.4.10                           generics_0.1.4
    ##  [3] SparseArray_1.12.0                    bitops_1.0-9
    ##  [5] lattice_0.22-9                        digest_0.6.39
    ##  [7] grid_4.6.0                            evaluate_1.0.5
    ##  [9] fastmap_1.2.0                         Matrix_1.7-5
    ## [11] jsonlite_2.0.0                        cigarillo_1.2.0
    ## [13] restfulr_0.0.16                       BiocManager_1.30.27
    ## [15] httr_1.4.8                            BSgenome_1.80.0
    ## [17] XML_3.99-0.23                         Biostrings_2.80.0
    ## [19] codetools_0.2-20                      jquerylib_0.1.4
    ## [21] abind_1.4-8                           cli_3.6.6
    ## [23] rlang_1.2.0                           crayon_1.5.3
    ## [25] XVector_0.52.0                        Biobase_2.72.0
    ## [27] DelayedArray_0.38.0                   cachem_1.1.0
    ## [29] yaml_2.3.12                           S4Arrays_1.12.0
    ## [31] tools_4.6.0                           parallel_4.6.0
    ## [33] BiocParallel_1.46.0                   Rhdf5lib_2.0.0
    ## [35] Rsamtools_2.28.0                      SummarizedExperiment_1.42.0
    ## [37] BiocGenerics_0.58.0                   curl_7.1.0
    ## [39] buildtools_1.0.0                      R6_2.6.1
    ## [41] matrixStats_1.5.0                     BiocIO_1.22.0
    ## [43] stats4_4.6.0                          lifecycle_1.0.5
    ## [45] rhdf5_2.56.0                          rtracklayer_1.72.0
    ## [47] Seqinfo_1.2.0                         S4Vectors_0.50.0
    ## [49] IRanges_2.46.0                        bslib_0.10.0
    ## [51] Rcpp_1.1.1-1.1                        xfun_0.57
    ## [53] GenomicRanges_1.64.0                  GenomicAlignments_1.48.0
    ## [55] MatrixGenerics_1.24.0                 sys_3.4.3
    ## [57] knitr_1.51                            rhdf5filters_1.22.0
    ## [59] rjson_0.2.23                          htmltools_0.5.9
    ## [61] rmarkdown_2.31                        maketools_1.3.2
    ## [63] compiler_4.6.0                        BSgenome.Dmelanogaster.UCSC.dm6_1.4.1
    ## [65] RCurl_1.98-1.18