---
title: "Introduction to HiContacts"
author: "Jacques Serizay"
date: "`r Sys.Date()`"
output: 
    BiocStyle::html_document
vignette: >
    %\VignetteIndexEntry{Introduction to HiContacts}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r, 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(HiCExperiment)
    library(HiContactsData)
    library(HiContacts)
})
```

# Citing HiContacts

```{r}
citation('HiContacts')
```

# Basics: importing `.(m)/cool` files as `HiCExperiment` objects

The `HiCExperiment` package provides classes and methods to import an .(m)cool 
file in R. The `HiContactsData` package gives access to a range of toy 
datasets stored by Bioconductor in the `ExperimentHub`. 

```{r}
library(dplyr)
library(ggplot2)
library(HiCExperiment)
library(HiContacts)
library(HiContactsData)
library(rtracklayer)
library(InteractionSet)
cool_file <- HiContactsData('yeast_wt', format = 'cool')
hic <- import(cool_file, format = 'cool')
hic
```

# Plotting matrices 

## Plot matrix heatmaps

The `plotMatrix` function takes a `HiCExperiment` object and plots it as a heatmap.  
Use the `use.scores` argument to specify which type of interaction scores to use 
in the contact maps (e.g. `count`, `balanced`, ...). By default, `plotMatrix()` 
looks for balanced scores. If they are not stored in the original `.(m)/cool` file, 
`plotMatrix()` simply takes the first scores available. 

```{r}
## Square matrix
plotMatrix(hic, use.scores = 'balanced', limits = c(-4, -1))

## Horizontal matrix
plotMatrix(
    refocus(hic, 'II'),
    use.scores = 'balanced', limits = c(-4, -1), 
    maxDistance = 200000
)
```

## Plot loops

Loops can be plotted on top of Hi-C matrices by providing a `GInteractions` 
object to the `loops` argument. 

*Note:* 
Loops in `.bedpe` format can be imported in R using the `import()` function, 
and converted into `GInteractions` with the 
`InteractionSet::makeGInteractionsFromGRangesPairs()` function. 

```{r}
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
loops <- system.file("extdata", 'S288C-loops.bedpe', package = 'HiCExperiment') |> 
    import() |> 
    makeGInteractionsFromGRangesPairs()
p <- import(mcool_file, format = 'mcool', focus = 'IV') |> 
    plotMatrix(loops = loops, limits = c(-4, -1), dpi = 120)
```

## Plot borders

```{r}
borders <- system.file("extdata", 'S288C-borders.bed', package = 'HiCExperiment') |> 
    import()
p <- import(mcool_file, format = 'mcool', focus = 'IV') |> 
    plotMatrix(loops = loops, borders = borders, limits = c(-4, -1), dpi = 120)
```

## Plot aggregated matrices over features

```{r}
aggr_centros <- HiContacts::aggregate(
    hic, targets = loops, BPPARAM = BiocParallel::SerialParam()
)
plotMatrix(
    aggr_centros, use.scores = 'detrended', limits = c(-1, 1), scale = 'linear', 
    cmap = bgrColors()
)
```

# Arithmetics

## Computing autocorrelated contact map 

```{r}
mcool_file <- HiContactsData('mESCs', format = 'mcool')
hic <- import(mcool_file, format = 'mcool', focus = 'chr2', resolution = 160000)
hic <- autocorrelate(hic)
scores(hic)
summary(scores(hic, 'autocorrelated'))
plotMatrix(hic, use.scores = 'autocorrelated', limits = c(-1, 1), scale = 'linear')
```

## Detrending contact map (map of scores over expected)

```{r}
hic <- import(mcool_file, format = 'mcool', focus = 'chr18:20000000-35000000', resolution = 40000)
detrended_hic <- detrend(hic)
patchwork::wrap_plots(
    plotMatrix(detrended_hic, use.scores = 'expected', scale = 'log10', limits = c(-3, -1), dpi = 120),
    plotMatrix(detrended_hic, use.scores = 'detrended', scale = 'linear', limits = c(-1, 1), dpi = 120)
)
```

## Summing two maps

```{r}
mcool_file_1 <- HiContactsData('yeast_eco1', format = 'mcool')
mcool_file_2 <- HiContactsData('yeast_wt', format = 'mcool')
hic_1 <- import(mcool_file_1, format = 'mcool', focus = 'II:1-300000', resolution = 2000)
hic_2 <- import(mcool_file_2, format = 'mcool', focus = 'II:1-300000', resolution = 2000)
merged_hic <- merge(hic_1, hic_2) 
hic_1
hic_2
merged_hic
```

## Computing ratio between two maps

```{r}
hic_1 <- import(mcool_file_1, format = 'mcool', focus = 'II', resolution = 2000)
hic_2 <- import(mcool_file_2, format = 'mcool', focus = 'II', resolution = 2000)
div_hic <- divide(hic_1, by = hic_2) 
div_hic
p <- patchwork::wrap_plots(
    plotMatrix(hic_1, use.scores = 'balanced', scale = 'log10', limits = c(-4, -1)),
    plotMatrix(hic_2, use.scores = 'balanced', scale = 'log10', limits = c(-4, -1)),
    plotMatrix(div_hic, use.scores = 'balanced.fc', scale = 'log2', limits = c(-2, 2), cmap = bwrColors())
)
```

## Despeckling (smoothing out) a contact map

```{r}
hic_1_despeckled <- despeckle(hic_1)
hic_1_despeckled5 <- despeckle(hic_1, focal.size = 5)
p <- patchwork::wrap_plots(
    plotMatrix(hic_1, use.scores = 'balanced', scale = 'log10', limits = c(-4, -1)),
    plotMatrix(hic_1_despeckled, use.scores = 'balanced.despeckled', scale = 'log10', limits = c(-4, -1)),
    plotMatrix(hic_1_despeckled5, use.scores = 'balanced.despeckled', scale = 'log10', limits = c(-4, -1))
)
```

# Mapping topological features 

## Chromosome compartments 

```{r}
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
hic <- import(mcool_file, format = 'mcool', resolution = 16000)

# - Get compartments
hic <- getCompartments(hic, chromosomes = c('XV', 'XVI'))
hic

# - Export compartments as bigwig and bed files
export(IRanges::coverage(metadata(hic)$eigens, weight = 'eigen'), 'compartments.bw')
export(
    topologicalFeatures(hic, 'compartments')[topologicalFeatures(hic, 'compartments')$compartment == 'A'], 
    'A-compartments.bed'
)
export(
    topologicalFeatures(hic, 'compartments')[topologicalFeatures(hic, 'compartments')$compartment == 'B'], 
    'B-compartments.bed'
)

# - Generate saddle plot
plotSaddle(hic)
```

## Diamond insulation score and chromatin domains borders

```{r}
# - Compute insulation score
hic <- refocus(hic, 'II:1-300000') |> 
    zoom(resolution = 1000) |> 
    getDiamondInsulation(window_size = 8000) |> 
    getBorders()
hic

# - Export insulation as bigwig track and borders as bed file
export(IRanges::coverage(metadata(hic)$insulation, weight = 'insulation'), 'insulation.bw')
export(topologicalFeatures(hic, 'borders'), 'borders.bed')
```

# Contact map analysis

## Virtual 4C

```{r}
mcool_file <- HiContactsData('mESCs', format = 'mcool')
hic <- import(mcool_file, format = 'mcool', focus = 'chr18:20000000-35000000', resolution = 40000)
v4C <- virtual4C(hic, viewpoint = GRanges('chr18:31000000-31050000'))
plot4C(v4C, ggplot2::aes(x = center, y = score))
```

## Cis-trans ratios

```{r}
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
hic <- import(mcool_file, format = 'mcool', resolution = 1000)
cisTransRatio(hic)
```

## P(s)

```{r}
# Without a pairs file
mcool_file <- HiContactsData('yeast_wt', format = 'mcool')
hic <- import(mcool_file, format = 'mcool', resolution = 1000)
ps <- distanceLaw(hic)
plotPs(ps, ggplot2::aes(x = binned_distance, y = norm_p))

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

# Comparing P(s) curves
c1 <- import(
    HiContactsData('yeast_wt', format = 'mcool'), 
    format = 'mcool',
    resolution = 1000, 
    pairsFile = HiContactsData('yeast_wt', format = 'pairs.gz')
)
c2 <- import(
    HiContactsData('yeast_eco1', format = 'mcool'), 
    format = 'mcool',
    resolution = 1000, 
    pairsFile = 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, ggplot2::aes(x = binned_distance, y = norm_p, group = sample, color = sample))
plotPsSlope(ps, ggplot2::aes(x = binned_distance, y = slope, group = sample, color = sample))
```

# Session info

```{r}
sessionInfo()
```