CXE1graphics.R

From Organic Design wiki
Revision as of 03:23, 14 February 2006 by Sven (talk | contribs) (Code for paper)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

library(lattice) library(grid)

  1. =============================== Setup ========================================
  2. ================ Read in esterase realtime RT-PCR data =======================

folder <- "/Volumes/HD2/Clinton/Data" filename <- "esteraseRTcompare4runs2.txt" dset <- read.table(file.path(folder,filename), header=T, sep="\t", as.is=T) dset$Time <- dset$Time / 24

  1. dset$Time <- dset$Time

dset$Treatment <- factor(dset$Treatment, levels=c("Sucrose","Esterase"), 

                                             labels=c("Sucrose","CXE1 dsRNA"))

dset$Run <- factor(dset$Run, levels=unique(dset$Run)[c(3,2,4,1)], 

                                   labels=c("Expt 4","Expt 3","Expt 2","Expt 1"))
  1. ==============================================================================


  1. ============================ Normalization ===================================
  1. Normalization of realtime RT-PCR

dset$RE <- dset$Quantity/dset$NormFactor

  1. ==============================================================================


  1. ============================ Diagnostics =====================================
  1. Details

dim(dset) summary(dset) table(dset$Time, dset$Treatment) table(dset$Run, dset$Time, dset$Treatment)

quartz() par(mfrow=c(1,3)) hist(dset$RE) plot(density(dset$RE)) plot(sort(dset$RE), pch=".") dev.off()

  1. ==============================================================================


  1. ========================= Trellis Graphs (draft) =============================

draft <- F

if( draft ) { plot(RE ~ Treatment, dset) }

if( draft ) { histogram(~ Quantity | Treatment*Time, data=dset) histogram(~ RE | Treatment*Time, data=dset) }

frac <- 0.8

  1. Graph Presentation Re versus Time in the Run:Time:Treatment stratum of variation

if( draft ) { quartz(width=frac * 210/25.4, height=frac* 297/25.4) xyplot(Quantity ~ Time | Treatment * Run, data=dset, aspect=1, 

           panel=function(x,y){
              panel.loess(x,y, span=1)
              panel.xyplot(x,y)
                                              })

quartz(width=frac * 210/25.4, height=frac* 297/25.4) xyplot(RE ~ Time | Treatment * Run, data=dset, aspect=1, 

           panel=function(x,y){
              panel.loess(x,y, span=1)
              panel.xyplot(x,y)
                                              })

dev.off() dev.off() }

xlab <- "Time post feeding (days)" figx <- 0.2 figy <- 0.97

plot1 <- xyplot(RE ~ Time | Treatment * Run, data=dset, aspect=1, groups = 1:3, col=1, between=list(x=0.2,y=0.2), 

               ylab="Relative Expression", xlab =  xlab,
               page = function(n) 
               grid.text(paste("A"),
                         x = figx, y = figy, gp=gpar(fontsize=20),
                         default.units = "npc", 
                         just = c("right", "bottom")),
               panel=function(x,y,...){
                 panel.loess(x,y, span=1.4)
                 panel.xyplot(x,y,...)
                 panel.grid(..., lty=3, lwd=0.3)
               })

if( draft ) { print(plot1) }

  1. suc <- dset[dset$Treatment=="Sucrose",]
  2. est <- dset[dset$Treatment=="CXE1 dsRNA",]
  3. combined <- cbind(suc[,c(1,2,4,5,12)], est[,c(1,2,4,5,12)])
  4. combined$ratio <- combined[,10]/combined[,5]
  1. ============================= Calculate ratio means ===========================
  1. 1) Take means over Time:Run:Treatment

colsToRemove <- c("Run","Time","Treatment","Sample.Name", "Detector.Name", 

                 "Reporter","Task","Well")

dsetMeans <- aggregate(dset[,-match(colsToRemove, names(dset))], list(Time=dset$Time, 

                                                       Run=dset$Run,
                                                       Treatment=dset$Treatment),
                      mean)

tmp <- dsetMeans[1:20, 1:2]

  1. Bug in here factor reversed after aggregate...

tmp$Run <- factor(tmp$Run, levels = unique(tmp$Run))

tmp$ratio <- dsetMeans$RE[21:40] / dsetMeans$RE[1:20]

  1. check...

run <- "Run 1" times <- 7 mean(dset$RE[dset$Time==times & dset$Run==run & dset$Treatment=="CXE1 dsRNA"]) / 

 mean(dset$RE[dset$Time==times & dset$Run==run & dset$Treatment=="Sucrose"])
  1. ===============================================================================

plot2 <- xyplot(ratio ~ Time | Run, data=tmp, layout=c(1,4), pch=4, 

               aspect=1, groups = 1:3, col=1,
               between=list(x=0.2,y=0.2),
               ylab=expression(over("CXE1 dsRNA","Sucrose")),
               xlab=xlab,
               page = function(n) grid.text(paste("B"), x = figx+0.15, y = figy+0.01, gp=gpar(fontsize=20),
                 default.units = "npc", just = c("right", "bottom")),
               panel=function(x,y, ...){ panel.loess(x,y, span=1.4)
                                         panel.xyplot(x,y,...)
                                         panel.abline(h=1, lty=2)
                                         panel.grid(..., lty=3, lwd=0.3) })

if( draft ) { print(plot2) }

dev.off()

  1. ==============================================================================
  1. ========================= Trellis Graphs (final) =============================

draft <- F directory <- "/Volumes/HD2/Clinton/Graphs" printType <- "postscript" dpi <- 1024 / ((4/5)*17) if(!printType=="postscript") { 

 trellis.par.set("background", list(col=0))
 trellis.par.set("strip.background", list(alpha=0, col="gray90"))

}


if( draft ) { 

 quartz(width=8,height=8)

} else { 

 if(printType=="postscript") {
   sizemm <- 80
   names(postscriptFonts()) # Available fonts
   trellis.device ("postscript", color=T, file=file.path(directory, "CXE1.eps"),width=sizemm, height=sizemm, horizontal=F) #, font="serif")
   trellis.par.set("strip.background", list(alpha=0, col=c("gray80")))
 } else {
   png(file.path(directory, "CXE1.eps"), width=8*dpi,height=8*dpi)
 }

}


  1. tweaky <- 0.985 # for png

tweaky <- 0.85 # for postscript print(plot1, split=c(1,1,2,1), position=c(0,0,1,1), more=T) print(plot2, split=c(2,1,2,1), position=c(-0.15,0+(1-tweaky),1,1*tweaky), more=F) 

 dev.off()
  1. Put graphs on one page for RE and ratio information
  1. ==============================================================================

REmeans <- tapply(dset$RE, paste(dset$Time, dset$Treatment, sep=":"), mean, na.rm=T) times <- as.numeric(substr(names(REmeans), 1,1)) treatment <- substr(names(REmeans), 3, nchar(REmeans))

tmp <- data.frame(Time=times, Treatment=treatment, RE = REmeans)

lerrorbars <- function(x, y, se, yl = y - se, yu = y + se, eps, ...) { 

   if(any(is.na(yl)) | any(is.na(yu)))
      ltext(x[is.na(yl) | is.na(yu)], y[is.na(yl) | is.na(yu)], "NA", adj = 1)
   lsegments(x, yl, x, yu, ...)
   lsegments(x - eps, yl, x + eps, yl, ...)
   lsegments(x - eps, yu, x + eps, yu, ...)

}

trellis.device("quartz") source("CXE1Anovas.R")

mykey <- simpleKey(text=levels(tmp$Treatment), space="top", pch=1:2, columns=2) mykey$points <- list(alpha=1, cex=0.8, col="black", font=1, pch=1:2)

xyplot(RE ~ Time, data=dset, col=as.numeric(dset$Treatment), key = simpleKey(c("foo", "fodda"), points = F))

trellis.device ("postscript", color=T, file = file.path(directory, "Esteraseinteraction.eps"), width=169, height=150, horizontal=F) plot3 <- xyplot(RE ~ Time, aspect=1, data=tmp, pch=as.numeric(tmp$Treatment), 

               ylab="Mean Relative Expression",

=xlab, 

               scales=list(x=list(at=0:7, labels=0:7)),
               panel=function(x,y, subscripts) {
                 panel.superpose(x,y, subscripts, groups=as.numeric(tmp$Treatment), pch=as.numeric(tmp$Treatment), col=1)
                 ltext(5, 0.45, "99% lsd")
                 lerrorbars(4, 0.45, CXElsd/2, eps = 0.1)   # lsd comes from CXE1Anovas.R
                 for( i in levels(tmp$Treatment)) {
                   panel.lines(x[tmp$Treatment==i],y[tmp$Treatment==i], col="black")
                 }
               }
               , key = mykey)

print(plot3) dev.off()

  1. LOOKS LIKE REDUNDANT CODE HERE...
  2. simpleKey(text=levels(tmp$Treatment), space="top", pch=1:2, columns=2

trellis.device("postscript") a <- xyplot(log2(RE) ~ Time, aspect=1, data=dset, col=as.numeric(dset$Treatment), ylim=log2(c(0.01,1)), xlim=-1:8) print(a, more=T) b <- xyplot(log2(RE) ~ Time, aspect=1, data=tmp, pch=as.numeric(tmp$Treatment), col="orange", cex=4, ylim=log2(c(0.01,1)), xlim=-1:8, 

           ylab="Relative expression", xlab=xlab,
           scales=list(x=list(at=0:7, labels=0:7)))

print(b, more=T) dev.off()

Retrieved from "https://organicdesign.nz/wiki/index.php?title=CXE1graphics.R&oldid=10388"