R/BarplotOnOff.R

In CellScore: Tool for Evaluation of Cell Identity from Transcription Profiles

# BarplotOnOff.R
#
# from version OnOff_v1.3.1.R

#' Pyramid barplot of on/off scores
#'
#' This function will generate a horizontal barplot of the OnOff scores of test
#' cell types, as defined by the eset$sub_cell_type1 column of the
#' input dataset. Note that if the cell types as provided in the second argument
#' (score data frame as produced by the function OnOff, are not matching the
#' phenotype of the input dataset, the function will return an error.
#' @param eset an ExpressionSet containing data matrices of normalized
#'   expression data, present/absent calls, a gene annotation data.frame and
#'   a phenotype data.frame.
#' @param group.score a data frame with cell type specific on/off scores as
#'   generated by the OnOff function.
#' @return This function returns a list of two objects, as follows:
#' %\describe{
#'   \item{GroupComparisonsForPlot}{an ordered data.frame of on/off scores,}
#'   \item{OnOffBarplotData}{a data frame of marker gain/loss and
#'   aditional features, used for making the plot.}
#'   %}
#' @keywords onoff score, barplot
#' @seealso \code{\link[CellScore]{OnOff}} for details on on/off score
#'   calculations, and \code{\link[hgu133plus2CellScore]{hgu133plus2CellScore}} for details
#'   on the specific expressionSet object that should be provided as an input.
#' @export
#' @importClassesFrom Biobase ExpressionSet
#' @importMethodsFrom Biobase pData
#' @importFrom graphics par barplot plot mtext legend
#' @examples
#' ## Load the expression set for the standard cell types
#' library(Biobase)
#' library(hgu133plus2CellScore) # eset.std
#'
#' ## Locate the external data files in the CellScore package
#' rdata.path <- system.file("extdata", "eset48.RData", package = "CellScore")
#' tsvdata.path <- system.file("extdata", "cell_change_test.tsv",
#'                             package = "CellScore")
#'
#' if (file.exists(rdata.path) && file.exists(tsvdata.path)) {
#'
#'    ## Load the expression set with normalized expressions of 48 test samples
#'    load(rdata.path)
#'
#'    ## Import the cell change info for the loaded test samples
#'    cell.change <- read.delim(file= tsvdata.path, sep="\t",
#'                              header=TRUE, stringsAsFactors=FALSE)
#'
#'    ## Combine the standards and the test data
#'    eset <- combine(eset.std, eset48)
#'
#'    ## Generate a marker list
#'    group.OnOff <- OnOff(eset, cell.change, out.put="marker.list")
#'
#'    ## Calculate on/off score for individual samples
#'    individ.OnOff <- OnOff(eset, cell.change, out.put="individual")
#'    
#'    ## Plot pyramid bar plot of on/off scores
#'    BarplotOnOff(eset, group.OnOff$scores)
#' }

BarplotOnOff <- function(eset, group.score) {

    ############################################################################
    ## PART 0. Check function arguments
    ############################################################################
    fun.main <- deparse(match.call()[[1]])
    .stopIfNotExpressionSet(eset, 'eset', fun.main)
    .stopIfNotDataFrame(group.score, 'group.score', fun.main)

    ############################################################################
    ## PART I. Filter samples according to phenoData including their scores
    ############################################################################
    pdata <- .filterPheno(pData(eset), fun.main, 'na')
    ## Select only rows that have group scores for derived cells-donor cell
    sel <- group.score$test %in% pdata$sub_cell_type1
    if (sum(sel) == 0) {
        stop(paste("Test cell type in the score data frame are not defined",
                    "exiting function", fun.main))
    }

    ## Table with specific group comparisons: use this for plotting
    mytable <- group.score[sel,]
    mytable <- mytable[order(mytable$OnOffScore, decreasing=FALSE),]

    ############################################################################
    ## PART II. Prepare data for plotting
    ############################################################################
    plot.me <-
        data.frame(position=rep(c("Below", "Above"), each=nrow(mytable)),
                   markers=c(mytable$loss.start.mkrs, mytable$gain.target.mkrs),
                   markername=rep(c("loss.start.mkrs", "gain.target.mkrs"),
                                  each=nrow(mytable)),
                   tags=c(mytable$test, mytable$test),
                   stringsAsFactors=FALSE)

    ## Transition names
    transition <- paste(mytable$start, mytable$target, sep="->")

    map.table = data.frame(position=c("Below", "Above"),
                           col=c("royalblue", "limegreen") ,
                           legend=c("loss of donor markers",
                                    "gain of target markers"),
                           abline=c(-1, 1),
                           axes=c(FALSE, TRUE),
                           stringsAsFactors=FALSE)

    ############################################################################
    ## PART III. Plot
    ############################################################################
    old.par <- par(no.readonly = TRUE)

    ## Get the y-values for bar position
    y.values <- barplot(height=plot.me$markers[plot.me$position == 'Below'],
                        plot=FALSE)
    ## Set up plot but don't draw anything
    cex.lab <- 0.85
    par(mar=c(5, 7, 4, 7) + 0.1)
    plot(c(-1, 1), c(0, max(y.values)+5), type = "n",
         main="Transition progression scored by on/off genes",
         xlab="fraction of markers",
         yaxt="n", ylab="",
         axes=FALSE,
         frame.plot=TRUE,
         cex.main=1.5,
         cex.lab=1.5,
         bty="L")

    lapply(seq_len(nrow(map.table)),
          function(i){
              instance <- map.table[i, ]
              h.value <- plot.me$markers[plot.me$position == instance$position]
              barplot(height=instance$abline*h.value,
                      add=TRUE, horiz=TRUE,
                      axes=instance$axes,
                      col=instance$col,
                      names.arg=transition,
                      las=2, cex.names=cex.lab, cex.axis=1.2)
              abline(v=instance$abline, col="red", lty=2)
          })
    ## Location of right-hand side axis labels
    mtext(side=4, mytable$test, at=y.values, las=2, cex=cex.lab)
    legend("top", legend=map.table$legend, fill=map.table$col, cex=1)

    ## Reset graphical parameters
    par(old.par)

    return(list(GroupComparisonsForPlot=mytable, OnOffBarplotData=plot.me))
}