R/plot.neighbourLD.r

In synbreed: Framework for the Analysis of Genomic Prediction Data using R

#' Plot neighbour linkage disequilibrium
#'
#' A function to visualize Linkage Disequilibrium estimates between adjacent
#' markers.
#'
#' The plot is similar to \code{plotGenMap} with the option \code{dense=TRUE},
#' but here the LD between adjacent markers is plotted along the chromosomes.
#'
#' @param LD object of class \code{LDmat}, i.e the output of function
#' \code{pairwiseLD} using argument \code{type="matrix"}.
#' @param gpData object of class \code{gpData} with object \code{map} or a
#' \code{data.frame} with columns 'chr' (specifying the chromosome of the
#' marker) and 'pos' (position of the marker within chromosome measured with
#' genetic or physical distances)
#' @param dense \code{logical}. Should density visualization for high-density
#' genetic maps be used?
#' @param nMarker \code{logical}. Print number of markers for each chromosome?
#' @param centr \code{numeric} vector. Positions for the centromeres in the
#' same order as chromosomes in \code{map}. If \code{"maize"}, centromere
#' positions of maize in Mbp are used.
#' @param file Optionally a path to a file where the plot is saved to
#' @param fileFormat \code{character}. At the moment two file formats are
#' supported: pdf and png. Default is \code{"pdf"}.
#' @param \dots further graphical arguments for function \code{plot}
#' @return Plot of neighbour LD along each chromosome. One chromosome is
#' displayed from the first to the last marker.
#' @author Theresa Albrecht and Hans-Juergen Auinger
#' @seealso \code{\link{plotGenMap}}, \code{\link{pairwiseLD}}
#' @keywords hplot
#' @examples
#'
#' \dontrun{
#' library(synbreedData)
#' data(maize)
#' maize2 <- codeGeno(maize)
#' LD <- pairwiseLD(maize2, chr = 1:10, type = "matrix")
#' plotNeighbourLD(LD, maize2, nMarker = FALSE)
#' }
#'
#' @export plotNeighbourLD
#' @importFrom grDevices dev.off pdf png
#' @importFrom methods is
#' @importFrom stats smooth
#' @importFrom graphics axis box image layout legend par plot points text polygon
#'
plotNeighbourLD <- function(LD, gpData, dense = FALSE, nMarker = TRUE, centr = NULL, file = NULL, fileFormat = "pdf", ...) {
  oldPar <- par()
  if (class(gpData) == "gpData") {
    gpData.unit <- gpData$info$map.unit
    map <- gpData$map
    class(map) <- "data.frame"
  }
  else {
    map.unit <- "unit"
  }
  chr <- unique(map$chr)
  chr <- chr[!is.na(chr)]
  map <- map[!is.na(map$chr), ]
  if (class(map$chr) != "factor") bord <- NULL else if (par()$bg == "transparent") bord <- "white" else bord <- "transparent"

  # centromere positions of maize
  if (!is.null(centr)) if (centr == "maize") centr <- c(133, 90, 95, 104.6, 105.5, 50, 55.3, 46.5, 68.8, 59.9)

  # norm pos
  if (!is.null(centr)) map$pos <- map$pos - centr[map$chr]

  # output in files
  if (!is.null(file)) {
    if (substr(file, nchar(file) - nchar(fileFormat) + 1, nchar(file)) != fileFormat | nchar(file) < 5) {
      file <- paste(file, ".", fileFormat, sep = "")
    }
    if (fileFormat == "pdf") {
      pdf(file)
    } else if (fileFormat == "png") {
      png(file)
    } else {
      stop("not supported file format choosen!")
    }
  }

  # initialize map
  layout(matrix(1:2, ncol = 2), widths = c(0.82, 0.18))

  # make an empty plot
  if (!is.null(centr)) {
    plot(map,
      type = "n", xaxt = "n", xlim = c(0.5, length(chr) + 0.5), border = bord,
      ylim = c(max(map$pos, na.rm = TRUE) * 1.1, min(map$pos, na.rm = TRUE)), axes = FALSE, ...
    )
  }
  else {
    plot(map,
      type = "n", xaxt = "n", xlim = c(0.5, length(chr) + 0.5), border = bord,
      ylim = c(max(map$pos, na.rm = TRUE) * 1.1, min(map$pos, na.rm = TRUE)), ...
    )
  }

  # x-axis
  axis(side = 1, at = seq(along = chr), labels = chr)
  # y-axis
  if (!is.null(centr)) {
    box()
    axis(side = 2, at = -seq(-round(max(map$pos, na.rm = TRUE), -2), round(max(map$pos, na.rm = TRUE), -2), by = 25), labels = abs(-seq(-round(max(map$pos, na.rm = TRUE), -2), round(max(map$pos, na.rm = TRUE), -2), by = 25)), las = 1)
  }

  # loop over chromosomes
  for (i in seq(along = chr)) {
    mNam <- rownames(map[map$chr == chr[i], ])
    LDm <- matrix(NA, ncol = 10, nrow = length(mNam))
    if (class(LD) == "LDdf") {
      for (j in 1:10) {
        LDm[1:(nrow(LDm) - j), j] <- LD[[i]]$r2[paste(LD[[i]]$marker1, LD[[i]]$marker2, sep = "") %in% paste(mNam[j:(length(mNam) - 1)], mNam[(1 + j):length(mNam)], sep = "")]
      }
    } else if (class(LD) == "LDmat") {
      for (j in 1:10) {
        LDm[1:(nrow(LDm) - j), j] <- diag(LD$LD[[i]][j:(nrow(LDm) - 1), (1 + j):nrow(LDm)])
      }
    } else {
      stop(paste("Wrong class of object", substitute(LD)))
    }

    n <- sum(map$chr == i, na.rm = TRUE)
    start <- min(map$pos[map$chr == chr[i]], na.rm = TRUE)
    end <- max(map$pos[map$chr == chr[i]], na.rm = TRUE)

    # display.brewer.pal(7, "Reds")s
    cols <- c("#FCBBA1", "#FC9272", "#FB6A4A", "#EF3B2C", "#CB181D", "#99000D")

    if (dense) { # calculating averaged LD
      # map positions
      mapPos <- smooth(map$pos[map$chr == chr[i]])
      # matrices with 10 rows
      avLD <- rowMeans(LDm[1:(nrow(LDm) - 11), ])
      avPos <- matrix(NA, ncol = 10, nrow = nrow(map[map$chr == i, ]))
      for (j in 1:10) {
        avPos[1:(nrow(LDm) - j), j] <- mapPos[1:(length(mapPos) - j)]
      }
      avPos <- rowMeans(avPos[1:(nrow(avPos) - 11), ])
    } # end of if,  smooth LD calculation
    else { # using LD directly
      # map posittions
      mapPos <- map$pos[map$chr == chr[i]]
      avLD <- LDm[-nrow(LDm), 1]
      avPos <- mapPos[-length(mapPos)]
    }
    LDPos <- data.frame(avLD, avPos)
    LDPos <- LDPos[!duplicated(LDPos$avPos), ]
    # LDPos <- LDPos[LDPos$avLD>=0.2, ]
    # visualisation of map
    image(seq(i - 0.35, i + 0.35, length = 20), LDPos$avPos, matrix(rep(LDPos$avLD, 20), nrow = 20, byrow = TRUE), col = cols, add = TRUE, zlim = c(0, 1))
    if (!is.null(centr)) {
      # centromere
      polygon(x = c(i - 0.4, i - 0.1, i - 0.1, i - 0.4, i - 0.4), y = c(-10, -1, 1, 10, -10), col = "white", border = "white")
      polygon(x = c(i + 0.4, i + 0.1, i + 0.1, i + 0.4, i + 0.4), y = c(-10, -1, 1, 10, -10), col = "white", border = "white")
    }
    if (nMarker) text(i, max(map$pos) * 1.05, sum(map$chr == chr[i], na.rm = TRUE))
  } # end chromosome loop

  # add legend
  par(mar = c(5, 1, 4, 3.8) + 0.1)
  image(seq(-0.4, 0.4, length = 20), seq(from = 0, to = 1, length = 6), matrix(rep(seq(from = 0, to = 1, length = 6), 20), nrow = 20, byrow = TRUE), main = expression(paste(r^{
    2
  })), col = cols, axes = FALSE, xlab = "")
  axis(side = 4, at = round(seq(from = 0, to = 1, length = 6), 4), las = 1)
  # close graphic device
  if (!is.null(file)) {
    dev.off()
  } else {
    oldPar$cin <- oldPar$cra <- oldPar$csi <- oldPar$cxy <- oldPar$din <- NULL
    par(oldPar)
  }
}