Hello,

I would like to process a set of mouse gene st v1.1 arrays in Oligo with the function fitProbeLevelModel(), in analogy of using the function fitPLM (from the library affyPLM). However, this won't work when i would like to summarize the probes into probesets using a 'custom CDF'. Any pointers on how to solve this would be appreciated (or is it a bug?). Please note that it works as expected when using the 'official' PdInfo package.

Thanks,

Guido

> #Preliminaries
> source("http://bioconductor.org/biocLite.R")
Bioconductor version 3.2 (BiocInstaller 1.20.1), ?biocLite for help
> biocLite("pd.mogene.1.1.st.v1")
> install.packages("http://mbni.org/customcdf/20.0.0/entrezg.download/pd.mogene11st.mm.entrezg_20.0.0.tar.gz")
>
> library(oligo)
> celFiles <- list.celfiles(full.names = TRUE)
> affy.data <- read.celfiles(celFiles, pkgname = "pd.mogene11st.mm.entrezg") #using remapped, up-to-date chip design/annotation.
Loading required package: pd.mogene11st.mm.entrezg
Loading required package: RSQLite
Loading required package: DBI
Platform design info loaded.
Reading in : ./G168_A5_1_WT_COLON_LFD_2.CEL
<<snip>>
>
> x.norm.plm <- fitProbeLevelModel(affy.data)
Error in sqliteSendQuery(con, statement, bind.data) :
  error in statement: no such table: corepm
>
> affy.data
GenericFeatureSet (storageMode: lockedEnvironment)
assayData: 1178100 features, 25 samples
  element names: exprs
protocolData
  rowNames: G168_A5_1_WT_COLON_LFD_2.CEL (25 total)
  varLabels: exprs dates
  varMetadata: labelDescription channel
phenoData
  rowNames: G168_A5_1_WT_COLON_LFD_2.CEL (25 total)
  varLabels: index
  varMetadata: labelDescription channel
featureData: none
experimentData: use 'experimentData(object)'
Annotation: pd.mogene11st.mm.entrezg
>
>
>
>#Works fine when using 'official' PdInfo package...
> affy.data <- read.celfiles(celFiles)
Loading required package: pd.mogene.1.1.st.v1
Platform design info loaded.
Reading in : ./G168_A5_1_WT_COLON_LFD_2.CEL
<<snip>>
> x.norm.plm <- fitProbeLevelModel(affy.data)
Background correcting... OK
Normalizing... OK
Summarizing... OK
Extracting...
  Estimates... OK
  StdErrors... OK
  Weights..... OK
  Residuals... OK
  Scale....... OK
>
>
> sessionInfo()
R version 3.2.2 Patched (2015-11-03 r69595)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets
[8] methods   base     

other attached packages:
 [1] pd.mogene.1.1.st.v1_3.14.1     pd.mogene11st.mm.entrezg_0.0.1
 [3] RSQLite_1.0.0                  DBI_0.3.1                     
 [5] oligo_1.34.2                   Biostrings_2.38.4             
 [7] XVector_0.10.0                 IRanges_2.4.8                 
 [9] S4Vectors_0.8.11               Biobase_2.30.0                
[11] oligoClasses_1.32.0            BiocGenerics_0.16.1           
[13] BiocInstaller_1.20.1          

loaded via a namespace (and not attached):
 [1] affxparser_1.42.0          splines_3.2.2             
 [3] GenomicRanges_1.22.4       zlibbioc_1.16.0           
 [5] bit_1.1-12                 foreach_1.4.3             
 [7] GenomeInfoDb_1.6.3         tools_3.2.2               
 [9] SummarizedExperiment_1.0.2 ff_2.2-13                 
[11] iterators_1.0.8            preprocessCore_1.32.0     
[13] affyio_1.40.0              codetools_0.2-14          
>

It's a bug. There is some code that builds a SQL query and I think the logic is based on an earlier iteration of what the various table names would be. I don't know how long it will take to get a fix, but it's easy enough for you to make changes on the fly to get it to work for now. Just do this:

> library(oligo)
> debug(getProbeInfo)
> x.norm.plm <- fitProbeLevelModel(affy.data)

Then step through the debugger until you see this:

Browse[3]>
debug: if (class(object) == "GenericFeatureSet") {
    probeTypeTable <- paste0(probeType, "feature")
    metaTable <- paste0(target, probeType)
    featureTable <- paste0("featureSet", substr(target, 4, 4))
    fields <- c(sprintf("%s.fid as fid", metaTable), sprintf("%s.man_fsetid as man_fsetid",
        featureTable), sprintf("%s.fsetid as fsetid", metaTable))
    fields <- paste(fields, collapse = ", ")
    sql <- paste("SELECT", fields, "FROM", metaTable, "INNER JOIN",
        featureTable, "USING(fsetid)")
    rm(fields, probeTypeTable, metaTable, featureTable)
} else {
    fields <- unique(c("fid", "man_fsetid", field))
    fields[fields == "fsetid"] <- paste(probeTable, "fsetid",
        sep = ".")
    if (isST | class(object) == "TilingFeatureSet")
        fields[fields == "man_fsetid"] <- paste(probeTable, "fsetid as man_fsetid",
            sep = ".")
    fields <- paste(fields, collapse = ", ")
    sql <- paste("SELECT", fields, "FROM", probeTable, "INNER JOIN featureSet",
        "USING(fsetid)")
    rm(fields)
}

Now you have to pay attention. Keep stepping through until you see this line:

Browse[3]>
debug: fields <- c(sprintf("%s.fid as fid", metaTable), sprintf("%s.man_fsetid as man_fsetid",
    featureTable), sprintf("%s.fsetid as fsetid", metaTable))

Before you hit Enter the next time, you need to do this:

Browse[3]> metaTable <- "mps1pm"
Browse[3]> featureTable <- "featureSet1"

To correct the metaTable and featureTable objects. Now if you hit Enter two more times, you should be able to inspect the SQL query:

Browse[3]>
debug: sql <- paste("SELECT", fields, "FROM", metaTable, "INNER JOIN",
    featureTable, "USING(fsetid)")
Browse[3]>
debug: rm(fields, probeTypeTable, metaTable, featureTable)
Browse[3]> sql
[1] "SELECT mps1pm.fid as fid, featureSet1.man_fsetid as man_fsetid, mps1pm.fsetid as fsetid FROM mps1pm INNER JOIN featureSet1 USING(fsetid)"

And if your SQL query looks like that, just hit 'c' to get the debugger to continue.

Hey Guido,

My bad - it's not a bug. It's a feature. Normally there are a set of 'meta' tables in a pdInfoPackage that map probes to probesets, and they are named things like 'core' 'full', 'extended' etc. But in the case of a GenericPDInfoPackage like what you get from MBNI, they are called 'mps1' -  'mpsX' and it's up to you to know what you want. Luckily for you, the MBNI summarization is just to a single level, so there is just an mps1 table. So here is how you do things:

> dat <- read.celfiles(list.celfiles(), pkgname = "pd.mogene20st.mm.entrezg")
Platform design info loaded.

> eset <- fitProbeLevelModel(dat, target = "mps1")
Background correcting... OK
Normalizing... OK
Summarizing... OK
Extracting...
  Estimates... OK
  StdErrors... OK
  Weights..... OK
  Residuals... OK
  Scale....... OK