################################################# ## ## ## Plot gene expression profiles from ## ## Milk time series expriments ## ## ## ## Anne de Jong ## ## may 28, 2010 ## ## ## ################################################# # set the work space #setwd("F:/Array_data/MG1363/Mortens milk project/SerieA data - Anne/GeneCyclePlot") setwd("/home/anne") # load functions library("stats") # --------------------------------------------------------------------------- # 1. import the raw data file milkdata <- read.table("/usr/milkts/R_SerieA_Signal_Average_inv.txt", row.names=1, header=1, sep="\t") # 2. import the growth curve growth_data <- read.table("/usr/milkts/growth_serieA.txt", header=1, sep="\t") # 3. import the growth curve pH_data <- read.table("/usr/milkts/pH_serieA.txt", header=1, sep="\t") # 4. Parse parameters # strange 'known!' thing with R is that is counts all the arguments, so inclusing the options :( # R command line: R --vanilla --slave --args SESSIONID dnaA oppA < Plot_Expression_profiles.R # also check if the column name exists countargs = 0 countvalidgenes = 0 genes <- c() for (e in commandArgs()) { countargs = countargs+1 if (countargs>5) { if (any(colnames(milkdata)== e)) { countvalidgenes = countvalidgenes + 1 genes[countvalidgenes] <- e } } } sessionfolder <- paste ("/tmp/milkts/",commandArgs()[5], sep="") filenamePNG = paste(sessionfolder,"/milkts-expression-profiles.png", sep= "") # --------------------------------------------------------------------------- # determine max expression maxexp = 1 for (gene in genes) { maxgene = max(milkdata[gene], na.rm = TRUE) if (maxgene>maxexp) {maxexp = maxgene} } # the output picture format png(file = filenamePNG, width = 1280, height = 1024, units = "px", pointsize = 12, bg = "white") # Set graph properties par(new=FALSE, mfrow = c(2, 1), col=1, lty=1, lwd=2, mar=c(2.2,4,2,5), usr=c(0,25,-2,2), pty="m") ylim = maxexp # plot expression profiles plot(milkdata[,3], col=0, type="o", ylim=c(0,ylim), xlim=c(1,24), ylab="Expression level", main=paste("Expression profile"), xlab =" samples" ) # Draw growth stages as blocks in the background polygon(c(7.5,11.5,11.5,7.5),c(-0.12,-0.12,maxexp+0.15,maxexp+0.15), col=8, border=0) polygon(c(13.5,19.5,19.5,13.5),c(-0.12,-0.12,maxexp+0.15,maxexp+0.15), col=8, border=0) # add position of time points ypos <- maxexp/25; pos <- 1; lab="T1"; x <- c(pos,pos); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 2.5; lab="T2"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 4.5; lab="T3"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 6.5; lab="T4"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 8.5; lab="T5"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 10.5; lab="T6"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 12.5; lab="T7"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 14.5; lab="T8"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 16.5; lab="T9"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 18.5; lab="T10"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 20.5; lab="T11"; x <- c(pos-0.5,pos+0.5); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); pos <- 22; lab="T12"; x <- c(pos,pos); y <- c(ypos, ypos); lines(x, y, type="o", col=1, pch="."); text(pos,ypos/1.8,labels = lab , col=4); # add PCA clusters ypos <- maxexp/1.5 ; pos <- 4 ; lab="1"; text(pos,ypos,labels = lab , cex=5, col=grey(0.5)); pos <- 9.5; lab="2"; text(pos,ypos,labels = lab , cex=5, col=grey(0.9)); pos <- 12.5; lab="3"; text(pos,ypos,labels = lab , cex=5, col=grey(0.5)); pos <- 16; lab="4"; text(pos,ypos,labels = lab , cex=5, col=grey(0.9)); pos <- 20.5; lab="5"; text(pos,ypos,labels = lab , cex=5, col=grey(0.5)); # Draw the gene expressions kleur = 0 xpos <- 23 ysteps <- maxexp / 25 palette(rainbow(countvalidgenes+1)); #palette(heat.colors(countvalidgenes)); #palette(gray(1:countvalidgenes)); for (gene in genes) { kleur = kleur+1 ; lines(milkdata[gene], type="o", ylim=c(0,3), col=kleur); ylabel = maxexp-kleur*ysteps; lines(c(xpos,xpos-0.4),c(ylabel,ylabel), col=kleur); text(xpos+0.1,maxexp-kleur*ysteps, labels = colnames(milkdata[gene]) , adj = c(0,0)); } palette("default") # plot growth curve and an empty expression profile graph plot(growth_data, log="y", type="o", col=0, main=paste("Growth curve"), ylab="CFU" , xlab="Time in minutes") polygon(c(275,375,375,275),c(36,36,2300,2300), col=8, border=0) polygon(c(450,900,900,450),c(36,36,2300,2300), col=8, border=0) lines(growth_data, type="o", col=1) # add position of time points pos <- 103; lab="T1"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 154; lab="T2"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 200; lab="T3"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 250; lab="T4"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 300; lab="T5"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 349; lab="T6"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 400; lab="T7"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 500; lab="T8"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 605; lab="T9"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 803; lab="T10"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 1002; lab="T11"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); pos <- 1116; lab="T12"; x <- c(pos,pos); y <- c(40, 55); lines(x, y, type="o", col=1, pch="."); text(pos,60,labels = lab , col=4); # add PCA clusters ypos <- 500 ; pos <- 140 ; lab="1"; text(pos,ypos,labels = lab , cex=5, col=grey(0.5)); pos <- 325; lab="2"; text(pos,ypos,labels = lab , cex=5, col=grey(0.9)); pos <- 410; lab="3"; text(pos,ypos,labels = lab , cex=5, col=grey(0.5)); pos <- 650; lab="4"; text(pos,ypos,labels = lab , cex=5, col=grey(0.9)); pos <- 1000; lab="5"; text(pos,ypos,labels = lab , cex=5, col=grey(0.5)); # add the pH curve to the current plot par(new=TRUE) plot(pH_data,xaxt="n",yaxt="n",xlab="",ylab="",pch=16, type="o", col=4) axis(4) text(1210, 5.5, "pH", srt = 90, xpd = TRUE, cex=1) text(30, 6.3, "pH", xpd = TRUE, col=4, cex=2) dev.off() # write the file that indicated end of the session filename = paste(sessionfolder,"/sessionstop", sep= "") write(genes, file=filename, append=FALSE) # write the result file resulthtml <- paste (" Expression profile

", sep="") filename = paste(sessionfolder,"/result.html", sep= "") write(resulthtml, file=filename, append=FALSE) # ---------------------------------------------------------------------------