Dear BioC users and developers: I'm working with Codelink plattform microarray data and i want to apply some of the knowledge currently available for other platforms as Affymetrix o cDNA spotted microarrays to deal with this data. I have been using Affymetrix data for a while and i used suscessfully the great affy package. I also worked with spoted cDNA data and use the also great package limma mainly. Currently this is what I have done: 1. Write a parser for exported text data from the codelink software. 2. I found convenient to create a new class for storing data from this platform because there are no unique identifiers as in Affymetrix and so an exprSet object was not appropiate. I think an RGlist/MAlist-like definition could do the job fine so I almost cut and paste the classes definition in the limma package an modified it a little (any concerns about using other's code is wellcome. I'm willing to release this code under LGPL licence) to make a currently named "Codelink" class. The class currently stores signal values, flag values, type of spot and codelink identifiers. Also I make some annotation packages for the human10k, humanwholegenome and ratwholegenome bioarrays but i have to test it deeper. A remember and older thread in the BioC list about modifiying the exprSet definition to deal with this problem. Is there any plan to change it so I can use a more standard approach better that create a new class? 3. Created/adapted some functions for plotting: Density plots, MAplot, Scatter plot that let see the position of diferent type of spots in codelink chips: that is FIDUCIAL, DISCOVERY, POSITIVE, etc. All that worked more or less. Data can be imported as RAW data (default) or normalized (Codelink median normalization) and as log2 values (default) or not. When I calculate log2 values, i initially decided to set negative values (obtained from subtracting background to signal sometimes) to somethign like 0.01. This works fine with normalize.quantiles and normalize.loess from affy package but create artifact point in a MAplot and a low intensity peak in a density plot. Then, I saw a post from Gordon Smyth (sorry dont have link) telling that limma make negative values NA prior to log2. I wanted to do that and modified the functions involved. Finaly this is when I found a mayor problem: 1) normalize.loess: I got to modify this function to allow NA values. The function compares each array whith the others so if you have 10 arrays then each array is compared and modified (after loess estimation) 10 times. If in each iteration I select the non NA values created when you calculate A or M then in each iteration you have a different subset of genes used in the normalization process. As a little example, if you have 4 arrays with 5 genes: chip1 chip2 chip3 chip4 gen1 1 10 NA 5 gen2 3 NA 40 6 gen3 4 NA NA 7 gen4 3 12 45 6 gen5 2 1 4 3 normalize.loess take a matrix as argument and do all pairwise comparisons: chip1-chip2: used gen1,gen4,gen5 to estimate loess curve. chip1-chip3: used gen2,gen4,gen5 to estimate loess curve. chip1-chip4: used all genes to estimate loess curve. chip2-chip3: used gen4,gen5 to estimate loess curve. chip2-chip4: used gen1,gen4,gen5 to estimate loess curve. chip3-chip4: used gen2,gen4,gen5 to estimate loess curve. The code modified/added to allow for NA is between #: (from normalize.loess in affy 1.5.8) ---------------------------------------------------------------------- -- normalize.loess <- function (mat, subset = sample(1:(dim(mat)[1]), min(c(5000, nrow(mat)))), epsilon = 10^-2, maxit = 1, log.it = TRUE, verbose = TRUE, span = 2/3, family.loess = "symmetric") { J <- dim(mat)[2] II <- dim(mat)[1] newData <- mat if log.it) { mat <- log2(mat) newData <- log2(newData) } change <- epsilon + 1 fs <- matrix(0, II, J) iter <- 0 w <- c(0, rep(1, length(subset)), 0) while (iter < maxit) { iter <- iter + 1 means <- matrix(0, II, J) for (j in 1:(J - 1)) { for (k in (j + 1):J) { y <- newData[, j] - newData[, k] x <- (newData[, j] + newData[, k])/2 ##################################### # Select genes that are not set to NA sel <- which(!is.na(as.character(y))) y <- y[sel] x <- x[sel] ##################################### index <- c(order(x)[1], subset, order(-x)[1]) xx <- x[index] yy <- y[index] aux <- loess(yy ~ xx, span = span, degree = 1, weights = w, family = family.loess) aux <- predict(aux, data.frame(xx = x))/J ##################################### # Apply correction to genes not NA. means[sel, j] <- means[sel, j] + aux means[sel, k] <- means[sel, k] - aux ##################################### if (verbose) cat("Done with", j, "vs", k, " in iteration ", iter, "\n") } } fs <- fs + means newData <- mat - fs change <- max(colMeans((means[subset, ])^2)) if (verbose) cat(iter, change, "\n") oldfs <- fs } if ((change > epsilon) & (maxit > 1)) warning(paste("No convergence after", maxit, "iterations.\n")) if log.it) { return(2^newData) } else return(newData) } ---------------------------------------------------------------------- -- I'm not an statitician so not sure if this could be a correct approach: Is this correct or I cannot do that? 2) With normalize.quantiles all seems to work fine... the function don't say nothing about NA values and the density plots seems correct (almost same density plot for all chips) but the MAplot have changed greatly with a great number of genes with greater values of M (in absolute value) that is, a wider cloud of points. I don't know why is this happening. I also tried to see at the C code for normalize.quantiles (at qnorm.c) but not found an easy answer. Hopefully I finally found that I can use normalizeQuantiles from limma package that allows for NA values. By the way, is this really necesary or can I assign a low value to negative ones?. I think it is better to use NA but... Something I'm going to do sonner is to import all the values obtained from the Codelink data (mean foreground and background, median foreground and background, etc) so it could be possible to calculate a signal value within R and then use a method that avoids negative values. But that could be the same in some cases to assign a small value (Or I'm wrong) I'm working with a linux box (Debian Sarge) with R 2.0.1 compiled from sources and BioC 1.5. Any comment on all that would be very appreciated. Thanks in advanced! D.