R/plot.haplen.R

In rehh: Searching for Footprints of Selection using 'Extended Haplotype Homozygosity' Based Tests

#'Plot the length of extended haplotypes around a focal marker
#'@description Plot the length of extended haplotype around a focal marker.
#'@param x an object of class \code{haplen} generated by \code{\link{calc_haplen}}.
#'@param allele if \code{NA} (default), haplotypes of all alleles are plotted,
#'otherwise for the specified alleles. Alleles must be specified by their
#'internal coding, i.e. '0' for ancestral resp. major allele, etc.
#'@param group_by_allele logical. If \code{TRUE} (default), group chromosomes
#'by their allele at the focal marker; alleles are ordered by their internal coding
#'unless parameter \code{alleles} is specified. If \code{FALSE}, haplotypes are drawn by
#'their order in the input file.
#'@param order_by_length if \code{TRUE}, chromosomes are ordered by their
#'shared haplotype length.
#'@param col color for each allele (as coded internally).
#'@param mrk.col color of the vertical line at the focal marker position.
#'@param lwd line width.
#'@param hap.names a vector containing the names of chromosomes.
#'@param cex.lab relative letter size of labels. See \code{\link[graphics]{par}}.
#'@param family.lab font family for labels. See \code{\link[graphics]{par}}.
#'@param offset.lab offset of labels. See \code{\link[graphics]{par}}.
#'@param pos.lab position of haplotype labels. Either \code{"left"} , \code{"right"} or \code{"both"}.
#'@param legend legend text.
#'@param legend.xy.coords if \code{"automatic"} (default) places legend either top left or top right;
#'if \code{"none"}, no legend is drawn; otherwise argument is passed to \code{\link[graphics]{legend}}.
#'@param ... other parameters to be passed to \code{\link[graphics]{plot.default}}.
#'@seealso \code{\link{calc_haplen}}, \code{\link{plot.furcation}}.
#'@examples #example haplohh object (280 haplotypes, 1424 SNPs)
#'#see ?haplohh_cgu_bta12 for details
#'data(haplohh_cgu_bta12)
#'#plotting length of extended haplotypes for both ancestral and derived allele
#'#of the marker "F1205400"
#'#which displays a strong signal of selection
#'f <- calc_furcation(haplohh_cgu_bta12, mrk = "F1205400")
#'h <- calc_haplen(f)
#'plot(h)
#'plot(h, hap.names = hap.names(haplohh_cgu_bta12), cex.lab = 0.3)
#'@export
#'@importFrom graphics lines text
plot.haplen <-
  function(x,
           allele = NA,
           group_by_allele = TRUE,
           order_by_length = FALSE,
           col = c("blue", "red", "violet", "orange"),
           mrk.col = "gray",
           lwd = 1,
           hap.names = NULL,
           cex.lab = 1.0,
           family.lab = "",
           offset.lab = 0.5,
           pos.lab = "left",
           legend = NA,
           legend.xy.coords = "automatic",
           ...) {
    ##check parameters
    if (!is.haplen(x)) {
      stop("The data is not a valid object of a haplen", call. = FALSE)
    }
    
    if (!is.null(hap.names)) {
      if (length(hap.names) != nrow(x$haplen)) {
        stop(
          "Number of specified haplotype names has to be equal to number of haplotypes.",
          call. = FALSE
        )
      }
    }
    
    if (!anyNA(allele) & !is.numeric(allele)) {
      stop("alleles have to be specified by integers 0,1,2,...", call. = FALSE)
    }
    
    dot_args <- list(...)
    
    if (!is.null(dot_args$xlim)) {
      if (!is.numeric(dot_args$xlim) |
          length(dot_args$xlim) != 2 |
          dot_args$xlim[2] < dot_args$xlim[1]) {
        stop("Incorrect specification of xlim.", call. = FALSE)
      }
      if (!((dot_args$xlim)[1] <= x$position &
            x$position <= (dot_args$xlim)[2])) {
        stop("Focal marker must lie in specified xlim.", call. = FALSE)
      }
    }
    
    if (!anyNA(allele)) {
      if (anyDuplicated(allele)) {
        stop("Repeated specification of the same allele.", call. = FALSE)
      }
      for (i in allele) {
        if (!(i %in% unique(x$haplen$ALLELE))) {
          stop(paste0("No allele ", i, " found."), call. = FALSE)
        }
      }
    } else{
      allele <- sort(unique(x$haplen$ALLELE))
    }
    
    if (!anyNA(legend.xy.coords)) {
      if (is.numeric(legend.xy.coords) & length(legend.xy.coords) == 2) {
        legend.x <- legend.xy.coords[1]
        legend.y <- legend.xy.coords[2]
      } else if (length(legend.xy.coords) == 1) {
        legend.x <- legend.xy.coords
        legend.y <- NULL
      } else{
        stop("Invalid legend coordinate specification.", call. = FALSE)
      }
    } else{
      legend.x <- "none"
    }
    
    # perform plot
    
    #remove NAs and non-specified alleles
    if (!is.null(hap.names)) {
      hap.names <- hap.names[x$haplen$ALLELE %in% allele]
    }
    haplen_subset <- x$haplen[x$haplen$ALLELE %in% allele, ]
    
    if (order_by_length) {
      ord <- order(haplen_subset$MAX - haplen_subset$MIN)
      haplen_subset <- haplen_subset[ord, ]
      hap.names <- hap.names[ord]
    }
    
    if (group_by_allele) {
      #sort in order of specified alleles ("radix" -> stable sort)
      ord <- order(vapply(haplen_subset$ALLELE, function(x) {
        which(x == allele)
      }, FUN.VALUE = 0L),
      method = "radix")
      haplen_subset <- haplen_subset[ord, ]
      hap.names <- hap.names[ord]
    }
    
    haplen_subset$MIN[is.na(haplen_subset$MIN)] <- x$position
    haplen_subset$MAX[is.na(haplen_subset$MAX)] <- x$position
    
    if (is.null(dot_args$xlim)) {
      dot_args$xlim <- c(min(haplen_subset$MIN),
                         max(haplen_subset$MAX))
    } else{
      haplen_subset$MIN <- pmax(haplen_subset$MIN, dot_args$xlim[1])
      haplen_subset$MAX <- pmin(haplen_subset$MAX, dot_args$xlim[2])
    }
    
    p <- floor(log(dot_args$xlim[2], 1000))
    ## only shrink big scales, but never magnify small ones (p<0)
    scale <- 1000 ** max(0, p)
    ## no unit if p < 0
    unit <- c("", "(bp)", "(kb)", "(Mb)", "(Gb)")[max(-1, p)  + 2]
    
    dot_args$xlim <- dot_args$xlim / scale
    
    dot_args$ylim <- c(0, 1)
    
    if (is.null(dot_args$xlab)) {
      dot_args$xlab <- paste("Position", unit)
    }
    
    if (is.null(dot_args$ylab)) {
      dot_args$ylab <- ""
    }
    
    if (is.null(dot_args$bty)) {
      dot_args$bty <- "n"
    }
    
    if (is.null(dot_args$yaxt)) {
      dot_args$yaxt <- "n"
    }
    
    if (is.null(dot_args$main)) {
      dot_args$main <-
        paste0("Haplotype length around '", x$mrk.name, "'")
    }
    
    ##perform plot
    description_colors <-
      get_description_colors(haplen_subset$DESCRIPTION, col)
    
    #invisible plot to get coordinate system
    do.call("plot", c(list(NULL),
                      dot_args))
    
    #dashed vertical line at focal mrk
    abline(v = x$position / scale,
           lty = 2,
           col = mrk.col)
    
    for (i in seq_len(nrow(haplen_subset))) {
      #draw haplo-lines
      lines(
        c(haplen_subset$MIN[i], haplen_subset$MAX[i]) / scale,
        rep(1 - (i - 0.5) / nrow(haplen_subset), 2),
        col = description_colors[i],
        lwd = lwd
      )
      
      #add labels
      if (!is.null(hap.names)) {
        if (pos.lab == "left" || pos.lab == "both") {
          text(
            x = dot_args$xlim[1],
            y = rep(1 - (i - 0.5) / nrow(haplen_subset), 2),
            labels = hap.names[i],
            cex = cex.lab,
            family = family.lab,
            offset = offset.lab,
            pos = 2,
            xpd = TRUE
          )
        }
        if (pos.lab == "right" || pos.lab == "both") {
          text(
            x = dot_args$xlim[2],
            y = rep(1 - (i - 0.5) / nrow(haplen_subset), 2),
            labels = hap.names[i],
            cex = cex.lab,
            family = family.lab,
            offset = offset.lab,
            pos = 4,
            xpd = TRUE
          )
        }
      }
    }
    if (legend.x != "none") {
      if (legend.x == "automatic") {
        picturemiddle <- mean(dot_args$xlim)
        legend.x <-
          ifelse(x$position / scale > picturemiddle, "topleft", "topright")
      } else if (is.numeric(legend.x)) {
        legend.x <- legend.x / scale
      }
      
      descriptions <- unique(haplen_subset$DESCRIPTION)
      
      if (is.na(legend)) {
        legend <- descriptions
      }
      
      legend(
        legend.x,
        legend.y,
        legend = legend,
        bty = dot_args$bty,
        col = get_description_colors(descriptions, col),
        lwd = lwd,
        xpd = TRUE
      )
    }
  }