R/plotCalls.R

In rohan-shah/magicCalling: Genotype calling for multi-parent RILs

#' @title Plot the results of a hierarchical Bayesian model
#' @description Plot the results of a hierarchical Bayesian model
#' @param x An object of class hierarchicalBayesianModel (generated by \code{\link[mpMap2]{fitClusterModel}}) or magicHeuristicHBC (generated by \code{\link[mpMap2]{runHeuristics}}).
#' @param allFounderNames Optional list of the names of the founder lines replicates within the genetic data. 
#' @param chainIndex The index of the MCMC chain to show. 
#' @param ... Extra arguments to these functions are ignored.
#' @details
#' This function plots the result of using the end-point of one of the Markov chains as a fitted model. Homozygote clusters are shown in red and green, and the heterozygote cluster in dark blue. 
#' Light blue represents points falling in the outlier cluster. Black points are unclassified. Ovals represent the covariance matrixces of the homozygote and heterozygote clusters. 
#' If Input allfounderNames is input, then the founding lines of the population are given a different shape, and shown in black. 
#' @method plot hierarchicalBayesianModel
#' @S3method plot hierarchicalBayesianModel
#' @rdname plot.model
#' @examples
#' data("eightWayExampleData", package="magicCalling")
#' data <- eightWayExampleData[[1]]
#' meanY <- mean(data[,2])
#' startingPoints <- list(
#'  rbind(c(0.5, meanY), c(0.5, meanY)),
#'      rbind(c(0.5, meanY), c(0.5, meanY)),
#'      rbind(c(0.25, meanY), c(0.5, meanY)),
#'      rbind(c(0.25, meanY), c(0.5, meanY)),
#'      rbind(c(0.75, meanY), c(0.5, meanY)),
#'      rbind(c(0.75, meanY), c(0.5, meanY)),
#'      rbind(c(0.8, meanY), c(0.2, meanY)),
#'      rbind(c(0.8, meanY), c(0.2, meanY))
#' )
#' result <- fitClusterModel(data, startingPoints, n.iter = 200, D_hom = diag(2)*4, V_hom = cbind(c(0.005, 0), c(0, 0.1))/3, n_hom = 30, D_err = diag(2), V_err = diag(2)*10/3, n_err = 300, V_het = diag(2)*0.025/3, n_het = 1500)
#' plot(result, chainIndex = 1)
#' heuristicResults <- runHeuristics(result, minHomozygoteSize = 200)
#' plot(heuristicResults, chainIndex = heuristicResults$chainIndex)
#' @export 
plot.hierarchicalBayesianModel <- function(x, allFounderNames, chainIndex, ...)
{
  heuristicResults <- runHeuristics(x, minHomozygoteSize = 200)
  if(missing(chainIndex) && length(heuristicResults$classification) == 1)
  {
    chainIndex <- 1
  }
  if(missing(chainIndex))
  {
    stop("Input chainIndex is required, if there are multiple Markov chains")
  }
  data <- heuristicResults$data
  ellipse1 <- as.data.frame(ellipse(center = heuristicResults$clusterMeans[[chainIndex]][1, ], shape = heuristicResults$covariances[[chainIndex]][1,,], radius=5, col = 2, center.pch = NULL, draw = FALSE))
  ellipse2 <- as.data.frame(ellipse(center = heuristicResults$clusterMeans[[chainIndex]][2, ], shape = heuristicResults$covariances[[chainIndex]][2,,], radius=5, col = 3, center.pch = NULL, draw = FALSE))
  ellipse3 <- as.data.frame(ellipse(center = heuristicResults$clusterMeans[[chainIndex]][3, ], shape = heuristicResults$covariances[[chainIndex]][3,,], radius=5, col = 4, center.pch = NULL, draw = FALSE))
  if(missing(allFounderNames))
  {
    data <- data.frame(x = data[,1], y = data[,2], colour = factor(heuristicResults$classification[[chainIndex]], levels = 1:5), row.names = rownames(data))
    basePlot <- ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y"), data = data) + ggplot2::geom_point(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), shape = 16)
    built <- ggplot2::ggplot_build(basePlot)
    y.range <- built$layout$panel_ranges[[1]]$y.range
    x.range <- built$layout$panel_ranges[[1]]$x.range
    basePlot + ggplot2::geom_path(data = ellipse1, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 2) + ggplot2::geom_path(data = ellipse2, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 3) + ggplot2::geom_path(data = ellipse3, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 4) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + ggplot2::scale_shape_identity() + coord_cartesian(xlim = x.range, ylim = y.range, expand = FALSE)
  }
  else
  {
    data <- data.frame(x = data[,1], y = data[,2], colour = factor(heuristicResults$classification[[chainIndex]], levels = 1:5), shape = 16, row.names = rownames(data))
    founderIndices <- na.omit(match(allFounderNames, rownames(data)))
    data[founderIndices, "shape"] <- 2
    data[founderIndices, "colour"] <- 1
    #Split out the founder indices subset, because these will have different plot symbols, which we want to appear *on top* of everything else.
    founderIndicesSubset <- data[data[,"shape"] == 2, ]
    basePlot <-  ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = data) + ggplot2::geom_point(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour", shape = "shape")) + ggplot2::scale_shape_identity()
    built <- ggplot2::ggplot_build(basePlot)
    y.range <- built$layout$panel_ranges[[1]]$y.range
    x.range <- built$layout$panel_ranges[[1]]$x.range
    basePlot + ggplot2::geom_path(data = ellipse1, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 2) + ggplot2::geom_path(data = ellipse2, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 3) + ggplot2::geom_path(data = ellipse3, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 4) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + ggplot2::geom_point(data = founderIndicesSubset) + coord_cartesian(xlim = x.range, ylim = y.range, expand = FALSE)
  }
}
#' @rdname plot.model
#' @method plot magicHeuristicHBC
#' @S3method plot magicHeuristicHBC
#' @export 
plot.magicHeuristicHBC <- function(x, allFounderNames, chainIndex, ...)
{
  if(missing(chainIndex) && length(x$classifications) == 1)
  {
    chainIndex <- 1
  }
  if(missing(chainIndex) && x$chainIndex == -1)
  {
    stop("Input chainIndex is required, if there are multiple Markov chains, and no chain was assessed as being good quality.")
  }
  if(missing(chainIndex) && x$chainIndex != -1) chainIndex <- x$chainIndex
  heuristicResults <- x
  data <- heuristicResults$data
  ellipse1 <- as.data.frame(ellipse(center = heuristicResults$clusterMeans[[chainIndex]][1, ], shape = heuristicResults$covariances[[chainIndex]][1,,], radius=5, col = 2, center.pch = NULL, draw = FALSE))
  ellipse2 <- as.data.frame(ellipse(center = heuristicResults$clusterMeans[[chainIndex]][2, ], shape = heuristicResults$covariances[[chainIndex]][2,,], radius=5, col = 3, center.pch = NULL, draw = FALSE))
  ellipse3 <- as.data.frame(ellipse(center = heuristicResults$clusterMeans[[chainIndex]][3, ], shape = heuristicResults$covariances[[chainIndex]][3,,], radius=5, col = 4, center.pch = NULL, draw = FALSE))
  if(missing(allFounderNames))
  {
    data <- data.frame(x = data[,1], y = data[,2], colour = factor(heuristicResults$classifications[[chainIndex]], levels = 1:5), row.names = rownames(data))
    basePlot <- ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = data) + ggplot2::geom_point(shape = 16) 
    built <- ggplot2::ggplot_build(basePlot)
    y.range <- built$layout$panel_ranges[[1]]$y.range
    x.range <- built$layout$panel_ranges[[1]]$x.range
    basePlot + ggplot2::geom_path(data = ellipse1, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 2) + ggplot2::geom_path(data = ellipse2, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 3) + ggplot2::geom_path(data = ellipse3, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 4) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + ggplot2::scale_shape_identity() + coord_cartesian(xlim = x.range, ylim = y.range, expand = FALSE)
  }
  else
  {
    data <- data.frame(x = data[,1], y = data[,2], colour = factor(heuristicResults$classifications[[chainIndex]], levels = 1:5), shape = 16, row.names = rownames(data))
    founderIndices <- na.omit(match(allFounderNames, rownames(data)))
    data[founderIndices, "shape"] <- 2
    data[founderIndices, "colour"] <- 1
    #Split out the founder indices subset, because these will have different plot symbols, which we want to appear *on top* of everything else.
    founderIndicesSubset <- data[data[,"shape"] == 2,]
    basePlot <- ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = data) + ggplot2::geom_point(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour", shape = "shape")) + ggplot2::scale_shape_identity()
    built <- ggplot2::ggplot_build(basePlot)
    y.range <- built$layout$panel_ranges[[1]]$y.range
    x.range <- built$layout$panel_ranges[[1]]$x.range
    basePlot + ggplot2::geom_path(data = ellipse1, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 2) + ggplot2::geom_path(data = ellipse2, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 3) + ggplot2::geom_path(data = ellipse3, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 4) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + ggplot2::geom_point(data = founderIndicesSubset) + ggplot2::scale_shape_identity() + coord_cartesian(xlim = x.range, ylim = y.range, expand = FALSE)
  }
}
#' @title Plot a manually called marker
#' @description Plot the results of manually calling a marker
#' @param x An object of class markerResult. 
#' @param allFounderNames Optional list of the names of the founder lines replicates within the genetic data. 
#' @param ... Extra arguments to these functions are ignored.
#' @S3method plot markerResult
#' @details
#' This function plots the results of calling a marker manually using interactiveCall. In these cases a marker might have been called using 
#' the hierarchical Bayesia approach, or using DBSCAN. 
#' @export 
plot.markerResult <- function(x, allFounderNames, ...)
{
  data <- x$data
  if(x$hasVariability)
  {
    if("dbscan" %in% names(x) && x$dbscan)
    {
    	if(missing(allFounderNames))
	    {
          data <- data.frame(x = data[,1], y = data[,2], colour = factor(x$classification + 1, levels = 1:(max(x$classification)+1)))
          ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = data) + ggplot2::geom_point(shape = 16) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE)
	    }
	    else
	    {
        founderIndices <- na.omit(match(allFounderNames, rownames(data)))
        data <- data.frame(x = data[,1], y = data[,2], colour = factor(x$classification+1, levels = 1:(max(x$classification)+1)), shape = 16, row.names = rownames(data))
	      data[founderIndices, "shape"] <- 2
	      data[founderIndices, "colour"] <- 1
        #Split out the founder indices subset, because these will have different plot symbols, which we want to appear *on top* of everything else.
        founderIndicesSubset <- data[data[,"shape"] == 2,]
        ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour", shape = "shape"), data = data) + ggplot2::geom_point() + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + ggplot2::geom_point(data = founderIndicesSubset) + ggplot2::scale_shape_identity()
	    }
    }
    else
    {
      ellipse1 <- as.data.frame(ellipse(center = x$clusterMeans[1, ], shape = x$covariances[1,,], radius=5, col = 2, center.pch = NULL, draw = FALSE))
      ellipse2 <- as.data.frame(ellipse(center = x$clusterMeans[2, ], shape = x$covariances[2,,], radius=5, col = 3, center.pch = NULL, draw = FALSE))
      ellipse3 <- as.data.frame(ellipse(center = x$clusterMeans[3, ], shape = x$covariances[3,,], radius=5, col = 4, center.pch = NULL, draw = FALSE))
      if(missing(allFounderNames))
      {
        data <- data.frame(x = data[,1], y = data[,2], colour = factor(x$classification, levels = 1:5), row.names = rownames(data))
        basePlot <- ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = data) + ggplot2::geom_point(shape = 16) 
        built <- ggplot2::ggplot_build(basePlot)
        y.range <- built$layout$panel_ranges[[1]]$y.range
        x.range <- built$layout$panel_ranges[[1]]$x.range
	      basePlot + ggplot2::geom_path(data = ellipse1, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 2) + ggplot2::geom_path(data = ellipse2, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 3) + ggplot2::geom_path(data = ellipse3, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 4) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + coord_cartesian(xlim = x.range, ylim = y.range, expand = FALSE)
      }
      else
      {
        data <- data.frame(x = data[,1], y = data[,2], colour = factor(x$classification, levels = 1:5), shape = 16, row.names = rownames(data))
        founderIndices <- na.omit(match(allFounderNames, rownames(data)))
	      data[founderIndices, "shape"] <- 2
	      data[founderIndices, "colour"] <- 1
	      #Split out the founder indices subset, because these will have different plot symbols, which we want to appear *on top* of everything else.
	      founderIndicesSubset <- data[data[,"shape"] == 2,]
        basePlot <- ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = data) + ggplot2::geom_point(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour", shape = "shape")) + ggplot2::scale_shape_identity()
        built <- ggplot2::ggplot_build(basePlot)
        y.range <- built$layout$panel_ranges[[1]]$y.range
        x.range <- built$layout$panel_ranges[[1]]$x.range
	      notFoundersIndicesSubset <- data[data[,"shape"] != 2,]
	      ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour"), data = notFoundersIndicesSubset) + ggplot2::geom_point(mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour", shape = "shape")) + ggplot2::scale_shape_identity() + ggplot2::geom_path(data = ellipse1, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 2) + ggplot2::geom_path(data = ellipse2, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 3) + ggplot2::geom_path(data = ellipse3, mapping = ggplot2::aes_string(x = "x", y = "y"), colour = 4) + ggplot2::theme_bw() + ggplot2::scale_color_manual(values = palette(), guide = FALSE) + ggplot2::geom_point(data = founderIndicesSubset, mapping = ggplot2::aes_string(x = "x", y = "y", colour = "colour", shape = "shape")) + coord_cartesian(xlim = x.range, ylim = y.range, expand = FALSE)
      }
    }
  }
  else
  {
    data <- data.frame(x = data[,1], y = data[,2], row.names = rownames(data))
    ggplot2::ggplot(mapping = ggplot2::aes_string(x = "x", y = "y"), data = data) + ggplot2::geom_point(colour = 1, shape = 16) + ggplot2::scale_shape_identity()
  }
}