[BioC] affy: Trouble Creating Custom Processing Methods for summary expression
Alexander C Cambon
accamb01 at louisville.edu
Sat Jun 11 15:42:58 CEST 2005
Thanks for your help, Ben. The changes you made below work for my data set too. It also helps considerably to load the MASS pacakge ( since I had forgotten to put in the code for the huber function that you had in your vignette).
Alex
>>> Ben Bolstad <bolstad at stat.berkeley.edu> 06/10/05 11:20 PM >>>
This is my code:
library(affy);library(affydata); library(MASS)
data(Dilution) # demonstration data
generateExprVal.method.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")
express.summary.stat.methods<-c(express.summary.stat.methods,"huber")
eset.bg.huber<-expresso(Dilution, bgcorrect.method="mas",
pmcorrect.method="pmonly", summary.method="huber")
exprs(eset.bg.huber)[1:10,]
My results:
> exprs(eset.bg.huber)[1:10,]
20A 20B 10A 10B
1000_at 435.75360 402.08952 411.59458 362.68261
1001_at 91.53807 85.69595 84.98978 77.48265
1002_f_at 176.19728 162.36399 146.33448 152.49209
1003_s_at 238.53642 193.41824 205.91676 209.08518
1004_at 176.02760 159.17141 164.96742 155.93032
1005_at 502.82416 526.42168 452.69040 465.98990
1006_at 76.96024 69.57741 73.88633 63.37821
1007_s_at 605.03563 591.35033 537.21003 546.18388
1008_f_at 545.60434 513.74078 616.34378 574.13054
1009_at 1959.30326 1986.06121 2118.56979 2007.63819
>
> >>> Ben Bolstad <bolstad at stat.berkeley.edu> 06/10/05 6:30 PM >>>
> Try renaming "computeExprVal.huber" to "generateExprVal.method.huber".
> That should make it work. I think summary methods for expresso need to
> have the "generateExprVal.method." prefix. I will fix the vignette
> accordingly.
>
> Ben
>
>
>
>
>
>
> On Fri, 2005-06-10 at 16:52 -0400, Alexander C Cambon wrote:
> > 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
> >
> > _______________________________________________
> > Bioconductor mailing list
> > Bioconductor at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/bioconductor
More information about the Bioconductor
mailing list