[BioC] Question

[Laurent Gautier]

On Tue, Jan 06, 2004 at 04:33:35PM -0500, Riki Kawaguchi wrote:
> Dear Bioconductor Administrators/Developers,
> 
> Hello. I am a beginner of Bioconductor affy package. I found a great
> potential of this software and I really appreciate people who made this. I
> was wondering if affy package has a function to do similar analysis to
> MAS5.0 pairwise comparison to generate signal log ratios of PMs (not PM-MM).
> Thank you so much for your help!
> 
> Best wishes,
> 
> Riki Kawaguchi
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor

As far as I know, there is no such function. However, an 'apply-like'
function was developped for this kind of work with
probe sets (objects of class ProbeSet)...
...and when time it will make its way to Bioconductor CVS repository
(it should be for the next release).

The code for it (with the doc file) will hopefully help you to achieve
what you want without too much effort. Check the example that
performs t-test for each probe in probe sets.


Hoping it helps,


L.


-- 
--------------------------------------------------------------
Laurent Gautier			CBS, Building 208, DTU
PhD. Student			DK-2800 Lyngby,Denmark	
tel: +45 45 25 24 89		http://www.cbs.dtu.dk/laurent
-------------- next part --------------
ppsetApply <- function(abatch, FUN, genenames=NULL, ...) {

  if (! inherits(abatch, "AffyBatch"))
    stop("abatch must be inheriting from class AffyBatch")

  if (! inherits(FUN, "function"))
    stop("FUN must be a function")
  
  cdfenv <- getCdfInfo(abatch)

  if (is.null(genenames))
    genenames <- ls(cdfenv)

  ##
  e1 <- new.env(parent = environment(FUN))
  multiassign(names(pData(abatch)), pData(abatch), env = e1)
  environment(FUN) <- e1
  
  ppset <- new("ProbeSet", pm=matrix(), mm=matrix())

  r <- vector("list", length=length(genenames))
  names(r) <- genenames
  
  for (i in seq(along=genenames)) {
    ## use multiget to get NA when genenames[i] not found
    probes.i <- multiget(genenames[i], envir = cdfenv)[[1]]
    if (all(is.na(probes.i)))
      next
    ppset at pm <- intensity(abatch)[probes.i[, 1], , drop=FALSE]
    ppset at mm <- intensity(abatch)[probes.i[, 2], , drop=FALSE]
    ppset at id <- genenames[i]
    r[[i]] <- FUN(ppset, ...)
  }

  return(r)
}

ppsetClusterApply <- function(abatch, FUN, genenames=NULL, ...) {

  if (! inherits(abatch, "AffyBatch"))
    stop("abatch must be inheriting from class AffyBatch")

  if (! inherits(FUN, "function"))
    stop("FUN must be a function")
  
  cdfenv <- getCdfInfo(abatch)

  if (is.null(genenames))
    genenames <- ls(cdfenv)

  ##
  e1 <- new.env(parent = environment(FUN))
  multiassign(names(pData(abatch)), pData(abatch), env = e1)
  environment(FUN) <- e1
  
  ppset <- new("ProbeSet", pm=matrix(), mm=matrix())

  r <- vector("list", length=length(genenames))
  names(r) <- genenames
  
  for (i in seq(along=genenames)) {
    ## use multiget to get NA when genenames[i] not found
    probes.i <- multiget(genenames[i], envir = cdfenv)[[1]]
    if (all(is.na(probes.i)))
      next
    ppset at pm <- intensity(abatch)[probes.i[, 1], , drop=FALSE]
    ppset at mm <- intensity(abatch)[probes.i[, 2], , drop=FALSE]
    ppset at id <- genenames[i]
    r[[i]] <- FUN(ppset, ...)
  }

  return(r)
}


ppset.ttest <- function(ppset, covariate, pmcorrect.fun = pmcorrect.pmonly, ...) {
  probes <- do.call("pmcorrect.fun", list(ppset))
  my.ttest <- function(x) {
    y <- split(x, get(covariate))
    t.test(y[[1]], y[[2]])$p.value
  }
  r <- apply(probes, 1, my.ttest)
  return(r)
}



# make.ppset.logitt <- function(abatch) {
  
#   ppset.logitt <- function(ppset, covariate, pmcorrect.fun = pmcorrect.pmonly, A, N) {
#     probes <- do.call("pmcorrect.fun", list(ppset))
#     probes.logit <- 
#     }

#   return(ppset.logitt)
# }

-------------- next part --------------
\name{ppsetApply}
\alias{ppsetApply}
\title{ Apply a function over the ProbeSets in an AffyBatch }
\description{
  Apply a function over the ProbeSets in an AffyBatch
}
\usage{
ppsetApply(abatch, FUN, genenames = NULL, ...)

ppset.ttest(ppset, covariate, pmcorrect.fun = pmcorrect.pmonly, ...)
}
\arguments{
  \item{abatch}{ An object inheriting from \code{AffyBatch}.}
  \item{FUN}{ A function working on a \code{ProbeSet} }
  \item{genenames}{ A list of Affymetrix probesets ids to work with. All
  probe set ids used when \code{NULL}.}
  \item{\dots}{ Optional parameters to the function \code{FUN} }
}
\details{
  
}
\value{
  Returns a \code{list} of objects, or values, as returned by the
  function \code{FUN}
  for each \code{ProbeSet} it processes.
}
\author{Laurent Gautier <laurent at cbs.dtu.dk>}
\seealso{\code{\link[affy]{ProbeSet-class}} }
\examples{
ppset.ttest <- function(ppset, covariate, pmcorrect.fun = pmcorrect.pmonly, ...) {
  probes <- do.call("pmcorrect.fun", list(ppset))
  my.ttest <- function(x) {
    y <- split(x, get(covariate))
    t.test(y[[1]], y[[2]])$p.value
  }
  r <- apply(probes, 1, my.ttest)
  return(r)
}

## craft a dataset
data(affybatch.example)
abatch <- merge(affybatch.example, affybatch.example)
intensity(abatch) <- jitter(intensity(abatch))
chip.variate <- c("a", "b", "a", "a", "b", "a", "a")
pData(abatch) <- data.frame(whatever = chip.variate)

## run a test over _all_ probes.
all.ttest <- ppsetApply(abatch, ppset.ttest, covariate="whatever")

}
\keyword{ manip }