CHiCANE supports built-in plotting capabilities using 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.

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.

# 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
        );
    }
}
#> Warning: package 'GenomicRanges' was built under R version 3.6.1
#> Warning: package 'BiocGenerics' was built under R version 3.6.1
#> Warning: package 'IRanges' was built under R version 3.6.1
#> Warning: package 'Biobase' was built under R version 3.6.1
#> Warning in plotTracks(list(ideogram.track, genome.axis.track, gene.track, : The
#> track chromosomes in 'trackList' differ. Setting all tracks to chromosome 'chr2'

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.

# 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
        );
    }
}
#> Warning in plotTracks(c(list(ideogram.track, genome.axis.track), read.tracks, :
#> The track chromosomes in 'trackList' differ. Setting all tracks to chromosome
#> 'chr2'

Visualising Interactions

2021-11-03

Single restriction fragment

Longer bait regions