CXE1graphics.R

Code snipits and programs written in R, S or S-PLUS library(lattice) library(grid)

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

showDiagnostics <- F draft <- F

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. Creating factors

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. Normalization of realtime RT-PCR

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

1. Diagnostics

if( showDiagnostics) {

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

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

}

1. Draft plots

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

1. RE versus Time in the Run:Time:Treatment stratum of variation

if( draft ) {

 xyplot(Quantity ~ Time | Treatment * Run, data=dset, aspect=1,
        panel=function(x,y){
          panel.loess(x,y, span=1)
          panel.xyplot(x,y)
        })

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

1. 2) -------------------------------- Figures --------------------------------- #

xlab <- "Time post feeding (days)" figx <- 0.2 figy <- 0.97 myformula <- RE ~ Time | Treatment*Run

plot1 <- xyplot(myformula, data=dset, aspect=1, col=1, groups = dset$Run, between=list(x=0.2,y=0.2),

               ylab="CXE1 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,subscripts,groups,...){
                 panel.loess(x,y, span=1.4)
                 panel.xyplot(x,y,...)
                 panel.grid(col=1, lty=3, lwd=0.3)
               })

if( draft ) { print(plot1) }

1. 3) ------------------------- Creating ratio data ---------------------------- #
2. Taking 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] tmp$Time <- as.numeric(as.vector(tmp$Time))

1. Check...

run <- "Expt 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. ===============================================================================

myformula <- ratio ~ Time |Run plot2 <- xyplot(myformula , 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.39)
                                         panel.xyplot(x,y,...)
                                         panel.abline(h=1, lty=2)
                                         panel.grid(..., lty=3, lwd=0.3) })

print(plot2) if( draft ) { print(plot2) } dev.off()

1. ------------------------- Final figure for paper ---------------------------- #

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

if( draft ) {

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

} 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")))

   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()
 } else {
   png(file.path(directory, "CXE1.png"), width=8*dpi,height=8*dpi)
   tweaky <- 0.985     # for png
   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)

1. CXElsd is from the ANOVA function!

plot3 <- xyplot(RE ~ Time, aspect=1, data=tmp, pch=as.numeric(tmp$Treatment), ylab="Mean Relative Expression", xlab=xlab,

               scales=list(x=list(at=0:7, labels=0:7)), key = mykey,
               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")
                 }
               })
               

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