[BioC] pa.calls-function in package "panp"
Samuel Wuest
wuests at tcd.ie
Mon Apr 25 13:35:07 CEST 2011
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
More information about the Bioconductor
mailing list