I have been attempting to obtain a list of differentially expressed genes from a pairwise comparison of microarray data. I am running into an error when calling the get.annotation function from package: "simpleaffy" because the object I call it on has no probe-sets in the calls slot of the PairComp object. This occurs even though the tt slot which includes t-test values between groups does seemingly have probe-sets. If anyone has any insights about this issue, I would really appreciate it!

> summary <- results.summary(pw.filtered, "hgu133plus2.db")
hgu133plus2.db  -cdfname1NULL
hgu133plus2  -cdfnameNULL
Error in mget(x, envir = get(paste(cdfname, "SYMBOL", sep = "")), ifnotfound = list(simpleaffy:::.if.probeset.not.found)) :
  only NA is currently supported for 'ifnotfound'
>pw.filtered

    3.811797e-04     1.307617e-02     2.166568e-04     8.291114e-04
        37547_at       37549_g_at         37892_at         37950_at
    3.818242e-03     7.050013e-03     7.331992e-05     4.211978e-03
        38037_at         38149_at         39248_at         39249_at
    5.177551e-05     3.071316e-02     1.483113e-03     1.732506e-02
        39402_at         41512_at         41660_at         45297_at
    5.762073e-03     1.370193e-02     5.776765e-03     7.799225e-04
        45714_at         50221_at         51226_at         53991_at
    4.207862e-04     2.384170e-03     8.849272e-03     6.814868e-05
        57715_at       64408_s_at         65472_at         65517_at
    3.240915e-03     2.841543e-02     1.042202e-02     1.073591e-03
AFFX-M27830_3_at
    2.282442e-02

Slot "calls":
<0 x 0 matrix>

Slot "group":
[1] "treatment"

Slot "members":
[1] "n" "a"

Slot "pData":
              sample treatment
GSM190178.CEL      1         n
GSM190179.CEL      2         n
GSM190180.CEL      3         n
GSM190185.CEL      4         a
GSM190186.CEL      5         a

Slot "calculated.from":
ExpressionSet (storageMode: lockedEnvironment)
assayData: 54675 features, 5 samples
  element names: exprs
protocolData: none
phenoData
  sampleNames: GSM190178.CEL GSM190179.CEL ... GSM190186.CEL (5 total)
  varLabels: sample treatment
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu133plus2

Best,

Matthew P. Grieshop

University of Wisconsin

Department of Biochemistry

This doesn't have anything to do with the PairComp object, or what it has in any of its slots. Instead, the issue arises from the mget() method as pertaining to a Bimap object.

The simpleaffy package is one of the oldest packages in Bioconductor, and was around back in the days that annotation packages were based on environments instead of the current SQLite databases. If we still used environments for the annotation packages, then this function would work just fine. But when you use mget() on a current annotation package, it dispatches a different function that no longer accepts anything but NA as the argument to ifnotfound.

You can get around this problem by specifying your own functions (untested!):

matts.get.annotation <- function(nam, cdfname) {
require(cdfname, character.only = TRUE)
cdfname <- get(cdfname)
out <- select(cdfname, nam, c("SYMBOL","GENENAME","ACCNUM","UNIGENE"))
out <- out[!duplicated(out[,1]),2:5]
out
}

matts.results.summary <- function(results, cdfname) {
res <- cbind(means(results), fc(results), sapply(fc(results),
        function(x) {
            if (x < 0) {
                -1 * 2^(-1 * x)
            }
            else {
                2^x
            }
        }), tt(results), matts.get.annotation(names(fc(results)), cdfname))
    cns <- colnames(res)
    cns[3] <- "log2(fc)"
    cns[4] <- "fc"
    cns[5] <- "tt"
    colnames(res) <- cns
    return(res)
}

And then summary <- matts.results.summary(pw.filtered, "hgu133plus2.db") should do what you expect.