[BioC] affy: Trouble Creating Custom Processing Methods for summary expression

[Alexander C Cambon]

I am trying to use the "Custom Processing Methods" available in package "affy" to create new expression methods. But for some reason I am getting "NA"'s for the expression values. I start out using 4 cel files with rae230a data. 

> x
AffyBatch object
size of arrays=602x602 features (11330 kb)
cdf=RAE230A (15923 affyids)
number of samples=4
number of genes=15923
annotation=rae230a

## And I proceed by using the example in the "Custom Processing Methods (How To)" vignette on page 4:

 > computeExprVal.huber <- function(probes) {
+  res <- apply(probes, 2, huber)
+ mu <- unlist(lapply(res, function(x) x$mu))
+ s <- unlist(lapply(res, function(x) x$s))
+ return(list(exprs = mu, se.exprs = s))
+ }

>generateExprSet.methods <- c(generateExprSet.methods, "huber")

###Then I use the expresso function... and get an error message (see below)

>eset.bg.huber<-expresso(x, bgcorrect.method="mas",  pmcorrect.method="pmonly",  
+ summary.method="huber")
background correction: mas 
PM/MM correction : pmonly 
expression values: huber 
background correcting...done.
normalizing...done.
Error in match.arg(summary.method, express.summary.stat.methods) : 
        'arg' should be one of avgdiff, liwong, mas, medianpolish, playerout

## So, after looking at the error, I decide to add (perhaps wrongly)  the line

> express.summary.stat.methods<-c(express.summary.stat.methods,"huber")

## and I rerun...

>eset.bg.trm<-expresso(x, bgcorrect.method="mas",  pmcorrect.method="pmonly",  
+  summary.method="huber")


## Now I get:

background correction: mas 
normalization: 
PM/MM correction : pmonly 
expression values: huber 
background correcting...done.
normalizing...done.
15923 ids to be processed
|                    |
|####################|

## This looks reassuring, and I start  thinking I am over the hump. However when I ask for the first few lines, I get all ##"NA"'s for expression ..

 > exprs(eset.bg.huber)[1:10,1:4]
           AK_0A.CEL AK_0B.CEL AK_4A.CEL AK_4B.CEL
1367452_at        NA        NA        NA        NA
1367453_at        NA        NA        NA        NA
1367454_at        NA        NA        NA        NA
1367455_at        NA        NA        NA        NA
1367456_at        NA        NA        NA        NA
1367457_at        NA        NA        NA        NA
1367458_at        NA        NA        NA        NA
1367459_at        NA        NA        NA        NA

###But when I do the same thing, using a "built-in" expresso method (using "x"), I get 

 > exprs(eset44)[1:10,1:4]
            AK_0A.CEL AK_0B.CEL AK_4A.CEL AK_4B.CEL
1367452_at  564.68667 400.93697 476.22935 396.00641
1367453_at  330.41054 281.98687 279.02867 305.24909
1367454_at  158.99028 165.78677 134.08273 151.77046
1367455_at  302.39252 223.56433 192.13621 204.74223
1367456_at  881.07080 931.81819 792.18972 583.09672
1367457_at  139.75233 142.67949 140.01944 153.02465
1367458_at  144.34063 154.16194 176.33189 190.89218
1367459_at 1010.72613 709.47988 886.92512 957.64834
1367460_at   48.19831  51.70533  59.43114  54.13907
1367461_at  213.03203 194.16510 210.99492 182.22481

Any ideas what I am doing wrong?

Alex



Alexander Cambon
Biostatistician
School of Public Health
Dept of Biostatistics and Bioinformatics
University of Louisville

R 2.1.0
Bioconductor 1.6
Windows XP
Dell Optiplex
1 GB RAM