## ---- eval = TRUE, echo=FALSE, results="hide", warning=FALSE------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>", 
    crop = NULL
)
suppressPackageStartupMessages({
    library(ggplot2)
    library(dplyr)
    library(GenomicRanges)
    library(HiContactsData)
    library(HiContacts)
})

## -----------------------------------------------------------------------------
showClass("Contacts")
contacts <- contacts_yeast()
contacts

## -----------------------------------------------------------------------------
citation('HiContacts')

## -----------------------------------------------------------------------------
library(HiContacts)
library(HiContactsData)
cool_file <- HiContactsData('yeast_wt', format = 'cool')
cool_file

## -----------------------------------------------------------------------------
range <- 'I:20000-80000' # range of interest
contacts <- Contacts(cool_file, focus = range)
summary(contacts)
contacts

## -----------------------------------------------------------------------------
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
range <- 'II:0-800000' # range of interest
lsCoolResolutions(mcool_file, verbose = TRUE)
# contacts <- Contacts(mcool_file, focus = range) # This would throw an error!
contacts <- Contacts(mcool_file, focus = range, res = 1000)
contacts

## -----------------------------------------------------------------------------
contacts <- Contacts(mcool_file, focus = 'II:0-500000 x II:100000-600000', res = 1000)
is_symmetrical(contacts)

## -----------------------------------------------------------------------------
fileName(contacts)
focus(contacts)
resolutions(contacts)
resolution(contacts)
interactions(contacts)
scores(contacts)
tail(scores(contacts, 1))
tail(scores(contacts, 'balanced'))
topologicalFeatures(contacts)
pairsFile(contacts)
metadata(contacts)

## -----------------------------------------------------------------------------
seqinfo(contacts) ## To recover the `Seqinfo` object from the `.(m)cool` file
bins(contacts) ## To bin the genome at the current resolution
regions(contacts) ## To extract unique regions of the contact matrix
anchors(contacts) ## To extract "first" and "second" anchors for each interaction

## -----------------------------------------------------------------------------
scores(contacts, 'random') <- runif(length(contacts))
scores(contacts)
tail(scores(contacts, 'random'))

## -----------------------------------------------------------------------------
topologicalFeatures(contacts, 'viewpoints') <- GRanges("II:300000-320000")
topologicalFeatures(contacts)
topologicalFeatures(contacts, 'viewpoints')

## -----------------------------------------------------------------------------
as(contacts, "GInteractions")
as(contacts, "ContactMatrix")
as(contacts, "matrix")[1:10, 1:10]

## -----------------------------------------------------------------------------
plotMatrix(contacts, use.scores = 'balanced', limits = c(-4, -1), dpi = 120)

## -----------------------------------------------------------------------------
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
loops <- system.file("extdata", 'S288C-loops.bedpe', package = 'HiContacts') |> 
    rtracklayer::import() |> 
    InteractionSet::makeGInteractionsFromGRangesPairs()
p <- Contacts(mcool_file, focus = 'IV', res = 1000) |> 
    plotMatrix(loops = loops, limits = c(-4, -1), dpi = 120)

## -----------------------------------------------------------------------------
borders <- system.file("extdata", 'S288C-borders.bed', package = 'HiContacts') |> 
    rtracklayer::import()
p <- Contacts(mcool_file, focus = 'IV', res = 1000) |> 
    plotMatrix(loops = loops, borders = borders, limits = c(-4, -1), dpi = 120)

## -----------------------------------------------------------------------------
contacts <- Contacts(mcool_file, focus = range, res = 2000)
contacts_subset <- contacts['II:20000-50000']
InteractionSet::boundingBox(interactions(contacts))
InteractionSet::boundingBox(interactions(contacts_subset))

## -----------------------------------------------------------------------------
mcool_file <- HiContactsData('mESCs', format = 'mcool')
contacts <- Contacts(mcool_file, focus = 'chr2', res = 160000)
autocorrelated <- autocorrelate(contacts, ignore_ndiags = 20)
scores(autocorrelated)
summary(scores(autocorrelated, 1))

## -----------------------------------------------------------------------------
mcool_file <- HiContactsData('mESCs', format = 'mcool')
contacts <- Contacts(mcool_file, focus = 'chr18:20000000-35000000', res = 40000)
detrended_contacts <- detrend(contacts)
cowplot::plot_grid(
    plotMatrix(detrended_contacts, use.scores = 'expected', dpi = 120),
    plotMatrix(detrended_contacts, use.scores = 'detrended', scale = 'linear', limits = c(-3, 3), dpi = 120)
)

## -----------------------------------------------------------------------------
mcool_file_1 <- HiContactsData('yeast_eco1', format = 'mcool')
mcool_file_2 <- HiContactsData('yeast_wt', format = 'mcool')
contacts_1 <- Contacts(mcool_file_1, focus = 'II', res = 2000)
contacts_1 <- contacts_1['II:1-300000']
contacts_2 <- Contacts(mcool_file_2, focus = 'II', res = 2000)
contacts_2 <- contacts_2['II:1-300000']
merged_contacts <- merge(contacts_1, contacts_2) 
contacts_1
contacts_2
merged_contacts

## -----------------------------------------------------------------------------
mcool_file_1 <- HiContactsData('yeast_eco1', format = 'mcool')
mcool_file_2 <- HiContactsData('yeast_wt', format = 'mcool')
contacts_1 <- Contacts(mcool_file_1, focus = 'II', res = 2000)
contacts_2 <- Contacts(mcool_file_2, focus = 'II', res = 2000)
div_contacts <- divide(contacts_1, by = contacts_2) 
div_contacts
p <- cowplot::plot_grid(
    plotMatrix(contacts_1, use.scores = 'balanced', scale = 'log10', limits = c(-4, -1)),
    plotMatrix(contacts_2, use.scores = 'balanced', scale = 'log10', limits = c(-4, -1)),
    plotMatrix(div_contacts, use.scores = 'ratio', scale = 'log2', limits = c(-2, 2), cmap = bwrColors())
)

## -----------------------------------------------------------------------------
mcool_file <- HiContactsData('mESCs', format = 'mcool')
contacts <- Contacts(mcool_file, focus = 'chr18:20000000-35000000', res = 40000)
v4C <- virtual4C(contacts, viewpoint = GRanges('chr18:31000000-31050000'))
plot4C(v4C, aes(x = center, y = score))

## -----------------------------------------------------------------------------
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
contacts <- Contacts(mcool_file, res = 1000)
cisTransRatio(contacts)

## -----------------------------------------------------------------------------
# Without a pairs file
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
contacts <- Contacts(mcool_file, res = 1000)
ps <- distanceLaw(contacts)
plotPs(ps, aes(x = binned_distance, y = norm_p))

# With a pairs file
pairsFile(contacts) <- HiContactsData('yeast_wt', format = 'pairs.gz')
ps <- distanceLaw(contacts)
plotPs(ps, aes(x = binned_distance, y = norm_p))
plotPsSlope(ps, aes(x = binned_distance, y = slope))

# Comparing P(s) curves
c1 <- Contacts(
    HiContactsData('yeast_wt', format = 'mcool'), 
    res = 1000, 
    pairs = HiContactsData('yeast_wt', format = 'pairs.gz')
)
c2 <- Contacts(
    HiContactsData('yeast_eco1', format = 'mcool'), 
    res = 1000, 
    pairs = HiContactsData('yeast_eco1', format = 'pairs.gz')
)
ps_1 <- distanceLaw(c1) |> mutate(sample = 'WT')
ps_2 <- distanceLaw(c2) |> mutate(sample = 'Eco1-AID')
ps <- rbind(ps_1, ps_2)
plotPs(ps, aes(x = binned_distance, y = norm_p, group = sample, color = sample))
plotPsSlope(ps, aes(x = binned_distance, y = slope, group = sample, color = sample))

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