## ---- message = FALSE------------------------------------------------------
library("HiCBricks")


## --------------------------------------------------------------------------

Bintable.path <- system.file("extdata",
"Bintable_40kb.txt", package = "HiCBricks")
Chromosomes <- "chr19"


## --------------------------------------------------------------------------

Path_to_cached_file <- CreateBrick(ChromNames = Chromosomes,
    BinTable = Bintable.path, bin.delim = " ",
    Output.Filename = file.path(tempdir(),"test.hdf"), exec = "cat",
remove.existing = TRUE)


## --------------------------------------------------------------------------

Test.mat <- matrix(runif(800*800),nrow = 800, ncol = 800)


## --------------------------------------------------------------------------

Matrix.file <- file.path(tempdir(),"Test_matrix.txt")
write.table(x = Test.mat, file = Matrix.file, sep = " ", quote = FALSE,
row.names = FALSE, col.names = FALSE)


## --------------------------------------------------------------------------

Brick.file <- Path_to_cached_file

Brick_load_matrix(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19",
    matrix.file = Matrix.file,
    delim = " ",
    exec = "cat",
    remove.prior = TRUE)


## --------------------------------------------------------------------------

Brick_load_cis_matrix_till_distance(Brick = Brick.file,
    chr = "chr19",
    matrix.file = Matrix.file,
    delim = " ",
    distance = 100,
    remove.prior = TRUE)


## ----eval = FALSE----------------------------------------------------------
#  require(curl)
#  Consortium.home = "https://data.4dnucleome.org/files-processed"
#  File = file.path(Consortium.home,"4DNFI7JNCNFB",
#      "@@download","4DNFI7JNCNFB.mcool")
#  curl_download(url = File,
#      destfile = file.path(temp.dir(),"H1-hESC-HiC-4DNFI7JNCNFB.mcool"))

## ----eval = FALSE----------------------------------------------------------
#  
#  mcoolName=file.path(temp.dir(),"H1-hESC-HiC-4DNFI7JNCNFB.mcool")
#  Brick_list_mcool_resolutions(mcool = mcoolName)
#  

## ----eval = FALSE----------------------------------------------------------
#  
#  Brick_mcool_normalisation_exists(mcool = mcoolName,
#      norm.factor = "Iterative-Correction",
#      binsize = 10000)
#  

## ----eval = FALSE----------------------------------------------------------
#  
#  Brick.name <- file.path(tempdir(),
#    "H1-hESC-HiC-4DNFI7JNCNFB-10000-ICE-normalised-chr1.brick")
#  
#  Output.brick <- CreateBrick_from_mcool(Brick = Brick.name,
#      mcool = mcoolName, binsize = 10000, chrs = "chr1", remove.existing = TRUE)
#  

## ----eval = FALSE----------------------------------------------------------
#  
#  Brick_load_data_from_mcool(Brick = Output.brick,
#      mcool = mcoolName,
#      chr1 = "chr1",
#      chr2 = "chr1",
#      binsize = 10000,
#      cooler.batch.size = 1000000,
#      matrix.chunk = 2000,
#      dont.look.for.chr2 = TRUE,
#      remove.prior = TRUE,
#      norm.factor = "Iterative-Correction")
#  

## --------------------------------------------------------------------------

Brick.file <- system.file("extdata", "test.hdf", package = "HiCBricks")

Brick_list_rangekeys(Brick.file)


## --------------------------------------------------------------------------

Brick_get_bintable(Brick.file)


## --------------------------------------------------------------------------

Brick_get_ranges(Brick = Brick.file,
    rangekey = "Bintable")


## --------------------------------------------------------------------------

Brick_get_ranges(Brick = Brick.file,
    rangekey = "Bintable",
    chr = "chr19")


## --------------------------------------------------------------------------

Brick_return_region_position(Brick = Brick.file,
    region = "chr19:5000000:10000000")



## --------------------------------------------------------------------------

Brick_fetch_range_index(Brick = Brick.file,
    chr = "chr19",
    start = 5000000,
    end = 10000000)



## --------------------------------------------------------------------------

Values <- Brick_get_values_by_distance(Brick = Brick.file,
    chr = "chr19",
    distance = 4)


## --------------------------------------------------------------------------

Failsafe_median_log10 <- function(x){
    x[is.nan(x) | is.infinite(x) | is.na(x)] <- 0
    return(median(log10(x+1)))
}

Brick_get_values_by_distance(Brick = Brick.file,
    chr = "chr19",
    distance = 4,
    FUN = Failsafe_median_log10)


## --------------------------------------------------------------------------

Failsafe_median_log10 <- function(x){
    x[is.nan(x) | is.infinite(x) | is.na(x)] <- 0
    return(median(log10(x+1)))
}

Brick_get_values_by_distance(Brick = Brick.file,
    chr = "chr19",
    distance = 4,
    constrain.region = "chr19:1:5000000",
    FUN = Failsafe_median_log10)


## --------------------------------------------------------------------------

Sub.matrix <- Brick_get_matrix_within_coords(Brick = Brick.file,
    x.coords="chr19:5000000:10000000",
    force = TRUE,
    y.coords = "chr19:5000000:10000000")

dim(Sub.matrix)


## --------------------------------------------------------------------------

x.axis <- 5000000/40000
y.axis <- 10000000/40000

Sub.matrix <- Brick_get_matrix(Brick = Brick.file,
    chr1 = "chr19", chr2 = "chr19",
    x.vector = c(x.axis:y.axis),
    y.vector = c(x.axis:y.axis))

dim(Sub.matrix)


## --------------------------------------------------------------------------

Coordinate <- c("chr19:1:40000","chr19:40001:80000")
Brick.file <- system.file("extdata",
    "test.hdf",
    package = "HiCBricks")
Test_Run <- Brick_fetch_row_vector(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19",
    by = "ranges",
    vector = Coordinate)


## --------------------------------------------------------------------------

Coordinate <- c(1,2)
Brick.file <- system.file("extdata",
    "test.hdf",
    package = "HiCBricks")
Test_Run <- Brick_fetch_row_vector(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19",
    by = "position",
    vector = Coordinate)


## --------------------------------------------------------------------------

Coordinate <- c(1,2)
Brick.file <- system.file("extdata",
    "test.hdf",
    package = "HiCBricks")
Test_Run <- Brick_fetch_row_vector(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19",
    by = "position",
    vector = Coordinate,
    regions = c("chr19:1:1000000", "chr19:40001:2000000"))


## --------------------------------------------------------------------------

Brick_list_matrix_mcols()


## --------------------------------------------------------------------------

Brick.file <- system.file("extdata",
    "test.hdf",
    package = "HiCBricks")
MCols.dat <- Brick_get_matrix_mcols(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19",
    what = "row.sums")
head(MCols.dat, 100)


## --------------------------------------------------------------------------

Brick_matrix_isdone(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------

Brick_matrix_issparse(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------

Brick_matrix_maxdist(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------

Brick_matrix_exists(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------

Brick_matrix_minmax(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------

Brick_matrix_dimensions(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------

Brick_matrix_filename(Brick = Brick.file,
    chr1 = "chr19",
    chr2 = "chr19")


## --------------------------------------------------------------------------
Brick.file <- system.file("extdata",
    "test.hdf", package = "HiCBricks")

Chromosome <- "chr19"
di_window <- 10
lookup_window <- 30
TAD_ranges <- Brick_local_score_differentiator(Brick = Brick.file,
    chrs = Chromosome,
    di.window = di_window,
    lookup.window = lookup_window,
    strict = TRUE,
    fill.gaps=TRUE,
    chunk.size = 500)

## --------------------------------------------------------------------------

Name <- paste("LSD",
    di_window,
    lookup_window,
    Chromosome,sep = "_")
Brick_add_ranges(Brick = Path_to_cached_file,
    ranges = TAD_ranges,
    rangekey = Name)


## --------------------------------------------------------------------------

Brick_list_rangekeys(Brick = Path_to_cached_file)

TAD_ranges <- Brick_get_ranges(Brick = Path_to_cached_file,
    rangekey = Name)


## ----plot, fig.cap = "A normal heatmap without any transformations", fig.small = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(), 
  "chr19-5MB-10MB-normal.pdf"),
    Bricks = Brick.file,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    palette = "Reds",
    width = 10,
    height = 11,
    return.object=TRUE)


## ----plot2, fig.cap = "A normal heatmap with colours computed in log10 scale", fig.small = TRUE----

Failsafe_log10 <- function(x){
    x[is.na(x) | is.nan(x) | is.infinite(x)] <- 0
    return(log10(x+1))
}

Brick_vizart_plot_heatmap(File = file.path(tempdir(), 
  "chr19-5MB-10MB-normal-colours-log10.pdf"),
    Bricks = Brick.file,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    legend.title = "Log10 Hi-C signal",
    palette = "Reds",
    width = 10,
    height = 11,
    return.object=TRUE)


## ----plot3, fig.cap = "A normal heatmap with colours computed in log10 scale after capping values to the 99th percentile", fig.small = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-normal-colours-log10-valuecap-99.pdf"),
    Bricks = Brick.file,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    legend.title = "Log10 Hi-C signal",
    palette = "Reds",
    width = 10,
    height = 11,
    return.object=TRUE)


## ----plot4, fig.cap = "Same heatmap as before with colours computed in log10 scale after capping values to the 99th percentile with 45 degree rotation", fig.small = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-normal-colours-log10-rotate.pdf"),
    Bricks = Brick.file,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    distance = 60,
    legend.title = "Log10 Hi-C signal",
    palette = "Reds",
    width = 10,
    height = 11,
    rotate = TRUE,
    return.object=TRUE)


## ----plot5, fig.cap = "Same heatmap as previous, but now the heatmaps are wider than they are taller", fig.wide = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-normal-colours-log10-rotate-2.pdf"),
    Bricks = Brick.file,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    distance = 60,
    legend.title = "Log10 Hi-C signal",
    palette = "Reds",
    width = 15,
    height = 5,
    rotate = TRUE,
    return.object=TRUE)


## ----plot6, fig.cap = "Normal rectangular heatmap with colours computed in the log scale after capping values to the 99th percentile with TAD calls", fig.small = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-normal-colours-log10-rotate-2-tads.pdf"),
    Bricks = Brick.file,
    tad.ranges = TAD_ranges,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    colours = "#230C0F",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    legend.title = "Log10 Hi-C signal",
    palette = "Reds",
    width = 10,
    height = 11,
    return.object=TRUE)

## ----plot7, fig.cap = "Normal rotated heatmap with colours computed in the log scale after capping values to the 99th percentile with TAD calls", fig.wide = TRUE----
Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-normal-colours-log10-rotate-3-tads.pdf"),
    Bricks = Brick.file,
    tad.ranges = TAD_ranges,
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    colours = "#230C0F",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    distance = 60,
    legend.title = "Log10 Hi-C signal",
    palette = "Reds",
    width = 15,
    height = 5,
    line.width = 0.8,
    cut.corners = TRUE,
    rotate = TRUE,
    return.object=TRUE)

## ----plot8, fig.cap = "A normal two sample heatmap with colours computed in log10 scale after capping values to the 99th percentile", fig.small = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-colours-log10-valuecap-99.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    legend.title = "Log10 Hi-C signal",
    palette = "YlGnBu",
    width = 10,
    height = 11,
    return.object=TRUE)


## ----plot9, fig.cap = "A normal two sample heatmap with colours computed in log10 scale on values until the 30th diagonal and capping these values to the 99th percentile.", fig.small = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-colours-log10-valuecap-99-2.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    legend.title = "Log10 Hi-C signal",
    palette = "RdGy",
    distance = 30,
    width = 10,
    height = 11,
    return.object=TRUE)


## ----plot10, fig.cap = "A rotated two sample heatmap with colours computed in log10 scale and capping these values to the 99th percentile.", fig.wide = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-colours-log10-valuecap-99-rotate.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    legend.title = "Log10 Hi-C signal",
    palette = "YlGnBu",
    distance = 30,
    width = 15,
    height = 4,
    rotate = TRUE,
    return.object=TRUE)


## --------------------------------------------------------------------------

Brick.file <- system.file("extdata",
    "test.hdf", package = "HiCBricks")

Chromosome <- "chr19"
di_windows <- c(5,10)
lookup_windows <- c(10, 20)
for (i in seq_along(di_windows)) {

    di_window <- di_windows[i]
    lookup_window <- lookup_windows[i]
    
    TAD_ranges <- Brick_local_score_differentiator(Brick = Brick.file,
        chrs = Chromosome,
        di.window = di_window,
        lookup.window = lookup_window,
        strict = TRUE,
        fill.gaps=TRUE,
        chunk.size = 500)
    
    Name <- paste("LSD",
        di_window,
        lookup_window,
        Chromosome,sep = "_")
    
    Brick_add_ranges(Brick = Path_to_cached_file, ranges = TAD_ranges,
        rangekey = Name)
}


## --------------------------------------------------------------------------

Chromosome <- "chr19"
di_windows <- c(5,10)
lookup_windows <- c(10, 20)
TADs_list <- list()

for (i in seq_along(di_windows)) {

    di_window <- di_windows[i]
    lookup_window <- lookup_windows[i]
    
    Name <- paste("LSD",
        di_window,
        lookup_window,
        Chromosome,sep = "_")
    
    TAD_ranges <- Brick_get_ranges(Brick = Path_to_cached_file, 
        rangekey = Name)
    # Map TADs to their Hi-C maps
    TAD_ranges$group <- i
    # Map TADs to a specific categorical value for the colours
    TAD_ranges$colour.group <- paste("LSD", di_window, lookup_window, 
        sep = "_")
    TADs_list[[Name]] <- TAD_ranges
}

TADs_ranges <- do.call(c, unlist(TADs_list, use.names = FALSE))


## ----plot11, fig.cap = "A normal two sample heatmap with colours computed in log10 scale and capping these values to the 97th percentile, with TAD borders", fig.small = TRUE----

Brick.file <- system.file("extdata",
    "test.hdf", package = "HiCBricks")

Colours <- c("#B4436C", "#F78154")
Colour.names <- unique(TADs_ranges$colour.group)

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-colours-log10-valuecap-99-tads.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000001:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.97,
    legend.title = "Log10 Hi-C signal",
    palette = "YlGnBu",
    tad.ranges = TADs_ranges,
    group.col = "group",
    tad.colour.col = "colour.group",
    colours = Colours,
    colours.names = Colour.names,
    distance = 30,
    width = 9,
    height = 11,
    return.object=TRUE)


## ----plot12, fig.cap = "A rotated two sample heatmap with colours computed in log10 scale and capping these values to the 97th percentile, with non-truncated TAD borders", fig.wide = TRUE----

Brick.file <- system.file("extdata",
    "test.hdf", package = "HiCBricks")

Colours <- c("#B4436C", "#F78154")
Colour.names <- unique(TADs_ranges$colour.group)

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-colours-log10-valuecap-99-rotate-tads.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.97,
    legend.title = "Log10 Hi-C signal",
    palette = "YlGnBu",
    tad.ranges = TADs_ranges,
    group.col = "group",
    tad.colour.col = "colour.group",
    colours = Colours,
    colours.names = Colour.names,
    distance = 30,
    width = 15,
    height = 4,
    rotate = TRUE,
    return.object=TRUE)


## ----plot13, fig.cap = "A rotated two sample heatmap with colours computed in log10 scale and capping these values to the 97th percentile, with truncated TAD borders", fig.wide = TRUE----

Brick.file <- system.file("extdata",
    "test.hdf", package = "HiCBricks")

Colours <- c("#B4436C", "#F78154")
Colour.names <- unique(TADs_ranges$colour.group)

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-colours-log10-valuecap-99-rotate-tads-2.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.97,
    legend.title = "Log10 Hi-C signal",
    palette = "YlGnBu",
    tad.ranges = TADs_ranges,
    group.col = "group",
    tad.colour.col = "colour.group",
    colours = Colours,
    colours.names = Colour.names,
    distance = 30,
    width = 15,
    height = 4,
    cut.corners = TRUE,
    rotate = TRUE,
    return.object=TRUE)


## ----plot14, fig.cap = "A rotated two sample heatmap with colours computed in log10 scale and capping these values to the 99th percentile.", fig.wide = TRUE----

Brick_vizart_plot_heatmap(File = file.path(tempdir(),
  "chr19-5MB-10MB-bipartite-final.pdf"),
    Bricks = c(Brick.file, Brick.file),
    x.coords = "chr19:5000000:10000000",
    y.coords = "chr19:5000001:10000000",
    FUN = Failsafe_log10,
    value.cap = 0.99,
    legend.title = "Log10 Hi-C signal",
    palette = "YlGnBu",
    tad.ranges = TADs_ranges,
    group.col = "group",
    tad.colour.col = "colour.group",
    colours = Colours,
    colours.names = Colour.names,
    distance = 30,
    width = 15,
    height = 4,
    legend.key.width = unit(10, "mm"), 
    legend.key.height = unit(5, "mm"),
    line.width = 1.2,
    cut.corners = TRUE,
    rotate = TRUE,
    return.object=TRUE)