[BioC] pa.calls-function in package "panp"

[Kyle R Covington]

I found a fix for this. All that is required to skip the error is to check if 
the x value in the cuttoff_fcn definition is NA.  I decided to return "A" in 
this case since the model can't tell me that the value is "P", but a more 
liberal person may want to use "P" as the default.

Complete code is below.

pa.calls2<-function (object = NULL, looseCutoff = 0.02, tightCutoff = 0.01, 
    verbose = FALSE) 
{
    if (is.null(object)) {
        cat("\nUSAGE:\n\tpa.calls(object, looseCutoff=0.02, tightCutoff=0.01, 
verbose = FALSE)\n\nINPUTS: \n\tobject - ExpressionSet, returned from 
expression-generating function, \n\t   such as expresso(), rma(), mas5(), 
gcrma(), etc. \n\tlooseCutoff - the larger P-value cutoff\n\ttightCutoff - the 
smaller, more strict P-value cutoff\n\tverbose - TRUE or FALSE\n\nOUTPUTS: 
\n\tReturns a list of two matrices, Pcalls and Pvals:\n\tPvals - a matrix of P-
values of same dimensions as exprs(input object). Each\n\t   datapoint is the P-
value for the probeset at the same x,y coordinates. \n\tPcalls - a matrix of 
Presence (P), Marginal (M), Absent (A) indicators\n\n")
        return()
    }
    if (require(affy) == FALSE) {
        stop("\npa.calls() requires BioConductor 'affy' package.\n  Currently, 
it is not installed. Please install and load, then\n  rerun pa.calls()\n\n")
    }
    if (verbose) {
        cat("\nInvoking function 'pa.calls'\n")
        cat("tightCutoff is ", tightCutoff, "\nlooseCutoff is ", 
            looseCutoff, "\n")
    }
    if (!is(object, "ExpressionSet")) {
        stop("\nAborting: object must be an ExpressionSet\n\n")
    }
    chip = annotation(object)
    if (chip == "hgu133b") {
        stop("\nHG-U133B is not currently supported. Supported chip types 
are\n\nHG-U133A and HG-U133 Plus 2.0\n\n")
    }
    if ((chip != "hgu133a") & (chip != "hgu133atag") & (chip != 
        "hgu133plus2")) {
        stop("\nAborting: chip type must be either HG-U133A or HG-U133 Plus 2.0 
\n\n")
    }
    if ((tightCutoff > 1) | (tightCutoff < 0) | (looseCutoff > 
        1) | (looseCutoff < 0)) {
        stop("\nAborting: cutoffs must be between 0.0 and 1.0\n\n")
    }
    if (tightCutoff > looseCutoff) {
        stop("\nAborting: tightCutoff must be lower than looseCutoff\n\n")
    }
    if (verbose) {
        cat("Outputs will be a list containing two matrices of same dimensions 
as input full dataset:\n  1. Pcalls - a matrix of P/A/M indicators: \n  2. Pvals 
- a matrix of P-values \n  'P': P-values <= tightCutoff ", 
            tightCutoff, "\n  'A': P-values > looseCutoff ", 
            looseCutoff, "\n  'M': P-values between ", tightCutoff, 
            " and ", looseCutoff, "\n")
    }
    if ((chip == "hgu133a") | (chip == "hgu133atag")) {
        data(NSMPnames.hgu133a)
        NSMPnames <- NSMPnames.hgu133a
        rm(NSMPnames.hgu133a, envir = globalenv())
    }
    else {
        data(NSMPnames.hgu133plus2)
        NSMPnames <- NSMPnames.hgu133plus2
        rm(NSMPnames.hgu133plus2, envir = globalenv())
    }
    AllExprs <- exprs(object)
    NegExprs <- AllExprs[NSMPnames, ]
    AllExprs <- as.matrix(AllExprs)
    NegExprs <- as.matrix(NegExprs)
    cutoff_fcn <- function(x) {
        if (is.na(x))                      ## fix by KRC
            return("A")
        if (x <= tightCutoff) 
            return("P")
        else if (x > looseCutoff) 
            return("A")
        else return("M")
    }
    len <- length(colnames(AllExprs))
    cat("\nProcessing", len, "chips: ")
    flush.console()
    for (chip in 1:len) {
        xNeg <- NegExprs[, chip]
        xNeg <- sort(xNeg, decreasing = TRUE)
        yNeg <- seq(0, 1, 1/(length(xNeg) - 1))
        interp <- approxfun(xNeg, yNeg, yleft = 1, yright = 0)
        PV <- sapply(AllExprs[, chip], interp)
        PC <- sapply(PV, cutoff_fcn)
        if (chip == 1) {
            Pvals <- PV
            Pcalls <- PC
        }
        else {
            Pvals <- cbind(Pvals, PV)
            Pcalls <- cbind(Pcalls, PC)
        }
        cat("#")
        flush.console()
    }
    Pvals <- as.matrix(Pvals)
    Pcalls <- as.matrix(Pcalls)
    colnames(Pvals) <- colnames(AllExprs)
    colnames(Pcalls) <- colnames(AllExprs)
    outlist <- list(Pcalls = Pcalls, Pvals = Pvals)
    cat("\nProcessing complete.\n\n")
    flush.console()
    myX <- NegExprs
    maxLen <- 0
    for (i in 1:len) {
        stringLen <- nchar(colnames(AllExprs)[i])
        if (stringLen > maxLen) 
            maxLen <- stringLen
    }
    maxLen <- maxLen - 6
    if (maxLen < 2) {
        maxLen = 2
    }
    for (i in 1:len) {
        myY <- seq(0, 1, 1/(length(myX[, i]) - 1))
        myX[, i] <- sort(myX[, i], decreasing = TRUE)
        revInterp <- approxfun(myY, myX[, i], yleft = 1, yright = 0)
        revTight <- revInterp(tightCutoff)
        revLoose <- revInterp(looseCutoff)
        if (i == 1) {
            cat("\nIntensities at cutoff P-values of ", looseCutoff, 
                " and ", tightCutoff, ":\n")
            cat("Array:")
            for (j in 1:maxLen) {
                cat(" ")
            }
            cat("value at", looseCutoff, "   value at", tightCutoff, 
                "\n")
        }
        cat(colnames(AllExprs)[i], "\t", format.pval(revLoose, 
            digits = 3), "\t\t", format.pval(revTight, digits = 3), 
            "\n")
    }
    cat("\n")
    cat("[NOTE: 'Collapsing to unique x values...' warning messages are 
benign.]\n\n")
    return(outlist)
}