CXE1graphics.R

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()