R/plotAberration.r

In igordot/copynumber: Segmentation of single- and multi-track copy number data by penalized least squares regression (with hg38 and mm10).

####################################################################
## Author: Gro Nilsen, Knut Liestøl and Ole Christian Lingjærde.
## Maintainer: Gro Nilsen <gronilse@ifi.uio.no>
## License: Artistic 2.0
## Part of the copynumber package
## Reference: Nilsen and Liestøl et al. (2012), BMC Genomics
####################################################################


# Function that plots aberrations given limits - by genome og chromosomes

## Input:
### segments: segmentation results from pcf or multipcf
### upper.lim, lower.lim: a vector with limit(s) to be applied for abberration calling
### pos.unit: unit used for positions (bp,kbp,mbp)
### chrom: a vector with chromosomes to be plotted. If specified the frequencies are plotted with one panel for each chromosome
### layout: number of columns and rows in plot
### ... : other optional plot parameters


## Required by: none

## Requires:
### numericChrom
### checkChrom
### genomeHeat
### chromosomeHeat
### is.multiseg
### checkSegments
### pullOutContent

# Main function for heatmap plotting:
plotAberration <- function(segments, thres.gain, thres.loss = -thres.gain, pos.unit = "bp", chrom = NULL, layout = c(1, 1), ...) {
  # If chrom is unspecified, the whole genome is plotted. Otherwise, selected chromosomes are plotted
  type <- ifelse(is.null(chrom), "genome", "bychrom")

  # Check input in segments:
  segments <- pullOutContent(res = segments, what = "segments")

  # Check segments (convert multi segments to unisegments etc)
  segments <- checkSegments(segments, type)

  # Check and if necessary modify chrom to be plotted
  chrom <- checkChrom(data = NULL, segments = segments, chrom)

  # Convert segments back to data frame (is returned as list in checkSegments)
  segments <- as.data.frame(segments)

  # Get sampleIds
  sampleID <- unique(segments[, 1])

  # Check pos.unit input:
  if (!pos.unit %in% c("bp", "kbp", "mbp")) {
    stop("pos.unit must be one of bp, kbp and mbp", call. = FALSE)
  }

  # Making sure number of upperlimits and lowerlimits are the same:
  upper.lim <- thres.gain
  lower.lim <- thres.loss
  nT <- min(length(upper.lim), length(lower.lim))
  upper.lim <- upper.lim[1:nT]
  lower.lim <- lower.lim[1:nT]

  # plot heatmap, either over genome or by chromosomes:
  switch(type,
    genome = genomeAberration(segments, upper.lim, lower.lim, pos.unit, sampleID, layout, ...),
    bychrom = chromosomeAberration(segments, upper.lim, lower.lim, pos.unit, sampleID, chrom, layout, ...)
  )
}



# Function that plots aberrations by genome

## Input:
### segments: segmentation data frame
### upper.lim, lower.lim: a vector with limit(s) to be applied for abberration calling
### pos.unit: unit used for positions (bp,kbp,mbp)
### sampleID: ids for the samples to be plotted
### layout: number of columns and rows in plot
### ... : other optional plot parameters


## Required by: plotHeatmap

## Requires:
### getHeatParameters
### getGlobal.xlim
### adjustSegPos
### framedim
### colorSetup
### getCol
### addChromlines
### chromPattern
### addArmlines

genomeAberration <- function(segments, upper.lim, lower.lim, pos.unit, sampleID, layout, ...) {
  nT <- length(upper.lim)
  nr <- layout[1]
  nc <- layout[2]
  rc <- nr * nc
  nSample <- length(sampleID)

  op <- getHeatParameters(type = "genome", nc = nc, nr = nr, nSample = nSample, upper.lim = upper.lim, lower.lim = lower.lim, ab = TRUE, ...)

  # Set global xlimits if not specified by user:
  if (is.null(op$xlim)) {
    op$xlim <- getGlobal.xlim(op = op, pos.unit = pos.unit, chrom = unique(segments[, 2]))
  }


  # Check if there should be more than one file/window with plot(s), and get file.name accordingly
  if (dev.cur() <= 1) { # to make Sweave work
    dev.new(width = op$plot.size[1], height = op$plot.size[2], record = TRUE)
  }

  # Initialize:
  row <- 1
  clm <- 1
  new <- FALSE

  # Division of plotting window:
  frames <- framedim(nr, nc)


  for (t in 1:nT) {
    # Frame dimensions for plot t:
    fig.t <- c(frames$left[clm], frames$right[clm], frames$bot[row], frames$top[row])
    par(fig = fig.t, new = new, oma = c(0, 0, 0.5, 0), mar = op$mar)

    # Empty plot with correct dimensions:
    plot(1, 1,
      type = "n", ylim = c(0, nSample), ylab = op$ylab, xlab = op$xlab, xlim = op$xlim,
      xaxs = "i", yaxt = "n", xaxt = "n", yaxs = "i", cex.lab = 0.9, mgp = op$mgp, main = ""
    )
    # Let background be white:
    rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = "white")
    # Add shifting white/grey pattern to backgroud to separate chromosomes:
    chromPattern(pos.unit, op)


    # main title for this plot
    title(main = op$main[t], line = op$main.line, cex.main = op$cex.main)


    # Plot heatmap for each sample:
    for (i in 1:nSample) {
      sample.segments <- segments[segments[, 1] == sampleID[i], ]

      # Adjust positions to be plotted along xaxis; i.e. get global positions, scale according to plot.unit, and get left and right pos for rectangles
      # to be plotted (either continuous or 1 probe long):
      x <- adjustSegPos(chrom = sample.segments[, 2], char.arms = sample.segments[, 3], start = sample.segments[, 4], stop = sample.segments[, 5], type = "genome", unit = pos.unit, op = op)
      xleft <- x$use.start
      xright <- x$use.stop

      # Find appropriate colour for each probe:
      heat.col <- rep("white", nrow(sample.segments))
      heat.col[sample.segments[, 2] %% 2 == 0] <- "grey95"
      heat.col[sample.segments[, 7] > upper.lim[t]] <- op$colors[2]
      heat.col[sample.segments[, 7] < lower.lim[t]] <- op$colors[1]
      ytop <- i - op$sep.samples
      ybottom <- i - (1 - op$sep.samples)

      # Plot rectangles with appropriate color for each probe:
      rect(xleft, ybottom, xright, ytop, col = heat.col, border = NA)

      # Add sampleid on yaxis
      if (op$sample.labels) {
        axis(side = 2, at = (ytop - (ytop - ybottom) / 2), labels = sampleID[i], line = op$sample.line, tcl = 0, cex.axis = op$sample.cex, las = 1, mgp = op$mgp, tick = FALSE)
      }
    } # endfor


    # Separate chromosomes and arms by vertical lines:
    op$chrom.lty <- 1
    addChromlines(chromosomes = segments[, 2], xaxis = "pos", unit = pos.unit, cex = op$cex.chrom, op = op)
    addArmlines(chromosomes = segments[, 2], xaxis = "pos", unit = pos.unit, cex = op$cex.chrom, op = op)

    # Box:
    abline(v = op$xlim)
    abline(h = c(0, nSample))

    # Get new page, or update column/row:
    if (t %% (nr * nc) == 0) {
      # Start new plot page (prompted by user)
      devAskNewPage(ask = TRUE)

      # Reset columns and row in layout:
      clm <- 1
      row <- 1
      new <- FALSE
    } else {
      # Update column and row index:
      if (clm < nc) {
        clm <- clm + 1
      } else {
        clm <- 1
        row <- row + 1
      } # endif
      new <- TRUE
    } # endif
  } # endfor
}

# Function that plots aberrations by chromosomes (each chrom in separate panel)

## Input:
### segments: segmentation data frame
### upper.lim, lower.lim: a vector with limits(s) to be applied for abberration calling
### pos.unit: unit used for positions (bp,kbp,mbp)
### sampleID: the IDs for samples to be plotted
### chrom: the chromosomes to be plotted
### layout: number of columns and rows in plot
### ... : other optional plot parameters


## Required by: plotHeatmap

## Requires:
### getHeatParameters
### adjustSegPos
### framedim
### plotIdeogram
### chromMax
### get.xticks


chromosomeAberration <- function(segments, upper.lim, lower.lim, pos.unit, sampleID, chrom, layout, ...) {
  nT <- length(upper.lim)
  nr <- layout[1]
  nc <- layout[2]
  nChrom <- length(chrom)
  nSample <- length(sampleID)


  # Default plot options:
  op <- getHeatParameters(type = "bychrom", nc = nc, nr = nr, nSample = nSample, upper.lim = upper.lim, lower.lim = lower.lim, chrom = chrom, ab = TRUE, ...)

  # Margins for entire plot in window:
  if (all(op$title == "")) {
    oma <- c(0, 0, 0, 0)
  } else {
    oma <- c(0, 0, 1, 0)
  }
  mar <- c(0.2, 0.2, 0.3, 0.2)


  # Divide the plotting window by the function "framedim":
  frames <- framedim(nr, nc)


  # make separate plots for each value of limits
  for (t in 1:nT) {
    # Start new window/file:
    if (dev.cur() <= 1) { # to make Sweave work
      dev.new(width = op$plot.size[1], height = op$plot.size[2], record = TRUE)
    }

    # Initialize row and column index:
    row <- 1
    clm <- 1
    new <- FALSE

    # Separate plots for each chromosome:
    for (c in 1:nChrom) {
      # Frame dimensions for plot c:
      fig.c <- c(frames$left[clm], frames$right[clm], frames$bot[row], frames$top[row])
      par(fig = fig.c, new = new, oma = oma, mar = mar)

      # Make list with frame dimensions:
      frame.c <- list(left = frames$left[clm], right = frames$right[clm], bot = frames$bot[row], top = frames$top[row])

      # Select relevant chromosome number
      k <- chrom[c]

      # Pick out indeces for this chromosome
      ind.c <- which(segments[, 2] == k)
      chrom.segments <- segments[ind.c, ]

      # Find maximum position for this chromosome (and scale according to plot unit)
      scale.fac <- convert.unit(unit1 = op$plot.unit, unit2 = pos.unit)
      xlim <- c(0, max(chrom.segments[, 5])) * scale.fac

      # Plot ideogram below heatmaps:
      if (op$plot.ideo) {
        # Ideogram frame:
        ideo.frame <- frame.c
        ideo.frame$top <- frame.c$bot + (frame.c$top - frame.c$bot) * op$ideo.frac

        par(fig = unlist(ideo.frame), new = new, mar = op$mar.i)
        # Plot ideogram and get maximum probe position in ideogram:
        plotIdeogram(chrom = k, cyto.text = op$cyto.text, cyto.data = op$assembly, cex = op$cex.cytotext, unit = op$plot.unit)

        # Get maximum position for this chromosome:
        xmaxI <- chromMax(chrom = k, cyto.data = op$assembly, pos.unit = op$plot.unit)
        xlim <- c(0, xmaxI)

        new <- TRUE
      }

      # Plot-dimensions:
      frame.c$bot <- frame.c$bot + (frame.c$top - frame.c$bot) * op$ideo.frac
      par(fig = unlist(frame.c), new = new, mar = op$mar)

      # Limits:
      if (!is.null(op$xlim)) {
        xlim <- op$xlim
      }

      # empty plot set up:
      plot(1, 1,
        type = "n", ylim = c(0, nSample), xlim = xlim, ylab = op$ylab, xlab = op$xlab,
        xaxs = "i", yaxt = "n", xaxt = "n", yaxs = "i", mgp = op$mgp, main = ""
      )

      # Let background be white:
      rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = "white")

      # main title for this plot
      title(main = op$main[c], line = op$main.line, cex.main = op$cex.main)

      # Plot heatmap for each sample:
      for (i in 1:nSample) {
        sample.segments <- chrom.segments[chrom.segments[, 1] == sampleID[i], ]
        # Adjust positions to be plotted along xaxis; i.e. scale according to plot.unit, and get left and right pos for rectangles to be plotted (either continuous
        # or 1 probe long):
        x <- adjustSegPos(chrom = sample.segments[, 2], char.arms = sample.segments[, 3], start = sample.segments[, 4], stop = sample.segments[, 5], type = "chromosome", unit = pos.unit, op = op)
        xleft <- x$use.start
        xright <- x$use.stop

        # Find appropriate colour for each probe:
        heat.col <- rep("white", nrow(sample.segments))
        # heat.col[sample.segments[,2]%%2==0] <- "grey95"
        heat.col[sample.segments[, 7] > upper.lim[t]] <- op$colors[2]
        heat.col[sample.segments[, 7] < lower.lim[t]] <- op$colors[1]
        ytop <- i - op$sep.samples
        ybottom <- i - (1 - op$sep.samples)

        rect(xleft, ybottom, xright, ytop, col = heat.col, border = NA)

        # Add sampleid on yaxis
        if (op$sample.labels) {
          axis(side = 2, at = (ytop - (ytop - ybottom) / 2), labels = sampleID[i], line = op$sample.line, tcl = 0, cex.axis = op$sample.cex, las = 1, mgp = op$mgp, tick = FALSE)
        }
      } # endfor

      if (!op$plot.ideo) {
        # Add xaxis:
        at.x <- get.xticks(xlim[1], xlim[2], unit = op$plot.unit, ideal.n = 6)
        axis(side = 1, tcl = -0.2, at = at.x, cex.axis = op$cex.axis, mgp = op$mgp)
        title(xlab = op$xlab, cex.lab = op$cex.lab, line = op$mgp[1])
      }


      # Add box around plot:
      abline(v = xlim)
      abline(h = c(0, nSample))



      # If page is full; start plotting on new page
      if (c %% (nr * nc) == 0 && c != nChrom) {
        # Add main title to page:
        title(op$title[t], outer = TRUE)

        devAskNewPage(ask = TRUE)

        # Reset columns and row in layout:
        clm <- 1
        row <- 1
        new <- FALSE
      } else {
        # Update column and row index:
        if (clm < nc) {
          clm <- clm + 1
        } else {
          clm <- 1
          row <- row + 1
        } # endif
        new <- TRUE
      } # endif
    } # endfor

    # Add main title to page:
    title(op$title[t], outer = TRUE)
    if (t != nT) {
      devAskNewPage(ask = TRUE)
    }
  } # endfor
} # endfunction