## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
    eval = TRUE,
    echo = TRUE,
    collapse = TRUE,
    comment = "#>"
)
library(HiCPotts)

## ----get-data-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)

## ----process-data-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
)

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

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

## ----fit-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
)

## ----relabel-example----------------------------------------------------------
oriented_chains <- relabel_hicpotts(fit$all_fits)
oriented_chains[[1L]]$relabel_rule
oriented_chains[[1L]]$noise_relationship_probability

## ----classify-example---------------------------------------------------------
classification <- classify_hicpotts(
    fit,
    data = hic,
    use = "selected",
    min_draws = 5L,
    reflect = "auto"
)

table(classification$classification)
head(classification[c("prob1", "prob2", "prob3", "probability_margin")])

## ----secondary-probability-example--------------------------------------------
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
)

## ----probability-summary-example----------------------------------------------
probability_summary <- summarise_hicpotts_probabilities(
    classification,
    ci_level = 0.95,
    include_hard_calls = TRUE,
    include_interaction_summary = TRUE
)

## ----diagnose-example---------------------------------------------------------
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")]

## ----allocation-example-------------------------------------------------------
allocation <- allocation_diagnostics(fit)
allocation$worst_cell_ess
allocation$max_mcse
allocation$between_chain_disagreement

## ----parameter-summary-example------------------------------------------------
# Display estimates and diagnostic indicators from the small worked example.
parameter_summary <- summarise_hicpotts_parameters(
    fit,
    x_vars = x_matrices,
    require_reliable = FALSE
)
head(parameter_summary[c("parameter", "estimate", "resolved")])

## ----posterior-predictive-example---------------------------------------------
ppc <- posterior_predictive_hicpotts(
    fit,
    x_vars = x_matrices,
    y = y_matrix,
    dist = "ZINB",
    n_rep = 5L,
    seed = 4921L
)
ppc$summary

## ----trace-example------------------------------------------------------------
plot_hicpotts_mcmc_by_component(
    fit$fits,
    index = 1L,
    burnin = 20L,
    plot_globals = TRUE,
    plot_size = TRUE,
    ask = FALSE
)

## ----map-example--------------------------------------------------------------
map <- plot_upper_prob_lower_count(
    results = classification,
    bin1_col = "start",
    bin2_col = "start.j.",
    prob_col = "prob2",
    count_col = "interactions",
    chr_label = "chr4",
    title = "Signal probability and observed counts",
    use_log_count = TRUE,
    symmetric_matrix = TRUE
)
print(map)

## ----block-example, eval=FALSE------------------------------------------------
# block_fit <- combine_hicpotts_blocks(
#     fits = fits_by_block,
#     data = data_by_block,
#     processed = processed_blocks,
#     what = "fit"
# )
# 
# block_classification <- classify_hicpotts(
#     block_fit,
#     min_draws = 100L,
#     reflect = "auto"
# )
# 
# block_parameters <- summarise_hicpotts_parameters(
#     block_fit,
#     pool_blocks = "none",
#     require_reliable = TRUE
# )
# 
# # Descriptive cross-block aggregation, not a joint shared-parameter fit.
# pooled_parameters <- summarise_hicpotts_parameters(
#     block_fit,
#     pool_blocks = "posterior_weighted",
#     require_reliable = FALSE,
#     pooling_draws = 100000L,
#     pooling_seed = 1L
# )

## ----complete-template, eval=FALSE--------------------------------------------
# library(HiCPotts)
# 
# input <- read.csv("contacts.csv")
# N <- length(unique(input$start))
# 
# processed <- process_data(
#     input,
#     N = N,
#     standardization_y = FALSE,
#     mirror = FALSE
# )
# y <- processed$y[[1L]]
# x_vars <- lapply(processed$x_vars, function(x) x[[1L]])
# 
# fit <- run_chain_betas(
#     N = N,
#     x_vars = x_vars,
#     y = y,
#     dist = "ZINB",
#     gamma_start = 0.3,
#     theta_start = 0.05,
#     size_start = rep(11, 3L),
#     iterations = 20000L,
#     burnin = 5000L,
#     n_chains = 4L,
#     seeds = c(101L, 202L, 303L, 404L),
#     robust = TRUE,
#     use_data_priors = TRUE,
#     mcse_stop = FALSE,
#     mc_cores = 4L
# )
# 
# classification <- classify_hicpotts(
#     fit,
#     data = input,
#     min_draws = 100L,
#     reflect = "auto"
# )
# 
# diagnostics <- diagnose_hicpotts_fit(
#     fit,
#     prob_result = classification
# )
# print(diagnostics$reliability_flags)
# 
# allocation <- allocation_diagnostics(fit)
# 
# parameters <- summarise_hicpotts_parameters(
#     fit,
#     x_vars = x_vars,
#     require_reliable = TRUE
# )
# 
# probability_summary <- summarise_hicpotts_probabilities(classification)
# 
# ppc <- posterior_predictive_hicpotts(
#     fit,
#     x_vars = x_vars,
#     y = y,
#     dist = "ZINB",
#     n_rep = 100L,
#     seed = 1L
# )
# 
# map <- plot_upper_prob_lower_count(
#     classification,
#     bin1_col = "start",
#     bin2_col = "start.j.",
#     prob_col = "prob2",
#     count_col = "interactions",
#     symmetric_matrix = TRUE
# )
# print(map)

## ----session-info, echo=FALSE-------------------------------------------------
sessionInfo()