The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Visualising Interactions

2021-11-03

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.

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.

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.