PBP1graphics.R

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

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

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

1. dset$Time <- factor(dset$Time)

dset$Treatment <- factor(dset$Treatment, levels = sort(unique(dset$Treatment)),

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

dset$Run <- factor(dset[,1], levels = sort(unique(dset$Run)), labels=c("Expt 2","Expt 1"))

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

1. ============================ Normalization ===================================

1. Normalization of realtime RT-PCR

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

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

n <- 5 xyplot(1 ~ 1,

       page = function(n) 
           grid.text("Some text to add...",
                     x = 0.5, y = 0.99, gp=gpar(fontsize=20),
                     default.units = "npc", 
                     just = c("right", "bottom")))

dev.off()

1. ========================= Trellis Graphs =====================================

draft <- F

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)
              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)
              panel.xyplot(x,y)
                                              })

dev.off() dev.off() }

xlab <- "Time post eclosion (days)" figx <- 0.2 figy <- 0.95

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,  scales=list(x=list(at=1:4, labels=1:4)),
               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=0.8)
                 panel.xyplot(x,y,...)
                 panel.grid(..., lty=3, lwd=0.3)
               })

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

1. Removing dubious rows of data

1. suc <- dset[dset$Treatment=="Sucrose",]
2. est <- dset[dset$Treatment=="PBP1 dsRNA",]
3. combined <- cbind(suc[,c(1,4,5,12)], est[,c(1,4,5,12)])

1. combined$Run <- factor(combined$Run, levels=levels(factor(combined$Run)), labels=c("Run 2","Run 1"))
2. names(combined)[c(4,8)] <-c("Sucrose.RE","PBP1 dsRNA")
3. combined$ratio <- combined[,"PBP1 dsRNA"]/combined[,"Sucrose.RE"]

1. ============================= Calculate ratio means ===========================

1. 1) Take means over Time:Run:Treatment

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

                 "Reporter","Task","Well.ID")

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

                                                       Run=dset$Run,
                                                       Treatment=dset$Treatment),
                      mean, na.rm=TRUE)

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

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

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

tmp$ratio <- dsetMeans$RE[9:16] / dsetMeans$RE[1:8]

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. ===============================================================================

plot2 <- xyplot(ratio ~ Time | Run, data=tmp, pch=4, layout=c(1,2), aspect=1, groups = 1:3, col=1, between=list(x=0.2,y=0.2), ylab=expression(over("PBP1 dsRNA","Sucrose")), xlab= xlab, scales=list(x=list(at=1:4, labels=1:4)),

               page = function(n) 
               grid.text(paste("B"),
                         x = figx+0.15, y = figy+0.02, gp=gpar(fontsize=20),
                         default.units = "npc", 
                         just = c("right", "bottom")),
               panel=function(x,y, ...){

1. panel.loess(x,y, span= 10, degree=1)

                   panel.lines(x,y, col="black")
                   panel.xyplot(x,y,...)
                   panel.abline(h=1, lty=2)
                   panel.grid(..., lty=3, lwd=0.3)
               })

if( draft ) { print(plot2) }

directory <- "/Volumes/HD2/Clinton/Graphs" printType <- "postscript" if(!printType=="postscript") {

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

} dpi <- 1024 / ((4/5)*17)

if( draft ) {

 quartz(width=8,height=8)

} else { if(printType=="postscript") {

 	trellis.device ("postscript", color=T, file = file.path(directory, "PBP1.eps"), width=169, height=169, horizontal=F) # font="serif"
 trellis.par.set("strip.background", list(alpha=0, col="gray80"))	

} else { png(file.path(directory, "PBP1.png"), width=8*dpi,height=8*dpi) } }

1. tweaky <- 0.79 # for png

tweaky <- 0.7 # 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.1,0+(1-tweaky),1,1*tweaky), more=F)

dev.off()

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

1. ======================== Interaction plot ====================================

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

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("PBP1Anovas.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)

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

               ylab="Mean Relative expression", xlab=xlab,
               scales=list(x=list(at=1:4, labels=1:4)),
               panel=function(x,y, subscripts){
                 panel.superpose(x,y, subscripts, groups=as.numeric(tmp$Treatment), pch=as.numeric(tmp$Treatment), col=1)
                 ltext(1.4, 0.4, "99% lsd")
                 lerrorbars(1, 0.4, se=PBPlsd/2, eps=0.05)
                 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. ==============================================================================