---
title: "Bayesian Analysis of Hi-C Interactions with HiCPotts"
author: "Itunu Godwin Osuntoki, Nicolae Radu Zabet"
date: "`r Sys.Date()`"
output:   
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Bayesian Analysis of Hi-C Interactions with HiCPotts}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# 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:

a.  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.


b.  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.

c.  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.

d.  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


```{r setup, include=FALSE}
library(BiocStyle)
knitr::opts_chunk$set(eval = TRUE, echo = TRUE, collapse = TRUE)
```

```{r install-cran, echo=TRUE, eval=FALSE}
install.packages(c("Rcpp", "RcppArmadillo", "parallel"))
```

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

```{r load-package, echo=TRUE, eval=TRUE}
# 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 \times 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

```{r simulate-data, include=FALSE, eval=TRUE}
set.seed(123)

N <- 10

bins <- seq(1, 20000, by = 2000)

data <- expand.grid(start = bins, start.j. = bins)

data$end.i. <- data$start + 1999

data$end <- data$start.j. + 1999

data$chrom <- "chr2L"

data$GC <- runif(nrow(data), 0.3, 0.7) # Simulated GC content

data$ACC <- runif(nrow(data), 0, 1) # Simulated accessibility

data$TES <- rpois(nrow(data), 2) # Simulated TE counts

data$interactions <- rpois(nrow(data), lambda = 5) # Simulated interaction counts
data <- as.data.frame(data)

head(data)
```

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.

```{r get_data, include=FALSE, eval=FALSE}
# Example (not run)

data <- prepare_data(
  file_path = "path/to/hic_file.h5",
  chr = "chr2L",
  start = 1,
  end = 20000,
  resolution = 2000,
  genome_package = "BSgenome.Dmelanogaster.UCSC.dm6",
  acc_wig = "path/to/dnase.wig",
  chain_file = "path/to/chain.file",
  te_granges = "path/to/te.gtf"
)
```

## 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.

```{r process_data, include=FALSE, eval=TRUE}
processed <- process_data(data, N = N, scale_max = 500, standardization_y = TRUE)

x_vars <- processed[["x_vars"]] # Single matrix for this example
y <- processed[["y"]]

str(x_vars) # Lists of distance, GC, TEs, ACC matrices
```

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.

```{r setup_MCMC, include=FALSE, eval=TRUE}
# Set MCMC parameters
gamma_prior <- 0.3

iterations <- 20

theta_start <- c(0.5)

# size_start <- c(1, 1, 1)

dist <- "ZIP"
```

```{r run_MCMC, include=FALSE, eval=TRUE}
# Run MCMC
results <- run_chain_betas(
  N = N,
  gamma_prior = gamma_prior,
  iterations = iterations,
  x_vars = x_vars,
  y = y,
  theta_start = theta_start,
  use_data_priors = TRUE,
  # user_fixed_priors = FALSE,
  # epsilon = NULL,
  distance_metric = "manhattan",
  dist = dist,
  # size_start = size_start,
  mc_cores = 1
)
```

```{r MCMC_results, include=FALSE, eval=TRUE}
# Extract results
chains <- results[["chains"]]

gamma <- results[["gamma"]]

theta <- results[["theta"]]

# size <- results[["size]]
```

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.

```{r compute_result, include=FALSE, eval=TRUE}
# Compute probabilities
probs <- compute_HMRFHiC_probabilities(
  data = data,
  chain_betas = chains,
  iterations = iterations,
  dist = "ZIP"
)

probs$prob1 # This produces matrices of probabilities for the first component, which translate to assigning interactions to biological states (e.g., active vs. inactive chromatin).
```

# Advanced Usage

## Custom Priors

To use user-specified priors instead of data-driven priors:
```{r User_prior, include=FALSE, eval=TRUE}
user_priors <- list(
  list(mean = c(0, 0, 0, 0, 0), sd = c(1, 1, 1, 1, 1)),
  list(mean = c(0, 0, 0, 0, 0), sd = c(1, 1, 1, 1, 1)),
  list(mean = c(0, 0, 0, 0, 0), sd = c(1, 1, 1, 1, 1))
)

# Run MCMC
results <- run_chain_betas(
  N = N,
  gamma_prior = gamma_prior,
  iterations = iterations,
  x_vars = x_vars,
  y = y,
  theta_start = theta_start,
  # use_data_priors = TRUE,
  user_fixed_priors = user_priors,
  # epsilon = NULL,
  distance_metric = "manhattan",
  dist = dist,
  # size_start = size_start,
  mc_cores = 1
)
```

## Parallel Processing

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

```{r multiple_matrices, include=FALSE, eval=FALSE}
y_list <- list(y, y) # Example with two matrices

results <- run_chain_betas(
  N = N,
  gamma_prior = gamma_prior,
  iterations = iterations,
  x_vars = x_vars,
  y = y_list,
  theta_start = theta_start,
  size_start = size_start,
  use_data_priors = TRUE,
  user_fixed_priors = NULL,
  epsilon = NULL,
  distance_metric = "manhattan",
  dist = dist,
  mc_cores = 2
)
```

## Other Distributions

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

```{r Additional_distributions, include=FALSE, eval=TRUE}
# Run MCMC
results <- run_chain_betas(
  N = N,
  gamma_prior = gamma_prior,
  iterations = iterations,
  x_vars = x_vars,
  y = y,
  theta_start = theta_start,
  # use_data_priors = TRUE,
  user_fixed_priors = user_priors,
  # epsilon = NULL,
  distance_metric = "manhattan",
  dist = "Poisson",
  # size_start = size_start,
  mc_cores = 1
)
```

# 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, echo = FALSE}
sessionInfo()
```
