## Add warning if less than 100 poihnts ## squeleton.pDensitites y plotSqueletonDensities are from ## limma, plotDensities, but to avoid repeating the same ## computations many times. squeleton.pDensities <- function(object, log.transform = FALSE, arrays = NULL, singlechannels = NULL, groups = NULL, col = NULL) { matDensities <- function(X) { densXY <- function(Z) { zd <- density(Z, na.rm = TRUE) x <- zd$x y <- zd$y cbind(x, y) } out <- apply(X, 2, densXY) outx <- out[(1:(nrow(out)/2)), ] outy <- out[(((nrow(out)/2) + 1):nrow(out)), ] list(X = outx, Y = outy) } if (is(object, "MAList")) { object <- RG.MA(object) } if (((is.null(arrays)) & (is.null(singlechannels)))) { arrays <- 1:(ncol(object$R)) x <- cbind(object$R, object$G) if (is.null(groups)) { groups <- c(length(arrays), length(arrays)) if (is.null(col)) colors <- rep(c("red", "green"), groups) if (!is.null(col)) { if (length(col) != 2) { print("number of groups=2 not equal to number of col") colors <- "black" } else { colors <- rep(col, groups) } } } else { if (!is.null(col)) { if (length(as.vector(table(groups))) != length(col)) { print("number of groups not equal to number of col") colors <- col } else { colors <- col[groups] } } else { print("warning no colors in col specified for the groups") colors <- "black" } } } else { if (!is.null(singlechannels)) { if (!is.null(arrays)) print("warning cant index using arrays AND singlechannels") x <- cbind(object$R, object$G)[, singlechannels] if (is.null(groups)) { groups <- c(length(intersect((1:ncol(object$R)), singlechannels)), length(intersect(((ncol(object$R) + 1):ncol(cbind(object$G, object$R))), singlechannels))) if (is.null(col)) colors <- rep(c("red", "green"), groups) if (!is.null(col)) { if (length(col) != 2) { print("number of groups=2 not equal to number of col") colors <- "black" } else { colors <- rep(col, groups) } } } else { if (!is.null(col)) { if (length(as.vector(table(groups))) != length(col)) { print("number of groups not equal to number of col") colors <- col } else { colors <- col[groups] } } else { print("warning no colors in col specified for the groups") colors <- "black" } } } else { if (!is.null(arrays)) { if (!is.null(singlechannels)) print("warning cant index using arrays AND singlechannels") x <- cbind(object$R[, arrays], object$G[, arrays]) if (is.null(groups)) { groups <- c(length(arrays), length(arrays)) if (is.null(col)) colors <- rep(c("red", "green"), groups) if (!is.null(col)) { if (length(col) != 2) { print("number of groups=2 not equal to number of col") colors <- "black" } else { colors <- rep(col, groups) } } } else { if (!is.null(col)) { if (length(as.vector(table(groups))) != length(col)) { print("number of groups not equal to number of col") colors <- "black" } else { colors <- col[groups] } } else { print("warning no colors in col specified for the groups") colors <- "black" } } } } } if (log.transform) x <- log(x, 2) dens.x <- matDensities(x) # matplot(dens.x$X, dens.x$Y, xlab = "Intensity", ylab = "Density", # main = "RG densities", type = "l", col = colors, lwd = 2, # lty = 1) return(dens.x) } plotSqueletonDensities <- function(dens.x, title, index = NULL) { if(!is.null(index)) { step <- dim(dens.x$X)[2]/2 index <- c(index, index + step) dens.x$X <- dens.x$X[, index] dens.x$Y <- dens.x$Y[, index] colors <- c("red", "blue") } else { colors <- rep(c("red", "blue"), c((dim(dens.x$X)[2])/2, (dim(dens.x$X)[2])/2)) } matplot(dens.x$X, dens.x$Y, xlab = "Intensity", ylab = "Density", main = title, type = "l", col = colors, lwd = 2, lty = 1) } ## The above work; just need to put them here; ## should be, in line 238 (aporx) of dnmad.R, a 5 row by 2 col plot, ## by rpow: hists ## densities (all, the array) ## BG ## FG ## M (unnorm and norm) plot.print.tip.lowess1 <- function (x, layout, norm = "n", image.id = 1, palette = rainbow(layout$ngrid.r * layout$ngrid.c), lty.palette = rep(1, layout$ngrid.r * layout$ngrid.c), bg.subtract = FALSE, ...) { ## from plot.print.tip.lowess, to eliminate warnings ## from ... if(bg.subtract) tmp <- ma.func(R = x$R[, image.id], G = x$G[, image.id], Rb = x$Rb[, image.id], Gb = x$Gb[, image.id], layout, norm = norm, pout = FALSE, ...) else tmp <- ma.func(R = x$R[, image.id], G = x$G[, image.id], Rb = NULL, Gb = NULL, layout, norm = norm, pout = FALSE, ...) plot(tmp$A, tmp$M, xlab = "A", ylab = "M", ...) npin <- layout$ngrid.r * layout$ngrid.c nspot <- layout$nspot.c * layout$nspot.r pin <- rep(1:npin, rep(nspot, npin)) for (i in 1:npin) { index <- pin == i tM <- tmp$M[index] tA <- tmp$A[index] ind2 <- is.na(tM) | is.na(tA) | is.infinite(tM) | is.infinite(tA) try(smoothnum <- lowess(tA[!ind2], tM[!ind2])) try(lines(approx(smoothnum), col = palette[i], lty = lty.palette[i])) } } plot.print.tip.lowess2 <- function (tmp, layout, image.id = 1, palette = rainbow(layout$ngrid.r * layout$ngrid.c), lty.palette = rep(1, layout$ngrid.r * layout$ngrid.c), ...) ## from plot.print.tip.lowess, but the data are ## already normalized { tmp <- tmp[, image.id] plot(tmp$A, tmp$M, xlab = "A", ylab = "M", ...) npin <- layout$ngrid.r * layout$ngrid.c nspot <- layout$nspot.c * layout$nspot.r pin <- rep(1:npin, rep(nspot, npin)) for (i in 1:npin) { index <- pin == i tM <- tmp$M[index] tA <- tmp$A[index] ind2 <- is.na(tM) | is.na(tA) | is.infinite(tM) | is.infinite(tA) try(smoothnum <- lowess(tA[!ind2], tM[!ind2])) try(lines(approx(smoothnum), col = palette[i], lty = lty.palette[i])) } } my.imageplot <- function (z, layout = list(ngrid.r = 12, ngrid.c = 4, nspot.r = 26, nspot.c = 26), low = NULL, high = NULL, ncolors = 123, zerocenter = NULL, zlim = NULL, ...) { ## from imageplot; to allow multiple plots in a page, and other minor ## things. gr <- layout$ngrid.r gc <- layout$ngrid.c sr <- layout$nspot.r sc <- layout$nspot.c if (is.null(gr) || is.null(gc) || is.null(sr) || is.null(sc)) stop("Layout needs to contain components ngrid.r, ngrid.c, nspot.r and spot.c") if (length(z) != gr * gc * sr * sc) stop("Number of image spots does not agree with layout dimensions") if (is.character(low)) low <- col2rgb(low)/255 if (is.character(high)) high <- col2rgb(high)/255 if (!is.null(low) && is.null(high)) high <- c(1, 1, 1) - low if (is.null(low) && !is.null(high)) low <- c(1, 1, 1) - high if (!is.null(zlim)) { z <- pmax(zlim[1], z) z <- pmin(zlim[2], z) } zr <- range(z, na.rm = TRUE) zmax <- max(abs(zr)) zmin <- zr[1] if (is.null(zerocenter)) zerocenter <- (zmin < 0) if (zerocenter) { if (is.null(low)) low <- c(0, 1, 0) if (is.null(high)) high <- c(1, 0, 0) if (is.null(zlim)) zlim <- c(-zmax, zmax) } else { if (is.null(low)) low <- c(1, 1, 1) if (is.null(high)) high <- c(0, 0, 1) if (is.null(zlim)) zlim <- c(zmin, zmax) } col <- rgb(seq(low[1], high[1], len = ncolors), seq(low[2], high[2], len = ncolors), seq(low[3], high[3], len = ncolors)) dim(z) <- c(sc, sr, gc, gr) z <- aperm(z, perm = c(2, 4, 1, 3)) dim(z) <- c(gr * sr, gc * sc) z1 <<- z # old.par <- par(no.readonly = TRUE) # on.exit(par(old.par)) # par(mar = rep(1, 4)) image(0:(gr * sr), 0:(gc * sc), z, zlim = zlim, col = col, axes = FALSE, ...) for (igrid in 0:gc) lines(c(0, gr * sr), rep(igrid * sc, 2)) for (igrid in 0:gr) lines(rep(igrid * sr, 2), c(0, gc * sc)) invisible() } my.imageplot2 <- function (z, layout = list(ngrid.r = 12, ngrid.c = 4, nspot.r = 26, nspot.c = 26), low = NULL, high = NULL, ncolors = 123, zerocenter = NULL, zlim = NULL, ...) { ## from imageplot; to allow multiple plots in a page, and other minor ## things. gr <- layout$ngrid.r gc <- layout$ngrid.c sr <- layout$nspot.r sc <- layout$nspot.c if (is.null(gr) || is.null(gc) || is.null(sr) || is.null(sc)) stop("Layout needs to contain components ngrid.r, ngrid.c, nspot.r and spot.c") if (length(z) != gr * gc * sr * sc) stop("Number of image spots does not agree with layout dimensions") if (is.character(low)) low <- col2rgb(low)/255 if (is.character(high)) high <- col2rgb(high)/255 if (!is.null(low) && is.null(high)) high <- c(1, 1, 1) - low if (is.null(low) && !is.null(high)) low <- c(1, 1, 1) - high if (!is.null(zlim)) { z <- pmax(zlim[1], z) z <- pmin(zlim[2], z) } zr <- range(z, na.rm = TRUE) zmax <- max(abs(zr)) zmin <- zr[1] if (is.null(zerocenter)) zerocenter <- (zmin < 0) if (zerocenter) { if (is.null(low)) low <- c(0, 1, 0) if (is.null(high)) high <- c(1, 0, 0) if (is.null(zlim)) zlim <- c(-zmax, zmax) } else { if (is.null(low)) low <- c(1, 1, 1) if (is.null(high)) high <- c(0, 0, 1) if (is.null(zlim)) zlim <- c(zmin, zmax) } col <- rgb(seq(low[1], high[1], len = ncolors), seq(low[2], high[2], len = ncolors), seq(low[3], high[3], len = ncolors)) dim(z) <- c(sc, sr, gc, gr) z <- aperm(z, perm = c(2, 4, 1, 3)) dim(z) <- c(gr * sr, gc * sc) z <- t(z) z <- z[,ncol(z):1] z2 <<- z # old.par <- par(no.readonly = TRUE) # on.exit(par(old.par)) # par(mar = rep(1, 4)) ## The following is to make the plots look nicer (spots are squares) horizontal.size <- gc * sc vertical.size <- gr * sr if(vertical.size < horizontal.size) { vert.span <- 0.8 * vertical.size/horizontal.size vl.inf <- (1 - vert.span)/2 vl.sup <- (1 + vert.span)/2 par(plt = c(0.1, 0.9, vl.inf, vl.sup)) } else if(vertical.size > horizontal.size) { horz.span <- 0.8 * horizontal.size/vertical.size hl.inf <- (1 - horz.span)/2 hl.sup <- (1 + horz.span)/2 par(plt = c(hl.inf, hl.sup, 0.1, 0.9)) } image(0:(gc * sc), 0:(gr * sr), z, zlim = zlim, col = col, axes = FALSE, xlab = "", ylab = "", ...) for (igrid in 0:gr) lines(c(0, gc * sc), rep(igrid * sr, 2)) for (igrid in 0:gc) lines(rep(igrid * sc, 2), c(0, gr * sr)) invisible() } ##diagnostics.for.limma <- function(d.unnorm, d.norm, outfile = NULL, ## png = FALSE) { ## ## plot in a file (or files, if png) my preferred diagnostics, ## ## with one plot per array (and the final boxplot for all arrays). ## n.arrays <- dim(d.unnorm$G)[2] ## names.arrays <- colnames(d.unnorm$G) ## if(!is.null(outfile)) #pdf(file = outfile, width = 12, height = 8) ## #postscript(file = outfile, width = 12, height = 8, horizontal = FALSE) ## if(png) ## png(file = paste(outfile, "%d.png", sep = ""), ## width = 1080, height = 720) ## else postscript(file = outfile, width = 12, height = 8, ## horizontal = FALSE) ## layout(matrix(c(1, 2, 5, 6, 3, 4, 7, 8), 2, 4, byrow = TRUE), ## widths = rep(7, 8)) ## par(mar = rep(1, 4)) ## par(oma = c(2, 2, 4, 2)) ## for(i in 1:n.arrays) { ## MA.raw <- MA.RG(d.unnorm[, i]) ## R.back <- log(d.unnorm$Rb[, i], 2) ## G.back <- log(d.unnorm$Gb[, i], 2) ## R.back[is.infinite(R.back)] <- NA ## G.back[is.infinite(G.back)] <- NA ## my.imageplot2(R.back, layout.array, low = "white", ## high = "red", main = "Red Bg") ## my.imageplot2(G.back, layout.array, low = "white", ## high = "green", main = "Green Bg") ## my.imageplot2(MA.raw$A, layout.array, low = "white", ## high = "blue", main = "Unnorm. A") ## my.imageplot2(MA.raw$M, layout.array, low = "green", ## high = "red", main = "Unnorm. M", ncolors = 10) ## par(mar = rep(2, 4)) ## plot.scale.box(MA.raw$M[, 1], layout.array, ## col = rainbow(layout.array$ngrid.r * ## layout.array$ngrid.c), ## main = "Unnormalized M") ## abline(h = 0, lty = 2) ## plot.scale.box(d.norm$M[, i], layout.array, ## col = rainbow(layout.array$ngrid.r * ## layout.array$ngrid.c), ## main = "Normalized M") ## abline(h = 0, lty = 2) ## par(mar = rep(2, 4)) ## plot.print.tip.lowess1(d.unnorm, layout.array, pch = ".", image.id = i, ## main = "Unnorm M-A", cex = 0.1) ## abline(h = 0, lty = 2) ## plot.print.tip.lowess2(d.norm, layout.array, pch = ".", image.id = i, ## main = "Norm M-A", cex = 0.1) ## abline(h = 0, lty = 2) ## mtext(paste("Array", names.arrays[i]), side = 3, ## outer = TRUE, line = 2, cex = 1.8) ## } ## layout(1) ## par(mfrow = c(1,1)) ## boxplot(d.norm$M ~ col(d.norm$M), names = names.arrays, ## main = "all arrays, normalized") ## abline(h = 0, lty = 2) ## if(!is.null(outfile)) dev.off() ##} MA.RG <- function (object, log.transform = TRUE, bg.subtract = FALSE) { ## Like the original, but: ### takes a parameter for background subtraction ### and always returns the whole object (i.e., does not eliminate columns) R <- object$R G <- object$G if (bg.subtract) { R <- R - object$Rb G <- G - object$Gb } if (log.transform) { R[R <= 0] <- NA G[G <= 0] <- NA R <- log(R, 2) G <- log(G, 2) } object$R <- object$G <- object$Rb <- object$Gb <- NULL object$M <- as.matrix(R - G) object$A <- as.matrix((R + G)/2) new("MAList", unclass(object)) } normalizeWithinArrays <- function (object, layout = object$printer, method = "printtiploess", weights = object$weights, span = 0.3, iterations = 4, controlspots = NULL, df = 5, robust = "M", bg.subtract = FALSE) {## From normalizeWithinArrays, in limma, but allowing for no points ## for loess and reporting a warning ## if(!bg.subtract) ## backgroundCorrect(object, method = "none") ## if (!is(object, "MAList")) object <- MA.RG(object, bg.subtract = bg.subtract) choices <- c("none", "median", "loess", "printtiploess", "composite", "robustspline") method <- match.arg(method, choices) if (method == "none") return(object) narrays <- ncol(object$M) ngenes <- nrow(object$M) if (method == "median") { for (j in 1:narrays) object$M[, j] <- object$M[, j] - median(object$M[, j], na.rm = TRUE) return(object) } switch(method, loess = { for (j in 1:narrays) { ### print(paste("We are at array", j)) y <- object$M[, j] x <- object$A[, j] w <- weights[, j] object$M[, j] <- loessFit(y, x, w, span = span, iterations = iterations)$residuals } }, printtiploess = { if (is.null(layout)) stop("Layout argument not specified") ngr <- layout$ngrid.r ngc <- layout$ngrid.c nspots <- layout$nspot.r * layout$nspot.c for (j in 1:narrays) { spots <- 1:nspots for (gridr in 1:ngr) for (gridc in 1:ngc) { y <- object$M[spots, j] x <- object$A[spots, j] w <- weights[spots, j] # print(paste("We are at array", j, "gridr", gridr, "gridc", gridc)) # print("******** summary y") # print(summary(y)) # print("******* summary x") # print(summary(x)) ## next two lines, and the try are my addition num.fill <- length(object$M[spots, j]) object$M[spots, j] <- rep(NA, num.fill) no.objects <- try(object$M[spots, j] <- loessFit(y, x, w, span = span, iterations = iterations)$residuals) if(class(no.objects) == "try-error") write.table( paste( "WARNING: in print-tip loess normalization, no samples in array", j, "grid row", gridr, "grid column", gridc), file = paste(datadir, "warning.file", sep =""), append = TRUE, row.names = FALSE, col.names = FALSE, quote = FALSE) else { sum.no.nas <- sum(!is.na(object$M[spots, j])) if(sum.no.nas < 100) { write.table( paste( "WARNING: in print-tip loess normalization, only", sum.no.nas, "samples in array", j, "grid row", gridr, "grid column", gridc), file = paste(datadir, "warning.file", sep =""), append = TRUE, row.names = FALSE, col.names = FALSE, quote = FALSE) } } spots <- spots + nspots } } }, composite = { if (is.null(layout)) stop("Layout argument not specified") ntips <- layout$ngrid.r * layout$ngrid.c nspots <- layout$nspot.r * layout$nspot.c for (j in 1:narrays) { y <- object$M[, j] x <- object$A[, j] w <- weights[, j] fit <- loess(y ~ x, weights = w, span = span, subset = controlspots, na.action = na.exclude, degree = 0, surface = "direct", family = "symmetric", trace.hat = "approximate", iterations = iterations) global <- predict(fit, newdata = x) alpha <- (rank(x) - 1)/sum(!is.na(x)) spots <- 1:nspots for (tip in 1:ntips) { y <- object$M[spots, j] x <- object$A[spots, j] w <- weights[spots, j] local <- loessFit(y, x, w, span = span, iterations = iterations)$fitted object$M[spots, j] <- object$M[spots, j] - alpha[spots] * global[spots] - (1 - alpha[spots]) * local spots <- spots + nspots } } }, robustspline = { if (is.null(layout)) stop("Layout argument not specified") for (j in 1:narrays) object$M[, j] <- normalizeRobustSpline(object$M[, j], object$A[, j], layout, df = df, method = robust) }) object }