################################################### ### chunk number 1: ged ################################################### library(dressCheck) if (!exists("c119")) data(c119) # pure RMA on CEL files, with trimming of sample names if (!exists("DrAsGiven")) data(DrAsGiven) # read of corrected platinum XLS to CSV # some names are not preserved between two images setdiff(sampleNames(c119), sampleNames(DrAsGiven)) # use the common ones okn = intersect(sampleNames(DrAsGiven), sampleNames(c119)) # the corrected XLS does not have all genes, so c119 needs trim c119r = c119[ featureNames(DrAsGiven), ] # now demonstrate allc1 = sapply(okn, function(i) cor(exprs(DrAsGiven)[,okn[1]], exprs(c119r[,i]))) png(file="corChk.png") par(mfrow=c(2,2)) remap1 = which.max(allc1) plot(exprs(DrAsGiven)[,okn[1]], exprs(c119r)[,okn[1]], xlab=paste("id", okn[1], "in RMA+SFR XLS"), ylab=paste("id", okn[1], "in pure RMA from CEL"), main="(a)") plot(exprs(DrAsGiven)[,okn[1]], exprs(c119r)[,names(remap1)], xlab=paste("id", okn[1], "in RMA+SFR XLS"), ylab=paste("id", names(remap1), "in pure RMA from CEL"), main="(b)") allc2 = sapply(okn, function(i) cor(exprs(DrAsGiven)[,okn[2]], exprs(c119r[,i]))) remap2 = which.max(allc2) plot(exprs(DrAsGiven)[,okn[2]], exprs(c119r)[,okn[2]], xlab=paste("id", okn[2], "in RMA+SFR XLS"), ylab=paste("id", okn[2], "in pure RMA from CEL"), main="(c)") plot(exprs(DrAsGiven)[,okn[2]], exprs(c119r)[,names(remap2)], xlab=paste("id", okn[2], "in RMA+SFR XLS"), ylab=paste("id", names(remap2), "in pure RMA from CEL"), main="(d)") dev.off() ################################################### ### chunk number 2: lk116 ################################################### library(chron) if (!exists("corrp116")) data(corrp116) dt = table(chron(corrp116$rundate)) cdt = chron(as.numeric(names(dt))) names(dt) = cdt dt ################################################### ### chunk number 3: lkte ################################################### library(limma) corrp116$cdate = factor(chron(corrp116$rundate)) des = model.matrix(~cdate, pData(corrp116)) if (!exists("f1")) f1 = lmFit(corrp116, des) ef1 = eBayes(f1) options(digits=4) tt = topTable(ef1,2:16)[,-c(2:16)] bigtt = topTable(ef1,2:16,n=1000)[,-c(2:16)] mm = max(bigtt[,5]) tops = tt[,1] par(mfrow=c(2,2)) plot(chron(corrp116$rundate), exprs(corrp116)[tops[1],], xlab="array rundate", ylab=tops[1]) plot(chron(corrp116$rundate), exprs(corrp116)[tops[2],], xlab="array rundate", ylab=tops[2]) plot(chron(corrp116$rundate), exprs(corrp116)[tops[3],], xlab="array rundate", ylab=tops[3]) plot(chron(corrp116$rundate), exprs(corrp116)[tops[4],], xlab="array rundate", ylab=tops[4])