## ----setup-opts, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
# # Install HiSpaR from Bioconductor
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("HiSpaR")
# 
# # HiSpaR depends on several packages for Hi-C data handling
# # These are automatically installed with HiSpaR:
# # - HiCExperiment: For representing Hi-C data
# # - Matrix: For sparse matrix handling
# 
# # Additionally, install example data package and HiContacts for data manipulation
# BiocManager::install("HiContactsData")
# BiocManager::install("HiContacts")
# # Optional: Install visualization and analysis packages
# install.packages("plotly")  # For interactive 3D visualization

## ----setup--------------------------------------------------------------------
library(HiSpaR)
library(HiContacts)
library(HiContactsData)
library(plotly)

# Load example data from HiCExperiment
cool_file <- CoolFile(HiContactsData::HiContactsData('yeast_wt', format = 'cool'))
hic <- import(cool_file, focus = "II:10000-100000")
# visualize the contact matrix
plotMatrix(hic)
# Load example contact matrix (from HiSpaR's built-in data)
data(su1_contact_mat)
cat("Loaded contact matrix dimensions:", dim(su1_contact_mat), "\n")

## ----example, eval=TRUE-------------------------------------------------------
# Run analysis on the example HiCExperiment object
result_yeast <- hispa_analyze(
  hic_experiment = hic,
  output_dir = tempdir(),
  mcmc_iterations = 1000,
  use_cluster_init = TRUE,
  mcmc_burn_in = 10,
  filter_quantile = 0.1,  # Filter out loci in the bottom 10% of contact counts
  verbose = TRUE
)
# Check the structure of results
cat("Result structure:\n")
str(result_yeast)
cat("Final positions (first 5 loci):\n")
print(head(result_yeast$positions, 5))
cat("Estimated parameters: beta0 =", result_yeast$beta0, ", beta1 =", result_yeast$beta1, "\n")

# Run analysis on the example contact matrix
result_su <- hispa_analyze(
  hic_experiment = su1_contact_mat,
  output_dir = tempdir(),
  mcmc_iterations = 2000,
  use_cluster_init = TRUE,
  mcmc_burn_in = 10,
  filter_quantile = 0.1,  # Filter out loci in the bottom 10% of contact counts
  verbose = TRUE
)

# Check the structure of results
cat("Result structure:\n")
str(result_su)
cat("Final positions (first 5 loci):\n")
print(head(result_su$positions, 5))
cat("Estimated parameters: beta0 =", result_su$beta0, ", beta1 =", result_su$beta1, "\n")

## ----visualization, eval=TRUE-------------------------------------------------
coords_yeast <- result_yeast$positions
# Create interactive 3D scatter plot
plot_ly(
  x = coords_yeast[,1],
  y = coords_yeast[,2],
  z = coords_yeast[,3],
  type = 'scatter3d',
  mode = 'markers+lines',
  marker = list(size = 5, color = ~seq_len(nrow(coords_yeast)), showscale = TRUE),
) %>%
  layout(
    title = "Inferred 3D Chromatin Structure",
    scene = list(
      xaxis = list(title = "X"),
      yaxis = list(title = "Y"),
      zaxis = list(title = "Z")
    )
  )
# Extract coordinates from results
coords_su <- result_su$positions

# Create interactive 3D scatter plot
plot_ly(
  x = coords_su[,1],
  y = coords_su[,2],
  z = coords_su[,3],
  type = 'scatter3d',
  mode = 'markers+lines',
  marker = list(size = 5, color = ~seq_len(nrow(coords_su)), showscale = TRUE),
) %>%
  layout(
    title = "Inferred 3D Chromatin Structure",
    scene = list(
      xaxis = list(title = "X"),
      yaxis = list(title = "Y"),
      zaxis = list(title = "Z")
    )
  )

## ----sessionInfo--------------------------------------------------------------
sessionInfo()