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.
