---
title: "Visualising Interactions"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{Visualising Interactions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dpi = 150,
  fig.width = 5,
  fig.height = 5,
  cache = TRUE,
  eval = TRUE
    );
```

CHiCANE supports built-in plotting capabilities using [Gviz](http://bioconductor.org/packages/Gviz/). For the built-in functionality, please refer to the documentation of `chiance::create.locus.plot()`. Below, we demonstrate (direct) use of Gviz API to visualise CHiCANE's interaction peaks.

```{r, cache = TRUE}
library(chicane);
data(bre80);

options(ucscChromosomeNames = FALSE);

bre80.interactions <- chicane(interactions = bre80);

```

# Single restriction fragment

To visualise the relationship between the number of reads and the significant interactions, we can look at a single restriction fragment. The HindIII fragment chr2:217035649-217042355 contains the breast cancer risk SNP rs13387042. By selecting fragment interactions that involve this fragment, we can plot the number of reads linking the fragment to all other fragments. 

```{r, message = FALSE}

# check dependency packages which must be installed
if ("GenomicInteractions" %in% installed.packages()[, "Package"] 
	&& "Gviz" %in% installed.packages()[, "Package"]) {

library(GenomicInteractions);
library(Gviz);

options(ucscChromosomeNames = FALSE);

locus <- list(chr = 'chr2', start = 217035649, end = 217042355);
locus.id <- paste0(locus$chr, ':', locus$start, '-', locus$end);

ideogram.track <- tryCatch({
	IdeogramTrack(genome = 'hg38', chromosome = locus$chr);
	}, error = function(e) {
	cat(e$message, '\n')
	});

genome.axis.track <- GenomeAxisTrack();

# get counts and interactions involving specific locus
# TO DO: need to switch this to generic other end in case of b2b
# add note about p-value/ q-value
locus.counts <- bre80[ target.id == locus.id | bait.id == locus.id ];
locus.interactions <- bre80.interactions[ p.value < 0.001 & (target.id == locus.id | bait.id == locus.id) ];


# create track for visualising read count
interaction.count.track <- DataTrack(
    range = GRanges(locus.counts$target.id),
    data = locus.counts$count,
    name = 'Reads'
    );

# create track for visualising significant interactions
interaction.object <- GenomicInteractions(
	anchor1 = GRanges( locus.interactions$bait.id ),
	anchor2 = GRanges( locus.interactions$target.id ),
	count = -log10(locus.interactions$p.value)
    );

interaction.track <- InteractionTrack(
    interaction.object,
    name = 'Chicane'
    );

gene.track <- GeneRegionTrack(
    system.file( 'extdata', 'gencode_2q35.gtf', package = 'chicane' ),
    chr = locus$chr,
    start = locus$start - 5e5,
    end = locus$end + 7.5e5,
    stacking = 'squish',
    stackHeight = 0.3,
    name = 'Genes'
    );

# plot
# handle Gviz - UCSC incompatibility with R-4.1.0
if (!is.null(ideogram.track)) {
    plotTracks(
        list(
            ideogram.track, 
            genome.axis.track, 
            gene.track,
            interaction.count.track,
            interaction.track
            ),
        sizes = c(0.4, 1, 1, 4, 3),
        type = 'histogram',
        transcriptAnnotation = 'symbol',
        collapseTranscripts = 'longest',
        col = NULL,
        from = locus$start - 6e5,
        to = locus$end + 8e5
        );
    }
}
```

Gene annotations rely on `GeneRegionTrack` objects. These can either be downloaded through BioMart, or loaded from GTF files. 


# Longer bait regions

It is also possible to visualise significant interactions for several baits simultaneously.

```{r, fig.height = 3}

# check dependency packages which must be installed
if ("GenomicInteractions" %in% installed.packages()[, "Package"] 
	&& "Gviz" %in% installed.packages()[, "Package"]) {

baits <- read.bed( system.file( 'extdata', '2q35.bed', package = 'chicane' ) );


sig.interactions <- bre80.interactions[ (bait.id %in% baits | target.id %in% baits) & p.value < 0.001 ];

# show baits
bait.track <- AnnotationTrack( GRanges(baits), name = 'Baits', fill = 'black');


read.tracks <- lapply(
    baits,
    function(bait, bre80, baits) {
        bait.counts <- bre80[ bait.id == bait & !(target.id %in% baits) ];

        bait.chr <- bait.counts$bait.chr[1];

        read.track <- DataTrack(
            range = GRanges(bait.counts$target.id),
            data = bait.counts$count,
            name = ' ', # empty string causes plotting to fail...
            chromosome = bait.chr,
            fill.histogram = 'darkgrey',
            col.histogram = 'darkgrey',
            cex.axis = 0
            );
        
        return(read.track);
    },
    bre80 = bre80,
    baits = baits
    );

# create track for visualising significant interactions
interaction.object <- GenomicInteractions(
	anchor1 = GRanges( sig.interactions$bait.id ),
	anchor2 = GRanges( sig.interactions$target.id ),
	count = -log10(sig.interactions$p.value)
    );

interaction.track <- InteractionTrack(
    interaction.object,
    name = 'Chicane'
    );

gene.track <- GeneRegionTrack(
    system.file( 'extdata', 'gencode_2q35.gtf', package = 'chicane' ),
    chr = locus$chr,
    start = 216150000,
    end = 217800000,
    stacking = 'squish',
    stackHeight = 0.3,
    name = 'Genes'
    );

# handle Gviz - UCSC incompatibility with R-4.1.0
if (!is.null(ideogram.track)) {
    plotTracks(
        c(
            list(
                ideogram.track, 
                genome.axis.track
            ),
            read.tracks,
            list(
                interaction.track,
                bait.track,
                gene.track
                )
            ),
        sizes = c(0.4, 1, rep(0.2, length(read.tracks)), 2, 0.5, 1),
        transcriptAnnotation = 'symbol',
        collapseTranscripts = 'longest',
        type = 'hist',
        col = NULL,
        from = 216150000,
        to = 217800000
        );
    }
}
```