--- /dev/null
+# Copyright (C) 1999 Anthony D. Long, Ecology & Evolutionary Biology
+# University of California, Irvine (tdlong@uci.edu)
+
+# Portions Copyright 2002,2003 Don Armstrong <don@donarmstrong.com>
+
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+
+# $Id: bayesian_analysis.R,v 1.5 2006/07/12 03:27:03 don Exp $
+
+
+microarray.bayesian.plot <-
+ function(data) {
+
+ ## png(file = paste(file,"_%03d.png",sep=""), width = 1600, height = 1200)
+ plotting.function <- getFunction("plot")
+ text.function <- getFunction("text")
+ use.grid <- 0
+ if (NROW(data)>1e4)
+ use.grid <- 1
+ if (use.grid) {
+ nbins <- 80
+ plotting.function <-
+ function(x,y,
+ xbins = nbins, shape = 1,
+ xbnds = range(x), ybnds = range(y),
+ xlab = NULL, ylab = NULL, IDs = FALSE,...) {
+ return(plot(hexbin(x,y,
+ xbins=xbins,
+ shape=shape,
+ xbnds=xbnds,
+ xlab=xlab,
+ ylab=ylab,
+ IDs=IDs),
+ colramp=BTY,
+ ...))
+ }
+ text.function <- function(x,y,text) {
+ grid.text(text,x=unit(x,"native"),y=unit(y,"native"))
+ }
+ }
+ plotting.function(data[, "controls.mean"]/1e3,
+ data[, "experiments.mean"]/1e3,
+ xlab = expression(paste("Mean of Control Expression (",10**3,")")),
+ ylab = expression(paste("Mean of Experimental Expression (",10**3,")")),
+ main = "Mean expression level comparison"
+ )
+ plotting.function(data[, "controls.mean.log"],
+ data[, "experiments.mean.log"],
+ xlab = "Mean of Controls",
+ ylab = "Mean of Experimentals",
+ main = "log transformed expression level comparison"
+ )
+ ## controls mean v controls sd
+ temp.y <- data[,"controls.sd"]/1000
+ temp.x <- data[, "controls.mean"]/1000
+ lmline <- lm(temp.y~temp.x)
+ temp <-
+ plotting.function(data[, "controls.mean"]/1e3,
+ data[, "controls.sd"]/1e3,
+ xlab = expression(paste("Mean of Control Expression (",10**3,")")),
+ ylab = expression(paste("SD of Control Expression (",10**3,")")),
+ main = "Control SD as a function of Control Mean"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ }
+ else {
+ abline(lmline)
+ }
+ xlim <- range(data[,"controls.mean"]/1e3)
+ ylim <- range(data[,"controls.sd"]/1e3)
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+ if (use.grid)
+ popViewport(1)
+
+ ## experiments mean v experiments sd
+ temp.y <- data[,"experiments.sd"]/1000
+ temp.x <- data[,"experiments.mean"]/1000
+ temp <-
+ plotting.function(data[, "experiments.mean"]/1e3,
+ data[, "experiments.sd"]/1e3,
+ xlab = expression(paste("Mean of Experimental Expression (",10**3,")")),
+ ylab = expression(paste( "SD of Experimental Expression (",10**3,")")),
+ main = "Experiment SD as a function of Experiment Mean"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ }
+ else {
+ abline(lmline)
+ }
+ xlim <- range(data[,"experiments.mean"])
+ ylim <- range(data[,"experiments.sd"])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+ if (use.grid)
+ popViewport(1)
+
+ ## controls log mean v controls log sd
+ lmline <- lm(data[,"controls.sd.log"]~data[, "controls.mean.log"])
+ temp <- plotting.function(data[, "controls.mean.log"],
+ data[, "controls.sd.log"],
+ xlab = "Mean of log(Control Expression)",
+ ylab = "SD of log(Control Expression)",
+ main = "log(Control SD) as a function of log(Control Mean)"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ }
+ else {
+ abline(lmline)
+ }
+ xlim <- range(data[,"controls.mean.log"])
+ ylim <- range(data[,"controls.sd.log"])
+ print(xlim)
+ print(ylim)
+ print(summary(lmline))
+ print(0.5*(xlim[2]-xlim[1])+xlim[1])
+ print(0.9*(ylim[2]-ylim[1])+ylim[1])
+ if (use.grid)
+ print(current.vpTree())
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+ if (use.grid)
+ popViewport(1)
+
+ ## experiments mean v experiments sd
+ lmline <- lm(data[,"experiments.sd.log"]~data[, "experiments.mean.log"])
+ temp <- plotting.function(data[, "experiments.mean.log"],
+ data[, "experiments.sd.log"],
+ xlab = "Mean of log(Experiment Expression)",
+ ylab = "SD of log(Experiment Expression)",
+ main = "log(Experiment SD) as a function of log(Experiment Mean)"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ }
+ else {
+ abline(lmline)
+ }
+ xlim <- range(data[,"experiments.mean.log"])
+ ylim <- range(data[,"experiments.sd.log"])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+ if (use.grid)
+ popViewport(1)
+
+ ## mean expression vs significance
+ mean.expression <- apply(data[,c("experiments.mean","controls.mean")],1,mean)
+ lmline <- lm(data[,"t.test.p"] ~
+ mean.expression)
+ temp <- plotting.function(mean.expression,
+ data[, "t.test.p"],
+ xlab = "Mean Expression",
+ ylab = "p value of T Test",
+ main = "Significance as a function of Mean Expression"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ panel.abline(h=0.05)
+ }
+ else {
+ abline(lmline)
+ abline(h=0.05)
+ }
+ xlim <- range(mean.expression)
+ ylim <- range(data[, "t.test.p"])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+
+ if (use.grid)
+ popViewport(1)
+
+ ## log transformed mean expression vs significance
+ mean.expression <- apply(data[,c("experiments.mean.log","controls.mean.log")],1,mean)
+ lmline <- lm(data[, "t.test.p.log"] ~
+ mean.expression)
+ temp <-
+ plotting.function(mean.expression,
+ data[, "t.test.p.log"],
+ xlab = "Mean Expression",
+ ylab = "p value of T Test log(Expression)",
+ main = "Significance of log(Expression) as a function of Mean Expression"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ panel.abline(h=0.05)
+ }
+ else {
+ abline(lmline)
+ abline(h=0.05)
+ }
+ xlim <- range(mean.expression)
+ ylim <- range(data[, "t.test.p.log"])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+ if (use.grid)
+ popViewport(1)
+
+
+
+ ## Fold change vs significance
+ lmline <- lm(data[abs(data[, "fold.change"]) < 20, "t.test.p"] ~
+ abs(data[abs(data[, "fold.change"]) < 20, "fold.change"]))
+ temp <- plotting.function(abs(data[abs(data[,"fold.change"]) < 20, "fold.change"]),
+ data[abs(data[,"fold.change"]) < 20, "t.test.p"],
+ xlab = "Fold Change of Expression",
+ ylab = "p value of T Test",
+ main = "Significance as a function of Fold Change"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ panel.abline(h=0.05)
+ }
+ else {
+ abline(lmline)
+ abline(h=0.05)
+ }
+ xlim <- range(abs(data[abs(data[,"fold.change"]) < 20, "fold.change"]))
+ ylim <- range(data[abs(data[,"fold.change"]) < 20, "t.test.p"])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+
+ if (use.grid)
+ popViewport(1)
+
+ ## log transformed fold change vs significance
+ print("log transformed fold change vs significance")
+ temp.data <- na.omit(cbind(data[abs(data[, "fold.change"]) < 20, "t.test.p.log"],
+ abs(data[abs(data[, "fold.change"]) < 20, "fold.change"])))
+ lmline <- lm(temp.data[,1]~temp.data[,2])
+ temp <-
+ plotting.function(temp.data[,2],temp.data[,1],
+ xlab = "Fold Change of Expression",
+ ylab = "p value of T Test log(Expression)",
+ main = "Significance of log(Expression) as a function of Fold Change"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ panel.abline(h=0.05)
+ }
+ else {
+ abline(lmline)
+ abline(h=0.05)
+ }
+ xlim <- range(temp.data[,2])
+ ylim <- range(temp.data[,1])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+ if (use.grid)
+ popViewport(1)
+
+ ## bayesian vs t test
+ temp.data <- log(na.omit(cbind(data[,"t.test.p"],
+ data[,"p.bayes"])))
+ lmline <- lm(temp.data[,1]~temp.data[,2])
+ temp <-
+ plotting.function(temp.data[,2],temp.data[,1],
+ xlab = "log(p value) of Bayesian",
+ ylab = "log(p value) of t test",
+ main = "Correspondance of t-test and bayesian p values"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ }
+ else {
+ abline(lmline)
+ }
+ xlim <- range(temp.data[,2])
+ ylim <- range(temp.data[,1])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+
+ if (use.grid)
+ popViewport(1)
+
+ ## log(bayesian) vs t test
+ print("log(bayesian) vs t test")
+ temp.data <- log(na.omit(cbind(data[,"t.test.p.log"],
+ data[,"p.log.bayes"])))
+ lmline <- lm(temp.data[,1]~temp.data[,2])
+ temp <-
+ plotting.function(temp.data[,2],temp.data[,1],
+ xlab = "log(p value) of Bayesian (log transformed)",
+ ylab = "log(p value) of t test (log transformed)",
+ main = "Correspondance of t-test and bayesian p values of log transformed values"
+ )
+ if (use.grid) {
+ pushHexport(temp$plot.vp)
+ panel.abline(lmline)
+ }
+ else {
+ abline(lmline)
+ }
+ xlim <- range(temp.data[,2])
+ ylim <- range(temp.data[,1])
+ text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
+ 0.9*(ylim[2]-ylim[1])+ylim[1],
+ substitute(r**2 == y,
+ list(y=format(summary(lmline)$r.squared,
+ digits=4)))
+ )
+
+ if (use.grid)
+ popViewport(1)
+
+
+ ## histograms
+ print("histograms")
+ plot(hist(data[,"t.test.p"],breaks=40,plot=FALSE),main="Histogram of T Test p values")
+ plot(hist(data[,"t.test.p.log"],breaks=40,plot=FALSE),main="Histogram of T Test p values for log transformed data")
+ plot(hist(data[,"p.bayes"],breaks=40,plot=FALSE),main="Histogram of Bayesian p values")
+ plot(hist(data[,"p.log.bayes"],breaks=40,plot=FALSE),main="Histogram of Bayesian p values for log transformed data")
+
+ }
+
projects / don_microarray.git / blobdiff