R/bayesian_analysis_plot.R

   1 #   Copyright (C) 1999 Anthony D. Long, Ecology & Evolutionary Biology
   2 #      University of California, Irvine (tdlong@uci.edu)
   3
   4 # Portions Copyright 2002,2003 Don Armstrong <don@donarmstrong.com>
   5
   6 # This program is free software; you can redistribute it and/or
   7 # modify it under the terms of the GNU General Public License
   8 # as published by the Free Software Foundation; either version 2
   9 # of the License, or (at your option) any later version.
  10 #
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 # $Id: bayesian_analysis.R,v 1.5 2006/07/12 03:27:03 don Exp $
  17
  18
  19 microarray.bayesian.plot <-
  20   function(data) {
  21
  22     ## png(file = paste(file,"_%03d.png",sep=""), width = 1600, height = 1200)
  23     plotting.function <- getFunction("plot")
  24     text.function <- getFunction("text")
  25     use.grid <- 0
  26     if (NROW(data)>1e4)
  27       use.grid <- 1
  28     if (use.grid) {
  29       nbins <- 80
  30       plotting.function <-
  31         function(x,y,
  32                  xbins = nbins, shape = 1,
  33                  xbnds = range(x), ybnds = range(y),
  34                  xlab = NULL, ylab = NULL, IDs = FALSE,...) {
  35           return(plot(hexbin(x,y,
  36                              xbins=xbins,
  37                              shape=shape,
  38                              xbnds=xbnds,
  39                              xlab=xlab,
  40                              ylab=ylab,
  41                              IDs=IDs),
  42                       colramp=BTY,
  43                       ...))
  44       }
  45       text.function <- function(x,y,text) {
  46         grid.text(text,x=unit(x,"native"),y=unit(y,"native"))
  47       }
  48     }
  49     plotting.function(data[,    "controls.mean"]/1e3,
  50                       data[, "experiments.mean"]/1e3,
  51                       xlab = expression(paste("Mean of Control Expression (",10**3,")")),
  52                       ylab = expression(paste("Mean of Experimental Expression (",10**3,")")),
  53                       main = "Mean expression level comparison"
  54                       )
  55     plotting.function(data[,    "controls.mean.log"],
  56                       data[, "experiments.mean.log"],
  57                       xlab = "Mean of Controls",
  58                       ylab = "Mean of Experimentals",
  59                       main = "log transformed expression level comparison"
  60                       )
  61     ## controls mean v controls sd
  62     temp.y <- data[,"controls.sd"]/1000
  63     temp.x <- data[, "controls.mean"]/1000
  64     lmline <- lm(temp.y~temp.x)
  65     temp <-
  66       plotting.function(data[, "controls.mean"]/1e3,
  67                         data[, "controls.sd"]/1e3,
  68                         xlab = expression(paste("Mean of Control Expression (",10**3,")")),
  69                         ylab = expression(paste("SD of Control Expression (",10**3,")")),
  70                         main = "Control SD as a function of Control Mean"
  71                         )
  72     if (use.grid) {
  73       pushHexport(temp$plot.vp)
  74       panel.abline(lmline)
  75     }
  76     else {
  77       abline(lmline)
  78     }
  79     xlim <- range(data[,"controls.mean"]/1e3)
  80     ylim <- range(data[,"controls.sd"]/1e3)
  81     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
  82          0.9*(ylim[2]-ylim[1])+ylim[1],
  83          substitute(r**2 == y,
  84                     list(y=format(summary(lmline)$r.squared,
  85                                     digits=4)))
  86          )
  87     if (use.grid)
  88       popViewport(1)
  89
  90     ## experiments mean v experiments sd
  91     temp.y <- data[,"experiments.sd"]/1000
  92     temp.x <- data[,"experiments.mean"]/1000
  93     temp <-
  94       plotting.function(data[, "experiments.mean"]/1e3,
  95                         data[, "experiments.sd"]/1e3,
  96                         xlab = expression(paste("Mean of Experimental Expression (",10**3,")")),
  97                         ylab = expression(paste(  "SD of Experimental Expression (",10**3,")")),
  98                         main = "Experiment SD as a function of Experiment Mean"
  99                         )
 100     if (use.grid) {
 101       pushHexport(temp$plot.vp)
 102       panel.abline(lmline)
 103     }
 104     else {
 105       abline(lmline)
 106     }
 107     xlim <- range(data[,"experiments.mean"])
 108     ylim <- range(data[,"experiments.sd"])
 109     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 110          0.9*(ylim[2]-ylim[1])+ylim[1],
 111          substitute(r**2 == y,
 112                     list(y=format(summary(lmline)$r.squared,
 113                                     digits=4)))
 114          )
 115     if (use.grid)
 116       popViewport(1)
 117
 118     ## controls log mean v controls log sd
 119     lmline <- lm(data[,"controls.sd.log"]~data[, "controls.mean.log"])
 120     temp <- plotting.function(data[, "controls.mean.log"],
 121                               data[, "controls.sd.log"],
 122                               xlab = "Mean of log(Control Expression)",
 123                               ylab = "SD of log(Control Expression)",
 124                               main = "log(Control SD) as a function of log(Control Mean)"
 125                               )
 126     if (use.grid) {
 127       pushHexport(temp$plot.vp)
 128       panel.abline(lmline)
 129     }
 130     else {
 131       abline(lmline)
 132     }
 133     xlim <- range(data[,"controls.mean.log"])
 134     ylim <- range(data[,"controls.sd.log"])
 135     print(xlim)
 136     print(ylim)
 137     print(summary(lmline))
 138     print(0.5*(xlim[2]-xlim[1])+xlim[1])
 139     print(0.9*(ylim[2]-ylim[1])+ylim[1])
 140     if (use.grid)
 141       print(current.vpTree())
 142     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 143          0.9*(ylim[2]-ylim[1])+ylim[1],
 144          substitute(r**2 == y,
 145                     list(y=format(summary(lmline)$r.squared,
 146                                     digits=4)))
 147          )
 148     if (use.grid)
 149       popViewport(1)
 150
 151     ## experiments mean v experiments sd
 152     lmline <- lm(data[,"experiments.sd.log"]~data[, "experiments.mean.log"])
 153     temp <- plotting.function(data[, "experiments.mean.log"],
 154                               data[, "experiments.sd.log"],
 155                               xlab = "Mean of log(Experiment Expression)",
 156                               ylab = "SD of log(Experiment Expression)",
 157                               main = "log(Experiment SD) as a function of log(Experiment Mean)"
 158                               )
 159     if (use.grid) {
 160       pushHexport(temp$plot.vp)
 161       panel.abline(lmline)
 162     }
 163     else {
 164       abline(lmline)
 165     }
 166     xlim <- range(data[,"experiments.mean.log"])
 167     ylim <- range(data[,"experiments.sd.log"])
 168     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 169          0.9*(ylim[2]-ylim[1])+ylim[1],
 170          substitute(r**2 == y,
 171                     list(y=format(summary(lmline)$r.squared,
 172                            digits=4)))
 173          )
 174     if (use.grid)
 175       popViewport(1)
 176
 177     ## mean expression vs significance
 178     mean.expression <- apply(data[,c("experiments.mean","controls.mean")],1,mean)
 179     lmline <- lm(data[,"t.test.p"] ~
 180                  mean.expression)
 181     temp <- plotting.function(mean.expression,
 182                               data[, "t.test.p"],
 183                               xlab = "Mean Expression",
 184                               ylab = "p value of T Test",
 185                               main = "Significance as a function of Mean Expression"
 186                               )
 187     if (use.grid) {
 188       pushHexport(temp$plot.vp)
 189       panel.abline(lmline)
 190       panel.abline(h=0.05)
 191     }
 192     else {
 193       abline(lmline)
 194       abline(h=0.05)
 195     }
 196     xlim <- range(mean.expression)
 197     ylim <- range(data[, "t.test.p"])
 198     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 199          0.9*(ylim[2]-ylim[1])+ylim[1],
 200          substitute(r**2 == y,
 201                     list(y=format(summary(lmline)$r.squared,
 202                                     digits=4)))
 203          )
 204
 205     if (use.grid)
 206       popViewport(1)
 207
 208     ## log transformed mean expression vs significance
 209     mean.expression <- apply(data[,c("experiments.mean.log","controls.mean.log")],1,mean)
 210     lmline <- lm(data[, "t.test.p.log"] ~
 211                  mean.expression)
 212     temp <-
 213       plotting.function(mean.expression,
 214                         data[, "t.test.p.log"],
 215                         xlab = "Mean Expression",
 216                         ylab = "p value of T Test log(Expression)",
 217                         main = "Significance of log(Expression) as a function of Mean Expression"
 218                         )
 219     if (use.grid) {
 220       pushHexport(temp$plot.vp)
 221       panel.abline(lmline)
 222       panel.abline(h=0.05)
 223     }
 224     else {
 225       abline(lmline)
 226       abline(h=0.05)
 227     }
 228     xlim <- range(mean.expression)
 229     ylim <- range(data[, "t.test.p.log"])
 230     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 231          0.9*(ylim[2]-ylim[1])+ylim[1],
 232          substitute(r**2 == y,
 233                     list(y=format(summary(lmline)$r.squared,
 234                                     digits=4)))
 235          )
 236     if (use.grid)
 237       popViewport(1)
 238
 239
 240
 241     ## Fold change vs significance
 242     lmline <- lm(data[abs(data[, "fold.change"]) < 20, "t.test.p"] ~
 243                  abs(data[abs(data[, "fold.change"]) < 20, "fold.change"]))
 244     temp <- plotting.function(abs(data[abs(data[,"fold.change"]) < 20, "fold.change"]),
 245                               data[abs(data[,"fold.change"]) < 20, "t.test.p"],
 246                               xlab = "Fold Change of Expression",
 247                               ylab = "p value of T Test",
 248                               main = "Significance as a function of Fold Change"
 249                               )
 250     if (use.grid) {
 251       pushHexport(temp$plot.vp)
 252       panel.abline(lmline)
 253       panel.abline(h=0.05)
 254     }
 255     else {
 256       abline(lmline)
 257       abline(h=0.05)
 258     }
 259     xlim <- range(abs(data[abs(data[,"fold.change"]) < 20, "fold.change"]))
 260     ylim <- range(data[abs(data[,"fold.change"]) < 20, "t.test.p"])
 261     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 262          0.9*(ylim[2]-ylim[1])+ylim[1],
 263          substitute(r**2 == y,
 264                     list(y=format(summary(lmline)$r.squared,
 265                                     digits=4)))
 266          )
 267
 268     if (use.grid)
 269       popViewport(1)
 270
 271     ## log transformed fold change vs significance
 272     print("log transformed fold change vs significance")
 273     temp.data <- na.omit(cbind(data[abs(data[, "fold.change"]) < 20, "t.test.p.log"],
 274                                abs(data[abs(data[, "fold.change"]) < 20, "fold.change"])))
 275     lmline <- lm(temp.data[,1]~temp.data[,2])
 276     temp <-
 277       plotting.function(temp.data[,2],temp.data[,1],
 278                         xlab = "Fold Change of Expression",
 279                         ylab = "p value of T Test log(Expression)",
 280                         main = "Significance of log(Expression) as a function of Fold Change"
 281                         )
 282     if (use.grid) {
 283       pushHexport(temp$plot.vp)
 284       panel.abline(lmline)
 285       panel.abline(h=0.05)
 286     }
 287     else {
 288       abline(lmline)
 289       abline(h=0.05)
 290     }
 291     xlim <- range(temp.data[,2])
 292     ylim <- range(temp.data[,1])
 293     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 294          0.9*(ylim[2]-ylim[1])+ylim[1],
 295          substitute(r**2 == y,
 296                     list(y=format(summary(lmline)$r.squared,
 297                                     digits=4)))
 298          )
 299     if (use.grid)
 300       popViewport(1)
 301
 302     ## bayesian vs t test
 303     temp.data <- log(na.omit(cbind(data[,"t.test.p"],
 304                                    data[,"p.bayes"])))
 305     lmline <- lm(temp.data[,1]~temp.data[,2])
 306     temp <-
 307       plotting.function(temp.data[,2],temp.data[,1],
 308                         xlab = "log(p value) of Bayesian",
 309                         ylab = "log(p value) of t test",
 310                         main = "Correspondance of t-test and bayesian p values"
 311                         )
 312     if (use.grid) {
 313       pushHexport(temp$plot.vp)
 314       panel.abline(lmline)
 315     }
 316     else {
 317       abline(lmline)
 318     }
 319     xlim <- range(temp.data[,2])
 320     ylim <- range(temp.data[,1])
 321     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 322          0.9*(ylim[2]-ylim[1])+ylim[1],
 323          substitute(r**2 == y,
 324                     list(y=format(summary(lmline)$r.squared,
 325                                     digits=4)))
 326          )
 327
 328     if (use.grid)
 329       popViewport(1)
 330
 331     ## log(bayesian) vs t test
 332     print("log(bayesian) vs t test")
 333     temp.data <- log(na.omit(cbind(data[,"t.test.p.log"],
 334                                    data[,"p.log.bayes"])))
 335     lmline <- lm(temp.data[,1]~temp.data[,2])
 336     temp <-
 337       plotting.function(temp.data[,2],temp.data[,1],
 338                         xlab = "log(p value) of Bayesian (log transformed)",
 339                         ylab = "log(p value) of t test (log transformed)",
 340                         main = "Correspondance of t-test and bayesian p values of log transformed values"
 341                         )
 342     if (use.grid) {
 343       pushHexport(temp$plot.vp)
 344       panel.abline(lmline)
 345     }
 346     else {
 347       abline(lmline)
 348     }
 349     xlim <- range(temp.data[,2])
 350     ylim <- range(temp.data[,1])
 351     text.function(0.5*(xlim[2]-xlim[1])+xlim[1],
 352          0.9*(ylim[2]-ylim[1])+ylim[1],
 353          substitute(r**2 == y,
 354                     list(y=format(summary(lmline)$r.squared,
 355                                     digits=4)))
 356          )
 357
 358     if (use.grid)
 359       popViewport(1)
 360
 361
 362     ## histograms
 363     print("histograms")
 364     plot(hist(data[,"t.test.p"],breaks=40,plot=FALSE),main="Histogram of T Test p values")
 365     plot(hist(data[,"t.test.p.log"],breaks=40,plot=FALSE),main="Histogram of T Test p values for log transformed data")
 366     plot(hist(data[,"p.bayes"],breaks=40,plot=FALSE),main="Histogram of Bayesian p values")
 367     plot(hist(data[,"p.log.bayes"],breaks=40,plot=FALSE),main="Histogram of Bayesian p values for log transformed data")
 368
 369   }
 370