R/bexy.R

In bexy: Visualize and Parse the Output of 'BeXY'

#-------------------------------------------------------------------------------
# Package BeXY
# Visualize and parse the output of 'BeXY'
# Developed by Madleina Caduff
#-------------------------------------------------------------------------------

library(TeachingDemos)
library(Ternary)

#-------------------------------------------------------------------------------
# Internal methods, not exported
#-------------------------------------------------------------------------------

#' Function to open a file and generate error message if it was not found
#' @param path A file path.
#' @param files A vector of character strings corresponding to file names found in the path.
#' @param patterns A vector of patterns to search for within 'files'.
#' @return A file connection.
#' @keywords internal
.openFile.bexy <- function(path, files, patterns){
  filename <- files
  for (p in 1:length(patterns)){
    filename <- filename[grepl(patterns[p], filename)]
  }

  if (length(filename) == 0){
    stop(paste0("No file containing the pattern '", paste0(patterns, collapse = ","), "' found in directory '", path, "'!"))
  }
  if (length(filename) > 1){
    stop(paste0("Multiple files containing the pattern '", paste0(patterns, collapse = ","), "' found in directory '", path, "'!"))
  }
  if (file.size(paste0(path, filename)) == 3){ # empty file (happens if z was not inferred)
    return(data.frame())
  }
  f <- read.table(paste0(path, filename), header = TRUE, check.names = FALSE)
  return(f)
}

#' Function to reformat column names
#' @param file A file whose header should be renamed.
#' @param pattern Where to split.
#' @return A new file with re-formatted column names.
#' @keywords internal
.renameHeader.bexy <- function(file, pattern){
  new <- sapply(names(file), function(x) return(strsplit(x, pattern)[[1]][2]))
  names(file) <- new
  return(file)
}

#' Function that returns all possible sex karyotypes used in BeXY
#' @return A vector of characters.
#' @keywords internal
.getAllPossibleSexKaryotypes.bexy <- function(){
  return(c("XY", "XX", "X0", "XXY", "XYY", "XXX", "XXYY"))
}

#' Function that returns all possible scaffold types used in BeXY
#' @return A vector of characters.
#' @keywords internal
.getAllPossibleScaffoldTypes.bexy <- function(){
  return(c("Autosome", "Y-linked", "X-linked", "Different"))
}

#' Function that returns all possible autosomal trisomy states used in BeXY
#' @param scaffoldName A vector of size one with the name of a scaffold.
#' @return A vector of characters.
#' @keywords internal
.getAllPossibleAutosomalTrisomy.bexy <- function(scaffoldName){
  return(c("diploid", paste0("trisomy_", scaffoldName)))
}

#' Function that returns the scaffold names of z
#' @param names_z A character vector with the column names of z
#' @param numSamples The number of samples
#' @return A vector of characters.
#' @keywords internal
.getScaffoldNames.bexy <- function(names_z, numSamples){
  all <- sapply(names_z, function(x) strsplit(strsplit(x, "z_")[[1]][2], "_")[[1]][1])
  uq <- all[seq(1, length(all), by = numSamples)]
  return(as.character(uq))
}

#' Function that returns the index of a given scaffold & checks if it is valid for autosomal trisomy
#' @param object A bexy object.
#' @param scaffoldName A vector of size one with the name of a scaffold.
#' @return A list containing the index of the scaffold within posterior_mode_z and counts.
#' @keywords internal
.getIndexScaffoldAutosomalTrisomy.bexy <- function(object, scaffoldName){
  # check if trisomy is allowed
  if (!object$hasZ){
    stop("BeXY was not run with autosomal trisomy!
         Please re-run BeXY using the argument '--triploidAutosomes A', where A are the autosome names where trisomy is allowed")
  }
  # check if it is only a single scaffold
  if (length(scaffoldName) > 1){
    stop("Please specify only a single scaffold at a time.")
  }

  # get all possible autosome names
  autosomes <- names(object$posterior_t)[getPosteriorModeScaffoldTypes(object) == "Autosome"]
  errMessage <- paste0("Available autosomes are: ", paste0(autosomes, collapse = ","))

  # get index of match in scaffold names
  ix_counts <- which(colnames(object$posterior_t) == as.character(scaffoldName))

  # check if this is an autosome
  if (getPosteriorModeScaffoldTypes(object)[ix_counts] != "Autosome"){
    stop(paste0("Scaffold '", scaffoldName, "' was not classified as an autosome. ", errMessage, "."))
  }

  # get index of match in z
  ix_z      <- which(row.names(object$posterior_mode_z) == as.character(scaffoldName))

  # check if there is a match
  if (length(ix_z) == 0){
    stop(paste0("Could not detect output of autosomal trisomy for scaffold '", scaffoldName, "'.
                Please re-run BeXY with the argument '--triploidAutosomes ", scaffoldName, "'."))
  }


  return(list(ix_z=ix_z, ix_counts=ix_counts))
}

#' Function that returns the index of a given scaffold
#' @param object A bexy object.
#' @param scaffoldName A string representing the name of a scaffold.
#' @param takeIfNotFound A string ("X-linked" or "Y-linked") that denotes the type of scaffold to use if scaffoldName could not be found
#' @return A single integer representing the index of the scaffold.
#' @keywords internal
.getIndexScaffold.bexy <- function(object, scaffoldName, takeIfNotFound){
  # check if it is only a single scaffold
  if (length(scaffoldName) > 1){
    stop("Please specify only a single scaffold at a time.")
  }
  # check if takeIfNotFound is either "X-linked" or "Y-linked"
  if (!(takeIfNotFound == "X-linked" | takeIfNotFound == "Y-linked")){
    stop("takeIfNotFound must be either X-linked or Y-linked.")
  }

  # check if scaffoldName is NULL -> find first scaffold of type 'takeIfNotFound'
  if (is.null(scaffoldName)){
    postMode <- getPosteriorModeScaffoldTypes(object)
    all <- names(postMode[postMode == takeIfNotFound])
    scaffoldName <- all[1]
    if (length(all) > 1){
      warning(paste0("Multiple scaffolds are ", takeIfNotFound, " according to their posterior mode, and only the first one (", scaffoldName, ") will be plotted. Use the argument 'scaffoldX' or 'scaffoldY' to specify other scaffolds for plotting."))
    }
    if (is.na(scaffoldName)){
      # no 'takeIfNotFound' found -> return NULL
      warning(paste0("No scaffolds are ", takeIfNotFound, " according to their posterior mode. Will use sample index for this axis."))
      return(NULL)
    }
  }

  # get index of match in scaffold names
  ix <- which(object$scaffoldNames == as.character(scaffoldName))
  if (length(ix) == 0){
    stop(paste0("Scaffold '", scaffoldName, "' does not exist."))
  }

  return(ix)
}

#' Function that transforms state posteriors of sex karyotype to ternary-format
#' @param object A bexy object.
#' @return A matrix with 3 rows: P(aneuploid), P(XX), P(XY).
#' @keywords internal
.getInputTernary.bexy <- function(object){
  return(rbind(colSums(object$posterior_s[3:nrow(object$posterior_s),]), object$posterior_s[2,], object$posterior_s[1,]))
}

#' Function to select samples that should be highlighted
#' @param object A bexy object.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @return A vector of integers corresponding to indices of samples to highlight
#' @keywords internal
.getSamplesToHighlight.bexy <- function(object, sampleNamesToHighlight){
  if (length(sampleNamesToHighlight) == 0){
    # automatically based on probabilities
    return(which(object$posterior_s[1,] < 0.9 & object$posterior_s[2,] < 0.9))
  }
  # get index of match
  return(which(names(object$posterior_s) %in% sampleNamesToHighlight))
}

#' Function to select samples that should be highlighted for autosomal trisomy
#' @param object A bexy object.
#' @param ix A list containing the index of a specific scaffold within z.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @return A vector of integers corresponding to indices of samples to highlight
#' @keywords internal
.getSamplesToHighlightAutosomalTrisomy.bexy <- function(object, ix, sampleNamesToHighlight){
  if (length(sampleNamesToHighlight) == 0){
    return(which(object$posterior_z[ix$ix_z,] > 0.9))
  }
  # get index of match
  return(which(colnames(object$posterior_z) %in% sampleNamesToHighlight))
}

#' Function to create a legend
#' @param x_coord The x-coordinate for legend.
#' @param y_coord The y-coordinate for legend.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @return No return value, called for side effects.
#' @keywords internal
.addLegend.bexy <- function(x_coord, y_coord, colors){
  legend(x_coord, y_coord,
         legend = c("Aneuploid", "XX", "XY"),
         col = colors,
         pch = rep(19, 3),
         bty = "n")
}

#' Function to create a legend for autosomal trisomy
#' @param x_coord The x-coordinate for legend.
#' @param y_coord The y-coordinate for legend.
#' @param scaffoldName The name of the scaffold to plot (has to be an autosome).
#' @param colors A vector of length two with the colors for triploid and diploid samples, respectively.
#' @return No return value, called for side effects.
#' @keywords internal
.addLegendAutosomalTrisomy.bexy <- function(x_coord, y_coord, scaffoldName, colors){
  legend(x_coord, y_coord,
         legend = c(paste0("Trisomy ", scaffoldName), "Euploid"),
         col = colors,
         pch = rep(19, 2),
         bty = "n")
}

#' Function to get the width of a string on the current plotting device in user coordinates
#' @param name A string.
#' @return The width of the input string on the current plotting device in user coordinates.
#' @keywords internal
.getStringWidth.bexy <- function(name){
  vec <- strsplit(name, split = "")[[1]]
  return(sum(sapply(vec, graphics::strwidth)))
}

#' Function to plot a ternary plot of the probabilities for being aneuploid, XX or XY
#' @param object A bexy object.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... More parameters passed to TernaryPlot()
#' @return The upper boundary of the plot.
#' @keywords internal
.plotTernary.bexy <- function(object,  colors, sampleNamesToHighlight, label, ...){
  # order: aneuploid (top), XX (right), XY (left)
  tern <- .getInputTernary.bexy(object)
  Ternary::TernaryPlot(point = "up", atip = "", btip = "", ctip = "",
              atip.rotate = 0, atip.pos = 3,
              grid.col = "lightgrey", grid.minor.col = "white", # no grid lines
              axis.labels = FALSE, axis.tick = FALSE, ...
  )
  # make nice labels for triangle nodes
  x_range <- Ternary::TernaryXRange()
  y_range <- Ternary::TernaryYRange()
  top <- y_range[2]+0.005
  text(x = mean(x_range), y = top, labels = "Aneuploid", font = 2, cex = 1)
  text(x = x_range[1]-0.03, y = y_range[1]+0.03, labels = "XY", font = 2, cex = 1)
  text(x = x_range[2]+0.03, y = y_range[1]+0.03, labels = "XX", font = 2, cex = 1)

  posterior_mode <- apply(tern, 2, which.max)
  col_per_sample <- colors[posterior_mode]

  Ternary::TernaryPoints(tern, pch = 19, col = col_per_sample, cex = 1)

  # Add labels
  highlight <- .getSamplesToHighlight.bexy(object, sampleNamesToHighlight)
  if (label & length(highlight) > 0){
    pointXY <- Ternary::TernaryToXY(as.matrix(tern[,highlight]))

    bottom <- pointXY[2,] < 0.1
    left <- pointXY[1,] < 0 & !bottom
    right <- pointXY[1,] >= 0 & !bottom

    segment_length_left <- 0.08
    segment_length_right <- 0.18

    # spread x
    xx <- pointXY[1,]
    xx[bottom] <- TeachingDemos::spread.labs(xx[bottom], 2*strheight('A'), maxiter=1000, min=0)
    if (any(left)){
      name_width <- sapply(names(object$posterior_s)[highlight][left], .getStringWidth.bexy)
      xx[left] <- xx[left] - name_width - graphics::strwidth('A') - segment_length_left # add a small extra space (A)
    }
    xx[right] <- xx[right] + 0.2

    # spread y
    yy <- pointXY[2,]
    yy[bottom] <- yy[bottom] - 0.1
    yy[left] <- TeachingDemos::spread.labs(yy[left], 2*strheight('A'), maxiter=1000, min=0)
    yy[right] <- TeachingDemos::spread.labs(yy[right], 2*strheight('A'), maxiter=1000, min=0)

    # define segment length
    x_segment <- pointXY[1,]
    x_segment[left] <- x_segment[left] - segment_length_left
    x_segment[right] <- x_segment[right] + segment_length_right

    # add label
    text(xx, yy, names(object$posterior_s)[highlight], adj = 0, xpd = TRUE, font = 3)
    segments(pointXY[1,], pointXY[2,],
             x_segment, yy)
  }
  return(list(top = top))
}

#' Function to plot a single 7-cell plot displaying posterior probabilities of karyotypes
#' @param object A bexy object.
#' @param ix The sample index.
#' @param label The sample name used as a label.
#' @param x_coords Left and right x-coordinates for 7-cell plot.
#' @param y_coords Bottom and top y-coordinates for 7-cell plot.
#' @param last Boolean indicating if this is the last sample to be plotted.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @return No return value, called for side effects.
#' @keywords internal
.plotTicTacToeWithCoordinates.bexy <- function(object, ix, label, x_coords, y_coords, last, colors){
  # draw grid
  width <- diff(x_coords)
  height <- diff(y_coords)
  width_third <- width / 3
  height_third <- height / 3

  xx <- c(x_coords[1], x_coords[1] + width_third, x_coords[1] + 2*width_third, x_coords[2])
  yy <- c(y_coords[1], y_coords[1] + height_third, y_coords[1] + 2*height_third, y_coords[2])

  # draw vertical lines
  lines(x = rep(xx[1], 2), y = y_coords)
  lines(x = rep(xx[2], 2), y = y_coords)
  lines(x = rep(xx[3], 2), y = y_coords)
  lines(x = rep(xx[4], 2), y = c(yy[1], yy[2]))

  # draw horizontal lines
  lines(x = x_coords, y = rep(y_coords[1], 2))
  lines(x = x_coords, y = rep(y_coords[1] + height_third, 2))
  lines(x = c(xx[1], xx[3]), y = rep(yy[3], 2))
  lines(x = c(xx[1], xx[3]), y = rep(yy[4], 2))

  # color according to probabilities
  createCol <- function(col, prob) { return(grDevices::adjustcolor(col, alpha.f = prob)) }

  rect(xleft = xx[1], ybottom = yy[1], xright = xx[2], ytop = yy[2], col = createCol(colors[1], object$posterior_s[3,ix])) # X
  rect(xleft = xx[2], ybottom = yy[1], xright = xx[3], ytop = yy[2], col = createCol(colors[2], object$posterior_s[2,ix])) # XX
  rect(xleft = xx[3], ybottom = yy[1], xright = xx[4], ytop = yy[2], col = createCol(colors[1], object$posterior_s[6,ix])) # XXX
  rect(xleft = xx[1], ybottom = yy[2], xright = xx[2], ytop = yy[3], col = createCol(colors[3], object$posterior_s[1,ix])) # XY
  rect(xleft = xx[2], ybottom = yy[2], xright = xx[3], ytop = yy[3], col = createCol(colors[1], object$posterior_s[4,ix])) # XXY
  rect(xleft = xx[1], ybottom = yy[3], xright = xx[2], ytop = yy[4], col = createCol(colors[1], object$posterior_s[5,ix])) # XYY
  rect(xleft = xx[2], ybottom = yy[3], xright = xx[3], ytop = yy[4], col = createCol(colors[1], object$posterior_s[7,ix])) # XXYY

  # add sample name
  text(x = x_coords[2] + width_third, y = yy[3] - height_third/2, labels = label, adj = 0, font = 3)

  # add Y labels
  text(x = x_coords[1] - width_third/2, y = yy[2] - height_third/2, labels = "0", srt = 0, adj = 1, cex = 0.8)
  text(x = x_coords[1] - width_third/2, y = yy[3] - height_third/2, labels = "Y", srt = 0, adj = 1, cex = 0.8)
  text(x = x_coords[1] - width_third/2, y = yy[4] - height_third/2, labels = "YY", srt = 0, adj = 1, cex = 0.8)

  # add X labels
  if (last){
    text(x = xx[1] + width_third/2, y = yy[1] - height_third/2 , labels = "X", srt = 90, adj = 1, cex = 0.8)
    text(x = xx[2] + width_third/2, y = yy[1] - height_third/2 , labels = "XX", srt = 90, adj = 1, cex = 0.8)
    text(x = xx[3] + width_third/2, y = yy[1] - height_third/2 , labels = "XXX", srt = 90, adj = 1, cex = 0.8)
  }
}

#' Function to plot ternary plot with 7-cell plots displaying posterior probabilities of karyotypes
#' @param object A bexy object.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param addSquares Boolean indicating whether 7-cell posterior probability square plots should be drawn.
#' @param ... Other parameters used for plotting
#' @return No return value, called for side effects.
#' @keywords internal
.plotTernaryAndTicTacToe.bexy <- function(object, colors, sampleNamesToHighlight, label, addSquares, ...){
  highlight <- .getSamplesToHighlight.bexy(object, sampleNamesToHighlight)

  if (!addSquares | !label | length(highlight) == 0){
    # plot ternary, no tic-tac-toe
    .plotTernary.bexy(object, colors, sampleNamesToHighlight, label, ...)
    .addLegend.bexy(par("usr")[1], par("usr")[4], colors)
  } else {
    oldPar <- par(no.readonly = TRUE)
    on.exit(par(oldPar))
    par(mar = c(3, 1, 1, 10), xpd = TRUE)

    # plot ternary
    top <- .plotTernary.bexy(object, colors, sampleNamesToHighlight, label, ...)

    # determine coordinates for tictactoes
    x_left <- par("usr")[2] + 0.01
    y_bottom <- par("usr")[3]
    y_top <- top$top

    height <- (y_top - y_bottom) / (length(highlight) + 1)
    height <- min(height, 0.1)
    dist <- height / (length(highlight) - 1)
    width <- height

    # add tictactoes
    cur_y_top <- y_top
    cur_y_bottom <- y_top - height
    names <- names(object$posterior_s)[highlight]
    numTicTacToes <- min(length(highlight), 11) # make sure we only plot as many tictactoe as there is space
    for (i in 1:numTicTacToes){
      x_coords <- c(x_left, x_left + width)
      y_coords <- c(cur_y_bottom, cur_y_top)

      .plotTicTacToeWithCoordinates.bexy(object = object,
                                         ix = highlight[i],
                                         label = names[i],
                                         x_coords = x_coords,
                                         y_coords = y_coords,
                                         last = (i == length(highlight)),
                                         colors = colors)

      cur_y_top <- cur_y_bottom - dist
      cur_y_bottom <- cur_y_top - height
    }
    .addLegend.bexy(par("usr")[1], par("usr")[4], colors)
  }
}

#' Function to plot X vs Y counts and color them according to BeXY posterior mode
#' @param object A bexy object.
#' @param ix_X An integer representing the index of the scaffold to be plotted on the x-axis.
#' @param ix_Y An integer representing the index of the scaffold to be plotted on the y-axis.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... Other parameters used for plotting
#' @return No return value, called for side effects.
#' @keywords internal
.plotXVsY.bexy <- function(object, ix_X, ix_Y, colors, sampleNamesToHighlight, label, ...){
  # get colors
  posterior_mode <- apply(.getInputTernary.bexy(object), 2, which.max)
  col_per_sample <- colors[posterior_mode]

  # get point types
  pchs <- rep(19, length(col_per_sample))

  # plot normalized counts
  n <- object$counts
  n$individual <- NULL
  n$sequencing_type <- NULL
  norm_counts <- sweep(n, 1, rowSums(n), "/")
  xx <- norm_counts[, ix_X] * 100
  yy <- norm_counts[, ix_Y] * 100

  # get axis labels
  getAxLabel <- function(ix){ return(paste0("Reads mapped to ", object$scaffoldNames[ix], " (%, ", getPosteriorModeScaffoldTypes(object)[ix], ")")) }
  ax_X <- getAxLabel(ix_X)
  ax_Y <- getAxLabel(ix_Y)

  # check if X or Y is NULL (does not exist) -> in this case, plot sample index on x-axis
  if (is.null(ix_X) | is.null(ix_Y)){
    if (is.null(ix_Y)){
      yy <- xx
      ax_Y <- ax_X
    }
    xx <- 1:length(yy)
    ax_X <- "Sample index"
  }
  plot(xx, yy, col = col_per_sample, pch = pchs, xlab = ax_X, ylab = ax_Y, ...)
  box()

  # identify samples to highlight
  highlight <- .getSamplesToHighlight.bexy(object, sampleNamesToHighlight)
  if (label & length(highlight) > 0){
    label_y <- yy[highlight]

    # add sample labels
    breaks <- hist(xx[highlight], plot = F)$breaks
    breaks <- unique(round(breaks, digits = 1))
    tmp_x <- sapply(xx[highlight], function(this) return(breaks[which.min(abs(breaks - this))]))
    for (i in unique(tmp_x)){
      label_y[tmp_x == i] <- TeachingDemos::spread.labs(label_y[tmp_x == i],
                                                        2*strheight('A'),
                                                        maxiter=1000,
                                                        min=0,
                                                        max = par("usr")[4] - 0.1*(par("usr")[4]-par("usr")[3]))
    }
    # overshoot x-axis on the right?
    label_x <- xx[highlight] + 0.2
    segment_x <- xx[highlight] + 0.15
    name_width <- sapply(names(object$posterior_s)[highlight], .getStringWidth.bexy)
    overshoot <- which(label_x + name_width > par("usr")[2])
    if (length(overshoot) > 0){
      label_y[overshoot] <- label_y[overshoot] + 2*strheight('A') # move up
      label_x[overshoot] <- par("usr")[2] - name_width[overshoot] - graphics::strwidth('A') # move left (with a bit extra: A)
      segment_x[overshoot] <- par("usr")[2] - graphics::strwidth('A')
    }

    # label!
    text(label_x, label_y,
         names(object$posterior_s)[highlight], adj = 0, xpd = T, font = 3)
    segments(xx[highlight], yy[highlight],
             segment_x, label_y )
  }

  # add legend
  .addLegend.bexy(par("usr")[1] + 0.75*(par("usr")[2] - par("usr")[1]),
                  par("usr")[4],
                  colors)
}

#' Function to plot counts of a specific autosome and color them according to BeXY posterior mode of autosomal trisomy
#' @param object A bexy object.
#' @param scaffoldName The name of the scaffold to plot (has to be an autosome)
#' @param plotCounts A boolean. If TRUE, the counts are plotted. If FALSE, the posterior probabilities are plotted.
#' @param colors A vector of length two with the colors for diploid and triploid samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... Other parameters used for plotting
#' @return No return value, called for side effects.
#' @keywords internal
.plotAutosomalTrisomy.bexy <- function(object, scaffoldName, plotCounts, colors, sampleNamesToHighlight, label, ...){
  ix <- .getIndexScaffoldAutosomalTrisomy.bexy(object, scaffoldName)

  # get colors
  posterior_mode <- object$posterior_mode_z[ix$ix_z,] + 1
  col_per_sample <- colors[posterior_mode]

  # get point types
  pchs <- rep(19, length(col_per_sample))


  if (plotCounts){  # plot normalized counts
    n <- object$counts
    n$individual <- NULL
    n$sequencing_type <- NULL
    norm_counts <- sweep(n, 1, rowSums(n), "/")
    xx <- 1:nrow(n)
    yy <- norm_counts[,ix$ix_counts] * 100
    plot(yy, col = col_per_sample, pch = pchs,
         xlab = "Sample",
         ylab = paste0("Reads mapped to scaffold ", scaffoldName, " (%)"))
    box()
  } else { # plot posterior probabilities
    xx <- 1:ncol(object$posterior_z)
    yy <- object$posterior_z[ix$ix_z,] * 100
    plot(yy, col = col_per_sample, pch = pchs,
         xlab = "Sample",
         ylab = paste0("Posterior probability of trisomy ", scaffoldName, " (%)"),
         ...)
  }


  # identify samples to highlight
  highlight <- .getSamplesToHighlightAutosomalTrisomy.bexy(object, ix, sampleNamesToHighlight)
  if (label & length(highlight) > 0){
    label_y <- yy[highlight]

    # add sample labels
    breaks <- hist(xx[highlight], plot = FALSE)$breaks
    breaks <- c(breaks[1] - diff(breaks[1:2]), breaks)
    breaks <- c(breaks, breaks[length(breaks)] + diff(breaks[1:2]))
    tmp_x <- sapply(xx[highlight], function(this) return(breaks[this >= breaks & this < (breaks+1)])) # ceiling(xx[highlight]*2) / 2
    for (i in unique(tmp_x)){
      label_y[tmp_x == i] <- TeachingDemos::spread.labs(label_y[tmp_x == i],
                                                        2*strheight('A'),
                                                        maxiter=1000,
                                                        min=0 )
    }
    text(xx[highlight]+length(xx)*0.05, label_y,
         names(object$posterior_s)[highlight], adj = 0, xpd = TRUE)
    segments(xx[highlight], yy[highlight],
             xx[highlight]+length(xx)*0.04, label_y )
  }

  # add legend
  .addLegendAutosomalTrisomy.bexy(par("usr")[1],
                                  par("usr")[3] + 0.2*(par("usr")[4] - par("usr")[3]),
                                  scaffoldName,
                                  colors)
}

#' Function to extract specific columns from a data frame or matrix
#' @param what A matrix or data frame
#' @param name The name of the parameter to extract.
#' @return A subset of the matrix or data frame containing only the columns whose names match the argument 'name'
#' @keywords internal
.getColumns.bexy <- function(what, name){
  return(what[,which(grepl(paste0("^", name), colnames(what)))])
}

#-------------------------------------------------------------------------------
# Constructor for BeXY
#-------------------------------------------------------------------------------

#' Accurate Bayesian inference of sex chromosome karyotypes and sex-linked scaffolds from low-depth sequencing data
#'
#' Visualization of results produced by BeXY
#' @param path The path where all the output files of BeXY are located. If NULL, the example files will be loaded.
#' @param readMCMCTrace If TRUE, read the full trace of the MCMC of all parameters (potentially slow).
#' @return An object of type bexy.
#' @export

bexy <- function(path = NULL, readMCMCTrace = FALSE){
  if (is.null(path)) {
    # use example files
    path <- paste0(system.file("extdata", package = "bexy"), "/")
    files <- dir(path)
  } else {
    if (!is.character(path)){ stop("Path should be of type 'character'!") }
    # get all files present in that directory
    files <- list.files(path)
  }

  if (length(files) == 0){ stop(paste0("Directory ", path, " is empty!")) }

  # read file containing state posterior probabilities of s, t and z
  posterior_s <- .openFile.bexy(path, files, "s_statePosteriors")
  posterior_t <- .openFile.bexy(path, files, "t_statePosteriors")
  posterior_z <- .openFile.bexy(path, files, "z_statePosteriors")

  # trim away s_ and t_ from header
  posterior_s <- .renameHeader.bexy(posterior_s, "s_")
  posterior_t <- .renameHeader.bexy(posterior_t, "t_")

  # calculate posterior mode for s, t and z
  posterior_mode_s <- apply(posterior_s, 2, which.max) - 1
  posterior_mode_t <- apply(posterior_t, 2, which.max) - 1
  posterior_mode_z <- apply(posterior_z, 2, which.max) - 1

  # re-format z as a matrix
  hasZ <- ncol(posterior_z) > 0
  if (hasZ){
    numScaffoldsForZ <- ncol(posterior_z) / length(posterior_s)
    scaffoldNamesZ <- .getScaffoldNames.bexy(colnames(posterior_z), length(posterior_s))
    posterior_z <- matrix(as.numeric(posterior_z[2,]), # just keep P(z = 1)
                          nrow = numScaffoldsForZ,
                          ncol = length(posterior_s),
                          byrow = TRUE,
                          dimnames = list(scaffoldNamesZ, names(posterior_s)))
    posterior_mode_z <- matrix(as.numeric(posterior_mode_z),
                               nrow = numScaffoldsForZ,
                               ncol = length(posterior_s),
                               byrow = TRUE,
                               dimnames = list(scaffoldNamesZ, names(posterior_s)))
  }

  # read file containing posterior mean and variance
  meanVar <- .openFile.bexy(path, files, "meanVar.txt")
  posterior_mean <- meanVar[1,]
  posterior_var <- meanVar[2,]

  # read counts
  counts <- .openFile.bexy(path, files, c("counts", "filtered.txt"))

  # read MCMC trace (if needed)
  trace <- NULL
  if (readMCMCTrace){
    trace <- .openFile.bexy(path, files, "trace.txt")
  }

  result <- list(posterior_s = posterior_s,
                 posterior_t = posterior_t,
                 hasZ = hasZ,
                 posterior_z = posterior_z,
                 posterior_mode_s = posterior_mode_s,
                 posterior_mode_t = posterior_mode_t,
                 posterior_mode_z = posterior_mode_z,
                 posterior_mean = posterior_mean,
                 posterior_var = posterior_var,
                 counts = counts,
                 trace = trace,
                 scaffoldNames = names(posterior_t))
  class(result) <- "bexy"

  return(result)
}

#-------------------------------------------------------------------------------
# Methods for printing
#-------------------------------------------------------------------------------

#' Printing a bexy object
#'
#' @param x A bexy object.
#' @param ... Additional parameters passed to print functions.
#' @return No return value, called for side effects.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' print(bex)
print.bexy <- function(x, ...){
  cat("BeXY posterior mode estimates:\n");
  # print scaffold types
  cat(" - scaffolds: ", sum(x$posterior_mode_t == 0), " autosomes, ",
                        sum(x$posterior_mode_t == 1), " Y-linked scaffolds, ",
                        sum(x$posterior_mode_t == 2), " X-linked scaffolds, ",
                        sum(x$posterior_mode_t == 3), " different scaffolds. ",
                        "\n", sep="");

  # print karyotypes
  string_karyo <- c()
  for (s in 1:length(.getAllPossibleSexKaryotypes.bexy())){
    string_karyo <- c(string_karyo, paste(sum(x$posterior_mode_s == (s-1)), .getAllPossibleSexKaryotypes.bexy()[s]))
  }
  string_karyo <- paste0(string_karyo, collapse = ", ")
  cat(" - sex karyotypes: ", string_karyo, ".\n", sep="");

  # print autosomal trisomies
  if (x$hasZ){
    for (c in 1:nrow(x$posterior_mode_z)){
      numSamplesWithTrisomy <- sum(x$posterior_mode_z[c,])
      if (numSamplesWithTrisomy > 0){
        cat(" - autosomal trisomy of scaffold ", row.names(x$posterior_mode_z)[c], ": ", numSamplesWithTrisomy, " individuals.\n", sep="");
      }
    }
  }

  invisible(x)
}

#' Summarizing a bexy object
#'
#' @param object A bexy object.
#' @param ... Additional parameters passed to summary functions.
#' @return No return value, called for side effects.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' print(bex)
summary.bexy <- function(object, ...){
  print.bexy(object, ...);
}

#-------------------------------------------------------------------------------
# Methods for parsing results
#-------------------------------------------------------------------------------

#' Getting the posterior mode for each sex karyotype
#'
#' @param object A bexy object.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All samples that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return A character vector containing the sex karyotype classification for each sample
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getPosteriorModeSexKaryotypes(bex)
getPosteriorModeSexKaryotypes <- function(object, threshold_certainty = 0.9){
  karyo <- .getAllPossibleSexKaryotypes.bexy()[object$posterior_mode_s + 1]
  certain <- apply(object$posterior_s, 2, function(x) any(x > threshold_certainty))
  karyo[!certain] <- "uncertain"
  names(karyo) <- names(object$posterior_mode_s)
  return(karyo)
}

#' Getting the posterior mode for each scaffold type
#'
#' @param object A bexy object.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All scaffolds that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return A character vector containing the scaffold type classification for each scaffold
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getPosteriorModeScaffoldTypes(bex)
getPosteriorModeScaffoldTypes <- function(object, threshold_certainty = 0.9){
  type <- .getAllPossibleScaffoldTypes.bexy()[object$posterior_mode_t + 1]
  certain <- apply(object$posterior_t, 2, function(x) any(x > threshold_certainty))
  type[!certain] <- "uncertain"
  names(type) <- names(object$posterior_mode_t)
  return(type)
}

#' Getting the posterior mode for each autosomal trisomy
#'
#' @param object A bexy object.
#' @param scaffoldName The scaffold name, must be an autosome.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All samples that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return A character vector containing the trisomy classification for each sample for the given scaffold
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getPosteriorModeAutosomalTrisomy(bex, "21")
getPosteriorModeAutosomalTrisomy <- function(object, scaffoldName, threshold_certainty = 0.9){
  ix <- .getIndexScaffoldAutosomalTrisomy.bexy(object, scaffoldName)
  classification <- .getAllPossibleAutosomalTrisomy.bexy(scaffoldName)[object$posterior_mode_z[ix$ix_z,] + 1]
  certain <- sapply(object$posterior_z[ix$ix_z,], function(x) return(x > threshold_certainty | x < (1 - threshold_certainty)))
  classification[!certain] <- "uncertain"
  names(classification) <- colnames(object$posterior_mode_z)
  return(classification)
}

#' Write a file with the posterior mode for each sex karyotype
#'
#' @param object A bexy object.
#' @param file The name of the output file.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All samples that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return No return value, called for side effects.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' tf <- tempfile(fileext = ".txt")
#' writePosteriorModeSexKaryotypes(bex, tf)
writePosteriorModeSexKaryotypes <- function(object, file, threshold_certainty = 0.9){
  karyo <- getPosteriorModeSexKaryotypes(object, threshold_certainty = threshold_certainty)
  karyo <- cbind(names(karyo), as.character(karyo))

  write.table(x = karyo, file = file, append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE, col.names = c("sample", "sex_karyotype"))
}

#' Write a file with the posterior mode for each scaffold type
#'
#' @param object A bexy object.
#' @param file The name of the output file.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All scaffolds that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return No return value, called for side effects.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' tf <- tempfile(fileext = ".txt")
#' writePosteriorModeScaffoldTypes(bex, tf)
writePosteriorModeScaffoldTypes <- function(object, file, threshold_certainty = 0.9){
  types <- getPosteriorModeScaffoldTypes(object, threshold_certainty = threshold_certainty)
  types <- cbind(names(types), as.character(types))

  write.table(x = types, file = file, append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE, col.names = c("scaffold", "type"))
}

#' Write a file with the posterior mode for each autosomal trisomy
#'
#' @param object A bexy object.
#' @param scaffoldName The scaffold name, must be an autosome.
#' @param file The name of the output file.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All samples that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return No return value, called for side effects.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' tf <- tempfile(fileext = ".txt")
#' writePosteriorModeAutosomalTrisomies(bex, "21", tf)
writePosteriorModeAutosomalTrisomies <- function(object, scaffoldName, file, threshold_certainty = 0.9){
  types <- getPosteriorModeAutosomalTrisomy(object, scaffoldName, threshold_certainty = threshold_certainty)
  types <- cbind(names(types), as.character(types))

  write.table(x = types, file = file, append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE, col.names = c("sample", paste0("autosomal_ploidy_", scaffoldName)))
}

#' Getting all samples classified as a specific sex karyotype
#'
#' @param object A bexy object.
#' @param karyotype One specific sex karyotype, can be XY, XX, X0, XXY, XYY, XXX or XXYY.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All samples that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return A character vector containing all sample names that are classified as the sex karyotype given by argument 'karyotype'.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getSamplesWithThisKaryotype(bex, "XX")
getSamplesWithThisKaryotype <- function(object, karyotype, threshold_certainty = 0.9){
  if (!(karyotype %in% .getAllPossibleSexKaryotypes.bexy())){
    stop(paste0("Unknown karyotype ", karyotype, "!"))
  }
  karyo <- getPosteriorModeSexKaryotypes(object, threshold_certainty = threshold_certainty)
  these <- karyo[karyo == karyotype]
  return(names(these))
}

#' Getting all scaffolds classified as a specific scaffold type
#'
#' @param object A bexy object.
#' @param type One specific scaffold type, can be Autosome, Y-linked, X-linked or Different.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All scaffolds that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return A character vector containing all scaffold names that are classified as the scaffold type given by argument 'type'.
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getScaffoldsWithThisType(bex, "X-linked")
getScaffoldsWithThisType <- function(object, type, threshold_certainty = 0.9){
  if (!(type %in% .getAllPossibleScaffoldTypes.bexy())){
    stop(paste0("Unknown type ", type, "! Possible types are: ", paste0(.getAllPossibleScaffoldTypes.bexy(), collapse = ", "), "."))
  }
  types <- getPosteriorModeScaffoldTypes(object, threshold_certainty = threshold_certainty)
  these <- types[types == type]
  return(names(these))
}

#' Getting all samples classified as a certain autosomal trisomy
#'
#' @param object A bexy object.
#' @param scaffoldName The scaffold name, must be an autosome.
#' @param threshold_certainty The threshold for certainty on the posterior probabilities. All samples that have a posterior probability less than this threshold are classified as 'uncertain'.
#' @return A character vector containing all sample names that are classified as having an autosomal trisomy of scaffold 'scaffoldName'
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getSamplesWithAutosomalTrisomy(bex, "21")
getSamplesWithAutosomalTrisomy <- function(object, scaffoldName, threshold_certainty = 0.9){
  classification <- getPosteriorModeAutosomalTrisomy(object, scaffoldName, threshold_certainty = threshold_certainty)
  these <- classification[classification != "diploid" & classification != "uncertain"]
  return(names(these))
}

#' Getting posterior mean rho (ploidy ratio parameter) for each scaffold
#'
#' @param object A bexy object.
#' @return A numeric vector containing the posterior mean of rho for each scaffold
#'
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' getPosteriorMeanRho(bex)
getPosteriorMeanRho <- function(object){
  rho <- .getColumns.bexy(object$posterior_mean, "rho")
  names(rho) <- names(object$posterior_mode_t)
  return(rho)
}

#-------------------------------------------------------------------------------
# Methods for plotting
#-------------------------------------------------------------------------------

#' Plotting a bexy object. To customize the plots, please call the respective plotting functions directly (e.g. plotTernary, plotCounts, plotBars...).
#'
#' @param x A bexy object.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plot(bex)
plot.bexy <- function(x, ...){
  plotTernary(x, ...)
  plotCounts(x, ...)
  plotBars(x, ...)
}

#' Plotting the ternary of a bexy object
#'
#' @param x A bexy object.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param addSquares Boolean indicating whether 7-cell posterior probability square plots should be drawn.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plotTernary(bex)
plotTernary <- function(x, colors = c("turquoise3", "darkorange", "royalblue4"), sampleNamesToHighlight = c(), label = TRUE, addSquares = TRUE, ...){
  if (length(colors) != 3){
    stop("Invalid colors! Should be of length 3 (aneuploid, XX and XY).")
  }
  .plotTernaryAndTicTacToe.bexy(x, colors, sampleNamesToHighlight, label, addSquares, ...)
}

#' Plotting the percentage of counts mapping to two scaffolds (by default a Y-linked scaffold against a X-linked scaffold), colored by sex karyotype
#'
#' @param x A bexy object.
#' @param scaffoldX A string with the name of the scaffold to be plotted on the x-axis. If NULL (default), a X-linked scaffold is used.
#' @param scaffoldY A string with the name of the scaffold to be plotted on the y-axis. If NULL (default), a Y-linked scaffold is used.
#' @param colors A vector of length three with the colors for aneuploid, XX and XY samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plotCounts(bex)
plotCounts <- function(x, scaffoldX = NULL, scaffoldY = NULL, colors = c("turquoise3", "darkorange", "royalblue4"), sampleNamesToHighlight = c(), label = TRUE, ...){
  if (length(colors) != 3){
    stop("Invalid colors! Should be of length 3 (aneuploid, XX and XY).")
  }

  # check if scaffoldX and scaffoldY exist
  ix_X <- .getIndexScaffold.bexy(x, scaffoldX, "X-linked")
  ix_Y <- .getIndexScaffold.bexy(x, scaffoldY, "Y-linked")
  if (is.null(ix_X) & is.null(ix_Y)){
    # neither X- or Y-linked scaffolds are found -> error
    stop("BeXY did not detect any X- or Y-linked scaffolds and it can thus not plot X vs Y.")
  }

  .plotXVsY.bexy(x, ix_X, ix_Y, colors, sampleNamesToHighlight, label, ...)
}

#' Plotting a barplot with the posterior probabilities of each sex karyotype per sample
#'
#' @param x A bexy object.
#' @param colors A vector of length 7 with the colors for XY, XX, X0, XXY, XYY, XXX and XXYY, respectively.
#' @param maxNumSamplesPerPlot How many samples to place next to each other into one barplot. If there are more samples, multiple barplots will be generated. This is to ensure that the labels are readable.
#' @param sortByKaryotype If TRUE, the samples (i.e. the bars) are ordered by the posterior mode of the sex karyotype.
#' @param sampleNames A vector of sample names that represents the order in which the bars are plotted. Overrides any argument given to sortByKaryotype.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plotBars(bex)
plotBars <- function(x, colors = c("lightblue", "lightpink", "khaki", "maroon1", "dodgerblue3", "red4", "navy"), maxNumSamplesPerPlot = 100, sortByKaryotype = TRUE, sampleNames = c(), ...){
  if (length(colors) != 7){
    stop("Invalid colors! Should be of length 7 (XY, XX, X0, XXY, XYY, XXX and XXYY).")
  }

  # sort samples by posterior mode karyotype
  s <- x$posterior_s

  if (length(sampleNames) > 0){ # match sample names
    matched <- match(sampleNames, names(s))
    if (any(is.na(matched))){
      stop("Sample(s) '", paste(sampleNames[which(is.na(matched))], collapse = ", "), "' could not be found!")
    }
    s <- s[,matched]
  } else if (sortByKaryotype){
    max_s <- apply(s, 2, which.max)
    sorted_max_s <- sort(max_s, index.return = T)
    s <- s[,sorted_max_s$ix]
  }


  # split into multiple plots, if necessary
  numPlots <- ceiling(ncol(s) / maxNumSamplesPerPlot)
  for (i in 1:numPlots){
    # get start and end of current plot
    start <- (i-1) * maxNumSamplesPerPlot + 1
    end <- min(ncol(s), start + maxNumSamplesPerPlot - 1)
    m <- s[,start:end]

    x <- barplot(as.matrix(m),
                 ylab = "Posterior probability",
                 las = 2,
                 col = colors,
                 space = 0,
                 border = "darkgray",
                 xaxt = "n",
                 ...)

    # add sample names (x axis)
    axis(side = 1,
         at = x,
         labels = names(m),
         tick = F,
         las = 2,
         line = -1,
         cex.axis = 0.5,
         ...)

    # add legend
    legend("top",
           inset=c(0, -0.1),
           legend = c("XY", "XX", "X0", "XXY", "XYY", "XXX", "XXYY"),
           fill = colors,
           border = NA,
           horiz = T,
           xpd = T,
           bty = "n",
           ...)
  }
}

#' Plotting the autosomal trisomies
#'
#' @param x A bexy object.
#' @param scaffoldName The name of the scaffold to plot (has to be an autosome)
#' @param colors A vector of length two with the colors for diploid and triploid samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plotAutosomalTrisomy(bex, "21")
plotAutosomalTrisomy <- function(x, scaffoldName, colors = c("deepskyblue4", "darkred"), sampleNamesToHighlight = c(), label = TRUE, ...){
  plotAutosomalTrisomyPosteriorProbabilities(x, scaffoldName, colors, sampleNamesToHighlight, label, ...)
  plotAutosomalTrisomyCounts(x, scaffoldName, colors, sampleNamesToHighlight, label, ...)
}

#' Plotting the autosomal trisomies: posterior probabilities
#'
#' @param x A bexy object.
#' @param scaffoldName The name of the scaffold to plot (has to be an autosome)
#' @param colors A vector of length two with the colors for diploid and triploid samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plotAutosomalTrisomyPosteriorProbabilities(bex, "21")
plotAutosomalTrisomyPosteriorProbabilities <- function(x, scaffoldName, colors = c("deepskyblue4", "darkred"), sampleNamesToHighlight = c(), label = TRUE, ...){
  if (length(colors) != 2){
    stop("Invalid colors! Should be of length 2 (diploid and triploid).")
  }
  # plot posterior probabilities
  .plotAutosomalTrisomy.bexy(object = x,
                             scaffoldName = scaffoldName,
                             plotCounts = FALSE,
                             colors = colors,
                             sampleNamesToHighlight = sampleNamesToHighlight,
                             label = label,
                             ...)
}

#' Plotting the autosomal trisomies: counts on that scaffold
#'
#' @param x A bexy object.
#' @param scaffoldName The name of the scaffold to plot (has to be an autosome)
#' @param colors A vector of length two with the colors for diploid and triploid samples, respectively.
#' @param sampleNamesToHighlight A vector of sample names that should be highlighted. If empty, samples are automatically highlighted based on the posterior probabilities.
#' @param label Boolean indicating whether samples should be labeled.
#' @param ... Other parameters used for plotting.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy()
#' plotAutosomalTrisomyCounts(bex, "21")
plotAutosomalTrisomyCounts <- function(x, scaffoldName, colors = c("deepskyblue4", "darkred"), sampleNamesToHighlight = c(), label = TRUE, ...){
  if (length(colors) != 2){
    stop("Invalid colors! Should be of length 2 (diploid and triploid).")
  }
  # plot counts
  .plotAutosomalTrisomy.bexy(object = x,
                             scaffoldName = scaffoldName,
                             plotCounts = TRUE,
                             colors = colors,
                             sampleNamesToHighlight = sampleNamesToHighlight,
                             label = label,
                             ...)
}

#' Plot the MCMC trace for a specific parameter
#'
#' @param object A bexy object.
#' @param parameterName The name of the parameter to plot (e.g. rho, logSigma, epsilon, s, t, a, f, ...)
#' @param maxNumPlots The maximum number of plots to plot, default 20.
#' @return No return value, called for side effects.
#' @export
#' @seealso \code{\link{bexy}}
#' @examples
#' bex <- bexy(readMCMCTrace = TRUE)
#' plotMCMCTrace(bex, "rho", maxNumPlots = 1)
plotMCMCTrace <- function(object, parameterName, maxNumPlots = 20){
  if (is.null(object$trace)){
    stop("Trace was not loaded in bexy object! Set 'readMCMCTrace = TRUE' when calling bexy().")
  }
  tr <- .getColumns.bexy(object$trace, parameterName)
  if (is.numeric(tr)){ # single parameter
    plot(tr, xlab = "MCMC iteration (thinned)", ylab = parameterName, main = parameterName)
  } else {
    numPlots <- min(ncol(tr), maxNumPlots)
    for (c in 1:numPlots){ # multiple parameters
      plot(tr[,c], xlab = "MCMC iteration (thinned)", ylab = names(tr)[c], main = names(tr)[c])
    }
  }
}