[BioC] Problem obtaining QC statistics from 'simpleaffy'

[Bio Sam]

Hi Greg, I had a similar problem when trying to run the affy qc stuff
on barley.  After some digging I found that the necessary entries were
not included for the barley chip in the packages data files
"qc.probes.tab", "alpha.tab", and "spikes.tab" which are located in
"/usr/local/lib/R/site-library/simpleaffy/data/" on my Linux system.
Adding entries for the barley chip did NOT fix the problem however, as
apparently this data isn't read from the data files once the plugin is
compiled.  To fix this problem, I re-wrote two functions which
basically override the original QCReport(), and qc.affy() functions.
I'm sure you could modify them for your own needs as well.  I stored
both functions in a file called "qc_barley.R" (shown below).  To use
them I run:
>source("qc_barley.R")
>QC_barley_Report(abatch)

It is sort of a quick-hack, and from Crispin Miller's post (mentioned
above) it looks like better things are on their way, but until then, I
hope this helps.

Sam Hunter
BCB - University of Idaho

Contents of my qc_barley.R:
#NOTE: qc.probenames, spike.probenames and bb, may not contain the
#correct probe names... but this will run
require("affyQCReport")

qc.barley <-
#> qc.affy
function (unnormalised, normalised = NULL, tau = 0.015, logged = TRUE,
    cdfn = cleancdfname(cdfName(unnormalised)))
{

    print(encodeString("1:normalizing"))

    if (is.null(normalised)) {
        getAllSpikeProbes("hgu133acdf")
        normalised <- call.exprs(unnormalised, "mas5")
    }

    print(encodeString("2:normalizing done"))

    x <- exprs(normalised)
    det <- detection.p.val(unnormalised, tau = tau, alpha1 = 0.05,
        alpha2 = 0.065)

    print(encodeString("3:detection.p.val"))

    dpv <- apply(det$call, 2, function(x) {
        x[x != "P"] <- 0
        x[x == "P"] <- 1
        x <- as.numeric(x)
        return(100 * sum(x)/length(x))
    })

    print(encodeString("4:apply(det$call)"))

    sfs <- normalised at description@preprocessing$sfs
    target <- normalised at description@preprocessing$tgt
    if (!logged) {
        x <- log2(x)
    }

    print(encodeString("5:stats"))

    bgsts <- .bg.stats(unnormalised)$zonebg
    meanbg <- apply(bgsts, 1, mean)
    minbg <- apply(bgsts, 1, min)
    maxbg <- apply(bgsts, 1, max)
    stdvbg <- sqrt(apply(bgsts, 1, var))
    qc.probenames <- c("AFFX-BioB-3_at", "AFFX-BioB-M_at",
"AFFX-BioB-5_at", "AFFX-TrpnX-3_at", "AFFX-TrpnX-M_at",
"AFFX-TrpnX-5_at")

    print(encodeString("6:end stats"))

    qc.probe.vals <- rbind(c(), (sapply(qc.probenames, function(y) {
        x[y, ]
    })))
    rownames(qc.probe.vals) <- colnames(x)
    spike.probenames <- c("AFFX-r2-Ec-bioB-3_at",
"AFFX-r2-Ec-bioC-3_at", "AFFX-r2-Ec-bioD-3_at", "AFFX-r2-P1-cre-3_at")
    spike.vals <- rbind(c(), (sapply(spike.probenames, function(y) {
        x[y, ]
    })))
    rownames(spike.vals) <- colnames(x)
    #bb <- getBioB(cdfn)
    bb <- c("AFFX-r2-Ec-bioB-3_at")
    if (!is.na(bb)) {
        biobcalls <- det$call[bb, ]
    }
    else {
        biobcalls <- NULL
    }
    return(new("QCStats", scale.factors = sfs, target = target,
        percent.present = dpv, average.background = meanbg,
minimum.background = minbg,
        maximum.background = maxbg, spikes = spike.vals, qc.probes =
qc.probe.vals,
        bioBCalls = biobcalls))
}


QC_barley_Report <-
function (object, file = "AffyQCReport.pdf", ...)
{
    pdf(file = file, width = 8, height = 11, onefile = TRUE)
    plot.window(c(1, 1), c(0, 1))
    plot.new()
    titlePage(object)
    signalDist(object)
    plot(qc.barley(object))
    borderQC1(object)
    borderQC2(object)
    correlationPlot(object)
    dev.off()
    return(TRUE)
}