R/plot-methods.R
In IoannisVardaxis/MACPET: Model based analysis for paired-end data
Defines functions
plot.PInter
plot.PIntra
plot.PSelf
plot.PSFit
Documented in
plot.PInter
plot.PIntra
plot.PSelf
plot.PSFit
#' @title plot methods for MACPET classes
#' @author Ioannis Vardaxis, \email{ioannis.vardaxis@@ntnu.no}
#'
#' @references
#' Vardaxis I, Drabløs F, Rye M and Lindqvist BH (2018). \emph{MACPET: Model-based Analysis for ChIA-PET}.
#' To be published.
#'
#' Vardaxis I, Drabløs F, Rye M and Lindqvist BH (2018). \emph{MACPET: Complete pipeline for ChIA-PET}.
#' To be published.
#'
#' @description Different plot methods for the classes in the
#' \code{\link{MACPET}} package.
#'
#' @param x An object of correct class used to create different plots.
#' @param ... further arguments to be passed in the plot functions.
#'
#' @seealso \code{\linkS4class{PSelf}},
#' \code{\linkS4class{PSFit}}
#' \code{\linkS4class{PInter}}, \code{\linkS4class{PIntra}},
#' \code{\linkS4class{GenomeMap}}
#'
#' @name plot
#' @include AllClasses.R
#' @importFrom plyr ddply . ldply
#' @importFrom utils methods
#' @importFrom gtools mixedsort
#' @importFrom S4Vectors metadata
NULL
# > NULL
#---PInter:
#' @rdname plot
#' @method plot PInter
#' @return For the \code{\linkS4class{PInter}} class:
#' A network plot.
#' Each node is a chromosome with size proportional to the total PETs of the
#' corresponding chromosome. Edges connect chromosomes which have common PETs,
#' where the thickness of an edge is proportional on the total number of PETs
#' connecting the two chromosomes.
#' @export
#'
#' @examples
#'
#' #load Inter-chromosomal data:
#' load(system.file('extdata', 'MACPET_pinterData.rda', package = 'MACPET'))
#' class(MACPET_pinterData)
#' requireNamespace('igraph')
#' plot(MACPET_pinterData)
plot.PInter = function(x, ...) {
# global variables for Rcheck:
V1 = NULL
#--------------------------
#-------check package:
if (!requireNamespace("igraph", quietly = TRUE)) {
stop("igraph needed for this function to work. Please install it.", call. = FALSE)
}
PETcounts = S4Vectors::metadata(x)$InteractionCounts
PETcounts = as.matrix(PETcounts)
nodes = PETcounts
nodes = as.data.frame(colSums(nodes))
nodes$sepnames1 = rownames(nodes)
colnames(nodes) = c("V1", "from")
nodes = nodes[, c("from", "V1")]
nodes$V1 = nodes$V1/max(nodes$V1)
nodes$from = as.character(nodes$from)
# network plot:
name.edges = colnames(PETcounts)
edges = split(PETcounts, rep(seq_len(ncol(PETcounts)), each = nrow(PETcounts)))
edges = plyr::ldply(seq_len(length(edges)), function(i, name.edges, edges) {
data.frame(from = name.edges[i], to = name.edges, V1 = edges[[i]])
}, name.edges = name.edges, edges = edges)
edges$V1 = edges$V1/max(edges$V1)
edges = subset(edges, V1 != 0)
edges$from = as.character(edges$from)
edges$to = as.character(edges$to)
net = igraph::graph_from_data_frame(d = edges, vertices = nodes, directed = TRUE)
res = igraph::plot.igraph(net, edge.color = "blue", vertex.color = "red", edge.arrow.size = 0,
vertex.size = igraph::V(net)$V1 * 10, edge.width = igraph::E(net)$V1 * 2,
main = "Inter Interaction Network Plot")
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PInter", plot.PInter)
#---Intra PETs:
#' @rdname plot
#' @method plot PIntra
#' @return For the \code{\linkS4class{PIntra}} class:
#' A bar-plot. Each bar
#' represents the total number of Intra-chromosomal PETs for each chromosome in
#' the data.
#' @export
#'
#' @examples
#'
#' #load Intra-chromosomal data:
#' load(system.file('extdata', 'MACPET_pintraData.rda', package = 'MACPET'))
#' class(MACPET_pintraData)
#' requireNamespace('ggplot2')
#' plot(MACPET_pintraData)
plot.PIntra = function(x, ...) {
# global variables for Rcheck:
Chrom = NULL
#--------------------------
if (!requireNamespace("ggplot2", quietly = TRUE)) {
stop("ggplot2 needed for this function to work. Please install it.", call. = FALSE)
}
x = S4Vectors::metadata(x)$InteractionCounts
x = data.frame(Chrom = rep(x$Chrom, x$Counts))
res = ggplot2::ggplot(x, ggplot2::aes(x = Chrom, fill = Chrom)) + ggplot2::geom_bar() +
ggplot2::xlab("Chromosome") + ggplot2::ylab("Intra count") + ggplot2::ggtitle("Intra-chromosomal PET counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PIntra", plot.PIntra)
#---Self PETs:
#' @rdname plot
#' @method plot PSelf
#' @return For the \code{\linkS4class{PSelf}} class:
#' A bar-plot. Each bar
#' represents the total number of Self-ligated PETs for each chromosome in
#' the data.
#' @export
#'
#' @examples
#'
#' #load Self-ligated data:
#' load(system.file('extdata', 'MACPET_pselfData.rda', package = 'MACPET'))
#' class(MACPET_pselfData)
#' requireNamespace('ggplot2')
#' plot(MACPET_pselfData)
plot.PSelf = function(x, ...) {
# global variables for Rcheck:
Chrom = NULL
#--------------------------
if (!requireNamespace("ggplot2", quietly = TRUE)) {
stop("ggplot2 needed for this function to work. Please install it.", call. = FALSE)
}
x = S4Vectors::metadata(x)$Self_info
x = data.frame(Chrom = rep(x$Chrom, x$PET.counts))
res = ggplot2::ggplot(x, ggplot2::aes(x = Chrom, fill = Chrom)) + ggplot2::geom_bar() +
ggplot2::xlab("Chromosome") + ggplot2::ylab("PET counts") + ggplot2::ggtitle("Self-ligated PET counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PSelf", plot.PSelf)
#----PSFit:
#' @rdname plot
#' @method plot PSFit
#' @return For the \code{\linkS4class{PSFit}} class:
#' Different plots depenting on the \code{kind} argument.
#' @param RegIndex an integer indicating which region to plot (1 means the biggest
#' in terms of total PETs.)
#' @param threshold The FDR cut-off when plotting the total significant peaks/interactions for each chromosome in the data.
#' @param kind A string with one of the following arguments. Note that if a region visualization is plotted, the vertical lines represent peak-summits.
#' \describe{
#' \item{\code{PETcounts}}{ For a bar-plot of the PET-counts in each chromosome.}
#' \item{\code{RegionCounts}}{ For a bar-plot for the region counts in each chromosome.}
#' \item{\code{PeakCounts}}{ For a bar-plot for the Peak-counts in each chromosome.}
#' \item{\code{RegionPETs}}{ For a ggplot for a visualization of the PETs in a region.}
#' \item{\code{RegionTags}}{ For a ggplot for a visualization of the tags in a region. The tags are classified by stream (upper/lower)}
#' \item{\code{PeakPETs}}{ For a ggplot for a visualization of the PETs in a region. The PETs are classified by the peak they belong to.}
#' \item{\code{PeakTags}}{ For a ggplot for a visualization of the tags in a region. The tags are classified by the peak they belong to.}
#' \item{\code{SigPETCounts}}{ For a bar-plot with the significant PET-counts in each chromosome.}
#' \item{\code{SigRegionCounts}}{For a bar-plot with the significant region-counts in each chromosome.}
#' \item{\code{SigPeakCounts}}{For a bar-plot with the significant peak-counts in each chromosome.}
#'
#' }
#'
#' @export
#'
#' @examples
#'
#' #load Self-ligated data:
#' load(system.file('extdata', 'MACPET_psfitData.rda', package = 'MACPET'))
#' class(MACPET_psfitData)
#' requireNamespace('ggplot2')
#' plot(MACPET_psfitData,kind='PETcounts')
#' plot(MACPET_psfitData,kind='PeakCounts')
#' plot(MACPET_psfitData,kind='PeakPETs',RegIndex=1)
#' plot(MACPET_psfitData,kind='PeakTags',RegIndex=1)
plot.PSFit = function(x, kind, RegIndex = NULL, threshold = NULL, ...) {
# global variables for Rcheck:
FDR = RegCount = X = Y = ymin = ymax = Dist = Tag = Stream = PeakID = NULL
#--------------------------
# check that package exists:
if (!requireNamespace("ggplot2", quietly = TRUE)) {
stop("ggplot2 needed for this function to work if
create.self_intra.image==T. Please install it.",
call. = FALSE)
}
if (!kind %in% c("PETcounts", "RegionCounts", "PeakCounts", "RegionPETs", "RegionTags",
"PeakPETs", "PeakTags", "SigPETCounts", "SigRegionCounts", "SigPeakCounts")) {
stop("kind has been given wrong value!", call. = FALSE)
}
if (!methods::is(threshold, "numeric"))
threshold = NULL
Peaks.Info = S4Vectors::metadata(x)$Peaks.Info
if (!is.null(threshold))
Peaks.Info = subset(Peaks.Info, FDR < threshold)
if (nrow(Peaks.Info) == 0)
stop("threshold too low, try a lower one!", call. = FALSE)
if (kind == "PETcounts") {
# plot PET count for chromosomes:
class(x) = "PSelf"
res = plot.PSelf(x = x)
} else if (kind == "SigPETCounts") {
# plot significant pet counts:
SigPETCounts = plyr::ddply(Peaks.Info, plyr::.(Chrom), function(y) sum(y$Pets))
SigPETCounts = data.frame(Chrom = rep(SigPETCounts$Chrom, SigPETCounts$V1))
res = ggplot2::ggplot(SigPETCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Significant PET counts") +
ggplot2::ggtitle("significant Self-ligated PET counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "RegionCounts") {
# plot region count for chromosomes:
RegionCounts = S4Vectors::metadata(x)$Self_info
RegionCounts = data.frame(Chrom = rep(RegionCounts$Chrom, RegionCounts$Region.counts))
res = ggplot2::ggplot(RegionCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Region counts") +
ggplot2::ggtitle("Self-ligated region counts by chromosome") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "SigRegionCounts") {
# plot significant region count for chromosomes:
SigRegionCounts = plyr::ddply(Peaks.Info, plyr::.(Chrom), function(y) nrow(y))
SigRegionCounts = data.frame(Chrom = rep(SigRegionCounts$Chrom, SigRegionCounts$V1))
res = ggplot2::ggplot(SigRegionCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Region counts") +
ggplot2::ggtitle("Significant Self-ligated region counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "PeakCounts") {
PeakCounts = S4Vectors::metadata(x)$Self_info
PeakCounts = data.frame(Chrom = rep(PeakCounts$Chrom, PeakCounts$Peak.counts))
res = ggplot2::ggplot(PeakCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Peak counts") +
ggplot2::ggtitle("Self-ligated peak counts by chromosome") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "SigPeakCounts") {
SigPeakCounts = plyr::ddply(Peaks.Info, plyr::.(Chrom), function(y) nrow(y))
SigPeakCounts = data.frame(Chrom = rep(SigPeakCounts$Chrom, SigPeakCounts$V1))
res = ggplot2::ggplot(SigPeakCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Peak counts") +
ggplot2::ggtitle("Significant Self-ligated peak counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind %in% c("RegionPETs", "RegionTags", "PeakPETs", "PeakTags")) {
# (with summits):
#----choose the region to plot:
Classification.Info = S4Vectors::metadata(x)$Classification.Info
xdf = as.data.frame(x)
Chrom = xdf$seqnames1
Chrom = as.character(Chrom)
Classification.Info$Chrom = Chrom[Classification.Info$MainIndex]
Classification.Info$RegCount = paste(Classification.Info$Region, "-", Classification.Info$Chrom,
sep = "")
MaxReg = table(Classification.Info$RegCount)
MaxReg = sort(MaxReg, decreasing = TRUE)
if (is.null(RegIndex) | !is.numeric(RegIndex)) {
# take max:
which.id.region = max(which(MaxReg == max(MaxReg)))
which.id.region = names(MaxReg)[which.id.region]
} else {
which.id.region = names(MaxReg)[min(RegIndex, length(MaxReg))]
}
# take subset which will be plotted:
Classification.Info = subset(Classification.Info, RegCount == which.id.region)
xsub = xdf[Classification.Info$MainIndex, ]
xsub$Peak.ID = 0
xsub$Peak.ID = Classification.Info$Peak.ID
xsub = xsub[, c("start1", "end1", "start2", "end2", "Peak.ID")]
# sort:
Tosort = which(xsub$start1 > xsub$start2)
if (length(Tosort) > 0) {
Start2new = xsub$start1[Tosort]
End2new = xsub$end1[Tosort]
strand2new = xsub$strand1[Tosort]
xsub$start1[Tosort] = xsub$start2[Tosort]
xsub$end1[Tosort] = xsub$end2[Tosort]
xsub$strand1[Tosort] = xsub$strand2[Tosort]
xsub$start2[Tosort] = Start2new
xsub$end2[Tosort] = End2new
xsub$strand2[Tosort] = strand2new
}
xsub$Dist = xsub$end2 - xsub$start1 + 1
# take the peak summits of the region:
Peaks.Info = S4Vectors::metadata(x)$Peaks.Info
Peaks.Info$RegCount = paste(Peaks.Info$Region, "-", Peaks.Info$Chrom, sep = "")
Peaks.Info = subset(Peaks.Info, RegCount == which.id.region)
if (nrow(Peaks.Info) == 0) {
Peaks.Info = NULL
}
#-----plot according to what is asked:
if (kind == "RegionPETs") {
# plot region PETs, no strand info.
Dens = data.frame(Y = (xsub$start1 + xsub$end2)/2, ymin = xsub$start1,
ymax = xsub$end2, X = xsub$Dist)
res = ggplot2::ggplot(Dens, ggplot2::aes(x = X, y = Y, ymin = ymin, ymax = ymax)) +
ggplot2::geom_errorbar(width = 35) + ggplot2::coord_flip() + ggplot2::ggtitle("Visualization of PETs in region") +
ggplot2::ylab("Midpoints of PETs") + ggplot2::xlab("PET sizes") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
# add peaks
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of PETs in region. Vertical lines correspond to peak-summits.")
}
} else if (kind == "RegionTags") {
# plot Region tags(with summits):
Dens = data.frame(Tag = c((xsub$start1 + xsub$end1)/2, (xsub$start2 +
xsub$end2)/2), ymin = c(xsub$start1, xsub$start2), ymax = c(xsub$end1,
xsub$end2), Dist = c(xsub$Dist, xsub$Dist), Stream = c(rep("Upper",
nrow(xsub)), rep("Lower", nrow(xsub))))
res = ggplot2::ggplot(Dens, ggplot2::aes(x = Dist, y = Tag, ymin = ymin,
ymax = ymax, color = factor(Stream))) + ggplot2::geom_errorbar(width = 10) +
ggplot2::coord_flip() + ggplot2::ggtitle("Visualization of Tags in region") +
ggplot2::ylab("Midpoints of Tags") + ggplot2::xlab("PET sizes") +
ggplot2::labs(color = "Stream") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of Tags in region. Vertical lines correspond to peak-summits.")
}
} else if (kind == "PeakPETs") {
# (with summits) plot region PETs, no strand info.
Dens = data.frame(Y = (xsub$start1 + xsub$end2)/2, ymin = xsub$start1,
ymax = xsub$end2, X = xsub$Dist, PeakID = xsub$Peak.ID)
res = ggplot2::ggplot(Dens, ggplot2::aes(x = X, y = Y, ymin = ymin, ymax = ymax,
color = factor(PeakID))) + ggplot2::geom_errorbar(width = 35) + ggplot2::coord_flip() +
ggplot2::ggtitle("Visualization of PETs in region") + ggplot2::ylab("Midpoints of PETs") +
ggplot2::xlab("PET sizes") + ggplot2::labs(color = "PeakID") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of PETs in region. Vertical lines correspond to peak-summits.")
}
} else if (kind == "PeakTags") {
# (with summits)
Dens = data.frame(Tag = c((xsub$start1 + xsub$end1)/2, (xsub$start2 +
xsub$end2)/2), ymin = c(xsub$start1, xsub$start2), ymax = c(xsub$end1,
xsub$end2), Dist = c(xsub$Dist, xsub$Dist), PeakID = c(xsub$Peak.ID,
xsub$Peak.ID))
res = ggplot2::ggplot(Dens, ggplot2::aes(x = Dist, y = Tag, ymin = ymin,
ymax = ymax, color = factor(PeakID))) + ggplot2::geom_errorbar(width = 10) +
ggplot2::coord_flip() + ggplot2::ggtitle("Visualization of Tags in region") +
ggplot2::ylab("Midpoints of Tags") + ggplot2::xlab("PET sizes") +
ggplot2::labs(color = "PeakID") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of Tags in region. Vertical lines correspond to peak-summits.")
}
}
}
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PSFit", plot.PSFit)
#---GenomeMap:
#' @rdname plot
#' @method plot GenomeMap
#' @return For the \code{\linkS4class{GenomeMap}} class:
#' A heatmap plot or different kinds of network plots showing the relations between
#' different chromosomes based on the interactions between them.
#' @param Type A string with one of the following values:
#' \code{'heatmap'}, \code{'network-random'},\code{'network-circle'},\code{'network-sphere'}.
#' Each one shows the relations between the chromosomes based on their interactions.
#' @export
#'
#' @examples
#'
#' #load Interactions data:
#' load(system.file('extdata', 'MACPET_GenomeMapData.rda', package = 'MACPET'))
#' class(MACPET_GenomeMapData)
#' requireNamespace('igraph')
#' requireNamespace('reshape2')
#' plot(MACPET_GenomeMapData,Type='network-circle')
plot.GenomeMap = function(x, Type, threshold = NULL, ...) {
# global variables for Rcheck:
seqnames1 = NULL
seqnames2 = NULL
value = NULL
V1 = NULL
#--------------------------
#-------check package:
if (!methods::is(Type, "character")) {
stop("Type variable has to be a character", call. = FALSE)
} else if (!Type %in% c("heatmap", "network-random", "network-circle", "network-sphere")) {
stop("Type variable has to be: 'heatmap', 'network-random', 'network-circle' or 'network-sphere'.",
call. = FALSE)
}
# take the data:
object = GetSignInteractions(object = x, threshold = threshold, ReturnedAs = "GInteractions")
#--------------------------
# take cases:
#--------------------------
# make the CountsMat:
object = as.data.frame(object)[, c("seqnames1", "seqnames2")]
CountsMat = table(object)
if (Type == "heatmap") {
# check package:
if (!requireNamespace("ggplot2", quietly = TRUE) | !requireNamespace("reshape2",
quietly = TRUE)) {
stop("ggplot2 and reshape2 are needed for this
function to work if heatmap==TRUE.
Please install it.",
call. = FALSE)
}
plot.data = reshape2::melt(as.matrix(CountsMat))
plot.data$value = plot.data$value/max(plot.data$value)
RES = ggplot2::ggplot(data = plot.data, ggplot2::aes(x = seqnames1, y = seqnames2,
fill = value)) + ggplot2::geom_tile(color = "red") + ggplot2::ggtitle("Heat-Map plot for Interactions") +
ggplot2::xlab("Chromosome") + ggplot2::ylab("Chromosome") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(RES)
} else if (Type %in% c("network-random", "network-circle", "network-sphere")) {
# check package
if (!requireNamespace("igraph", quietly = TRUE)) {
stop("igraph needed for this function to work. Please install it.", call. = FALSE)
}
# make data:
nodes = as.data.frame(colSums(CountsMat))
nodes$sepnames1 = rownames(nodes)
colnames(nodes) = c("V1", "from")
nodes = nodes[, c("from", "V1")]
nodes$V1 = nodes$V1/max(nodes$V1)
nodes$from = as.character(nodes$from)
# network plot:
name.edges = colnames(CountsMat)
edges = split(CountsMat, rep(seq_len(ncol(CountsMat)), each = nrow(CountsMat)))
edges = plyr::ldply(seq_len(length(edges)), function(i, name.edges, edges) {
data.frame(from = name.edges[i], to = name.edges, V1 = edges[[i]])
}, name.edges = name.edges, edges = edges)
edges$V1 = edges$V1/max(edges$V1)
edges = subset(edges, V1 != 0)
edges$from = as.character(edges$from)
edges$to = as.character(edges$to)
net = igraph::graph_from_data_frame(d = edges, vertices = nodes, directed = FALSE)
if (Type == "network-random") {
igraph::plot.igraph(net, edge.color = "blue", vertex.color = "red", main = "Interactions Network Plot")
} else if (Type == "network-circle") {
igraph::plot.igraph(net, layout = igraph::layout_in_circle, edge.curved = 0.3,
edge.color = "blue", vertex.color = "red")
} else if (Type == "network-sphere") {
igraph::plot.igraph(net, layout = igraph::layout_on_sphere(net), edge.color = "blue",
vertex.color = "red")
}
}
}
#defining the method as S4 method:
setMethod("plot", "GenomeMap", plot.GenomeMap)
IoannisVardaxis/MACPET documentation built on Aug. 9, 2019, 12:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.PInter plot.PIntra plot.PSelf plot.PSFit
Documented in plot.PInter plot.PIntra plot.PSelf plot.PSFit
#' @title plot methods for MACPET classes
#' @author Ioannis Vardaxis, \email{ioannis.vardaxis@@ntnu.no}
#'
#' @references
#' Vardaxis I, Drabløs F, Rye M and Lindqvist BH (2018). \emph{MACPET: Model-based Analysis for ChIA-PET}.
#' To be published.
#'
#' Vardaxis I, Drabløs F, Rye M and Lindqvist BH (2018). \emph{MACPET: Complete pipeline for ChIA-PET}.
#' To be published.
#'
#' @description Different plot methods for the classes in the
#' \code{\link{MACPET}} package.
#'
#' @param x An object of correct class used to create different plots.
#' @param ... further arguments to be passed in the plot functions.
#'
#' @seealso \code{\linkS4class{PSelf}},
#' \code{\linkS4class{PSFit}}
#' \code{\linkS4class{PInter}}, \code{\linkS4class{PIntra}},
#' \code{\linkS4class{GenomeMap}}
#'
#' @name plot
#' @include AllClasses.R
#' @importFrom plyr ddply . ldply
#' @importFrom utils methods
#' @importFrom gtools mixedsort
#' @importFrom S4Vectors metadata
NULL
# > NULL
#---PInter:
#' @rdname plot
#' @method plot PInter
#' @return For the \code{\linkS4class{PInter}} class:
#' A network plot.
#' Each node is a chromosome with size proportional to the total PETs of the
#' corresponding chromosome. Edges connect chromosomes which have common PETs,
#' where the thickness of an edge is proportional on the total number of PETs
#' connecting the two chromosomes.
#' @export
#'
#' @examples
#'
#' #load Inter-chromosomal data:
#' load(system.file('extdata', 'MACPET_pinterData.rda', package = 'MACPET'))
#' class(MACPET_pinterData)
#' requireNamespace('igraph')
#' plot(MACPET_pinterData)
plot.PInter = function(x, ...) {
# global variables for Rcheck:
V1 = NULL
#--------------------------
#-------check package:
if (!requireNamespace("igraph", quietly = TRUE)) {
stop("igraph needed for this function to work. Please install it.", call. = FALSE)
}
PETcounts = S4Vectors::metadata(x)$InteractionCounts
PETcounts = as.matrix(PETcounts)
nodes = PETcounts
nodes = as.data.frame(colSums(nodes))
nodes$sepnames1 = rownames(nodes)
colnames(nodes) = c("V1", "from")
nodes = nodes[, c("from", "V1")]
nodes$V1 = nodes$V1/max(nodes$V1)
nodes$from = as.character(nodes$from)
# network plot:
name.edges = colnames(PETcounts)
edges = split(PETcounts, rep(seq_len(ncol(PETcounts)), each = nrow(PETcounts)))
edges = plyr::ldply(seq_len(length(edges)), function(i, name.edges, edges) {
data.frame(from = name.edges[i], to = name.edges, V1 = edges[[i]])
}, name.edges = name.edges, edges = edges)
edges$V1 = edges$V1/max(edges$V1)
edges = subset(edges, V1 != 0)
edges$from = as.character(edges$from)
edges$to = as.character(edges$to)
net = igraph::graph_from_data_frame(d = edges, vertices = nodes, directed = TRUE)
res = igraph::plot.igraph(net, edge.color = "blue", vertex.color = "red", edge.arrow.size = 0,
vertex.size = igraph::V(net)$V1 * 10, edge.width = igraph::E(net)$V1 * 2,
main = "Inter Interaction Network Plot")
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PInter", plot.PInter)
#---Intra PETs:
#' @rdname plot
#' @method plot PIntra
#' @return For the \code{\linkS4class{PIntra}} class:
#' A bar-plot. Each bar
#' represents the total number of Intra-chromosomal PETs for each chromosome in
#' the data.
#' @export
#'
#' @examples
#'
#' #load Intra-chromosomal data:
#' load(system.file('extdata', 'MACPET_pintraData.rda', package = 'MACPET'))
#' class(MACPET_pintraData)
#' requireNamespace('ggplot2')
#' plot(MACPET_pintraData)
plot.PIntra = function(x, ...) {
# global variables for Rcheck:
Chrom = NULL
#--------------------------
if (!requireNamespace("ggplot2", quietly = TRUE)) {
stop("ggplot2 needed for this function to work. Please install it.", call. = FALSE)
}
x = S4Vectors::metadata(x)$InteractionCounts
x = data.frame(Chrom = rep(x$Chrom, x$Counts))
res = ggplot2::ggplot(x, ggplot2::aes(x = Chrom, fill = Chrom)) + ggplot2::geom_bar() +
ggplot2::xlab("Chromosome") + ggplot2::ylab("Intra count") + ggplot2::ggtitle("Intra-chromosomal PET counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PIntra", plot.PIntra)
#---Self PETs:
#' @rdname plot
#' @method plot PSelf
#' @return For the \code{\linkS4class{PSelf}} class:
#' A bar-plot. Each bar
#' represents the total number of Self-ligated PETs for each chromosome in
#' the data.
#' @export
#'
#' @examples
#'
#' #load Self-ligated data:
#' load(system.file('extdata', 'MACPET_pselfData.rda', package = 'MACPET'))
#' class(MACPET_pselfData)
#' requireNamespace('ggplot2')
#' plot(MACPET_pselfData)
plot.PSelf = function(x, ...) {
# global variables for Rcheck:
Chrom = NULL
#--------------------------
if (!requireNamespace("ggplot2", quietly = TRUE)) {
stop("ggplot2 needed for this function to work. Please install it.", call. = FALSE)
}
x = S4Vectors::metadata(x)$Self_info
x = data.frame(Chrom = rep(x$Chrom, x$PET.counts))
res = ggplot2::ggplot(x, ggplot2::aes(x = Chrom, fill = Chrom)) + ggplot2::geom_bar() +
ggplot2::xlab("Chromosome") + ggplot2::ylab("PET counts") + ggplot2::ggtitle("Self-ligated PET counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PSelf", plot.PSelf)
#----PSFit:
#' @rdname plot
#' @method plot PSFit
#' @return For the \code{\linkS4class{PSFit}} class:
#' Different plots depenting on the \code{kind} argument.
#' @param RegIndex an integer indicating which region to plot (1 means the biggest
#' in terms of total PETs.)
#' @param threshold The FDR cut-off when plotting the total significant peaks/interactions for each chromosome in the data.
#' @param kind A string with one of the following arguments. Note that if a region visualization is plotted, the vertical lines represent peak-summits.
#' \describe{
#' \item{\code{PETcounts}}{ For a bar-plot of the PET-counts in each chromosome.}
#' \item{\code{RegionCounts}}{ For a bar-plot for the region counts in each chromosome.}
#' \item{\code{PeakCounts}}{ For a bar-plot for the Peak-counts in each chromosome.}
#' \item{\code{RegionPETs}}{ For a ggplot for a visualization of the PETs in a region.}
#' \item{\code{RegionTags}}{ For a ggplot for a visualization of the tags in a region. The tags are classified by stream (upper/lower)}
#' \item{\code{PeakPETs}}{ For a ggplot for a visualization of the PETs in a region. The PETs are classified by the peak they belong to.}
#' \item{\code{PeakTags}}{ For a ggplot for a visualization of the tags in a region. The tags are classified by the peak they belong to.}
#' \item{\code{SigPETCounts}}{ For a bar-plot with the significant PET-counts in each chromosome.}
#' \item{\code{SigRegionCounts}}{For a bar-plot with the significant region-counts in each chromosome.}
#' \item{\code{SigPeakCounts}}{For a bar-plot with the significant peak-counts in each chromosome.}
#'
#' }
#'
#' @export
#'
#' @examples
#'
#' #load Self-ligated data:
#' load(system.file('extdata', 'MACPET_psfitData.rda', package = 'MACPET'))
#' class(MACPET_psfitData)
#' requireNamespace('ggplot2')
#' plot(MACPET_psfitData,kind='PETcounts')
#' plot(MACPET_psfitData,kind='PeakCounts')
#' plot(MACPET_psfitData,kind='PeakPETs',RegIndex=1)
#' plot(MACPET_psfitData,kind='PeakTags',RegIndex=1)
plot.PSFit = function(x, kind, RegIndex = NULL, threshold = NULL, ...) {
# global variables for Rcheck:
FDR = RegCount = X = Y = ymin = ymax = Dist = Tag = Stream = PeakID = NULL
#--------------------------
# check that package exists:
if (!requireNamespace("ggplot2", quietly = TRUE)) {
stop("ggplot2 needed for this function to work if
create.self_intra.image==T. Please install it.",
call. = FALSE)
}
if (!kind %in% c("PETcounts", "RegionCounts", "PeakCounts", "RegionPETs", "RegionTags",
"PeakPETs", "PeakTags", "SigPETCounts", "SigRegionCounts", "SigPeakCounts")) {
stop("kind has been given wrong value!", call. = FALSE)
}
if (!methods::is(threshold, "numeric"))
threshold = NULL
Peaks.Info = S4Vectors::metadata(x)$Peaks.Info
if (!is.null(threshold))
Peaks.Info = subset(Peaks.Info, FDR < threshold)
if (nrow(Peaks.Info) == 0)
stop("threshold too low, try a lower one!", call. = FALSE)
if (kind == "PETcounts") {
# plot PET count for chromosomes:
class(x) = "PSelf"
res = plot.PSelf(x = x)
} else if (kind == "SigPETCounts") {
# plot significant pet counts:
SigPETCounts = plyr::ddply(Peaks.Info, plyr::.(Chrom), function(y) sum(y$Pets))
SigPETCounts = data.frame(Chrom = rep(SigPETCounts$Chrom, SigPETCounts$V1))
res = ggplot2::ggplot(SigPETCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Significant PET counts") +
ggplot2::ggtitle("significant Self-ligated PET counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "RegionCounts") {
# plot region count for chromosomes:
RegionCounts = S4Vectors::metadata(x)$Self_info
RegionCounts = data.frame(Chrom = rep(RegionCounts$Chrom, RegionCounts$Region.counts))
res = ggplot2::ggplot(RegionCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Region counts") +
ggplot2::ggtitle("Self-ligated region counts by chromosome") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "SigRegionCounts") {
# plot significant region count for chromosomes:
SigRegionCounts = plyr::ddply(Peaks.Info, plyr::.(Chrom), function(y) nrow(y))
SigRegionCounts = data.frame(Chrom = rep(SigRegionCounts$Chrom, SigRegionCounts$V1))
res = ggplot2::ggplot(SigRegionCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Region counts") +
ggplot2::ggtitle("Significant Self-ligated region counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "PeakCounts") {
PeakCounts = S4Vectors::metadata(x)$Self_info
PeakCounts = data.frame(Chrom = rep(PeakCounts$Chrom, PeakCounts$Peak.counts))
res = ggplot2::ggplot(PeakCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Peak counts") +
ggplot2::ggtitle("Self-ligated peak counts by chromosome") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind == "SigPeakCounts") {
SigPeakCounts = plyr::ddply(Peaks.Info, plyr::.(Chrom), function(y) nrow(y))
SigPeakCounts = data.frame(Chrom = rep(SigPeakCounts$Chrom, SigPeakCounts$V1))
res = ggplot2::ggplot(SigPeakCounts, ggplot2::aes(x = Chrom, fill = Chrom)) +
ggplot2::geom_bar() + ggplot2::xlab("Chromosome") + ggplot2::ylab("Peak counts") +
ggplot2::ggtitle("Significant Self-ligated peak counts by chromosome") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
} else if (kind %in% c("RegionPETs", "RegionTags", "PeakPETs", "PeakTags")) {
# (with summits):
#----choose the region to plot:
Classification.Info = S4Vectors::metadata(x)$Classification.Info
xdf = as.data.frame(x)
Chrom = xdf$seqnames1
Chrom = as.character(Chrom)
Classification.Info$Chrom = Chrom[Classification.Info$MainIndex]
Classification.Info$RegCount = paste(Classification.Info$Region, "-", Classification.Info$Chrom,
sep = "")
MaxReg = table(Classification.Info$RegCount)
MaxReg = sort(MaxReg, decreasing = TRUE)
if (is.null(RegIndex) | !is.numeric(RegIndex)) {
# take max:
which.id.region = max(which(MaxReg == max(MaxReg)))
which.id.region = names(MaxReg)[which.id.region]
} else {
which.id.region = names(MaxReg)[min(RegIndex, length(MaxReg))]
}
# take subset which will be plotted:
Classification.Info = subset(Classification.Info, RegCount == which.id.region)
xsub = xdf[Classification.Info$MainIndex, ]
xsub$Peak.ID = 0
xsub$Peak.ID = Classification.Info$Peak.ID
xsub = xsub[, c("start1", "end1", "start2", "end2", "Peak.ID")]
# sort:
Tosort = which(xsub$start1 > xsub$start2)
if (length(Tosort) > 0) {
Start2new = xsub$start1[Tosort]
End2new = xsub$end1[Tosort]
strand2new = xsub$strand1[Tosort]
xsub$start1[Tosort] = xsub$start2[Tosort]
xsub$end1[Tosort] = xsub$end2[Tosort]
xsub$strand1[Tosort] = xsub$strand2[Tosort]
xsub$start2[Tosort] = Start2new
xsub$end2[Tosort] = End2new
xsub$strand2[Tosort] = strand2new
}
xsub$Dist = xsub$end2 - xsub$start1 + 1
# take the peak summits of the region:
Peaks.Info = S4Vectors::metadata(x)$Peaks.Info
Peaks.Info$RegCount = paste(Peaks.Info$Region, "-", Peaks.Info$Chrom, sep = "")
Peaks.Info = subset(Peaks.Info, RegCount == which.id.region)
if (nrow(Peaks.Info) == 0) {
Peaks.Info = NULL
}
#-----plot according to what is asked:
if (kind == "RegionPETs") {
# plot region PETs, no strand info.
Dens = data.frame(Y = (xsub$start1 + xsub$end2)/2, ymin = xsub$start1,
ymax = xsub$end2, X = xsub$Dist)
res = ggplot2::ggplot(Dens, ggplot2::aes(x = X, y = Y, ymin = ymin, ymax = ymax)) +
ggplot2::geom_errorbar(width = 35) + ggplot2::coord_flip() + ggplot2::ggtitle("Visualization of PETs in region") +
ggplot2::ylab("Midpoints of PETs") + ggplot2::xlab("PET sizes") +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
# add peaks
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of PETs in region. Vertical lines correspond to peak-summits.")
}
} else if (kind == "RegionTags") {
# plot Region tags(with summits):
Dens = data.frame(Tag = c((xsub$start1 + xsub$end1)/2, (xsub$start2 +
xsub$end2)/2), ymin = c(xsub$start1, xsub$start2), ymax = c(xsub$end1,
xsub$end2), Dist = c(xsub$Dist, xsub$Dist), Stream = c(rep("Upper",
nrow(xsub)), rep("Lower", nrow(xsub))))
res = ggplot2::ggplot(Dens, ggplot2::aes(x = Dist, y = Tag, ymin = ymin,
ymax = ymax, color = factor(Stream))) + ggplot2::geom_errorbar(width = 10) +
ggplot2::coord_flip() + ggplot2::ggtitle("Visualization of Tags in region") +
ggplot2::ylab("Midpoints of Tags") + ggplot2::xlab("PET sizes") +
ggplot2::labs(color = "Stream") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of Tags in region. Vertical lines correspond to peak-summits.")
}
} else if (kind == "PeakPETs") {
# (with summits) plot region PETs, no strand info.
Dens = data.frame(Y = (xsub$start1 + xsub$end2)/2, ymin = xsub$start1,
ymax = xsub$end2, X = xsub$Dist, PeakID = xsub$Peak.ID)
res = ggplot2::ggplot(Dens, ggplot2::aes(x = X, y = Y, ymin = ymin, ymax = ymax,
color = factor(PeakID))) + ggplot2::geom_errorbar(width = 35) + ggplot2::coord_flip() +
ggplot2::ggtitle("Visualization of PETs in region") + ggplot2::ylab("Midpoints of PETs") +
ggplot2::xlab("PET sizes") + ggplot2::labs(color = "PeakID") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of PETs in region. Vertical lines correspond to peak-summits.")
}
} else if (kind == "PeakTags") {
# (with summits)
Dens = data.frame(Tag = c((xsub$start1 + xsub$end1)/2, (xsub$start2 +
xsub$end2)/2), ymin = c(xsub$start1, xsub$start2), ymax = c(xsub$end1,
xsub$end2), Dist = c(xsub$Dist, xsub$Dist), PeakID = c(xsub$Peak.ID,
xsub$Peak.ID))
res = ggplot2::ggplot(Dens, ggplot2::aes(x = Dist, y = Tag, ymin = ymin,
ymax = ymax, color = factor(PeakID))) + ggplot2::geom_errorbar(width = 10) +
ggplot2::coord_flip() + ggplot2::ggtitle("Visualization of Tags in region") +
ggplot2::ylab("Midpoints of Tags") + ggplot2::xlab("PET sizes") +
ggplot2::labs(color = "PeakID") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
if (!is.null(Peaks.Info)) {
res = res + ggplot2::geom_hline(yintercept = Peaks.Info$Peak.Summit,
color = "red", linetype = "dashed") + ggplot2::ggtitle("Visualization of Tags in region. Vertical lines correspond to peak-summits.")
}
}
}
return(res)
}
#defining the method as S4 method:
setMethod("plot", "PSFit", plot.PSFit)
#---GenomeMap:
#' @rdname plot
#' @method plot GenomeMap
#' @return For the \code{\linkS4class{GenomeMap}} class:
#' A heatmap plot or different kinds of network plots showing the relations between
#' different chromosomes based on the interactions between them.
#' @param Type A string with one of the following values:
#' \code{'heatmap'}, \code{'network-random'},\code{'network-circle'},\code{'network-sphere'}.
#' Each one shows the relations between the chromosomes based on their interactions.
#' @export
#'
#' @examples
#'
#' #load Interactions data:
#' load(system.file('extdata', 'MACPET_GenomeMapData.rda', package = 'MACPET'))
#' class(MACPET_GenomeMapData)
#' requireNamespace('igraph')
#' requireNamespace('reshape2')
#' plot(MACPET_GenomeMapData,Type='network-circle')
plot.GenomeMap = function(x, Type, threshold = NULL, ...) {
# global variables for Rcheck:
seqnames1 = NULL
seqnames2 = NULL
value = NULL
V1 = NULL
#--------------------------
#-------check package:
if (!methods::is(Type, "character")) {
stop("Type variable has to be a character", call. = FALSE)
} else if (!Type %in% c("heatmap", "network-random", "network-circle", "network-sphere")) {
stop("Type variable has to be: 'heatmap', 'network-random', 'network-circle' or 'network-sphere'.",
call. = FALSE)
}
# take the data:
object = GetSignInteractions(object = x, threshold = threshold, ReturnedAs = "GInteractions")
#--------------------------
# take cases:
#--------------------------
# make the CountsMat:
object = as.data.frame(object)[, c("seqnames1", "seqnames2")]
CountsMat = table(object)
if (Type == "heatmap") {
# check package:
if (!requireNamespace("ggplot2", quietly = TRUE) | !requireNamespace("reshape2",
quietly = TRUE)) {
stop("ggplot2 and reshape2 are needed for this
function to work if heatmap==TRUE.
Please install it.",
call. = FALSE)
}
plot.data = reshape2::melt(as.matrix(CountsMat))
plot.data$value = plot.data$value/max(plot.data$value)
RES = ggplot2::ggplot(data = plot.data, ggplot2::aes(x = seqnames1, y = seqnames2,
fill = value)) + ggplot2::geom_tile(color = "red") + ggplot2::ggtitle("Heat-Map plot for Interactions") +
ggplot2::xlab("Chromosome") + ggplot2::ylab("Chromosome") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(RES)
} else if (Type %in% c("network-random", "network-circle", "network-sphere")) {
# check package
if (!requireNamespace("igraph", quietly = TRUE)) {
stop("igraph needed for this function to work. Please install it.", call. = FALSE)
}
# make data:
nodes = as.data.frame(colSums(CountsMat))
nodes$sepnames1 = rownames(nodes)
colnames(nodes) = c("V1", "from")
nodes = nodes[, c("from", "V1")]
nodes$V1 = nodes$V1/max(nodes$V1)
nodes$from = as.character(nodes$from)
# network plot:
name.edges = colnames(CountsMat)
edges = split(CountsMat, rep(seq_len(ncol(CountsMat)), each = nrow(CountsMat)))
edges = plyr::ldply(seq_len(length(edges)), function(i, name.edges, edges) {
data.frame(from = name.edges[i], to = name.edges, V1 = edges[[i]])
}, name.edges = name.edges, edges = edges)
edges$V1 = edges$V1/max(edges$V1)
edges = subset(edges, V1 != 0)
edges$from = as.character(edges$from)
edges$to = as.character(edges$to)
net = igraph::graph_from_data_frame(d = edges, vertices = nodes, directed = FALSE)
if (Type == "network-random") {
igraph::plot.igraph(net, edge.color = "blue", vertex.color = "red", main = "Interactions Network Plot")
} else if (Type == "network-circle") {
igraph::plot.igraph(net, layout = igraph::layout_in_circle, edge.curved = 0.3,
edge.color = "blue", vertex.color = "red")
} else if (Type == "network-sphere") {
igraph::plot.igraph(net, layout = igraph::layout_on_sphere(net), edge.color = "blue",
vertex.color = "red")
}
}
}
#defining the method as S4 method:
setMethod("plot", "GenomeMap", plot.GenomeMap)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.