R/barplot-methods.R
In pappewaio/AllelicImbalance: Investigates Allele Specific Expression
Defines functions
barplotLatticeFraction
Documented in
barplotLatticeFraction
#' barplot ASEset objects
#'
#' Generates barplots for ASEset objects. Two levels of plotting detail are
#' provided: a detailed barplot of read counts by allele useful for fewer
#' samples and SNPs, and a less detailed barplot of the fraction of imbalance,
#' useful for more samples and SNPs.
#'
#' \code{filter.pValue.fraction} is intended to remove p-value annotation with
#' very large difference in frequency, which could just be a sequencing
#' mistake. This is to avoid p-values like 1e-235 or similar.
#'
#' \code{sampleColour}User specified colours, either given as named colours
#' ('red', 'blue', etc) or as hexadecimal code. Can be either length 1 for all
#' samples, or else of a length corresponding to the number of samples for
#' individual colouring.
#'
#' @name ASEset-barplot
#' @rdname ASEset-barplot
#' @aliases ASEset-barplot barplot barplot,ASEset-method
#' @docType methods
#' @param height An \code{ASEset} object
#' @param type 'count' or 'fraction'
#' @param sampleColour.top User specified colours for top fraction
#' @param sampleColour.bot User specified colours for bottom fraction
#' @param legend Display legend
#' @param pValue Display p-value
#' @param strand four options, '+', '-', 'both' or '*'
#' @param testValue if set, a matrix or vector with user p-values
#' @param testValue2 if set, a matrix or vector with user p-values
#' @param OrgDb an OrgDb object which provides annotation
#' @param TxDb a TxDb object which provides annotation
#' @param annotationType select one or more from
#' 'gene','exon','transcript','cds'.
#' @param main text to use as main label
#' @param ylim set plot y-axis limit
#' @param yaxis wheter the y-axis is to be displayed or not
#' @param xaxis wheter the x-axis is to be displayed or not
#' @param ylab showing labels for the tic marks
#' @param ylab.text ylab text
#' @param xlab showing labels for the tic marks
#' @param xlab.text xlab text
#' @param legend.colnames gives colnames to the legend matrix
#' @param las.ylab orientation of ylab text
#' @param las.xlab orientation of xlab text
#' @param cex.main set main label size (max 2)
#' @param cex.pValue set pValue label size
#' @param cex.ylab set ylab label size
#' @param cex.xlab set xlab label size
#' @param cex.legend set legend label size
#' @param cex.annotation size of annotation text
#' @param ypos.annotation relative ypos for annotation text
#' @param annotation.interspace space between annotation text
#' @param legend.interspace set legend space between fills and text
#' @param add \code{boolean} indicates if a new device should be started
#' @param lowerLeftCorner integer that is only useful when \code{add}=TRUE
#' @param size Used internally by locationplot. Rescales each small barplot
#' window
#' @param addHorizontalLine adds a horizontal line that marks the default
#' fraction of 0.5 - 0.5
#' @param add.frame \code{boolean} to give the new plot a frame or not
#' @param filter.pValue.fraction \code{numeric} between 0 and 1 that filter
#' away pValues where the main allele has this frequency.
#' @param legend.fill.size size of the fill/boxes in the legend (default:NULL)
#' @param verbose Makes function more talkative
#' @param top.fraction.criteria 'maxcount', 'ref' or 'phase'
#' @param ... for simpler generics when extending function
#' @author Jesper R. Gadin, Lasse Folkersen
#' @seealso \itemize{ \item The \code{\link{ASEset}} class which the barplot
#' function can be called up on. }
#' @keywords barplot
#' @examples
#'
#' data(ASEset)
#' barplot(ASEset[1])
#'
#' @importFrom graphics plot
#' @importFrom graphics barplot
#' @exportMethod barplot
NULL
#' @rdname ASEset-barplot
setGeneric("barplot")
#' @rdname ASEset-barplot
setMethod("barplot", signature(height = "ASEset"), function(height, type = "count",
sampleColour.top = NULL, sampleColour.bot = NULL, legend = TRUE, pValue = TRUE, strand = "*", testValue = NULL,
testValue2 = NULL, OrgDb = NULL, TxDb = NULL, annotationType = c("gene", "exon",
"transcript"), main = NULL, ylim = NULL, yaxis = TRUE, xaxis = FALSE, ylab = TRUE,
ylab.text=NULL, xlab.text="samples",
xlab = TRUE, legend.colnames = "", las.ylab = 1, las.xlab = 2, cex.main = 1,
cex.pValue = 0.7, cex.ylab = 0.7, cex.xlab = 0.7, cex.legend = 0.6, add = FALSE,
lowerLeftCorner = c(0, 0), size = c(1, 1), addHorizontalLine = 0.5, add.frame = TRUE,
filter.pValue.fraction = 0.99, legend.fill.size=1, legend.interspace=1, verbose = FALSE,
top.fraction.criteria="maxcount", cex.annotation=0.7, ypos.annotation=0, annotation.interspace=1,
...) {
# catch useful graphical parameters that can be used to later add onto plot. This
# list will be retireved by using 'glst <- barplot(x)'
graphParamList <- list()
x <- height
# check type
okPlotTypes <- c("fraction", "count")
if (!type %in% okPlotTypes)
stop(paste("type can't be '", type, "' - it should be one of these: ", paste(okPlotTypes,
collapse = ", "), sep = ""))
# check verbose
if (class(verbose)[1] != "logical")
stop(paste("verbose should be of class logical, not", class(verbose)[1]))
if (length(verbose) != 1)
stop(paste("verbose should be of length 1, not", length(verbose)))
# check strand
okStrandTypes <- c("both", "+", "-", "*")
if (!strand %in% okStrandTypes)
stop(paste("strand can't be '", strand, "' - it should be one of these: ",
paste(okStrandTypes, collapse = ", "), sep = ""))
# check if strand type is present in object
# there is a check at the alleleCountsList aswell.
# And since the new ASEset version the strands will always be present
#
#if (strand == "+") {
# el <- "countsPlus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else if (strand == "-") {
# el <- "countsMinus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else if (strand == "*") {
# el <- "countsUnknown"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else if (strand == "both") {
# el <- "countsPlus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
# el <- "countsMinus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else {
# stop("unknown strand option")
#}
# check legend
if (!is.logical(legend)) {
stop("legend has to be logical, TRUE or FALSE")
}
if (!length(legend) == 1) {
stop("legend has to be of length 1")
}
# check pValue
if (!is.logical(pValue)) {
stop("pValue has to be logical, TRUE or FALSE")
}
if (!length(pValue) == 1) {
stop("pValue has to be of length 1")
}
# check ylim
if (!is.null(ylim)) {
}
# check testValue
if (!is.null(testValue)) {
if (!class(testValue)[1] == "numeric" & !class(testValue)[1] == "matrix") {
stop("wrong class for testValue")
}
if (!class(testValue)[1] == "matrix") {
if (!sum(dim(testValue) == rev(dim(x))) == 2) {
stop("wrong dimensions of testValue")
}
}
if (!class(testValue)[1] == "numeric") {
if (!length(testValue) == ncol(x)) {
stop("wrong dimensions of testValue")
}
}
}
# check testValue2
if (!is.null(testValue2)) {
if (!class(testValue2)[1] == "numeric" & !class(testValue2)[1] == "matrix") {
stop("wrong class for testValue")
}
if (!class(testValue2)[1] == "matrix") {
if (!sum(dim(testValue2) == rev(dim(x))) == 2) {
stop("wrong dimensions of testValue")
}
}
if (!class(testValue2)[1] == "numeric") {
if (!length(testValue2) == ncol(x)) {
stop("wrong dimensions of testValue")
}
}
}
# check testValue combinations
if (!strand == "both" & !is.null(testValue2)) {
stop("testValue2 is supposed to only be \n\t\t\tspecified when using strand='both', and then be the values for the minus strand.")
}
# check ylim
if (!is.null(ylim)) {
if (type == "fraction") {
warning("an ylim argument given to barplot with type fraction is ignored")
ylim <- NULL
} else {
if (class(ylim)[1] != "numeric")
stop(paste("if given, ylim must be of class numeric, not", class(ylim)[1]))
if (length(ylim) != 2)
stop(paste("if given, ylim must be of length 2, not", length(ylim)))
if (ylim[2] < ylim[1])
stop(paste("if ylim is given, ylim[2] must be larger than ylim[1]"))
}
}
# set easy accessor to samples
samples <- colnames(x)
# set easy accessor to snps
snps <- rownames(x)
# check sampleColour
if (is.null(sampleColour.top)) {
sampleColour.top <- rep("darkolivegreen4", length(samples))
}
if (is.null(sampleColour.bot)) {
sampleColour.bot <- rep("brown2", length(samples))
}
# if (is.null(sampleColour)) {
# sampleColour <- rep("dodgerblue", length(samples))
# } else {
# if (class(sampleColour) != "character")
# stop(paste("if given, sampleColour should be of class character, not",
# class(sampleColour)))
# if (length(sampleColour) == 1) {
# sampleColour <- rep(sampleColour, length(samples))
# } else {
# if (length(sampleColour) != length(samples))
# stop(paste("if given, sampleColour should be of length 1 or length",
# length(samples), "(corresponding to the number of samples)"))
# }
# hexadecimals <- sampleColour[grep("^#[0-9A-F]{6}", sampleColour)]
# missing <- sampleColour[!(sampleColour %in% colors() | sampleColour %in%
# hexadecimals)]
# if (length(missing) > 0) {
# stop(paste("The following", length(unique(missing)), "colour(s) from sampleColour were not recognised:",
# paste(sort(unique(missing)), collapse = ", ")))
# }
# }
# check OrgDb
if (!is.null(OrgDb)) {
if (!class(OrgDb)[1] == "OrgDb") {
stop("class of OrgDb has to be an OrgDb class")
}
}
if (!is.null(TxDb)) {
if (!class(TxDb)[1] == "TxDb") {
stop("class of TxDb has to be a TxDb class")
}
}
# check annotationType
okAnnotationTypes <- c("gene", "exon", "cds", "transcript")
if (!sum(annotationType %in% okAnnotationTypes) == length(annotationType)) {
stop(paste("type can't be '", annotationType, "' - it should be one of these: ",
paste(okAnnotationTypes, collapse = ", "), sep = ""))
}
# check main
if (!is.null(main)) {
if (!class(main)[1] == "character") {
stop("main has to be of class character")
}
}
# check yaxis
if (!is.logical(yaxis)) {
stop("yaxis has to be logical, TRUE or FALSE")
}
if (!length(yaxis) == 1) {
stop("yaxis has to be of length 1")
}
# check xaxis
if (!is.logical(xaxis)) {
stop("xaxis has to be logical, TRUE or FALSE")
}
if (!length(xaxis) == 1) {
stop("xaxis has to be of length 1")
}
# check ylab
if (!is.logical(ylab)) {
stop("ylab has to be logical, TRUE or FALSE")
}
if (!length(ylab) == 1) {
stop("ylab has to be of length 1")
}
# check xlab
if (!is.logical(xlab)) {
stop("xlab has to be logical, TRUE or FALSE")
}
if (!length(xlab) == 1) {
stop("xlab has to be of length 1")
}
# check legend.colnames
if (!class(legend.colnames)[1] == "character") {
stop("legend.colnames has to be of class character")
}
# check las.ylab
if (!class(las.ylab)[1] == "numeric" & !class(las.ylab)[1] == "integer") {
stop("las.ylab has to be of class numeric or integer")
}
if (!las.ylab == 1 & !las.ylab == 2) {
stop("las.ylab has to be a number, either 1 or 2")
}
# check las.xlab
if (!class(las.xlab)[1] == "numeric" & !class(las.xlab)[1] == "integer") {
stop("las.xlab has to be of class numeric or integer")
}
if (!las.xlab == 1 & !las.xlab == 2) {
stop("las.xlab has to be a number, either 1 or 2")
}
# check cex all of them
if (!class(cex.main)[1] == "numeric") {
if (cex.main <= 0 | !length(cex.main) == 1)
stop("cex.main has to be of class numeric, and have length 1, and a value above 0 ")
}
if (!class(cex.pValue)[1] == "numeric") {
if (cex.pValue > 1 | cex.pValue <= 0 | !length(cex.pValue) == 1)
stop("cex.pValue has to be of class numeric, and have length 1, and a value above 0")
}
if (!class(cex.ylab)[1] == "numeric") {
if (cex.ylab <= 0 | !length(cex.ylab) == 1)
stop("cex.ylab has to be of class numeric, and have length 1, and a value above 0 ")
}
if (!class(cex.xlab)[1] == "numeric") {
if (cex.xlab <= 0 | !length(cex.xlab) == 1)
stop("cex.xlab has to be of class numeric, and have length 1, and a value above 0 ")
}
if (!class(cex.legend)[1] == "numeric") {
if (cex.legend <= 0 | !length(cex.legend) == 1)
stop("cex.legend has to be of class numeric, and have length 1, and a value above 0 ")
}
# check add
if (!is.logical(add)) {
stop("add has to be logical, TRUE or FALSE")
}
if (!length(add) == 1) {
stop("add has to be of length 1")
}
# check lowerLeftCorner
if (!class(lowerLeftCorner)[1] == "numeric" & !class(lowerLeftCorner)[1] == "integer") {
stop("lowerLeftCorner has to be of class numeric or integer")
}
if (!length(lowerLeftCorner) == 2) {
stop("lowerLeftCorner has be of length 2")
}
# check size
if (!class(size)[1] == "numeric" & !class(size)[1] == "integer") {
stop("size has to be of class numeric or integer")
}
if (!length(size) == 2) {
stop("size has be of length 2")
}
# check addHorizontalLine
if (!is.null(addHorizontalLine)) {
if (!class(addHorizontalLine)[1] == "numeric") {
stop("addHorizontalLine has to be either NULL or of class numeric")
}
if (!length(addHorizontalLine) == 1) {
stop("addHorizontalLine has be of length 1")
}
if (addHorizontalLine > 1 | addHorizontalLine <= 0) {
stop("addHorizontalLine has be numeric between 0 and 1")
}
}
# check add.frame
if (!is.logical(add.frame)) {
stop("add.frame has to be logical, TRUE or FALSE")
}
if (!length(add.frame) == 1) {
stop("add.frame has to be of length 1")
}
# check filter.pValue.fraction=0.99
if (!class(filter.pValue.fraction)[1] == "numeric") {
stop("filter.pValue.fraction has to be of class numeric")
}
if (!length(filter.pValue.fraction) == 1) {
stop("filter.pValue.fraction has to be of length 1")
}
if (filter.pValue.fraction < 0 | filter.pValue.fraction > 1) {
stop("filter.pValue.fraction has to be a value between 0 and 1")
}
# make annotation dataframes
if ((!is.null(OrgDb) | !is.null(TxDb)) & (strand == "+" |
strand == "*" | strand == "both")) {
annDfPlus <- getAnnotationDataFrame(rowRanges(x), strand = "+", annotationType = annotationType,
OrgDb = OrgDb, TxDb = TxDb, verbose = verbose)
} else {
annDfPlus <- NULL
}
if ((!is.null(OrgDb) | !is.null(TxDb)) & (strand == "-" |
strand == "*" | strand == "both")) {
annDfMinus <- getAnnotationDataFrame(rowRanges(x), strand = "-", annotationType = annotationType,
OrgDb = OrgDb, TxDb = TxDb, verbose = verbose)
} else {
annDfMinus <- NULL
}
# creating a background from foreground colour - basically just turning down the
# saturation a little compared to the foreground color except for very
# low-saturated fg's which are made more saturated, and very dark fg's which are
# made lighter
#hsvColoursBg <- hsvColoursFg <- rgb2hsv(col2rgb(sampleColour))
#veryLowSaturation <- hsvColoursFg["s", ] < 0.3
#veryDark <- hsvColoursFg["v", ] < 0.2
#remaining <- !veryLowSaturation & !veryDark
#hsvColoursBg["s", veryLowSaturation] <- hsvColoursFg["s", veryLowSaturation] +
# 0.3
#hsvColoursBg["v", veryDark] <- hsvColoursFg["v", veryDark] + 0.2
#hsvColoursBg["s", remaining] <- hsvColoursFg["s", remaining] - 0.4
#fgCol <- hsv(hsvColoursFg["h", ], hsvColoursFg["s", ], hsvColoursFg["v", ])
#bgCol <- hsv(hsvColoursBg["h", ], hsvColoursBg["s", ], hsvColoursBg["v", ])
#fgCol <- hsv(rgb2hsv(col2rgb(sampleColour.bot))["h",])
#bgCol <- hsv(rgb2hsv(col2rgb(sampleColour.top))["h",])
fgCol <- sampleColour.bot
bgCol <- sampleColour.top
if (type == "count") {
# the intention of this is to be able to create barplots which show both the
# amount of counts and the distribution between alleles this is the least
# digested format and consequently produces the most cluttered plots. It is
# therefore intended for subsets, few snps, few samples etc in these smaller
# subsets it is useful to get information on higher expression values at certain
# points of genes etc.
if (strand == "+" | strand == "-" | strand == "*" ) {
for (i in 1:length(snps)) {
# getting revelant data
snp <- snps[i]
tmp <- t(alleleCounts(x, strand)[[snp]][samples, , drop = FALSE])
if ((ncol(tmp) * nrow(tmp)) == sum(tmp == 0)) {
zeroRows <- TRUE
} else {
zeroRows <- FALSE
}
# alleles are
n <- names(sort(apply(tmp, 1, sum), decreasing = TRUE)[1:2])
# theor total count
tc <- sort(apply(tmp, 1, sum), decreasing = TRUE)[1:2]
# extracting only major and minor allele in correctly transposed format and
# ordering colours accordingly
height <- tmp[match(n, rownames(tmp)), , drop = FALSE]
# change all zero rows to X names
rownames(height)[apply(height, 1, sum) == 0] <- "X"
# warn if the two remaining alleles have too many counts
majorAndMinorFraction <- sum(height)/sum(tmp)
if (verbose & majorAndMinorFraction < 0.9) {
message(paste("Snp", snp, "was possible tri-allelic, but only two most frequent alleles were plotted. Counts:\n",
paste(paste(names(apply(tmp, 1, sum)), apply(tmp, 1, sum), sep = "="),
collapse = ", ")))
}
## make bars section
extraspaceOnTop <- 1.15
if (is.null(ylim)) {
if (!zeroRows) {
ylimHere <- c(0, max(abs(range(height, na.rm = TRUE))) * extraspaceOnTop)
} else {
ylimHere <- c(0, 10)
}
} else {
ylimHere <- ylim
}
# init frame
if (!add) {
size <- c(1, 1)
lowerLeftCorner <- c(0, 0)
plot.default(NULL, xlim = c(0, 1), ylim = c(0, 1), ylab = "", xlab = "",
xaxt = "n", yaxt = "n", frame.plot = FALSE)
}
# set so the ylimHere is at least 10.
if (ylimHere[2] < 10) {
ylimHere[2] <- 10
}
# make tics and labels prettier
ylabels <- pretty(ylimHere[1]:ylimHere[2])
if (strand == "both") {
ylimHere[1] <- -max(abs(ylabels))
ylimHere[2] <- max(abs(ylabels))
} else {
ylimHere[1] <- min(ylabels)
ylimHere[2] <- max(ylabels)
}
# all x positions
xPointsAll <- seq(lowerLeftCorner[1], lowerLeftCorner[1] + size[1],
length.out = length(height) * 1.5 + 2)[2:(length(height) * 1.5 +
1)]
widthsOfRectangles <- rep(size[1]/(length(height) * 1.5 + 1), length(height)/2)
# normalize toPlot-values to a range of c(0,size[2])
heightsOfRectangles <- ((height - ylimHere[1])/(ylimHere[2] - ylimHere[1]) *
size[2])
# take every third element (start at 2)
xPoints1 <- xPointsAll[seq(1, length(height) * 1.5, by = 3)]
yPoints <- lowerLeftCorner[2] + heightsOfRectangles[1, ]/2
# alleleColors1 <- rep(fgCol[1],9)
alleleColors1 <- fgCol
symbols(x = xPoints1, y = yPoints, bg = alleleColors1, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectangles[1, ])), add = TRUE, inches = FALSE)
# take every third element (start at 3)
xPoints2 <- xPointsAll[seq(2, length(height) * 1.5, by = 3)]
yPoints <- lowerLeftCorner[2] + heightsOfRectangles[2, ]/2
alleleColors2 <- rep(bgCol[1], 9)
alleleColors2 <- bgCol
widthsOfRectangles <- rep(size[1]/(length(height) * 1.5 + 1), length(height)/2)
symbols(x = xPoints2, y = yPoints, bg = alleleColors2, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectangles[2, ])), add = TRUE, inches = FALSE)
## make axis and labels section
# y-axis
yat <- (ylabels + abs(min(ylabels)))/(sum(abs(min(ylabels)) + max(ylabels)))
graphParamList[["yat"]] <- yat
if (yaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1]), c(lowerLeftCorner[2],
lowerLeftCorner[2] + 1))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x - size[1] * 0.01), c(y, y))
}, rep(lowerLeftCorner[1], length(yat)), (lowerLeftCorner[2] +
yat)))
}
if (ylab) {
if (las.ylab == 2) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.ylab)
}
if (las.ylab == 1) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.ylab)
}
}
# x-axis labels
xat <- xPoints1 + (xPoints2 - xPoints1)/2
graphParamList[["xat"]] <- xat
if (xaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1] + 1), c(lowerLeftCorner[2],
lowerLeftCorner[2]))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x), c(y, y - size[1] * 0.01))
}, (lowerLeftCorner[2] + xat), rep(lowerLeftCorner[1], length(xat))))
}
if (xlab) {
if (las.xlab == 2) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.xlab)
}
if (las.xlab == 1) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.xlab)
}
}
if (is.null(main)) {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = paste(snp), cex = cex.main,
adj = 0.5, xpd = TRUE)
} else {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = main[i], cex = cex.main, adj = 0.5)
}
if (pValue) {
# check pValue slot
if (!is.null(testValue)) {
if (class(testValue)[1] == "matrix") {
text(xat, apply(heightsOfRectangles, 2, max), signif(testValue[,
i], 1), cex = cex.pValue, adj = c(0.5, 0.01))
} else {
text(xat, apply(heightsOfRectangles, 2, max), signif(testValue,
1), cex = cex.pValue, adj = c(0.5, 0.01))
}
} else {
# use chisq by default
TFf <- t(fraction(x[i], strand = strand) < filter.pValue.fraction)
ct <- signif(chisq.test(x[i], strand), 1)
ct[!TFf] <- NA
text(xat, apply(heightsOfRectangles, 2, max), ct, cex = cex.pValue,
adj = c(0.5, 0.01))
}
}
if (!is.null(OrgDb) | !is.null(TxDb)) {
annotationBarplot(strand,i, lowerLeftCorner, annDfPlus, annDfMinus,
cex=cex.annotation, ypos=ypos.annotation,
interspace=annotation.interspace)
}
if (legend) {
# any strand
o <- apply(height, 1, sum)
# remove zeros count alleles
TFz <- rownames(height) == "X"
o <- o[!TFz]
textOver <- rownames(height)[!TFz]
# set another xLegendPos and yLegendPos for '*'.
if ((strand == "*" ) & (!is.null(OrgDb) |
!is.null(TxDb))) {
xLegendPos <- 0.5
yLegendPos <- 1.00
} else {
xLegendPos <- 0.94
yLegendPos <- 1.00
}
legendBarplot(lowerLeftCorner, size, rownames=textOver, colnames=legend.colnames,
boxsize=legend.fill.size, boxspace=legend.interspace, fgCol, bgCol,
ylegendPos=1, xlegendPos=0.96, cex=cex.legend)
}
#set y and x labs
if(is.null(ylab.text)){
ylab.text <- "reads"
}
title(ylab=ylab.text, xlab=xlab.text)
}
} else if (strand == "both") {
for (i in 1:length(snps)) {
if (verbose) {
cat("plotting snp", i, "\n")
}
snp <- snps[i]
# getting revelant data forward strand
tmpPlus <- t(alleleCounts(x, "+")[[snp]][samples, ])
if ((ncol(tmpPlus) * nrow(tmpPlus)) == sum(tmpPlus == 0)) {
zeroRowsPlus <- TRUE
} else {
zeroRowsPlus <- FALSE
}
tmpMinus <- t(alleleCounts(x, "-")[[snp]][samples, ])
if ((ncol(tmpMinus) * nrow(tmpMinus)) == sum(tmpMinus == 0)) {
zeroRowsMinus <- TRUE
} else {
zeroRowsMinus <- FALSE
}
# check most common alleles on + and - (should be same)
allVar <- data.frame(A = c(NA, NA), T = c(NA, NA), G = c(NA, NA),
C = c(NA, NA), del = c(NA, NA))
o <- apply(tmpPlus, 1, sum)
u <- apply(tmpMinus, 1, sum)
# to get their total count
allVar[1, match(names(o), colnames(allVar))] <- o
allVar[2, match(names(u), colnames(allVar))] <- u
# save to tri-allelic warner
nAll <- names(sort(apply(allVar, 2, sum), decreasing = TRUE))
tcAll <- sort(apply(allVar, 2, sum), decreasing = TRUE)
# set colnames to X for non present alleles
xNames <- colnames(allVar)[apply(allVar, 2, sum) == 0]
colnames(allVar)[colnames(allVar) %in% xNames] <- "X"
xNames <- rownames(tmpPlus)[apply(tmpPlus, 1, sum) == 0]
rownames(tmpPlus)[rownames(tmpPlus) %in% xNames] <- "X"
xNames <- rownames(tmpMinus)[apply(tmpMinus, 1, sum) == 0]
rownames(tmpMinus)[rownames(tmpMinus) %in% xNames] <- "X"
# alleles are
n <- names(sort(apply(allVar, 2, sum), decreasing = TRUE))[1:2]
# theor total count
tc <- sort(apply(allVar, 2, sum), decreasing = TRUE)[1:2]
# extracting only major and minor allele in correctly transposed format and
# ordering colours accordingly
over <- tmpPlus[match(n, rownames(tmpPlus)), ]
under <- tmpMinus[match(n, rownames(tmpMinus)), ]
# change all NAs to X names and zero values
over[is.na(rownames(over)), ] <- 0
rownames(over)[is.na(rownames(over))] <- "X"
under[is.na(rownames(under)), ] <- 0
rownames(under)[is.na(rownames(under))] <- "X"
# Plus Warn if the two remaining alleles have too many counts
if (!zeroRowsPlus) {
majorAndMinorFraction <- sum(allVar[1, ][c(n[1], n[2])])/sum(allVar[1,
])
if (verbose & majorAndMinorFraction < 0.9) {
message(paste("Snp", snp, "was possible tri-allelic, but only two most frequent alleles were plotted. Counts:\n",
paste(paste(nAll, tcAll, sep = "="), collapse = ", ")))
}
}
# Minus calculating major and minor allele, warn if the two remaining alleles
# have too many counts
if (!zeroRowsMinus) {
majorAndMinorFraction <- sum(allVar[2, ][c(n[1], n[2])])/sum(allVar[2,
])
if (verbose & majorAndMinorFraction < 0.9) {
message(paste("Snp", snp, "was possible tri-allelic, but only two most frequent alleles were plotted. Counts:\n",
paste(paste(nAll, tcAll, sep = "="), collapse = ", ")))
}
}
# set under as minus
under <- -under
# plotting function
if (!add) {
# init frame
plot.default(NULL, xlim = c(0, 1), ylim = c(0, 1), ylab = "", xlab = "",
xaxt = "n", yaxt = "n", , frame.plot = FALSE)
}
heightPlus <- over
heightMinus <- under
if (is.null(ylim)) {
if (zeroRowsPlus) {
ymax <- 5
} else {
ymax <- max(abs(range(heightPlus, na.rm = TRUE)))
}
if (zeroRowsMinus) {
ymin <- -5
} else {
ymin <- min(heightMinus)
}
ylimHere <- c(ymin, ymax)
} else {
ylimHere <- ylim
}
# increase ylimHere before making pretty
extraspaceOnTopAndBottom <- 0.15
extraRealDistance <- (dist(ylimHere) * extraspaceOnTopAndBottom)
ylimHere[1] <- ylimHere[1] - extraRealDistance
ylimHere[2] <- ylimHere[2] + extraRealDistance
# set lowest max and min values
if (ylimHere[2] < 10) {
ylimHere[2] <- 10
}
if (abs(ylimHere[1]) < 10) {
ylimHere[1] <- -10
}
# make more even limits (make more pretty)
ylabels <- pretty(ylimHere[1]:ylimHere[2])
ylimHere[1] <- min(ylabels)
ylimHere[2] <- max(ylabels)
# breakpoint Plus Minus bars
middleLine <- (1/(ylimHere[2] - ylimHere[1])) * abs(ylimHere[1])
# make a xPointsAll framework
xPointsAll <- seq(lowerLeftCorner[1], lowerLeftCorner[1] + size[1],
length.out = length(heightPlus) * 1.5 + 2)[2:(length(heightPlus) *
1.5 + 1)]
widthsOfRectangles <- rep(size[1]/(length(heightPlus) * 1.5 + 1),
length(heightPlus)/2)
# normalize toPlot-values to a range of c(0,size[2])
heightsOfRectanglesPlus <- round(((heightPlus - ylimHere[1])/(ylimHere[2] -
ylimHere[1]) * size[2]) - middleLine, 5)
heightsOfRectanglesMinus <- round(((-heightMinus - ylimHere[1])/(ylimHere[2] -
ylimHere[1]) * size[2]) - middleLine, 5)
# take every third element (start at 2)
xPoints1 <- xPointsAll[seq(2, length(heightPlus) * 1.5, by = 3)]
yPointsPlus <- lowerLeftCorner[2] + middleLine + heightsOfRectanglesPlus[1,
]/2
barMinusEdge <- matrix(rep(NA, length(heightsOfRectanglesMinus)),
nrow = 2) #is used later to plot p-value
# yPointMinus is a little tricky and needs a more complicated equation
yPointsMinus1 <- lowerLeftCorner[2] + heightsOfRectanglesMinus[1,
]/2
yPointsMinus2 <- 1 - yPointsMinus1
yPointsMinus3 <- yPointsMinus2 - (1 - middleLine)
barMinusEdge[1, ] <- (1 - (lowerLeftCorner[2] + heightsOfRectanglesMinus[1,
])) - (1 - middleLine)
yPointsMinus <- yPointsMinus3
alleleColors1 <- fgCol
symbols(x = xPoints1, y = yPointsPlus, bg = alleleColors1, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesPlus[1, ])), add = TRUE,
inches = FALSE)
symbols(x = xPoints1, y = yPointsMinus, bg = alleleColors1, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesMinus[1, ])), add = TRUE,
inches = FALSE)
# take every third element (start at 3)
xPoints2 <- xPointsAll[seq(3, length(heightPlus) * 1.5, by = 3)]
yPointsPlus <- lowerLeftCorner[2] + middleLine + heightsOfRectanglesPlus[2,
]/2
# yPointMinus is a little tricky and needs a more complicated equation
yPointsMinus1 <- lowerLeftCorner[2] + heightsOfRectanglesMinus[2,
]/2
yPointsMinus2 <- 1 - yPointsMinus1
yPointsMinus3 <- yPointsMinus2 - (1 - middleLine)
barMinusEdge[2, ] <- (1 - (lowerLeftCorner[2] + heightsOfRectanglesMinus[2,
])) - (1 - middleLine)
yPointsMinus <- yPointsMinus3
# alleleColors2 <- rep('#7F5959',9)
alleleColors2 <- bgCol
symbols(x = xPoints2, y = yPointsPlus, bg = alleleColors2, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesPlus[2, ])), add = TRUE,
inches = FALSE)
symbols(x = xPoints2, y = yPointsMinus, bg = alleleColors2, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesMinus[2, ])), add = TRUE,
inches = FALSE)
# y-axis
yat <- (ylabels + abs(min(ylabels)))/(sum(abs(min(ylabels)) + max(ylabels)))
graphParamList[["yat"]] <- yat
if (yaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1]), c(lowerLeftCorner[2],
lowerLeftCorner[2] + 1))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x - size[1] * 0.01), c(y, y))
}, rep(lowerLeftCorner[1], length(yat)), (lowerLeftCorner[2] +
yat)))
}
if (ylab) {
if (las.ylab == 2) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.ylab)
}
if (las.ylab == 1) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.ylab)
}
}
# x-axis labels
xat <- xPoints1 + (xPoints2 - xPoints1)/2
graphParamList[["xat"]] <- xat
if (xaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1] + 1), c(lowerLeftCorner[2],
lowerLeftCorner[2]))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x), c(y, y - size[1] * 0.01))
}, (lowerLeftCorner[2] + xat), rep(lowerLeftCorner[1], length(xat))))
}
if (xlab) {
if (las.xlab == 2) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.xlab)
}
if (las.xlab == 1) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.xlab)
}
}
if (is.null(main)) {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = paste(snp), cex = cex.main,
adj = 0.5, xpd = TRUE)
} else {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = main[i], cex = cex.main, adj = 0.5)
}
if (pValue) {
# check pValue slot
if (!is.null(testValue)) {
if (class(testValue)[1] == "matrix") {
text(xat, apply(middleLine + heightsOfRectanglesPlus, 2, max),
signif(testValue[, i], 1), cex = cex.pValue, adj = c(0.4,
-0.3))
} else {
text(xat, apply(middleLine + heightsOfRectanglesPlus, 2, max),
signif(testValue, 1), cex = cex.pValue, adj = c(0.4, -0.3))
}
} else {
# use chisq by default
TFf <- t(fraction(x[i], strand = "+") < filter.pValue.fraction)
ct <- signif(chisq.test(x[i], "+"), 1)
ct[!TFf] <- NA
text(xat, apply(middleLine + heightsOfRectanglesPlus, 2, max),
ct, cex = cex.pValue, adj = c(0.4, -0.3))
}
if (!is.null(testValue2)) {
if (class(testValue)[1] == "matrix") {
text(xat, apply(barMinusEdge, 2, min), signif(testValue2[,
i], 1), cex = cex.pValue, adj = c(0.4, 1, 1))
} else {
text(xat, apply(barMinusEdge, 2, min), signif(testValue2, 1),
cex = cex.pValue, adj = c(0.4, 1, 1))
}
} else {
# use chisq by default
TFf <- t(fraction(x[i], strand = "-") < filter.pValue.fraction)
ct <- signif(chisq.test(x[i], "-"), 1)
ct[!TFf] <- NA
text(xat, apply(barMinusEdge, 2, min), ct, cex = cex.pValue,
adj = c(0.4, 1, 1))
}
}
if (!is.null(OrgDb) | !is.null(TxDb)) {
annotationBarplot(strand,i, lowerLeftCorner, annDfPlus, annDfMinus,
cex=cex.annotation, ypos=ypos.annotation,
interspace=annotation.interspace)
}
if (legend) {
# plus strand
TFz <- rownames(over) == "X"
if (!sum(!TFz) == 0) {
textOver <- paste(rownames(over)[!TFz], apply(over, 1, sum)[!TFz])
#remove the total count for that nucleotide
textOver <- unlist(lapply(strsplit(textOver," "),function(x){x[1]}))
legendBarplot(lowerLeftCorner, size, rownames=textOver, colnames=legend.colnames,
boxsize=legend.fill.size, boxspace=legend.interspace, fgCol, bgCol,
ylegendPos=1, xlegendPos=0.96, cex=cex.legend)
}
}
}
#set y and x labs
if(is.null(ylab.text)){
ylab.text <- "reads"
}
title(ylab=ylab.text, xlab=xlab.text)
}
} else if (type == "fraction") {
# a plot type that shows the allelic imbalance, irrespective of absolute read
# counts at each Snp. this produce plots that are more easily overviewed than
# those from 'plot_reads_by_allele', but of course omits some detail information
# from these plots. Intended to use for a medium number of snps and samples
if (strand == "both") {
stop("strand must be '+', '-' or '*' for type='fraction' ")
}
for (i in 1:length(snps)) {
# getting revelant data
snp <- snps[i]
fgColHere <- fgCol
#tmp <- alleleCounts(x, strand)[[snp]][samples, , drop = FALSE]
df <- fractionPlotDf(x, snp, strand=strand, top.fraction.criteria=top.fraction.criteria)
# check for zeroRows
if (sum(df$na=="no")==0) {
zeroRows <- TRUE
} else {
zeroRows <- FALSE
}
# setting no-count samples to colour grey90 and fraction 0
#fgColHere[df$na[seq(1,nrow(df),by=2)+1]=="yes"] <- "grey90"
bgCol[df$na[seq(1,nrow(df),by=2)+1]=="yes"] <- "grey90"
#fraction[is.nan(fraction)] <- 0
# if add=FALSE then make a new device
if (!add) {
plot.default(NULL, xlim = c(0, 1), ylim = c(0, 1), ylab = "", xlab = "",
xaxt = "n", yaxt = "n")
}
# ylim is always (0,1) and usergiven ylims are reset to this with a warning (in
# start checkups).
ylim <- c(0, 1)
# find the centers of rectangles
xPoints <- seq(lowerLeftCorner[1], lowerLeftCorner[1] + size[1], length.out = nrow(df)/2 +
2)[2:(nrow(df)/2 + 1)]
widthsOfRectangles <- rep(size[1]/(nrow(df)/2 + 1), nrow(df)/2)
# normalize fraction-values to a range of c(0,size[2])
heightsOfRectangles <- ((df$values[seq(1,nrow(df),by=2)] - ylim[1])/(ylim[2] - ylim[1])) * size[2]
# find the Y-coordinate center points for rectangles
yPoints <- lowerLeftCorner[2] + heightsOfRectangles/2
# if bgCol is given, start by painting background
if (!is.null(bgCol)) {
yPointsBg <- rep(lowerLeftCorner[2] + size[2]/2, length(xPoints))
rectanglesBg <- matrix(ncol = 2, c(widthsOfRectangles, rep(size[2],
length(widthsOfRectangles))))
symbols(x = xPoints, y = yPointsBg, bg = bgCol, rectangles = rectanglesBg,
add = TRUE, inches = FALSE)
}
symbols(x = xPoints, y = yPoints, bg = fgColHere, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectangles)), add = TRUE, inches = FALSE)
# add a frame to the plot
if (add.frame) {
symbols(x = lowerLeftCorner[1] + size[1]/2, y = lowerLeftCorner[2] +
size[2]/2, rectangles = matrix(ncol = 2, c(size[1], size[2])),
add = TRUE, inches = FALSE)
}
# add a tick markers (easier for fraction-plots because they are restricted to
# 0-100%
marks <- c(0.25, 0.5, 0.75)
names(marks) <- c("25%", "50%", "75%")
for (markName in names(marks)) {
mark <- marks[markName]
yBase <- lowerLeftCorner[2] + ((mark - ylim[1])/(ylim[2] - ylim[1])) *
size[2]
xBase <- xPoints[1] - widthsOfRectangles[1]/2
lines(x = c(xBase - size[1] * 0.01, xBase + size[1] * 0.01), y = c(yBase,
yBase))
text(x = xBase - size[1] * 0.015, y = yBase, label = markName, adj = 1,
cex = 0.7)
}
# writing the name of each Snp on top of the plot(within frame)
text(x = 0.5, y = 1.02, label = snp, xpd = TRUE, cex = cex.main)
# add horizontal line at 0.5
if (!is.null(addHorizontalLine)) {
lines(x = c(lowerLeftCorner[1], lowerLeftCorner[1] + size[1]), y = c(lowerLeftCorner[2] +
size[2] * addHorizontalLine, lowerLeftCorner[2] + size[2] * addHorizontalLine))
mapBiasHere <- mapBias(x)[[snp]]
#if (!all(mapBiasHere %in% c(0.5, 0)) & !zeroRows) {
# # Only activate if anything is different from 0.5 or 0
#
# for (j in 1:length(fraction)) { #fraction doesnt exist (check df$values instead)
# biasfraction <- mapBiasHere[j, majorAllele]/(mapBiasHere[j, majorAllele] +
# mapBiasHere[j, minorAllele])
# yPoint <- lowerLeftCorner[2] + ((biasfraction - ylim[1])/(ylim[2] -
# ylim[1])) * size[2]
# lines(x = c(xPoints[j] - widthsOfRectangles[j]/2, xPoints[j] +
# widthsOfRectangles[j]/2), y = c(yPoint, yPoint), lty = 2)
# }
#}
}
#set y and x labs
if(is.null(ylab.text)){
ylab.text <- "fraction"
}
title(ylab=ylab.text, xlab=xlab.text)
}
}
# return graphical paramlist (does not work for fraction type yet)
invisible(graphParamList)
})
#' gbarplot ASEset objects
#'
#' Generates gbarplots for ASEset objects. Two levels of plotting detail are
#' provided: a detailed gbarplot of read counts by allele useful for fewer
#' samples and SNPs, and a less detailed gbarplot of the fraction of imbalance,
#' useful for more samples and SNPs.
#'
#' This function serves the same purpose as the normal barplot, but with
#' trellis graphics using lattice, to be able to integrate well with Gviz track
#' functionality.
#'
#' @name ASEset-gbarplot
#' @aliases ASEset-gbarplot gbarplot gbarplot,ASEset-method
#' @docType methods
#' @param x An \code{ASEset} object
#' @param type 'count' or 'fraction'
#' @param strand four options, '+', '-', 'both' or '*'
#' @param verbose Makes function more talkative
#' @param ... for simpler generics when extending function
#' @author Jesper R. Gadin
#' @seealso \itemize{ \item The \code{\link{ASEset}} class which the gbarplot
#' function can be called up on. \item The \code{\link{barplot}} non trellis
#' barplot. }
#' @keywords gbarplot
#' @examples
#'
#' data(ASEset)
#' gbarplot(ASEset[1])
#'
#' @exportMethod gbarplot
setGeneric("gbarplot", function(x, type = "count", strand = "*",
verbose = FALSE, ...) {
standardGeneric("gbarplot")
})
setMethod("gbarplot", signature(x = "ASEset"), function(x, type = "count", strand = "*",
usePhase=FALSE, verbose = FALSE, ...) {
if (length(list(...)) == 0) {
e <- new.env(hash = TRUE)
} else {
e <- list2env(list(...))
}
#print(ls(envir=e))
if (!exists("mainvec", envir = e, inherits = FALSE)) {
e$mainvec <- rep("",nrow(x))
}
if (!exists("cex.mainvec", envir = e, inherits = FALSE)) {
e$cex.mainvec <- 1
}
if (!exists("ylab", envir = e, inherits = FALSE)) {
e$ylab <- "reads"
}
if (!exists("xlab", envir = e, inherits = FALSE)) {
e$xlab <- "samples"
}
if (!exists("middleLine", envir = e, inherits = FALSE)) {
e$middleLine <- TRUE
}
if (!exists("deAnnoPlot", envir = e, inherits = FALSE)) {
e$deAnnoPlot <- FALSE
}
if (!exists("top.fraction.criteria", envir = e, inherits = FALSE)) {
e$top.fraction.criteria <- "maxcount"
}
for (i in 1:nrow(x)) {
name <- rownames(x)[i]
if (type == "fraction") {
b <- barplotLatticeFraction(
identifier = name,
x=x,
astrand=strand,
ids=rownames(x),
ylab=e$ylab,
xlab=e$xlab,
mainvec=e$mainvec,
middleLine=e$middleLine,
deAnnoPlot=e$deAnnoPlot,
cex.mainvec=e$cex.mainvec,
top.fraction.criteria=e$top.fraction.criteria)
} else if (type == "count") {
b <- barplotLatticeCounts(
identifier = name,
x=x,
astrand=strand,
ids=rownames(x),
ylab=e$ylab,
xlab=e$xlab,
mainvec=e$mainvec)
} else {
stop("type has to be fraction or count")
}
b
}
})
barplotLatticeFraction <- function(identifier, ...) {
if (length(list(...)) == 0) {
e <- new.env(hash = TRUE)
} else {
e <- list2env(list(...))
}
e$ids <- unlist(e$ids)
e$strand <- e$astrand
df <- fractionPlotDf(e$x, identifier, strand=e$strand, top.fraction.criteria=e$top.fraction.criteria)
# Grah params
my_cols <- c("green", "red")
# set default values
parset <- list(superpose.polygon=list(col=my_cols))
scales = list(rot = c(90, 0))
if (e$deAnnoPlot) {
parset <- list(
layout.widths = list(
left.padding = 0,
axis.left = 0,
ylab.axis.padding = 0,
right.padding = 0,
axis.right = 0
),
layout.heights = list(
top.padding = 0.1,
between = 0.1,
xlab.top= 0.1,
axis.top = 0,
main=1.1,
main.key.padding=1,
axis.xlab.padding = 1,
bottom.padding = 1,
axis.bottom = 0.3
),
superpose.polygon=list(col=my_cols)
)
scales = list(y = list(at = NULL, labels = NULL),
x = list(labels = rep("",ncol(e$x)))
)
#, rot = c(90, 0))
}
if(!e$middleLine) {
b <- barchart(values ~ sample, group = groups, data = df, origin = 0,
stack = TRUE, scales = scales,
main = list(label=e$main, cex=e$cex.mainvec),
ylab = e$ylab, xlab = e$xlab,
par.settings = parset)
}else {
# if middle line at 0.5
b <- barchart(values ~ sample, group = groups, data = df, origin = 0,
stack = TRUE, scales = scales,
main = list(label=e$main, cex=e$cex.mainvec),
ylab = e$ylab, xlab = e$xlab,
par.settings = parset, panel=function(x, y, ...) {
panel.barchart(x, y, ...)
panel.abline(h=0.5, lty=1)
} )
}
b
}
pappewaio/AllelicImbalance documentation built on April 11, 2020, 2:58 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions barplotLatticeFraction
Documented in barplotLatticeFraction
#' barplot ASEset objects
#'
#' Generates barplots for ASEset objects. Two levels of plotting detail are
#' provided: a detailed barplot of read counts by allele useful for fewer
#' samples and SNPs, and a less detailed barplot of the fraction of imbalance,
#' useful for more samples and SNPs.
#'
#' \code{filter.pValue.fraction} is intended to remove p-value annotation with
#' very large difference in frequency, which could just be a sequencing
#' mistake. This is to avoid p-values like 1e-235 or similar.
#'
#' \code{sampleColour}User specified colours, either given as named colours
#' ('red', 'blue', etc) or as hexadecimal code. Can be either length 1 for all
#' samples, or else of a length corresponding to the number of samples for
#' individual colouring.
#'
#' @name ASEset-barplot
#' @rdname ASEset-barplot
#' @aliases ASEset-barplot barplot barplot,ASEset-method
#' @docType methods
#' @param height An \code{ASEset} object
#' @param type 'count' or 'fraction'
#' @param sampleColour.top User specified colours for top fraction
#' @param sampleColour.bot User specified colours for bottom fraction
#' @param legend Display legend
#' @param pValue Display p-value
#' @param strand four options, '+', '-', 'both' or '*'
#' @param testValue if set, a matrix or vector with user p-values
#' @param testValue2 if set, a matrix or vector with user p-values
#' @param OrgDb an OrgDb object which provides annotation
#' @param TxDb a TxDb object which provides annotation
#' @param annotationType select one or more from
#' 'gene','exon','transcript','cds'.
#' @param main text to use as main label
#' @param ylim set plot y-axis limit
#' @param yaxis wheter the y-axis is to be displayed or not
#' @param xaxis wheter the x-axis is to be displayed or not
#' @param ylab showing labels for the tic marks
#' @param ylab.text ylab text
#' @param xlab showing labels for the tic marks
#' @param xlab.text xlab text
#' @param legend.colnames gives colnames to the legend matrix
#' @param las.ylab orientation of ylab text
#' @param las.xlab orientation of xlab text
#' @param cex.main set main label size (max 2)
#' @param cex.pValue set pValue label size
#' @param cex.ylab set ylab label size
#' @param cex.xlab set xlab label size
#' @param cex.legend set legend label size
#' @param cex.annotation size of annotation text
#' @param ypos.annotation relative ypos for annotation text
#' @param annotation.interspace space between annotation text
#' @param legend.interspace set legend space between fills and text
#' @param add \code{boolean} indicates if a new device should be started
#' @param lowerLeftCorner integer that is only useful when \code{add}=TRUE
#' @param size Used internally by locationplot. Rescales each small barplot
#' window
#' @param addHorizontalLine adds a horizontal line that marks the default
#' fraction of 0.5 - 0.5
#' @param add.frame \code{boolean} to give the new plot a frame or not
#' @param filter.pValue.fraction \code{numeric} between 0 and 1 that filter
#' away pValues where the main allele has this frequency.
#' @param legend.fill.size size of the fill/boxes in the legend (default:NULL)
#' @param verbose Makes function more talkative
#' @param top.fraction.criteria 'maxcount', 'ref' or 'phase'
#' @param ... for simpler generics when extending function
#' @author Jesper R. Gadin, Lasse Folkersen
#' @seealso \itemize{ \item The \code{\link{ASEset}} class which the barplot
#' function can be called up on. }
#' @keywords barplot
#' @examples
#'
#' data(ASEset)
#' barplot(ASEset[1])
#'
#' @importFrom graphics plot
#' @importFrom graphics barplot
#' @exportMethod barplot
NULL
#' @rdname ASEset-barplot
setGeneric("barplot")
#' @rdname ASEset-barplot
setMethod("barplot", signature(height = "ASEset"), function(height, type = "count",
sampleColour.top = NULL, sampleColour.bot = NULL, legend = TRUE, pValue = TRUE, strand = "*", testValue = NULL,
testValue2 = NULL, OrgDb = NULL, TxDb = NULL, annotationType = c("gene", "exon",
"transcript"), main = NULL, ylim = NULL, yaxis = TRUE, xaxis = FALSE, ylab = TRUE,
ylab.text=NULL, xlab.text="samples",
xlab = TRUE, legend.colnames = "", las.ylab = 1, las.xlab = 2, cex.main = 1,
cex.pValue = 0.7, cex.ylab = 0.7, cex.xlab = 0.7, cex.legend = 0.6, add = FALSE,
lowerLeftCorner = c(0, 0), size = c(1, 1), addHorizontalLine = 0.5, add.frame = TRUE,
filter.pValue.fraction = 0.99, legend.fill.size=1, legend.interspace=1, verbose = FALSE,
top.fraction.criteria="maxcount", cex.annotation=0.7, ypos.annotation=0, annotation.interspace=1,
...) {
# catch useful graphical parameters that can be used to later add onto plot. This
# list will be retireved by using 'glst <- barplot(x)'
graphParamList <- list()
x <- height
# check type
okPlotTypes <- c("fraction", "count")
if (!type %in% okPlotTypes)
stop(paste("type can't be '", type, "' - it should be one of these: ", paste(okPlotTypes,
collapse = ", "), sep = ""))
# check verbose
if (class(verbose)[1] != "logical")
stop(paste("verbose should be of class logical, not", class(verbose)[1]))
if (length(verbose) != 1)
stop(paste("verbose should be of length 1, not", length(verbose)))
# check strand
okStrandTypes <- c("both", "+", "-", "*")
if (!strand %in% okStrandTypes)
stop(paste("strand can't be '", strand, "' - it should be one of these: ",
paste(okStrandTypes, collapse = ", "), sep = ""))
# check if strand type is present in object
# there is a check at the alleleCountsList aswell.
# And since the new ASEset version the strands will always be present
#
#if (strand == "+") {
# el <- "countsPlus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else if (strand == "-") {
# el <- "countsMinus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else if (strand == "*") {
# el <- "countsUnknown"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else if (strand == "both") {
# el <- "countsPlus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
# el <- "countsMinus"
# if (!(el %in% assayNames(x))) {
# stop("strand is not present as assay in ASEset object")
# }
#} else {
# stop("unknown strand option")
#}
# check legend
if (!is.logical(legend)) {
stop("legend has to be logical, TRUE or FALSE")
}
if (!length(legend) == 1) {
stop("legend has to be of length 1")
}
# check pValue
if (!is.logical(pValue)) {
stop("pValue has to be logical, TRUE or FALSE")
}
if (!length(pValue) == 1) {
stop("pValue has to be of length 1")
}
# check ylim
if (!is.null(ylim)) {
}
# check testValue
if (!is.null(testValue)) {
if (!class(testValue)[1] == "numeric" & !class(testValue)[1] == "matrix") {
stop("wrong class for testValue")
}
if (!class(testValue)[1] == "matrix") {
if (!sum(dim(testValue) == rev(dim(x))) == 2) {
stop("wrong dimensions of testValue")
}
}
if (!class(testValue)[1] == "numeric") {
if (!length(testValue) == ncol(x)) {
stop("wrong dimensions of testValue")
}
}
}
# check testValue2
if (!is.null(testValue2)) {
if (!class(testValue2)[1] == "numeric" & !class(testValue2)[1] == "matrix") {
stop("wrong class for testValue")
}
if (!class(testValue2)[1] == "matrix") {
if (!sum(dim(testValue2) == rev(dim(x))) == 2) {
stop("wrong dimensions of testValue")
}
}
if (!class(testValue2)[1] == "numeric") {
if (!length(testValue2) == ncol(x)) {
stop("wrong dimensions of testValue")
}
}
}
# check testValue combinations
if (!strand == "both" & !is.null(testValue2)) {
stop("testValue2 is supposed to only be \n\t\t\tspecified when using strand='both', and then be the values for the minus strand.")
}
# check ylim
if (!is.null(ylim)) {
if (type == "fraction") {
warning("an ylim argument given to barplot with type fraction is ignored")
ylim <- NULL
} else {
if (class(ylim)[1] != "numeric")
stop(paste("if given, ylim must be of class numeric, not", class(ylim)[1]))
if (length(ylim) != 2)
stop(paste("if given, ylim must be of length 2, not", length(ylim)))
if (ylim[2] < ylim[1])
stop(paste("if ylim is given, ylim[2] must be larger than ylim[1]"))
}
}
# set easy accessor to samples
samples <- colnames(x)
# set easy accessor to snps
snps <- rownames(x)
# check sampleColour
if (is.null(sampleColour.top)) {
sampleColour.top <- rep("darkolivegreen4", length(samples))
}
if (is.null(sampleColour.bot)) {
sampleColour.bot <- rep("brown2", length(samples))
}
# if (is.null(sampleColour)) {
# sampleColour <- rep("dodgerblue", length(samples))
# } else {
# if (class(sampleColour) != "character")
# stop(paste("if given, sampleColour should be of class character, not",
# class(sampleColour)))
# if (length(sampleColour) == 1) {
# sampleColour <- rep(sampleColour, length(samples))
# } else {
# if (length(sampleColour) != length(samples))
# stop(paste("if given, sampleColour should be of length 1 or length",
# length(samples), "(corresponding to the number of samples)"))
# }
# hexadecimals <- sampleColour[grep("^#[0-9A-F]{6}", sampleColour)]
# missing <- sampleColour[!(sampleColour %in% colors() | sampleColour %in%
# hexadecimals)]
# if (length(missing) > 0) {
# stop(paste("The following", length(unique(missing)), "colour(s) from sampleColour were not recognised:",
# paste(sort(unique(missing)), collapse = ", ")))
# }
# }
# check OrgDb
if (!is.null(OrgDb)) {
if (!class(OrgDb)[1] == "OrgDb") {
stop("class of OrgDb has to be an OrgDb class")
}
}
if (!is.null(TxDb)) {
if (!class(TxDb)[1] == "TxDb") {
stop("class of TxDb has to be a TxDb class")
}
}
# check annotationType
okAnnotationTypes <- c("gene", "exon", "cds", "transcript")
if (!sum(annotationType %in% okAnnotationTypes) == length(annotationType)) {
stop(paste("type can't be '", annotationType, "' - it should be one of these: ",
paste(okAnnotationTypes, collapse = ", "), sep = ""))
}
# check main
if (!is.null(main)) {
if (!class(main)[1] == "character") {
stop("main has to be of class character")
}
}
# check yaxis
if (!is.logical(yaxis)) {
stop("yaxis has to be logical, TRUE or FALSE")
}
if (!length(yaxis) == 1) {
stop("yaxis has to be of length 1")
}
# check xaxis
if (!is.logical(xaxis)) {
stop("xaxis has to be logical, TRUE or FALSE")
}
if (!length(xaxis) == 1) {
stop("xaxis has to be of length 1")
}
# check ylab
if (!is.logical(ylab)) {
stop("ylab has to be logical, TRUE or FALSE")
}
if (!length(ylab) == 1) {
stop("ylab has to be of length 1")
}
# check xlab
if (!is.logical(xlab)) {
stop("xlab has to be logical, TRUE or FALSE")
}
if (!length(xlab) == 1) {
stop("xlab has to be of length 1")
}
# check legend.colnames
if (!class(legend.colnames)[1] == "character") {
stop("legend.colnames has to be of class character")
}
# check las.ylab
if (!class(las.ylab)[1] == "numeric" & !class(las.ylab)[1] == "integer") {
stop("las.ylab has to be of class numeric or integer")
}
if (!las.ylab == 1 & !las.ylab == 2) {
stop("las.ylab has to be a number, either 1 or 2")
}
# check las.xlab
if (!class(las.xlab)[1] == "numeric" & !class(las.xlab)[1] == "integer") {
stop("las.xlab has to be of class numeric or integer")
}
if (!las.xlab == 1 & !las.xlab == 2) {
stop("las.xlab has to be a number, either 1 or 2")
}
# check cex all of them
if (!class(cex.main)[1] == "numeric") {
if (cex.main <= 0 | !length(cex.main) == 1)
stop("cex.main has to be of class numeric, and have length 1, and a value above 0 ")
}
if (!class(cex.pValue)[1] == "numeric") {
if (cex.pValue > 1 | cex.pValue <= 0 | !length(cex.pValue) == 1)
stop("cex.pValue has to be of class numeric, and have length 1, and a value above 0")
}
if (!class(cex.ylab)[1] == "numeric") {
if (cex.ylab <= 0 | !length(cex.ylab) == 1)
stop("cex.ylab has to be of class numeric, and have length 1, and a value above 0 ")
}
if (!class(cex.xlab)[1] == "numeric") {
if (cex.xlab <= 0 | !length(cex.xlab) == 1)
stop("cex.xlab has to be of class numeric, and have length 1, and a value above 0 ")
}
if (!class(cex.legend)[1] == "numeric") {
if (cex.legend <= 0 | !length(cex.legend) == 1)
stop("cex.legend has to be of class numeric, and have length 1, and a value above 0 ")
}
# check add
if (!is.logical(add)) {
stop("add has to be logical, TRUE or FALSE")
}
if (!length(add) == 1) {
stop("add has to be of length 1")
}
# check lowerLeftCorner
if (!class(lowerLeftCorner)[1] == "numeric" & !class(lowerLeftCorner)[1] == "integer") {
stop("lowerLeftCorner has to be of class numeric or integer")
}
if (!length(lowerLeftCorner) == 2) {
stop("lowerLeftCorner has be of length 2")
}
# check size
if (!class(size)[1] == "numeric" & !class(size)[1] == "integer") {
stop("size has to be of class numeric or integer")
}
if (!length(size) == 2) {
stop("size has be of length 2")
}
# check addHorizontalLine
if (!is.null(addHorizontalLine)) {
if (!class(addHorizontalLine)[1] == "numeric") {
stop("addHorizontalLine has to be either NULL or of class numeric")
}
if (!length(addHorizontalLine) == 1) {
stop("addHorizontalLine has be of length 1")
}
if (addHorizontalLine > 1 | addHorizontalLine <= 0) {
stop("addHorizontalLine has be numeric between 0 and 1")
}
}
# check add.frame
if (!is.logical(add.frame)) {
stop("add.frame has to be logical, TRUE or FALSE")
}
if (!length(add.frame) == 1) {
stop("add.frame has to be of length 1")
}
# check filter.pValue.fraction=0.99
if (!class(filter.pValue.fraction)[1] == "numeric") {
stop("filter.pValue.fraction has to be of class numeric")
}
if (!length(filter.pValue.fraction) == 1) {
stop("filter.pValue.fraction has to be of length 1")
}
if (filter.pValue.fraction < 0 | filter.pValue.fraction > 1) {
stop("filter.pValue.fraction has to be a value between 0 and 1")
}
# make annotation dataframes
if ((!is.null(OrgDb) | !is.null(TxDb)) & (strand == "+" |
strand == "*" | strand == "both")) {
annDfPlus <- getAnnotationDataFrame(rowRanges(x), strand = "+", annotationType = annotationType,
OrgDb = OrgDb, TxDb = TxDb, verbose = verbose)
} else {
annDfPlus <- NULL
}
if ((!is.null(OrgDb) | !is.null(TxDb)) & (strand == "-" |
strand == "*" | strand == "both")) {
annDfMinus <- getAnnotationDataFrame(rowRanges(x), strand = "-", annotationType = annotationType,
OrgDb = OrgDb, TxDb = TxDb, verbose = verbose)
} else {
annDfMinus <- NULL
}
# creating a background from foreground colour - basically just turning down the
# saturation a little compared to the foreground color except for very
# low-saturated fg's which are made more saturated, and very dark fg's which are
# made lighter
#hsvColoursBg <- hsvColoursFg <- rgb2hsv(col2rgb(sampleColour))
#veryLowSaturation <- hsvColoursFg["s", ] < 0.3
#veryDark <- hsvColoursFg["v", ] < 0.2
#remaining <- !veryLowSaturation & !veryDark
#hsvColoursBg["s", veryLowSaturation] <- hsvColoursFg["s", veryLowSaturation] +
# 0.3
#hsvColoursBg["v", veryDark] <- hsvColoursFg["v", veryDark] + 0.2
#hsvColoursBg["s", remaining] <- hsvColoursFg["s", remaining] - 0.4
#fgCol <- hsv(hsvColoursFg["h", ], hsvColoursFg["s", ], hsvColoursFg["v", ])
#bgCol <- hsv(hsvColoursBg["h", ], hsvColoursBg["s", ], hsvColoursBg["v", ])
#fgCol <- hsv(rgb2hsv(col2rgb(sampleColour.bot))["h",])
#bgCol <- hsv(rgb2hsv(col2rgb(sampleColour.top))["h",])
fgCol <- sampleColour.bot
bgCol <- sampleColour.top
if (type == "count") {
# the intention of this is to be able to create barplots which show both the
# amount of counts and the distribution between alleles this is the least
# digested format and consequently produces the most cluttered plots. It is
# therefore intended for subsets, few snps, few samples etc in these smaller
# subsets it is useful to get information on higher expression values at certain
# points of genes etc.
if (strand == "+" | strand == "-" | strand == "*" ) {
for (i in 1:length(snps)) {
# getting revelant data
snp <- snps[i]
tmp <- t(alleleCounts(x, strand)[[snp]][samples, , drop = FALSE])
if ((ncol(tmp) * nrow(tmp)) == sum(tmp == 0)) {
zeroRows <- TRUE
} else {
zeroRows <- FALSE
}
# alleles are
n <- names(sort(apply(tmp, 1, sum), decreasing = TRUE)[1:2])
# theor total count
tc <- sort(apply(tmp, 1, sum), decreasing = TRUE)[1:2]
# extracting only major and minor allele in correctly transposed format and
# ordering colours accordingly
height <- tmp[match(n, rownames(tmp)), , drop = FALSE]
# change all zero rows to X names
rownames(height)[apply(height, 1, sum) == 0] <- "X"
# warn if the two remaining alleles have too many counts
majorAndMinorFraction <- sum(height)/sum(tmp)
if (verbose & majorAndMinorFraction < 0.9) {
message(paste("Snp", snp, "was possible tri-allelic, but only two most frequent alleles were plotted. Counts:\n",
paste(paste(names(apply(tmp, 1, sum)), apply(tmp, 1, sum), sep = "="),
collapse = ", ")))
}
## make bars section
extraspaceOnTop <- 1.15
if (is.null(ylim)) {
if (!zeroRows) {
ylimHere <- c(0, max(abs(range(height, na.rm = TRUE))) * extraspaceOnTop)
} else {
ylimHere <- c(0, 10)
}
} else {
ylimHere <- ylim
}
# init frame
if (!add) {
size <- c(1, 1)
lowerLeftCorner <- c(0, 0)
plot.default(NULL, xlim = c(0, 1), ylim = c(0, 1), ylab = "", xlab = "",
xaxt = "n", yaxt = "n", frame.plot = FALSE)
}
# set so the ylimHere is at least 10.
if (ylimHere[2] < 10) {
ylimHere[2] <- 10
}
# make tics and labels prettier
ylabels <- pretty(ylimHere[1]:ylimHere[2])
if (strand == "both") {
ylimHere[1] <- -max(abs(ylabels))
ylimHere[2] <- max(abs(ylabels))
} else {
ylimHere[1] <- min(ylabels)
ylimHere[2] <- max(ylabels)
}
# all x positions
xPointsAll <- seq(lowerLeftCorner[1], lowerLeftCorner[1] + size[1],
length.out = length(height) * 1.5 + 2)[2:(length(height) * 1.5 +
1)]
widthsOfRectangles <- rep(size[1]/(length(height) * 1.5 + 1), length(height)/2)
# normalize toPlot-values to a range of c(0,size[2])
heightsOfRectangles <- ((height - ylimHere[1])/(ylimHere[2] - ylimHere[1]) *
size[2])
# take every third element (start at 2)
xPoints1 <- xPointsAll[seq(1, length(height) * 1.5, by = 3)]
yPoints <- lowerLeftCorner[2] + heightsOfRectangles[1, ]/2
# alleleColors1 <- rep(fgCol[1],9)
alleleColors1 <- fgCol
symbols(x = xPoints1, y = yPoints, bg = alleleColors1, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectangles[1, ])), add = TRUE, inches = FALSE)
# take every third element (start at 3)
xPoints2 <- xPointsAll[seq(2, length(height) * 1.5, by = 3)]
yPoints <- lowerLeftCorner[2] + heightsOfRectangles[2, ]/2
alleleColors2 <- rep(bgCol[1], 9)
alleleColors2 <- bgCol
widthsOfRectangles <- rep(size[1]/(length(height) * 1.5 + 1), length(height)/2)
symbols(x = xPoints2, y = yPoints, bg = alleleColors2, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectangles[2, ])), add = TRUE, inches = FALSE)
## make axis and labels section
# y-axis
yat <- (ylabels + abs(min(ylabels)))/(sum(abs(min(ylabels)) + max(ylabels)))
graphParamList[["yat"]] <- yat
if (yaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1]), c(lowerLeftCorner[2],
lowerLeftCorner[2] + 1))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x - size[1] * 0.01), c(y, y))
}, rep(lowerLeftCorner[1], length(yat)), (lowerLeftCorner[2] +
yat)))
}
if (ylab) {
if (las.ylab == 2) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.ylab)
}
if (las.ylab == 1) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.ylab)
}
}
# x-axis labels
xat <- xPoints1 + (xPoints2 - xPoints1)/2
graphParamList[["xat"]] <- xat
if (xaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1] + 1), c(lowerLeftCorner[2],
lowerLeftCorner[2]))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x), c(y, y - size[1] * 0.01))
}, (lowerLeftCorner[2] + xat), rep(lowerLeftCorner[1], length(xat))))
}
if (xlab) {
if (las.xlab == 2) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.xlab)
}
if (las.xlab == 1) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.xlab)
}
}
if (is.null(main)) {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = paste(snp), cex = cex.main,
adj = 0.5, xpd = TRUE)
} else {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = main[i], cex = cex.main, adj = 0.5)
}
if (pValue) {
# check pValue slot
if (!is.null(testValue)) {
if (class(testValue)[1] == "matrix") {
text(xat, apply(heightsOfRectangles, 2, max), signif(testValue[,
i], 1), cex = cex.pValue, adj = c(0.5, 0.01))
} else {
text(xat, apply(heightsOfRectangles, 2, max), signif(testValue,
1), cex = cex.pValue, adj = c(0.5, 0.01))
}
} else {
# use chisq by default
TFf <- t(fraction(x[i], strand = strand) < filter.pValue.fraction)
ct <- signif(chisq.test(x[i], strand), 1)
ct[!TFf] <- NA
text(xat, apply(heightsOfRectangles, 2, max), ct, cex = cex.pValue,
adj = c(0.5, 0.01))
}
}
if (!is.null(OrgDb) | !is.null(TxDb)) {
annotationBarplot(strand,i, lowerLeftCorner, annDfPlus, annDfMinus,
cex=cex.annotation, ypos=ypos.annotation,
interspace=annotation.interspace)
}
if (legend) {
# any strand
o <- apply(height, 1, sum)
# remove zeros count alleles
TFz <- rownames(height) == "X"
o <- o[!TFz]
textOver <- rownames(height)[!TFz]
# set another xLegendPos and yLegendPos for '*'.
if ((strand == "*" ) & (!is.null(OrgDb) |
!is.null(TxDb))) {
xLegendPos <- 0.5
yLegendPos <- 1.00
} else {
xLegendPos <- 0.94
yLegendPos <- 1.00
}
legendBarplot(lowerLeftCorner, size, rownames=textOver, colnames=legend.colnames,
boxsize=legend.fill.size, boxspace=legend.interspace, fgCol, bgCol,
ylegendPos=1, xlegendPos=0.96, cex=cex.legend)
}
#set y and x labs
if(is.null(ylab.text)){
ylab.text <- "reads"
}
title(ylab=ylab.text, xlab=xlab.text)
}
} else if (strand == "both") {
for (i in 1:length(snps)) {
if (verbose) {
cat("plotting snp", i, "\n")
}
snp <- snps[i]
# getting revelant data forward strand
tmpPlus <- t(alleleCounts(x, "+")[[snp]][samples, ])
if ((ncol(tmpPlus) * nrow(tmpPlus)) == sum(tmpPlus == 0)) {
zeroRowsPlus <- TRUE
} else {
zeroRowsPlus <- FALSE
}
tmpMinus <- t(alleleCounts(x, "-")[[snp]][samples, ])
if ((ncol(tmpMinus) * nrow(tmpMinus)) == sum(tmpMinus == 0)) {
zeroRowsMinus <- TRUE
} else {
zeroRowsMinus <- FALSE
}
# check most common alleles on + and - (should be same)
allVar <- data.frame(A = c(NA, NA), T = c(NA, NA), G = c(NA, NA),
C = c(NA, NA), del = c(NA, NA))
o <- apply(tmpPlus, 1, sum)
u <- apply(tmpMinus, 1, sum)
# to get their total count
allVar[1, match(names(o), colnames(allVar))] <- o
allVar[2, match(names(u), colnames(allVar))] <- u
# save to tri-allelic warner
nAll <- names(sort(apply(allVar, 2, sum), decreasing = TRUE))
tcAll <- sort(apply(allVar, 2, sum), decreasing = TRUE)
# set colnames to X for non present alleles
xNames <- colnames(allVar)[apply(allVar, 2, sum) == 0]
colnames(allVar)[colnames(allVar) %in% xNames] <- "X"
xNames <- rownames(tmpPlus)[apply(tmpPlus, 1, sum) == 0]
rownames(tmpPlus)[rownames(tmpPlus) %in% xNames] <- "X"
xNames <- rownames(tmpMinus)[apply(tmpMinus, 1, sum) == 0]
rownames(tmpMinus)[rownames(tmpMinus) %in% xNames] <- "X"
# alleles are
n <- names(sort(apply(allVar, 2, sum), decreasing = TRUE))[1:2]
# theor total count
tc <- sort(apply(allVar, 2, sum), decreasing = TRUE)[1:2]
# extracting only major and minor allele in correctly transposed format and
# ordering colours accordingly
over <- tmpPlus[match(n, rownames(tmpPlus)), ]
under <- tmpMinus[match(n, rownames(tmpMinus)), ]
# change all NAs to X names and zero values
over[is.na(rownames(over)), ] <- 0
rownames(over)[is.na(rownames(over))] <- "X"
under[is.na(rownames(under)), ] <- 0
rownames(under)[is.na(rownames(under))] <- "X"
# Plus Warn if the two remaining alleles have too many counts
if (!zeroRowsPlus) {
majorAndMinorFraction <- sum(allVar[1, ][c(n[1], n[2])])/sum(allVar[1,
])
if (verbose & majorAndMinorFraction < 0.9) {
message(paste("Snp", snp, "was possible tri-allelic, but only two most frequent alleles were plotted. Counts:\n",
paste(paste(nAll, tcAll, sep = "="), collapse = ", ")))
}
}
# Minus calculating major and minor allele, warn if the two remaining alleles
# have too many counts
if (!zeroRowsMinus) {
majorAndMinorFraction <- sum(allVar[2, ][c(n[1], n[2])])/sum(allVar[2,
])
if (verbose & majorAndMinorFraction < 0.9) {
message(paste("Snp", snp, "was possible tri-allelic, but only two most frequent alleles were plotted. Counts:\n",
paste(paste(nAll, tcAll, sep = "="), collapse = ", ")))
}
}
# set under as minus
under <- -under
# plotting function
if (!add) {
# init frame
plot.default(NULL, xlim = c(0, 1), ylim = c(0, 1), ylab = "", xlab = "",
xaxt = "n", yaxt = "n", , frame.plot = FALSE)
}
heightPlus <- over
heightMinus <- under
if (is.null(ylim)) {
if (zeroRowsPlus) {
ymax <- 5
} else {
ymax <- max(abs(range(heightPlus, na.rm = TRUE)))
}
if (zeroRowsMinus) {
ymin <- -5
} else {
ymin <- min(heightMinus)
}
ylimHere <- c(ymin, ymax)
} else {
ylimHere <- ylim
}
# increase ylimHere before making pretty
extraspaceOnTopAndBottom <- 0.15
extraRealDistance <- (dist(ylimHere) * extraspaceOnTopAndBottom)
ylimHere[1] <- ylimHere[1] - extraRealDistance
ylimHere[2] <- ylimHere[2] + extraRealDistance
# set lowest max and min values
if (ylimHere[2] < 10) {
ylimHere[2] <- 10
}
if (abs(ylimHere[1]) < 10) {
ylimHere[1] <- -10
}
# make more even limits (make more pretty)
ylabels <- pretty(ylimHere[1]:ylimHere[2])
ylimHere[1] <- min(ylabels)
ylimHere[2] <- max(ylabels)
# breakpoint Plus Minus bars
middleLine <- (1/(ylimHere[2] - ylimHere[1])) * abs(ylimHere[1])
# make a xPointsAll framework
xPointsAll <- seq(lowerLeftCorner[1], lowerLeftCorner[1] + size[1],
length.out = length(heightPlus) * 1.5 + 2)[2:(length(heightPlus) *
1.5 + 1)]
widthsOfRectangles <- rep(size[1]/(length(heightPlus) * 1.5 + 1),
length(heightPlus)/2)
# normalize toPlot-values to a range of c(0,size[2])
heightsOfRectanglesPlus <- round(((heightPlus - ylimHere[1])/(ylimHere[2] -
ylimHere[1]) * size[2]) - middleLine, 5)
heightsOfRectanglesMinus <- round(((-heightMinus - ylimHere[1])/(ylimHere[2] -
ylimHere[1]) * size[2]) - middleLine, 5)
# take every third element (start at 2)
xPoints1 <- xPointsAll[seq(2, length(heightPlus) * 1.5, by = 3)]
yPointsPlus <- lowerLeftCorner[2] + middleLine + heightsOfRectanglesPlus[1,
]/2
barMinusEdge <- matrix(rep(NA, length(heightsOfRectanglesMinus)),
nrow = 2) #is used later to plot p-value
# yPointMinus is a little tricky and needs a more complicated equation
yPointsMinus1 <- lowerLeftCorner[2] + heightsOfRectanglesMinus[1,
]/2
yPointsMinus2 <- 1 - yPointsMinus1
yPointsMinus3 <- yPointsMinus2 - (1 - middleLine)
barMinusEdge[1, ] <- (1 - (lowerLeftCorner[2] + heightsOfRectanglesMinus[1,
])) - (1 - middleLine)
yPointsMinus <- yPointsMinus3
alleleColors1 <- fgCol
symbols(x = xPoints1, y = yPointsPlus, bg = alleleColors1, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesPlus[1, ])), add = TRUE,
inches = FALSE)
symbols(x = xPoints1, y = yPointsMinus, bg = alleleColors1, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesMinus[1, ])), add = TRUE,
inches = FALSE)
# take every third element (start at 3)
xPoints2 <- xPointsAll[seq(3, length(heightPlus) * 1.5, by = 3)]
yPointsPlus <- lowerLeftCorner[2] + middleLine + heightsOfRectanglesPlus[2,
]/2
# yPointMinus is a little tricky and needs a more complicated equation
yPointsMinus1 <- lowerLeftCorner[2] + heightsOfRectanglesMinus[2,
]/2
yPointsMinus2 <- 1 - yPointsMinus1
yPointsMinus3 <- yPointsMinus2 - (1 - middleLine)
barMinusEdge[2, ] <- (1 - (lowerLeftCorner[2] + heightsOfRectanglesMinus[2,
])) - (1 - middleLine)
yPointsMinus <- yPointsMinus3
# alleleColors2 <- rep('#7F5959',9)
alleleColors2 <- bgCol
symbols(x = xPoints2, y = yPointsPlus, bg = alleleColors2, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesPlus[2, ])), add = TRUE,
inches = FALSE)
symbols(x = xPoints2, y = yPointsMinus, bg = alleleColors2, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectanglesMinus[2, ])), add = TRUE,
inches = FALSE)
# y-axis
yat <- (ylabels + abs(min(ylabels)))/(sum(abs(min(ylabels)) + max(ylabels)))
graphParamList[["yat"]] <- yat
if (yaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1]), c(lowerLeftCorner[2],
lowerLeftCorner[2] + 1))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x - size[1] * 0.01), c(y, y))
}, rep(lowerLeftCorner[1], length(yat)), (lowerLeftCorner[2] +
yat)))
}
if (ylab) {
if (las.ylab == 2) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.ylab)
}
if (las.ylab == 1) {
text(lowerLeftCorner[1] - size[1] * 0.02, yat, labels = abs(ylabels),
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.ylab)
}
}
# x-axis labels
xat <- xPoints1 + (xPoints2 - xPoints1)/2
graphParamList[["xat"]] <- xat
if (xaxis) {
# vertical line
lines(c(lowerLeftCorner[1], lowerLeftCorner[1] + 1), c(lowerLeftCorner[2],
lowerLeftCorner[2]))
# horizontal lines
invisible(mapply(FUN = function(x, y) {
lines(c(x, x), c(y, y - size[1] * 0.01))
}, (lowerLeftCorner[2] + xat), rep(lowerLeftCorner[1], length(xat))))
}
if (xlab) {
if (las.xlab == 2) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 90, adj = c(1, 0.5), cex = cex.xlab)
}
if (las.xlab == 1) {
text(xat, lowerLeftCorner[2] - size[2] * 0.02, labels = samples,
xpd = TRUE, srt = 0, adj = c(1, 0.5), cex = cex.xlab)
}
}
if (is.null(main)) {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = paste(snp), cex = cex.main,
adj = 0.5, xpd = TRUE)
} else {
text(x = lowerLeftCorner[1] + size[1] * 0.5, y = lowerLeftCorner[2] +
size[2] * 1.05, labels = main[i], cex = cex.main, adj = 0.5)
}
if (pValue) {
# check pValue slot
if (!is.null(testValue)) {
if (class(testValue)[1] == "matrix") {
text(xat, apply(middleLine + heightsOfRectanglesPlus, 2, max),
signif(testValue[, i], 1), cex = cex.pValue, adj = c(0.4,
-0.3))
} else {
text(xat, apply(middleLine + heightsOfRectanglesPlus, 2, max),
signif(testValue, 1), cex = cex.pValue, adj = c(0.4, -0.3))
}
} else {
# use chisq by default
TFf <- t(fraction(x[i], strand = "+") < filter.pValue.fraction)
ct <- signif(chisq.test(x[i], "+"), 1)
ct[!TFf] <- NA
text(xat, apply(middleLine + heightsOfRectanglesPlus, 2, max),
ct, cex = cex.pValue, adj = c(0.4, -0.3))
}
if (!is.null(testValue2)) {
if (class(testValue)[1] == "matrix") {
text(xat, apply(barMinusEdge, 2, min), signif(testValue2[,
i], 1), cex = cex.pValue, adj = c(0.4, 1, 1))
} else {
text(xat, apply(barMinusEdge, 2, min), signif(testValue2, 1),
cex = cex.pValue, adj = c(0.4, 1, 1))
}
} else {
# use chisq by default
TFf <- t(fraction(x[i], strand = "-") < filter.pValue.fraction)
ct <- signif(chisq.test(x[i], "-"), 1)
ct[!TFf] <- NA
text(xat, apply(barMinusEdge, 2, min), ct, cex = cex.pValue,
adj = c(0.4, 1, 1))
}
}
if (!is.null(OrgDb) | !is.null(TxDb)) {
annotationBarplot(strand,i, lowerLeftCorner, annDfPlus, annDfMinus,
cex=cex.annotation, ypos=ypos.annotation,
interspace=annotation.interspace)
}
if (legend) {
# plus strand
TFz <- rownames(over) == "X"
if (!sum(!TFz) == 0) {
textOver <- paste(rownames(over)[!TFz], apply(over, 1, sum)[!TFz])
#remove the total count for that nucleotide
textOver <- unlist(lapply(strsplit(textOver," "),function(x){x[1]}))
legendBarplot(lowerLeftCorner, size, rownames=textOver, colnames=legend.colnames,
boxsize=legend.fill.size, boxspace=legend.interspace, fgCol, bgCol,
ylegendPos=1, xlegendPos=0.96, cex=cex.legend)
}
}
}
#set y and x labs
if(is.null(ylab.text)){
ylab.text <- "reads"
}
title(ylab=ylab.text, xlab=xlab.text)
}
} else if (type == "fraction") {
# a plot type that shows the allelic imbalance, irrespective of absolute read
# counts at each Snp. this produce plots that are more easily overviewed than
# those from 'plot_reads_by_allele', but of course omits some detail information
# from these plots. Intended to use for a medium number of snps and samples
if (strand == "both") {
stop("strand must be '+', '-' or '*' for type='fraction' ")
}
for (i in 1:length(snps)) {
# getting revelant data
snp <- snps[i]
fgColHere <- fgCol
#tmp <- alleleCounts(x, strand)[[snp]][samples, , drop = FALSE]
df <- fractionPlotDf(x, snp, strand=strand, top.fraction.criteria=top.fraction.criteria)
# check for zeroRows
if (sum(df$na=="no")==0) {
zeroRows <- TRUE
} else {
zeroRows <- FALSE
}
# setting no-count samples to colour grey90 and fraction 0
#fgColHere[df$na[seq(1,nrow(df),by=2)+1]=="yes"] <- "grey90"
bgCol[df$na[seq(1,nrow(df),by=2)+1]=="yes"] <- "grey90"
#fraction[is.nan(fraction)] <- 0
# if add=FALSE then make a new device
if (!add) {
plot.default(NULL, xlim = c(0, 1), ylim = c(0, 1), ylab = "", xlab = "",
xaxt = "n", yaxt = "n")
}
# ylim is always (0,1) and usergiven ylims are reset to this with a warning (in
# start checkups).
ylim <- c(0, 1)
# find the centers of rectangles
xPoints <- seq(lowerLeftCorner[1], lowerLeftCorner[1] + size[1], length.out = nrow(df)/2 +
2)[2:(nrow(df)/2 + 1)]
widthsOfRectangles <- rep(size[1]/(nrow(df)/2 + 1), nrow(df)/2)
# normalize fraction-values to a range of c(0,size[2])
heightsOfRectangles <- ((df$values[seq(1,nrow(df),by=2)] - ylim[1])/(ylim[2] - ylim[1])) * size[2]
# find the Y-coordinate center points for rectangles
yPoints <- lowerLeftCorner[2] + heightsOfRectangles/2
# if bgCol is given, start by painting background
if (!is.null(bgCol)) {
yPointsBg <- rep(lowerLeftCorner[2] + size[2]/2, length(xPoints))
rectanglesBg <- matrix(ncol = 2, c(widthsOfRectangles, rep(size[2],
length(widthsOfRectangles))))
symbols(x = xPoints, y = yPointsBg, bg = bgCol, rectangles = rectanglesBg,
add = TRUE, inches = FALSE)
}
symbols(x = xPoints, y = yPoints, bg = fgColHere, rectangles = matrix(ncol = 2,
c(widthsOfRectangles, heightsOfRectangles)), add = TRUE, inches = FALSE)
# add a frame to the plot
if (add.frame) {
symbols(x = lowerLeftCorner[1] + size[1]/2, y = lowerLeftCorner[2] +
size[2]/2, rectangles = matrix(ncol = 2, c(size[1], size[2])),
add = TRUE, inches = FALSE)
}
# add a tick markers (easier for fraction-plots because they are restricted to
# 0-100%
marks <- c(0.25, 0.5, 0.75)
names(marks) <- c("25%", "50%", "75%")
for (markName in names(marks)) {
mark <- marks[markName]
yBase <- lowerLeftCorner[2] + ((mark - ylim[1])/(ylim[2] - ylim[1])) *
size[2]
xBase <- xPoints[1] - widthsOfRectangles[1]/2
lines(x = c(xBase - size[1] * 0.01, xBase + size[1] * 0.01), y = c(yBase,
yBase))
text(x = xBase - size[1] * 0.015, y = yBase, label = markName, adj = 1,
cex = 0.7)
}
# writing the name of each Snp on top of the plot(within frame)
text(x = 0.5, y = 1.02, label = snp, xpd = TRUE, cex = cex.main)
# add horizontal line at 0.5
if (!is.null(addHorizontalLine)) {
lines(x = c(lowerLeftCorner[1], lowerLeftCorner[1] + size[1]), y = c(lowerLeftCorner[2] +
size[2] * addHorizontalLine, lowerLeftCorner[2] + size[2] * addHorizontalLine))
mapBiasHere <- mapBias(x)[[snp]]
#if (!all(mapBiasHere %in% c(0.5, 0)) & !zeroRows) {
# # Only activate if anything is different from 0.5 or 0
#
# for (j in 1:length(fraction)) { #fraction doesnt exist (check df$values instead)
# biasfraction <- mapBiasHere[j, majorAllele]/(mapBiasHere[j, majorAllele] +
# mapBiasHere[j, minorAllele])
# yPoint <- lowerLeftCorner[2] + ((biasfraction - ylim[1])/(ylim[2] -
# ylim[1])) * size[2]
# lines(x = c(xPoints[j] - widthsOfRectangles[j]/2, xPoints[j] +
# widthsOfRectangles[j]/2), y = c(yPoint, yPoint), lty = 2)
# }
#}
}
#set y and x labs
if(is.null(ylab.text)){
ylab.text <- "fraction"
}
title(ylab=ylab.text, xlab=xlab.text)
}
}
# return graphical paramlist (does not work for fraction type yet)
invisible(graphParamList)
})
#' gbarplot ASEset objects
#'
#' Generates gbarplots for ASEset objects. Two levels of plotting detail are
#' provided: a detailed gbarplot of read counts by allele useful for fewer
#' samples and SNPs, and a less detailed gbarplot of the fraction of imbalance,
#' useful for more samples and SNPs.
#'
#' This function serves the same purpose as the normal barplot, but with
#' trellis graphics using lattice, to be able to integrate well with Gviz track
#' functionality.
#'
#' @name ASEset-gbarplot
#' @aliases ASEset-gbarplot gbarplot gbarplot,ASEset-method
#' @docType methods
#' @param x An \code{ASEset} object
#' @param type 'count' or 'fraction'
#' @param strand four options, '+', '-', 'both' or '*'
#' @param verbose Makes function more talkative
#' @param ... for simpler generics when extending function
#' @author Jesper R. Gadin
#' @seealso \itemize{ \item The \code{\link{ASEset}} class which the gbarplot
#' function can be called up on. \item The \code{\link{barplot}} non trellis
#' barplot. }
#' @keywords gbarplot
#' @examples
#'
#' data(ASEset)
#' gbarplot(ASEset[1])
#'
#' @exportMethod gbarplot
setGeneric("gbarplot", function(x, type = "count", strand = "*",
verbose = FALSE, ...) {
standardGeneric("gbarplot")
})
setMethod("gbarplot", signature(x = "ASEset"), function(x, type = "count", strand = "*",
usePhase=FALSE, verbose = FALSE, ...) {
if (length(list(...)) == 0) {
e <- new.env(hash = TRUE)
} else {
e <- list2env(list(...))
}
#print(ls(envir=e))
if (!exists("mainvec", envir = e, inherits = FALSE)) {
e$mainvec <- rep("",nrow(x))
}
if (!exists("cex.mainvec", envir = e, inherits = FALSE)) {
e$cex.mainvec <- 1
}
if (!exists("ylab", envir = e, inherits = FALSE)) {
e$ylab <- "reads"
}
if (!exists("xlab", envir = e, inherits = FALSE)) {
e$xlab <- "samples"
}
if (!exists("middleLine", envir = e, inherits = FALSE)) {
e$middleLine <- TRUE
}
if (!exists("deAnnoPlot", envir = e, inherits = FALSE)) {
e$deAnnoPlot <- FALSE
}
if (!exists("top.fraction.criteria", envir = e, inherits = FALSE)) {
e$top.fraction.criteria <- "maxcount"
}
for (i in 1:nrow(x)) {
name <- rownames(x)[i]
if (type == "fraction") {
b <- barplotLatticeFraction(
identifier = name,
x=x,
astrand=strand,
ids=rownames(x),
ylab=e$ylab,
xlab=e$xlab,
mainvec=e$mainvec,
middleLine=e$middleLine,
deAnnoPlot=e$deAnnoPlot,
cex.mainvec=e$cex.mainvec,
top.fraction.criteria=e$top.fraction.criteria)
} else if (type == "count") {
b <- barplotLatticeCounts(
identifier = name,
x=x,
astrand=strand,
ids=rownames(x),
ylab=e$ylab,
xlab=e$xlab,
mainvec=e$mainvec)
} else {
stop("type has to be fraction or count")
}
b
}
})
barplotLatticeFraction <- function(identifier, ...) {
if (length(list(...)) == 0) {
e <- new.env(hash = TRUE)
} else {
e <- list2env(list(...))
}
e$ids <- unlist(e$ids)
e$strand <- e$astrand
df <- fractionPlotDf(e$x, identifier, strand=e$strand, top.fraction.criteria=e$top.fraction.criteria)
# Grah params
my_cols <- c("green", "red")
# set default values
parset <- list(superpose.polygon=list(col=my_cols))
scales = list(rot = c(90, 0))
if (e$deAnnoPlot) {
parset <- list(
layout.widths = list(
left.padding = 0,
axis.left = 0,
ylab.axis.padding = 0,
right.padding = 0,
axis.right = 0
),
layout.heights = list(
top.padding = 0.1,
between = 0.1,
xlab.top= 0.1,
axis.top = 0,
main=1.1,
main.key.padding=1,
axis.xlab.padding = 1,
bottom.padding = 1,
axis.bottom = 0.3
),
superpose.polygon=list(col=my_cols)
)
scales = list(y = list(at = NULL, labels = NULL),
x = list(labels = rep("",ncol(e$x)))
)
#, rot = c(90, 0))
}
if(!e$middleLine) {
b <- barchart(values ~ sample, group = groups, data = df, origin = 0,
stack = TRUE, scales = scales,
main = list(label=e$main, cex=e$cex.mainvec),
ylab = e$ylab, xlab = e$xlab,
par.settings = parset)
}else {
# if middle line at 0.5
b <- barchart(values ~ sample, group = groups, data = df, origin = 0,
stack = TRUE, scales = scales,
main = list(label=e$main, cex=e$cex.mainvec),
ylab = e$ylab, xlab = e$xlab,
par.settings = parset, panel=function(x, y, ...) {
panel.barchart(x, y, ...)
panel.abline(h=0.5, lty=1)
} )
}
b
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.