Lifecourse.R

Code snipits and programs written in R, S or S-PLUS 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, ...)

}

foo <- read.table("/Volumes/HD2/Clinton/Data/lifecoursedataset.txt", sep="\t", header=T, as.is=T)

names(foo) library(lattice) library(grid)

1. foo$Time <- factor(foo$Time, levels=unique(foo$Time))

foo$Sex <- factor(foo$Sex) foo$Body.Part <- factor(fooBody.Part)

Larvae <- foo[foo$Lifestage=="Larvae",] Larvae$Time <- factor(Larvae$Time, levels=unique(Larvae$Time), labels=c("Eggs","1*","2*","3*","4*","5*"))

xlab <- "Lifestage" directory <- "/Volumes/HD2/Clinton/Graphs" figx <- 0.2 figy <- 0.97 plot1 <- xyplot(RE ~ Time|Detector.Name, data=Larvae, aspect=1, between=list(x=0.3), groups = Std.Err, layout=c(2,1),

      ylab="Mean 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.xyplot(x,y, col=1)
        panel.lines(x,y, col=1)

1. panel.loess(x,y, span=0.8)

                                       # 3 replicates 
        se <- sqrt(mean(groups[subscripts]^2))
        lerrorbars(as.numeric(x)[2],0.075, se = se, eps=0.05)
        ltext(as.numeric(x)[2]+1.5, 0.075, "Pooled SEM")

1. lerrorbars(as.numeric(x),y, se = groups[subscripts]/2, eps=0.05)

      })

Pupae <- foo[foo$Lifestage=="Pupae",] Pupae$Time <- factor(Pupae$Time, levels=unique(Pupae$Time)) Pupae$Sex <- factor(Pupae$Sex, levels=unique(Pupae$Sex))

xlab <- "Time post pupation (days)" plot2 <- xyplot(RE ~ Time|Detector.Name*Sex, data=Pupae, ylim = c(-0.1,0.6), aspect=1,between=list(x=0.3, y=0.3),

      ylab="Mean Relative Expression", xlab = xlab,
                page = function(n) 
               grid.text(paste("B"),
                         x = figx, y = figy, gp=gpar(fontsize=20),
                         default.units = "npc", 
                         just = c("right", "bottom")),

groups = Std.Err,

panel=function(x,y, subscripts, groups) {
        panel.xyplot(x,y, col=1)
        panel.lines(x,y, col=1)

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

                                              # 3 replicates 
        se <- sqrt(mean(groups[subscripts]^2))
        lerrorbars(as.numeric(x)[1],0.15, se = se, eps=0.05)
        ltext(as.numeric(x)[1]+1.1, 0.15, "Pooled SEM")

1. lerrorbars(as.numeric(x),y, se = groups[subscripts]/2, eps=0.1)

      })

Adult <- foo[foo$Lifestage=="Adult",] Adult$Time <- factor(Adult$Time, levels=unique(Adult$Time)) Adult$Sex <- factor(Adult$Sex, levels=unique(Adult$Sex)) Adult$Body.Part <- factor(Adult$Body.Part, levels=unique(Adult$Body.Part))

xlab <- "Time post eclosion (days)"

mykey <- simpleKey(text=levels(Adult$Sex), space="top", pch=c(1,4), columns=2) mykey$points <- list(alpha=1, cex=0.8, col="black", font=1, pch=c(1, 4))

plot3 <- xyplot(RE ~ Time|Detector.Name*Body.Part, data=Adult, groups=Sex, between=list(x=0.3, y=0.3), ylab="Mean Relative Expression", xlab=xlab,

               page = function(n) 
               grid.text(paste("C"),
                         x = figx, y = figy, gp=gpar(fontsize=20),
                         default.units = "npc", 
                         just = c("right", "bottom")),
               aspect=1,
             panel=function(x,y, subscripts, groups) {
       #        panel.superpose(x,y, subscripts, groups)
               panel.superpose(x,y, subscripts, groups, pch=c(1,4), col=1)
                                                      # 3 replicates 
               se <- sqrt(mean(Adult[subscripts, "Std.Err"]^2))
               lerrorbars(as.numeric(x)[1],0.2, se = se, eps=0.05)
               ltext(as.numeric(x)[1] + 0.9, 0.2, "Pooled SEM")

               for( i in levels(Adult$Sex)) {
                 panel.lines(x[Adult$Sex[subscripts]==i],y[Adult$Sex[subscripts]==i], col="black")
               }

1. lerrorbars(as.numeric(x), y, se = Adult[subscripts, "Std.Err"]/2, eps=0.1)

             },
      key=mykey)

tweaky <- 0.85

1. print(plot2, split=c(1,2,1,2), position=c(0,0,1.1,1.1), more=T)

1. print(plot1, split=c(1,1,2,1), position=c(0,0,1,1), more=T)
2. print(plot2, split=c(2,1,2,1), position=c(-0.15,0+(1-tweaky),1,1*tweaky), more=F)

trellis.device ("postscript", color=T, file = file.path(directory, "lifestage1.eps"), width=100, height=100, horizontal=F) trellis.par.set("strip.background", list(alpha=0, col=c("gray80"))) print(plot1, position=c(0.1,0.1, 0.9,0.9)) dev.off()

trellis.device ("postscript", color=T, file = file.path(directory, "lifestage2.eps"), width=169, height=150, horizontal=F) trellis.par.set("strip.background", list(alpha=0, col=c("gray80"))) print(plot2) dev.off() trellis.device ("postscript", color=T, file = file.path(directory, "lifestage3.eps"), width=169, height=150, horizontal=F)

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

print(plot3) dev.off()

1. plot(Larvae$RE, Larvae$Std.Err)
2. plot(Pupae$RE, Pupae$Std.Err)
3. plot(Adult$RE, Adult$Std.Err)