Entering edit mode
Kyle R Covington
▴
10
@kyle-r-covington-4607
Last seen 10.2 years ago
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)
}
