[BioC] tiling array normalization

[Wolfgang Huber]

Dear Lana,

> I have implemented your normalization code from your vignette with some of our data.  I would like to 
> demonstrate your normalization results by looking at genes RPN2 and SER33.  After calling your
> function comparisonPlot, I am not seeing the "ps" file. 

Where did you get the idea that you should see a "ps" file? Neither the 
man page of the function nor any other place I am aware of suggest that 
it should produce one. It will create a plot on the current graphics 
device, using grid graphics. Note that for any R graphics going e.g. to 
a postscript device, you will need to call "dev.off()" before using the 
postscript file.

 > Would you be able to help figure out why I can't
> get the postscript file? 

Please send a code example that I (and others) can run. Otherwise I will 
not be able to help. The first line in your code example produces:

 > VY=readCel2eSet(fns,path=celpath)
Error in unique(c("AsIs", oldClass(x))) : object "fns" not found

> (I am not familiar with the function of "[" for lapply.)

lapply(dat,"[",sel) takes list dat and produces a new list where each 
element has been subset using indices "sel". This functionality of R is 
quite unrelated to writing postscript files.

  Best wishes
  Wolfgang

> R version 2.4.0 (2006-10-03) 
> i686-redhat-linux-gnu 
> 
> locale:
> LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
> 
> attached base packages:
>  [1] "splines"   "grid"      "tools"     "methods"   "stats"     "graphics" 
>  [7] "grDevices" "utils"     "datasets"  "base"     
> 
> other attached packages:
>  davidTiling           GO  tilingArray       pixmap  geneplotter     annotate 
>      "1.2.1"     "1.10.0"     "1.12.0"      "0.4-5"     "1.12.0"     "1.12.1" 
>   genefilter     survival          vsn  strucchange     sandwich          zoo 
>     "1.12.0"       "2.29"     "1.12.0"      "1.3-1"      "2.0-1"      "1.2-2" 
> RColorBrewer         affy       affyio      Biobase 
>      "0.2-3"     "1.12.2"      "1.2.0"     "1.12.2" 
> 
> VY=readCel2eSet(fns,path=celpath)
> library(davidTiling)
> data("probeAnno")
> whPM = PMindex(probeAnno)
> whBG = BGindex(probeAnno)
> length(whPM)
> length(whBG)
> all(whBG %in% whPM)
> VY$nucleicAcid <- c("DNA","total RNA","total RNA","total RNA","total RNA")
> isDNA = VY$nucleicAcid %in% "DNA"
> isRNA = VY$nucleicAcid %in% "total RNA"
> pfn = sprintf("assessNorm-normalize%d.pdf", seq(along = which(isRNA)))
> xn2 = normalizeByReference(VY[,isRNA], VY[,isDNA], 
>   pm=whPM, background=whBG, plotFileNames=pfn)
> #  found the plotFiles in the directory
> xn1 = normalizeByReference(VY[,isRNA],VY[,isDNA],pm=whPM,background=whBG,cutoffQuantile=0)
> sta = probeAnno$"9.-.start"
> end = probeAnno$"9.-.end"
> ind=probeAnno$"9.-.index"
> dat = vector(mode = "list", length = 5)
> dat[[1]] = log2(exprs(VY)[ind, which(isDNA)[1]])
> dat[[2]] = log2(exprs(VY)[ind, which(isRNA)[1]])
> dat[[3]] = dat[[2]] - dat[[1]]
> dat[[4]] = exprs(xn1)[ind, 1]
> dat[[5]] = exprs(xn2)[ind, 1]
> dat[[6]] = exprs(xn2)[ind, 1]
> for (j in 3:length(dat)) dat[[j]] = dat[[j]] - quantile(dat[[j]],0.05, na.rm = TRUE)
> names(dat) = letters[seq(along=dat)]
> sel = (sta >= 216600 & end <= 227000)
> ysc = sapply(dat, function(py) quantile(py, probs = c(0,1),na.rm=TRUE))
> ysc[,3:6] = c(-3,8)
> anno = data.frame(start=c(217860, 221078),end =c(220297,222487),name=I(c("RPN2","SER33")))
> ticks = c(217,223,224,225,226)
> comparisonPlot((sta+end)[sel]/2,lapply(dat,"[",sel),yscale=ysc,anno=anno,ticks=ticks,cex=0.2)
> 
> Lana Schaffer
> Biostatistics/Informatics
> The Scripps Research Institute
> DNA Array Core Facility
> La Jolla, CA 92037
> (858) 784-2263
> (858) 784-2994
> schaffer at scripps.edu
> 
> 
> 


-- 
------------------------------------------------------------------
Wolfgang Huber  EBI/EMBL  Cambridge UK  http://www.ebi.ac.uk/huber