pa.calls-function in package "panp"
1
0
Entering edit mode
@kyle-r-covington-4607
Last seen 9.6 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) }
• 867 views
ADD COMMENT
0
Entering edit mode
Samuel Wuest ▴ 330
@samuel-wuest-2821
Last seen 9.6 years ago
Hi Kyle, as far as my original post is concerned: there was a small bug in the approxfun-function used by pa.calls in R 2.11.X, and the problem has been fixed after R 2.11.1. I am using R 2.12 now, and there are no further problems. See also: https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=14377 If you are still using R 2.11, then an easy fix is to round the expression values in your ExprSet, as posted earlier (for example exprs(eset) <- round(exprs(eset), digits=6) or anything like this) Hope this helps, best, Sam On 20 April 2011 15:53, Kyle R Covington <kyle at="" red-r.org=""> wrote: > 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) > } > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor > > -- ----------------------------------------------------- Samuel Wuest Smurfit Institute of Genetics Trinity College Dublin Dublin 2, Ireland Phone: +353-1-896 2444 Web: http://www.tcd.ie/Genetics/wellmer-2/index.html Email: wuests at tcd.ie
ADD COMMENT

Login before adding your answer.

Traffic: 796 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6