In js2264/cooleR: Analysing cool files in R with HiContacts

knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>", 
    crop = NULL
)
suppressPackageStartupMessages({
    library(ggplot2)
    library(dplyr)
    library(GenomicRanges)
    library(HiCExperiment)
    library(HiContactsData)
    library(HiContacts)
})

Citing HiContacts

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.

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.

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

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

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

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

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)

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

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

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

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

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

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

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

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

P(s)

# 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

sessionInfo()

js2264/cooleR documentation built on Jan. 29, 2024, 10:47 p.m.