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