Dear Bioconductor community,

recenltly i started analyzing a dataset which has annotation: "primeview", which is a perfect-match-only(PM) array. Thus, i cant use filtering on present/absent calls with mas5calls function from the mas5algorithm, to filter based on absent/present calls. Is there another alternative methodology based on intensity but not on variance for filtering prior to statistical inference, regarding PM only arrays ? I found a paper by Wu, Z. and R. A. Irizarry (2005) about half-price, but i couldnt find any vignette or tutorial about the code or anything.

My second and also important question refers to a methodology i used both on HG-U133a & also HGU-133aplus2 microarray platform datasets, after normalization with GCRMA(as they both have PM & MM probes. More specifically, for instance(about the hgu133a)

eset <- gcrma(cancer_data) # cancer_data my affybatch object

mas5.eset <- mas5calls(cancer_data)

call.matrix <- exprs(mas5.eset)
AbsentCalls.index <- 
( ifelse(call.matrix == 'A',1,0) | ifelse(call.matrix == 'M',1,0) )
AbsentCalls.index <- which(AbsentCalls.index == TRUE)
exprs(eset)[AbsentCalls.index] <- NA
filtered.eset <- eset

My basic consern is (because im new in r) if with the above code i "mark" correctly as NA absent & present calls, because i use a different normalization algorithm than mas5 ??

dear svlachavas,

concerning your first question about filtering on intensity, what you could do is to plot the density of the probe set expression values after pre-processing and use the mode (i.e. highest peak of the distribution) as a value for a intensity based filter. example:

eset <- rma( abatch )
dens <- density( exprs( eset )[ , 1 ] )
mode.x <- dens$x[ which.max( dens$y ) ]
## plot it
plot( dens, xlab=expression( log[2]~expression ) )
abline( v=mode.x, col="red" )

you could then use this value to define which probe sets are very low or not at all expressed (assuming that less than ~50% of genes are expressed in the analyzed tissue).

AbsentOrLow.idx <- which( exprs( eset )[ , 1 ] < mode.x )

Regarding your second question: I don't think it is wise to combine data from different pre-processing approaches. To start with, I don't think you should use the mas5 pre-processing (plenty of papers out there) at all and you should also be careful with the absent call, as about 30% of MM probes yield higher signals then their PM counterpart [Naef and Magnasco, Phys. Rev. 2003, vol. 68 (1)]. If at all, I would use the approach from above, specifically so for the GCRMA pre-processing as you usually see there a bi-modal distribution of probe set expression values, with the lower peak most likely corresponding to non- or very low expressed genes.

hope this helps,

cheers, jo

The xps package has a function dabg.call(), which implements Affy's 'Detection above background' measure, which is what you should probably use if you want something like the P/A calls that you used to be able to get with the PM/MM arrays.

The method you use in your first post does something that I wouldn't recommend. I have a hard time coming up with a good rationale for taking data and converting to NA. This just makes downstream processing of the data more difficult, and may cause things to happen that you didn't expect. Instead, you should think of your data in terms of probesets over all samples, rather than individual probesets per sample.

What I mean by that is you should be looking for probesets where you have little or no evidence for expression over all the samples, and excluding those probesets from consideration, rather than picking and choosing probesets by subject that you don't like. As an example, say you have three groups of 5 subjects (treated, control, vehicle), and you find that probeset X is below detection in the control and vehicle but really high in treated. If you convert all the probeset data to NA for the control and vehicle, you now cannot make any comparisons, and you have effectively thrown that probeset out for no good reason.

Instead, I would recommend doing something like

ind <- rowSums(call.matrix == "P") > 4

eset <- eset[ind,]

In which case you are unbiasedly choosing probesets that appear to be present in at least 5 samples

The error you see with simpleaffy is saying that there isn't a set of QC parameters in the simpleaffy package for this array. See section 5.4 in the simpleaffy vignette (which you should be reading anyway, if you are using that package) http://bioconductor.org/packages/release/bioc/vignettes/simpleaffy/inst/doc/simpleAffy.pdf

You will likely get other errors with simpleaffy, as it was originally designed for PM/MM arrays, and IIRC uses MM measures as part of its processing. If so you might consider looking at arrayQualityMetrics.