CXE1graphics.R
library(lattice) library(grid)
- =============================== Setup ========================================
- ================ 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
- 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"))
- ==============================================================================
- ============================ Normalization ===================================
- Normalization of realtime RT-PCR
dset$RE <- dset$Quantity/dset$NormFactor
- ==============================================================================
- ============================ Diagnostics =====================================
- 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()
- ==============================================================================
- ========================= 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
- 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 myformula <- RE ~ Time | Treatment plot1 <- xyplot(myformula, data=dset, aspect=1, groups = 1:3, col=1, 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,...){
panel.loess(x,y, span=1.4)
panel.xyplot(x,y,...)
panel.grid(..., lty=3, lwd=0.3)
})
if( draft ) { print(plot1) }
- suc <- dset[dset$Treatment=="Sucrose",]
- est <- dset[dset$Treatment=="CXE1 dsRNA",]
- combined <- cbind(suc[,c(1,2,4,5,12)], est[,c(1,2,4,5,12)])
- combined$ratio <- combined[,10]/combined[,5]
- ============================= Calculate ratio means ===========================
- 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]
- 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]
- 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"])
- ===============================================================================
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()
- ==============================================================================
- ========================= 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)
}
}
- 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)
- print(plot2, split=c(2,1,2,1), position=c(-0.15,0+(1-tweaky),1,1*tweaky), more=F)
dev.off()
- Put graphs on one page for RE and ratio information
- ==============================================================================
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)
- 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()
- LOOKS LIKE REDUNDANT CODE HERE...
- 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()
