[BioC] normalize.AffyBatch.vsn with reference vsn object

[R Tagett]

Hi Wolfgang,
Thanks for your response. I found the problem thanks to looking into your suggestion. 
It is not because of the subsampling argument in the vsn2 call. It is because the affybatch and the reference fit sent to normalize.AffyBatch.vsn both contain all pm and mm probes, but inside the function, the abatch is automatically subsetted to use only pm probes, but the reference fit is NOT subsetted as well. So I just removed the "[ind, ]" and it works.
Below is the function. You can see the call to vsn2 uses "intensity(abatch)[ind, ]":
normalize.AffyBatch.vsn
function (abatch, reference, strata = NULL, subsample = if (nrow(exprs(abatch)) > 
    30000L) 30000L else 0L, subset, log2scale = TRUE, log2asymp = FALSE, 
    ...) 
{
    if (is.na(log2scale) || is.na(log2asymp) || (log2scale && 
        log2asymp)) 
        stop("'log2asymp' and 'log2scale' must not both be TRUE, and not be NA.")
    ind = indexProbes(abatch, "pm")
    if (!missing(subset)) 
        ind = ind[subset]
    ind = unlist(ind)
    vsn2res = vsn2(intensity(abatch)[ind, ], reference = reference, 
        returnData = FALSE, subsample = subsample, ...)
    description(abatch)@preprocessing = c(description(abatch)@preprocessing, 
        list(vsnReference = vsn2res))
    trsfx = predict(vsn2res, newdata = intensity(abatch), log2scale = log2scale)
    if (log2asymp) 
        trsfx = (trsfx - log(2))/log(2)
    intensity(abatch) <- 2^trsfx
    return(abatch)
}




 

----- Original Message -----
From: "Wolfgang Huber" <whuber at embl.de>
To: bioconductor at r-project.org
Sent: Tuesday, May 8, 2012 11:20:35 AM
Subject: Re: [BioC] normalize.AffyBatch.vsn with reference vsn object

Dear Rebecca

thanks. It's been a while since I have seen anyone use this function in 
this mode of use!
I haven't tested this answer, but I suspect that the problem is related 
to the subsampling that vsn2 does by default for AffyBatch objects. What 
happens if in your first call to vsn, you do:

  fit<- vsn2(ab, subsample=0)

Best wishes
	Wolfgang


R Tagett scripsit 05/08/2012 02:33 AM:
> Hello,
> I am trying to normalize an affybatch (U133plus2) against a reference set (all U133plus2) using
> normalize.AffyBatch.vsn(abatch, reference, ...)
> according to the vignette where the reference is a vsn object from a previous fit
> and abatch is an affybatch, but the objects do not seem compatible.
> I hope that you can help.
> A simplified script is shown.
> Thank you,
> Rebecca T
>
>> library(affy)
>> library(vsn)
>> ab<- read.affybatch(filenames= c("GSM467031.CEL.gz","GSM467033.CEL.gz","GSM467035.CEL.gz")   )
>> fit<- vsn2(ab)
> vsn2: 1354896 x 3 matrix (1 stratum). Please use 'meanSdPlot' to verify the fit.
>>
>> # affybatch with one cel in it
>> ab1<- read.affybatch(filenames=c( "GSM467037.CEL.gz"   ))
>
>> normAb<- normalize.AffyBatch.vsn(ab1, reference=fit)
> Error in vsnMatrix(as.matrix(x, ncol = 1), reference, strata, ...) :
>    'nrow(reference)' must be equal to 'nrow(x)'.
>>
>>
>> sessionInfo()
> R version 2.14.0 (2011-10-31)
> Platform: x86_64-redhat-linux-gnu (64-bit)
>
> locale:
>   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>   [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>   [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>   [7] LC_PAPER=C                 LC_NAME=C
>   [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] hgu133plus2cdf_2.9.1  AnnotationDbi_1.16.19 vsn_3.22.0
> [4] affy_1.32.1           Biobase_2.14.0
>
> loaded via a namespace (and not attached):
>   [1] affyio_1.22.0         BiocInstaller_1.2.1   DBI_0.2-5
>   [4] grid_2.14.0           IRanges_1.12.6        lattice_0.20-6
>   [7] limma_3.10.3          preprocessCore_1.16.0 RSQLite_0.11.1
> [10] tools_2.14.0          zlibbioc_1.0.1
>>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor


-- 
Best wishes
	Wolfgang

Wolfgang Huber
EMBL
http://www.embl.de/research/units/genome_biology/huber

_______________________________________________
Bioconductor mailing list
Bioconductor at r-project.org
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor