Entering edit mode
Alexander C Cambon
▴
40
@alexander-c-cambon-1251
Last seen 10.2 years ago
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