R/ChipEffectFile.R
In HenrikBengtsson/aroma.affymetrix: Analysis of Large Affymetrix Microarray Data Sets
###########################################################################/**
# @RdocClass ChipEffectFile
#
# @title "The ChipEffectFile class"
#
# \description{
# @classhierarchy
#
# This class represents estimates of chip effects in the probe-level models.
# }
#
# @synopsis
#
# \arguments{
# \item{...}{Arguments passed to @see "ParameterCelFile".}
# \item{probeModel}{The specific type of model, e.g. \code{"pm"}.}
# }
#
# \section{Fields and Methods}{
# @allmethods "public"
# }
#
# @author "HB"
#
# \seealso{
# An object of this class is typically obtained through the
# \code{getChipEffectSet()} method for the @see "ProbeLevelModel" class.
# An object of this class is typically part of a @see "ChipEffectSet".
# }
#*/###########################################################################
setConstructorS3("ChipEffectFile", function(..., probeModel=c("pm")) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'probeModel':
probeModel <- match.arg(probeModel)
this <- extend(ParameterCelFile(...), "ChipEffectFile",
"cached:.firstCells" = NULL,
probeModel = probeModel
)
setEncodeFunction(this, function(groupData, ...) {
theta <- .subset2(groupData, "theta")
stdvs <- .subset2(groupData, "sdTheta")
outliers <- .subset2(groupData, "thetaOutliers")
pixels <- NULL
if (!is.null(outliers))
pixels <- -as.integer(outliers)
res <- list()
if (!is.null(theta))
res$intensities <- theta
if (!is.null(stdvs))
res$stdvs <- stdvs
if (!is.null(pixels))
res$pixels <- pixels
res
})
## Not safe, at least not what I know of. /HB 2007-08-17
## setEncodeFunction(this, function(groupData, ...) {
## groupData[[3]] <- -as.integer(.subset2(groupData, 3))
## names(groupData) <- c("intensities", "stdvs", "pixels")
## groupData
## })
setDecodeFunction(this, function(groupData, ...) {
res <- list()
if (!is.null(groupData$intensities))
res$theta <- groupData$intensities
if (!is.null(groupData$stdvs))
res$sdTheta <- groupData$stdvs
if (!is.null(groupData$pixels))
res$thetaOutliers <- as.logical(-groupData$pixels)
res
})
# Parse attributes (all subclasses must call this in the constructor).
setAttributesByTags(this)
this
})
setMethodS3("as.character", "ChipEffectFile", function(x, ...) {
# To please R CMD check
this <- x
s <- NextMethod("as.character")
s <- c(s, sprintf("Parameters: %s", getParametersAsString(this)))
s
}, protected=TRUE)
setMethodS3("getParameters", "ChipEffectFile", function(this, ...) {
params <- NextMethod("getParameters")
params$probeModel <- this$probeModel
params
}, protected=TRUE)
setMethodS3("createParamCdf", "ChipEffectFile", function(static, sourceCdf, ..., verbose=FALSE) {
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
verbose && enter(verbose, "Creating CDF for chip effects")
verbose && cat(verbose, "Source chip type: ", getChipType(sourceCdf))
verbose && cat(verbose, "Source CDF: ", getPathname(sourceCdf))
# Search for existing monocell CDF
for (sep in c(",", "-")) {
chipType <- paste(getChipType(sourceCdf), "monocell", sep=sep)
verbose && cat(verbose, "Looking for chip type: ", chipType)
pathname <- AffymetrixCdfFile$findByChipType(chipType)
if (!is.null(pathname)) {
verbose && cat(verbose, "Found: ", pathname)
break
}
}
if (is.null(pathname)) {
verbose && cat(verbose, "Pathname: Not found!")
verbose && cat(verbose, "Will create CDF for the chip-effect files from the original CDF. NOTE: This will take several minutes or more!")
verbose && enter(verbose, "Creating CDF")
cdf <- createMonocellCdf(sourceCdf, verbose=less(verbose))
verbose && exit(verbose)
} else {
verbose && cat(verbose, "Pathname: ", pathname)
cdf <- AffymetrixCdfFile$fromFile(pathname)
}
verbose && exit(verbose)
cdf
}, static=TRUE, private=TRUE)
setMethodS3("readUnits", "ChipEffectFile", function(this, units=NULL, cdf=NULL, ..., force=FALSE, cache=FALSE, verbose=FALSE) {
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
# Check for cached data
key <- list(method="readUnits", class=class(this)[1],
pathname=getPathname(this),
cdf=cdf, units=units, ...)
if (getOption(aromaSettings, "devel/useCacheKeyInterface", FALSE)) {
key <- getCacheKey(this, method="readUnits", pathname=getPathname(this), cdf=cdf, units=units, ...)
}
id <- getChecksum(key)
res <- this$.readUnitsCache[[id]]
if (!force && !is.null(res)) {
verbose && cat(verbose, "readUnits.ChipEffectFile(): Returning cached data")
return(res)
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Retrieve the data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if (is.null(cdf)) {
cdf <- getCellIndices(this, units=units, verbose=less(verbose))
}
# Note that the actually call to the decoding is done in readUnits()
# of the superclass.
res <- NextMethod("readUnits", cdf=cdf, force=force, verbose=less(verbose))
# Store read units in cache?
if (cache) {
verbose && cat(verbose, "readUnits.ChipEffectFile(): Updating cache")
this$.readUnitsCache <- list()
this$.readUnitsCache[[id]] <- res
}
res
})
###########################################################################/**
# @RdocMethod getCellIndices
#
# @title "Retrieves tree list of cell indices for a set of units"
#
# \description{
# @get "title" from the associated CDF.
# }
#
# @synopsis
#
# \arguments{
# \item{...}{Additional arguments passed to \code{getCellIndices()}
# of @see "AffymetrixCdfFile".}
# \item{.cache}{Ignored.}
# }
#
# \value{
# Returns a @list structure, where each element corresponds to a unit.
# If argument \code{unlist=TRUE} is passed, an @integer @vector is returned.
# }
#
# \seealso{
# @seeclass
# }
#
# @keyword internal
#*/###########################################################################
setMethodS3("getCellIndices", "ChipEffectFile", function(this, ..., .cache=TRUE) {
cdf <- getCdf(this)
getCellIndices(cdf, ...)
}, protected=TRUE)
setMethodS3("updateUnits", "ChipEffectFile", function(this, units=NULL, cdf=NULL, data, ...) {
if (is.null(cdf))
cdf <- getCellIndices(this, units=units)
# Note that the actually call to the encoding is done in updateUnits()
# of the superclass.
NextMethod("updateUnits", cdf=cdf, data=data)
}, private=TRUE)
setMethodS3("findUnitsTodo", "ChipEffectFile", function(this, units=NULL, ..., force=FALSE, verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Identifying non-fitted units in chip-effect file")
verbose && cat(verbose, "Pathname: ", getPathname(this))
idxs <- NULL
if (is.null(units)) {
# Look up chip-type and parameter specific but data set independent data
cdf <- getCdf(this)
chipType <- getChipType(cdf)
key <- list(method="findUnitsTodo", class=class(this)[1],
chipType=chipType, params=getParameters(this))
dirs <- c("aroma.affymetrix", chipType)
if (!force) {
idxs <- loadCache(key, dirs=dirs)
if (!is.null(idxs))
verbose && cat(verbose, "Found indices cached on file")
}
}
if (is.null(idxs)) {
verbose && enter(verbose, "Identifying CDF units")
units0 <- units
if (is.null(units)) {
cdf <- getCdf(this)
units <- seq_len(nbrOfUnits(cdf))
}
nbrOfUnits <- length(units)
idxs <- lapplyInChunks(units, function(unitsChunk, ...) {
verbose && enter(verbose, "Reading CDF cell indices")
idxsChunk <- getCellIndices(this, units=unitsChunk, force=TRUE, verbose=less(verbose))
# Save memory 100% -> 94%
names(idxsChunk) <- NULL
verbose && exit(verbose)
verbose && enter(verbose, "Extracting first CDF group for each unit")
# Save memory 100% -> 94% -> 5% (all the nested named lists costs!)
idxsChunk <- lapply(idxsChunk, FUN=function(unit) {
groups <- .subset2(unit, "groups")
fields <- .subset2(groups, 1)
.subset2(fields, 1)
})
verbose && exit(verbose)
gc <- gc()
idxsChunk
}, chunkSize=100e3, useNames=FALSE, verbose=verbose)
# Not needed anymore
units <- units0
# Not needed anymore
units0 <- NULL
idxs <- unlist(idxs, use.names=FALSE)
gc <- gc()
verbose && print(verbose, gc)
# Verify correctness
if (length(idxs) != nbrOfUnits) {
throw("Internal error: Expected ", nbrOfUnits, " cell indices, but got ", length(idxs), ".")
}
if (is.null(units)) {
verbose && enter(verbose, "Saving to file cache")
saveCache(idxs, key=key, dirs=dirs)
verbose && exit(verbose)
}
verbose && exit(verbose)
}
# Read one cell from each unit
verbose && enter(verbose, "Reading data for these ", length(idxs), " cells")
value <- .readCel(getPathname(this), indices=idxs, readIntensities=FALSE,
readStdvs=TRUE, readPixels=FALSE)$stdvs
verbose && exit(verbose)
# Identify units for which the stdvs <= 0.
value <- which(value <= 0)
if (!is.null(units))
value <- units[value]
verbose && cat(verbose, "Looking for stdvs <= 0 indicating non-estimated units:")
verbose && str(verbose, value)
verbose && exit(verbose)
value
})
###########################################################################/**
# @RdocMethod getUnitGroupCellMap
# @aliasmethod getCellMap
#
# @title "Gets a (unit, group, cell) index map"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{units}{The units for which the map should be returned.
# If @NULL, all units are considered.}
# \item{force}{If @TRUE, cached cell indices are ignored.}
# \item{...}{Not used.}
# \item{verbose}{See @see "R.utils::Verbose".}
# }
#
# \value{
# Returns a @data.frame with @integer columns \code{unit}, \code{group},
# and \code{cell}.
# }
#
# \examples{\dontrun{
# unit group cell
# 1 104 1 335
# 2 104 2 336
# 3 105 1 337
# 4 105 2 338
# 5 105 3 339
# 6 105 4 340
# }}
#
# \seealso{
# @seeclass
# }
#*/###########################################################################
setMethodS3("getUnitGroupCellMap", "ChipEffectFile", function(this, units=NULL, force=FALSE, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'units':
if (inherits(units, "UnitGroupCellMap")) {
return(units)
} else if (is.null(units)) {
} else if (is.list(units)) {
} else {
units <- Arguments$getIndices(units)
}
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Retrieving unit-to-cell map")
# Special case: requesting zero units?
if (length(units) == 0 && !is.null(units)) {
map <- data.frame(unit=integer(0), group=integer(0), cell=integer(0))
class(map) <- c("UnitGroupCellMap", class(map))
verbose && exit(verbose)
return(map)
}
# Get the CDF
cdf <- getCdf(this)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Check cache
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# useFileCache <- (is.null(units) || (!is.list(units) && length(units) > 10000))
useFileCache <- (is.null(units) || (!is.list(units) && length(units) > 100))
if (useFileCache) {
chipType <- getChipType(cdf)
# Look up chip-type and parameter specific but data set independent data
key <- list(method="getUnitGroupCellMap", class=class(this)[1],
chipType=chipType, params=getParameters(this),
units=units)
dirs <- c("aroma.affymetrix", chipType)
if (!force) {
map <- loadCache(key, dirs=dirs)
if (!is.null(map)) {
verbose && cat(verbose, "Found (unit,group,cell) map cached on file")
verbose && exit(verbose)
return(map)
}
}
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Main
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Is 'units' already a CDF list?
if (is.list(units)) {
# No fancy validation for now.
cells <- units
units <- indexOf(cdf, names=names(units))
if (any(is.na(units))) {
throw("Argument 'units' is of unknown structure.")
}
verbose && enter(verbose, "Argument 'cells' is already a CDF cell-index structure")
# Get the unit names
unitNames <- names(cells)
# Get the number of groups per unit
unitSizes <- lapply(cells, FUN=function(unit) {
length(.subset2(unit, "groups"))
})
unitSizes <- unlist(unitSizes, use.names=FALSE)
cells <- unlist(cells, use.names=FALSE)
} else {
verbose && enter(verbose, "Retrieving cell indices for specified units")
if (is.null(units))
units <- seq_len(nbrOfUnits(cdf))
chunks <- splitInChunks(units, chunkSize=100e3)
nbrOfChunks <- length(chunks)
nbrOfUnits <- length(units)
# Store unit names and unit sizes
unitNames <- vector("character", nbrOfUnits)
unitSizes <- vector("integer", nbrOfUnits)
# Store cell indices first chunk-by-chunk, then as a vector.
cells <- vector("list", nbrOfChunks)
offset <- 0
for (kk in seq_len(nbrOfChunks)) {
verbose && printf(verbose, "Chunk #%d of %d\n", kk, length(chunks))
chunk <- chunks[[kk]]
chunks[[kk]] <- NA
cells0 <- getCellIndices(this, units=chunk, force=force, .cache=FALSE, verbose=less(verbose))
idxs <- offset + seq_along(chunk)
offset <- offset + length(chunk)
# Not needed anymore
chunk <- NULL
unitNames[idxs] <- names(cells0)
names(cells0) <- NULL
unitSizes0 <- lapply(cells0, FUN=function(unit) {
length(.subset2(unit, "groups"))
})
unitSizes[idxs] <- unlist(unitSizes0, use.names=FALSE)
# Not needed anymore
unitSizes0 <- NULL
cells[[kk]] <- unlist(cells0, use.names=FALSE)
# Not needed anymore
cells0 <- idxs <- NULL
}
# Not needed anymore
chunks <- NULL
gc <- gc()
verbose && print(verbose, gc)
}
cells <- unlist(cells, use.names=FALSE)
gc <- gc()
verbose && exit(verbose)
verbose && enter(verbose, "Creating return data frame")
uUnitSizes <- sort(unique(unitSizes))
verbose && cat(verbose, "Unique number of groups per unit: ",
paste(uUnitSizes, collapse=","))
verbose && cat(verbose, "Number of units: ", length(unitNames))
if (is.null(units))
units <- seq_len(nbrOfUnits(cdf))
# The following is too slow:
# groups <- sapply(unitSizes, FUN=function(n) seq_len(n))
# Instead, updated size by size
verbose && printf(verbose, "Allocating matrix of size %dx%d.\n",
max(uUnitSizes), length(unitNames))
naValue <- NA_integer_
units2 <- groups <- matrix(naValue, nrow=max(uUnitSizes), ncol=length(unitNames))
for (size in uUnitSizes) {
# Identify units with a certain number of groups
cc <- which(unitSizes == size)
# Assign group indices to the groups
seq <- seq_len(size)
groups[seq,cc] <- seq
# Assign unit indices to ditto
units2[seq,cc] <- rep(units[cc], each=size)
# Not needed anymore
seq <- NULL
gc <- gc()
}
keep <- !is.na(groups)
groups <- groups[keep]
units2 <- units2[keep]
# Not needed anymore
keep <- NULL
gc <- gc()
map <- data.frame(unit=units2, group=groups, cell=cells)
verbose && exit(verbose)
verbose && exit(verbose)
class(map) <- c("UnitGroupCellMap", class(map))
if (useFileCache) {
verbose && enter(verbose, "Saving to file cache")
saveCache(map, key=key, dirs=dirs)
verbose && exit(verbose)
}
map
}, private=TRUE)
setMethodS3("getUnitGroupCellChromosomePositionMap", "ChipEffectFile", function(this, units=NULL, chromosomes=NULL, orderByPosition=TRUE, ..., force=FALSE, verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
cdf <- getCdf(this)
# Argument 'units':
ugcMap <- NULL
if (is.null(units)) {
} else if (isUnitGroupCellMap(units)) {
ugcMap <- units
units <- ugcMap[,"unit"]
}
if (!is.null(units)) {
units <- Arguments$getIndices(units, range=c(1, nbrOfUnits(cdf)))
}
units0 <- units
# Get the genome position information
gi <- getGenomeInformation(cdf)
# Argument 'chromosomes':
if (!is.null(chromosomes)) {
allChromosomes <- getChromosomes(gi)
unknown <- chromosomes[!(chromosomes %in% allChromosomes)]
if (length(unknown) > 0) {
throw("Argument 'chromosomes' contains unknown values: ",
paste(unknown, collapse=", "))
}
}
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Getting (unit, group, cell, chromosome, position) map")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Check for cached results
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Look for results in file cache
verbose && enter(verbose, "Checking cache")
chipType <- getChipType(cdf)
key <- list(method="getUnitGroupCellChromosomePositionMap",
class=class(this)[1],
chipType=chipType, units=units, ugcMap=ugcMap,
chromosomes=chromosomes, orderByPosition=orderByPosition)
dirs <- c("aroma.affymetrix", chipType)
if (!force) {
map <- loadCache(key=key, dirs=dirs)
if (!is.null(map)) {
verbose && cat(verbose, "Found cached results")
verbose && exit(verbose)
return(map)
}
}
# Select by chromosome(s)?
if (!is.null(chromosomes)) {
verbose && cat(verbose, "Units:")
verbose && str(verbose, units)
verbose && cat(verbose, "Subset by chromosomes:")
verbose && str(verbose, chromosomes)
units <- getUnitsOnChromosomes(gi, chromosomes)
verbose && cat(verbose, "Units:")
verbose && str(verbose, units)
if (!is.null(units0)) {
units <- intersect(units, units0)
}
}
verbose && cat(verbose, "Units:")
verbose && str(verbose, units)
# Get the (unit, group, cell) map?
if (!isUnitGroupCellMap(ugcMap)) {
ugcMap <- getUnitGroupCellMap(this, units=units, force=force, verbose=less(verbose, 10))
verbose && cat(verbose, "(unit, group, cell) map:")
verbose && str(verbose, ugcMap)
}
# Get the (chromosome, position) map
cpMap <- getData(gi, units=ugcMap[,"unit"], force=force, verbose=less(verbose, 10))
verbose && cat(verbose, "(chromosome, position) map:")
verbose && str(verbose, cpMap)
# Sanity check
stopifnot(nrow(ugcMap) == nrow(cpMap))
# Merge the two maps
map <- cbind(ugcMap, cpMap)
# Not needed anymore
ugcMap <- cpMap <- NULL
if (orderByPosition) {
o <- with(map, order(chromosome, physicalPosition))
map <- map[o,,drop=FALSE]
# Not needed anymore
o <- NULL
verbose && cat(verbose, "Reordered by genomic position")
}
rownames(map) <- NULL
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Save to cache
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Save only results > 50kB
if (object.size(map) > 50e3) {
saveCache(map, key=key, dirs=dirs)
verbose && cat(verbose, "Saved to file cache")
}
verbose && exit(verbose)
map
}, private=TRUE)
setMethodS3("getDataFlat", "ChipEffectFile", function(this, units=NULL, fields=c("theta", "sdTheta", "outliers"), ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Retrieving data as a flat data frame")
# Get unit-to-cell map
suppressWarnings({
map <- getUnitGroupCellMap(this, units=units, ..., verbose=less(verbose))
})
verbose && enter(verbose, "Reading data fields")
celFields <- c(theta="intensities", sdTheta="stdvs", outliers="pixels")
suppressWarnings({
data <- getData(this, indices=map[,"cell"], fields=celFields[fields])
})
rownames(data) <- seq_len(nrow(data)); # Work around?!? /HB 2006-11-28
# Decode
names <- colnames(data)
names <- gsub("intensities", "theta", names)
names <- gsub("stdvs", "sdTheta", names)
names <- gsub("pixels", "outliers", names)
colnames(data) <- names
verbose && str(verbose, data)
if ("outliers" %in% names) {
data[,"outliers"] <- as.logical(-data[,"outliers"])
}
verbose && exit(verbose)
len <- sapply(data, FUN=length)
keep <- (len == nrow(map))
data <- data[keep]
data <- as.data.frame(data)
data <- cbind(map, data)
verbose && exit(verbose)
data
}, private=TRUE)
setMethodS3("updateDataFlat", "ChipEffectFile", function(this, data, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'data':
names <- colnames(data)
namesStr <- paste(names, collapse=", ")
if (!"cell" %in% names)
throw("Argument 'data' must contain a column 'cell': ", namesStr)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose2 <- -as.integer(verbose)-2
verbose && enter(verbose, "Storing flat data to file")
# Encode?
if ("outliers" %in% names) {
data[,"outliers"] <- -as.integer(data[,"outliers"])
}
names <- gsub("theta", "intensities", names)
names <- gsub("sdTheta", "stdvs", names)
names <- gsub("outliers", "pixels", names)
colnames(data) <- names
verbose && enter(verbose, "Updating file")
indices <- data[,"cell"]
keep <- (names %in% c("intensities", "stdvs", "pixels"))
data <- data[,keep]
pathname <- getPathname(this)
pathname <- Arguments$getWritablePathname(pathname)
.updateCel(pathname, indices=indices, data, verbose=verbose2)
verbose && exit(verbose)
verbose && exit(verbose)
invisible(data)
}, private=TRUE)
setMethodS3("mergeGroups", "ChipEffectFile", function(this, fcn, fields=c("theta", "sdTheta"), ..., pathname, overwrite=FALSE, verbose=FALSE) {
# Argument 'fcn':
if (!is.function(fcn)) {
throw("Argument 'fcn' is not a function: ", class(fcn)[1])
}
# Argument 'pathname':
pathname <- Arguments$getWritablePathname(pathname, mustNotExist=!overwrite)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Merging groups")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Test 'fcn':
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && enter(verbose, "Testing merge function")
for (size in 1:10) {
y <- matrix(1000+1:(size*4), nrow=size)
yOut <- fcn(y)
if (!identical(dim(yOut), dim(y))) {
throw("Function 'fcn' must not change the dimension of the data: ",
paste(dim(yOut), collapse="x"), " != ",
paste(dim(y), collapse="x"))
}
}
verbose && exit(verbose)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Get flat (unit, group, cell) map
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
map <- getUnitGroupCellMap(this, verbose=less(verbose))
verbose && str(verbose, map)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Merge data for each unit size separately (in reverse order!)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Get all possible unit sizes (number of groups per unit)
uSizes <- sort(unique(data[,"group"]))
verbose && cat(verbose, "Different number of groups per unit identified:")
verbose && print(verbose, uSizes)
data <- getDataFlat(this, units=map, ..., verbose=less(verbose))
verbose && str(verbose, data)
for (size in rev(uSizes)) {
verbose && enter(verbose, "Unit size ", size)
# Identify the units of that size
idxs <- which(data[,"group"] == size)
unitsS <- data[idxs, "unit"]
# Get the subset of the data for such units
idxs <- which(data[,"unit"] %in% unitsS)
for (field in fields) {
# Extract signals as a matrix where each column is one unit
y <- data[idxs, field]
y <- matrix(y, nrow=size)
verbose && str(verbose, y)
y <- fcn(y)
verbose && str(verbose, y)
# Update data table
data[idxs, field] <- as.vector(y)
}
# Get the cells where to store the merged data
verbose && exit(verbose)
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Copy CEL file and update the copy
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && enter(verbose, "Storing merged data")
verbose && cat(verbose, "Pathname: ", pathname)
# Create CEL file to store results, if missing
verbose && enter(verbose, "Creating CEL file for results, if missing")
cfN <- createFrom(this, filename=pathname, path=NULL, methods="create", clear=TRUE, verbose=less(verbose))
verbose && print(verbose, cfN)
verbose && exit(verbose)
verbose && enter(verbose, "Writing merged data")
updateDataFlat(cfN, data=data, verbose=less(verbose))
verbose && exit(verbose)
verbose && exit(verbose)
verbose && exit(verbose)
cfN
}, protected=TRUE)
setMethodS3("extractMatrix", "ChipEffectFile", function(this, ..., field=c("theta", "sdTheta")) {
# Argument 'field':
field <- match.arg(field)
NextMethod("extractMatrix", field=field)
})
HenrikBengtsson/aroma.affymetrix documentation built on Feb. 20, 2024, 9:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
###########################################################################/**
# @RdocClass ChipEffectFile
#
# @title "The ChipEffectFile class"
#
# \description{
# @classhierarchy
#
# This class represents estimates of chip effects in the probe-level models.
# }
#
# @synopsis
#
# \arguments{
# \item{...}{Arguments passed to @see "ParameterCelFile".}
# \item{probeModel}{The specific type of model, e.g. \code{"pm"}.}
# }
#
# \section{Fields and Methods}{
# @allmethods "public"
# }
#
# @author "HB"
#
# \seealso{
# An object of this class is typically obtained through the
# \code{getChipEffectSet()} method for the @see "ProbeLevelModel" class.
# An object of this class is typically part of a @see "ChipEffectSet".
# }
#*/###########################################################################
setConstructorS3("ChipEffectFile", function(..., probeModel=c("pm")) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'probeModel':
probeModel <- match.arg(probeModel)
this <- extend(ParameterCelFile(...), "ChipEffectFile",
"cached:.firstCells" = NULL,
probeModel = probeModel
)
setEncodeFunction(this, function(groupData, ...) {
theta <- .subset2(groupData, "theta")
stdvs <- .subset2(groupData, "sdTheta")
outliers <- .subset2(groupData, "thetaOutliers")
pixels <- NULL
if (!is.null(outliers))
pixels <- -as.integer(outliers)
res <- list()
if (!is.null(theta))
res$intensities <- theta
if (!is.null(stdvs))
res$stdvs <- stdvs
if (!is.null(pixels))
res$pixels <- pixels
res
})
## Not safe, at least not what I know of. /HB 2007-08-17
## setEncodeFunction(this, function(groupData, ...) {
## groupData[[3]] <- -as.integer(.subset2(groupData, 3))
## names(groupData) <- c("intensities", "stdvs", "pixels")
## groupData
## })
setDecodeFunction(this, function(groupData, ...) {
res <- list()
if (!is.null(groupData$intensities))
res$theta <- groupData$intensities
if (!is.null(groupData$stdvs))
res$sdTheta <- groupData$stdvs
if (!is.null(groupData$pixels))
res$thetaOutliers <- as.logical(-groupData$pixels)
res
})
# Parse attributes (all subclasses must call this in the constructor).
setAttributesByTags(this)
this
})
setMethodS3("as.character", "ChipEffectFile", function(x, ...) {
# To please R CMD check
this <- x
s <- NextMethod("as.character")
s <- c(s, sprintf("Parameters: %s", getParametersAsString(this)))
s
}, protected=TRUE)
setMethodS3("getParameters", "ChipEffectFile", function(this, ...) {
params <- NextMethod("getParameters")
params$probeModel <- this$probeModel
params
}, protected=TRUE)
setMethodS3("createParamCdf", "ChipEffectFile", function(static, sourceCdf, ..., verbose=FALSE) {
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
verbose && enter(verbose, "Creating CDF for chip effects")
verbose && cat(verbose, "Source chip type: ", getChipType(sourceCdf))
verbose && cat(verbose, "Source CDF: ", getPathname(sourceCdf))
# Search for existing monocell CDF
for (sep in c(",", "-")) {
chipType <- paste(getChipType(sourceCdf), "monocell", sep=sep)
verbose && cat(verbose, "Looking for chip type: ", chipType)
pathname <- AffymetrixCdfFile$findByChipType(chipType)
if (!is.null(pathname)) {
verbose && cat(verbose, "Found: ", pathname)
break
}
}
if (is.null(pathname)) {
verbose && cat(verbose, "Pathname: Not found!")
verbose && cat(verbose, "Will create CDF for the chip-effect files from the original CDF. NOTE: This will take several minutes or more!")
verbose && enter(verbose, "Creating CDF")
cdf <- createMonocellCdf(sourceCdf, verbose=less(verbose))
verbose && exit(verbose)
} else {
verbose && cat(verbose, "Pathname: ", pathname)
cdf <- AffymetrixCdfFile$fromFile(pathname)
}
verbose && exit(verbose)
cdf
}, static=TRUE, private=TRUE)
setMethodS3("readUnits", "ChipEffectFile", function(this, units=NULL, cdf=NULL, ..., force=FALSE, cache=FALSE, verbose=FALSE) {
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
# Check for cached data
key <- list(method="readUnits", class=class(this)[1],
pathname=getPathname(this),
cdf=cdf, units=units, ...)
if (getOption(aromaSettings, "devel/useCacheKeyInterface", FALSE)) {
key <- getCacheKey(this, method="readUnits", pathname=getPathname(this), cdf=cdf, units=units, ...)
}
id <- getChecksum(key)
res <- this$.readUnitsCache[[id]]
if (!force && !is.null(res)) {
verbose && cat(verbose, "readUnits.ChipEffectFile(): Returning cached data")
return(res)
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Retrieve the data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if (is.null(cdf)) {
cdf <- getCellIndices(this, units=units, verbose=less(verbose))
}
# Note that the actually call to the decoding is done in readUnits()
# of the superclass.
res <- NextMethod("readUnits", cdf=cdf, force=force, verbose=less(verbose))
# Store read units in cache?
if (cache) {
verbose && cat(verbose, "readUnits.ChipEffectFile(): Updating cache")
this$.readUnitsCache <- list()
this$.readUnitsCache[[id]] <- res
}
res
})
###########################################################################/**
# @RdocMethod getCellIndices
#
# @title "Retrieves tree list of cell indices for a set of units"
#
# \description{
# @get "title" from the associated CDF.
# }
#
# @synopsis
#
# \arguments{
# \item{...}{Additional arguments passed to \code{getCellIndices()}
# of @see "AffymetrixCdfFile".}
# \item{.cache}{Ignored.}
# }
#
# \value{
# Returns a @list structure, where each element corresponds to a unit.
# If argument \code{unlist=TRUE} is passed, an @integer @vector is returned.
# }
#
# \seealso{
# @seeclass
# }
#
# @keyword internal
#*/###########################################################################
setMethodS3("getCellIndices", "ChipEffectFile", function(this, ..., .cache=TRUE) {
cdf <- getCdf(this)
getCellIndices(cdf, ...)
}, protected=TRUE)
setMethodS3("updateUnits", "ChipEffectFile", function(this, units=NULL, cdf=NULL, data, ...) {
if (is.null(cdf))
cdf <- getCellIndices(this, units=units)
# Note that the actually call to the encoding is done in updateUnits()
# of the superclass.
NextMethod("updateUnits", cdf=cdf, data=data)
}, private=TRUE)
setMethodS3("findUnitsTodo", "ChipEffectFile", function(this, units=NULL, ..., force=FALSE, verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Identifying non-fitted units in chip-effect file")
verbose && cat(verbose, "Pathname: ", getPathname(this))
idxs <- NULL
if (is.null(units)) {
# Look up chip-type and parameter specific but data set independent data
cdf <- getCdf(this)
chipType <- getChipType(cdf)
key <- list(method="findUnitsTodo", class=class(this)[1],
chipType=chipType, params=getParameters(this))
dirs <- c("aroma.affymetrix", chipType)
if (!force) {
idxs <- loadCache(key, dirs=dirs)
if (!is.null(idxs))
verbose && cat(verbose, "Found indices cached on file")
}
}
if (is.null(idxs)) {
verbose && enter(verbose, "Identifying CDF units")
units0 <- units
if (is.null(units)) {
cdf <- getCdf(this)
units <- seq_len(nbrOfUnits(cdf))
}
nbrOfUnits <- length(units)
idxs <- lapplyInChunks(units, function(unitsChunk, ...) {
verbose && enter(verbose, "Reading CDF cell indices")
idxsChunk <- getCellIndices(this, units=unitsChunk, force=TRUE, verbose=less(verbose))
# Save memory 100% -> 94%
names(idxsChunk) <- NULL
verbose && exit(verbose)
verbose && enter(verbose, "Extracting first CDF group for each unit")
# Save memory 100% -> 94% -> 5% (all the nested named lists costs!)
idxsChunk <- lapply(idxsChunk, FUN=function(unit) {
groups <- .subset2(unit, "groups")
fields <- .subset2(groups, 1)
.subset2(fields, 1)
})
verbose && exit(verbose)
gc <- gc()
idxsChunk
}, chunkSize=100e3, useNames=FALSE, verbose=verbose)
# Not needed anymore
units <- units0
# Not needed anymore
units0 <- NULL
idxs <- unlist(idxs, use.names=FALSE)
gc <- gc()
verbose && print(verbose, gc)
# Verify correctness
if (length(idxs) != nbrOfUnits) {
throw("Internal error: Expected ", nbrOfUnits, " cell indices, but got ", length(idxs), ".")
}
if (is.null(units)) {
verbose && enter(verbose, "Saving to file cache")
saveCache(idxs, key=key, dirs=dirs)
verbose && exit(verbose)
}
verbose && exit(verbose)
}
# Read one cell from each unit
verbose && enter(verbose, "Reading data for these ", length(idxs), " cells")
value <- .readCel(getPathname(this), indices=idxs, readIntensities=FALSE,
readStdvs=TRUE, readPixels=FALSE)$stdvs
verbose && exit(verbose)
# Identify units for which the stdvs <= 0.
value <- which(value <= 0)
if (!is.null(units))
value <- units[value]
verbose && cat(verbose, "Looking for stdvs <= 0 indicating non-estimated units:")
verbose && str(verbose, value)
verbose && exit(verbose)
value
})
###########################################################################/**
# @RdocMethod getUnitGroupCellMap
# @aliasmethod getCellMap
#
# @title "Gets a (unit, group, cell) index map"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{units}{The units for which the map should be returned.
# If @NULL, all units are considered.}
# \item{force}{If @TRUE, cached cell indices are ignored.}
# \item{...}{Not used.}
# \item{verbose}{See @see "R.utils::Verbose".}
# }
#
# \value{
# Returns a @data.frame with @integer columns \code{unit}, \code{group},
# and \code{cell}.
# }
#
# \examples{\dontrun{
# unit group cell
# 1 104 1 335
# 2 104 2 336
# 3 105 1 337
# 4 105 2 338
# 5 105 3 339
# 6 105 4 340
# }}
#
# \seealso{
# @seeclass
# }
#*/###########################################################################
setMethodS3("getUnitGroupCellMap", "ChipEffectFile", function(this, units=NULL, force=FALSE, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'units':
if (inherits(units, "UnitGroupCellMap")) {
return(units)
} else if (is.null(units)) {
} else if (is.list(units)) {
} else {
units <- Arguments$getIndices(units)
}
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Retrieving unit-to-cell map")
# Special case: requesting zero units?
if (length(units) == 0 && !is.null(units)) {
map <- data.frame(unit=integer(0), group=integer(0), cell=integer(0))
class(map) <- c("UnitGroupCellMap", class(map))
verbose && exit(verbose)
return(map)
}
# Get the CDF
cdf <- getCdf(this)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Check cache
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# useFileCache <- (is.null(units) || (!is.list(units) && length(units) > 10000))
useFileCache <- (is.null(units) || (!is.list(units) && length(units) > 100))
if (useFileCache) {
chipType <- getChipType(cdf)
# Look up chip-type and parameter specific but data set independent data
key <- list(method="getUnitGroupCellMap", class=class(this)[1],
chipType=chipType, params=getParameters(this),
units=units)
dirs <- c("aroma.affymetrix", chipType)
if (!force) {
map <- loadCache(key, dirs=dirs)
if (!is.null(map)) {
verbose && cat(verbose, "Found (unit,group,cell) map cached on file")
verbose && exit(verbose)
return(map)
}
}
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Main
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Is 'units' already a CDF list?
if (is.list(units)) {
# No fancy validation for now.
cells <- units
units <- indexOf(cdf, names=names(units))
if (any(is.na(units))) {
throw("Argument 'units' is of unknown structure.")
}
verbose && enter(verbose, "Argument 'cells' is already a CDF cell-index structure")
# Get the unit names
unitNames <- names(cells)
# Get the number of groups per unit
unitSizes <- lapply(cells, FUN=function(unit) {
length(.subset2(unit, "groups"))
})
unitSizes <- unlist(unitSizes, use.names=FALSE)
cells <- unlist(cells, use.names=FALSE)
} else {
verbose && enter(verbose, "Retrieving cell indices for specified units")
if (is.null(units))
units <- seq_len(nbrOfUnits(cdf))
chunks <- splitInChunks(units, chunkSize=100e3)
nbrOfChunks <- length(chunks)
nbrOfUnits <- length(units)
# Store unit names and unit sizes
unitNames <- vector("character", nbrOfUnits)
unitSizes <- vector("integer", nbrOfUnits)
# Store cell indices first chunk-by-chunk, then as a vector.
cells <- vector("list", nbrOfChunks)
offset <- 0
for (kk in seq_len(nbrOfChunks)) {
verbose && printf(verbose, "Chunk #%d of %d\n", kk, length(chunks))
chunk <- chunks[[kk]]
chunks[[kk]] <- NA
cells0 <- getCellIndices(this, units=chunk, force=force, .cache=FALSE, verbose=less(verbose))
idxs <- offset + seq_along(chunk)
offset <- offset + length(chunk)
# Not needed anymore
chunk <- NULL
unitNames[idxs] <- names(cells0)
names(cells0) <- NULL
unitSizes0 <- lapply(cells0, FUN=function(unit) {
length(.subset2(unit, "groups"))
})
unitSizes[idxs] <- unlist(unitSizes0, use.names=FALSE)
# Not needed anymore
unitSizes0 <- NULL
cells[[kk]] <- unlist(cells0, use.names=FALSE)
# Not needed anymore
cells0 <- idxs <- NULL
}
# Not needed anymore
chunks <- NULL
gc <- gc()
verbose && print(verbose, gc)
}
cells <- unlist(cells, use.names=FALSE)
gc <- gc()
verbose && exit(verbose)
verbose && enter(verbose, "Creating return data frame")
uUnitSizes <- sort(unique(unitSizes))
verbose && cat(verbose, "Unique number of groups per unit: ",
paste(uUnitSizes, collapse=","))
verbose && cat(verbose, "Number of units: ", length(unitNames))
if (is.null(units))
units <- seq_len(nbrOfUnits(cdf))
# The following is too slow:
# groups <- sapply(unitSizes, FUN=function(n) seq_len(n))
# Instead, updated size by size
verbose && printf(verbose, "Allocating matrix of size %dx%d.\n",
max(uUnitSizes), length(unitNames))
naValue <- NA_integer_
units2 <- groups <- matrix(naValue, nrow=max(uUnitSizes), ncol=length(unitNames))
for (size in uUnitSizes) {
# Identify units with a certain number of groups
cc <- which(unitSizes == size)
# Assign group indices to the groups
seq <- seq_len(size)
groups[seq,cc] <- seq
# Assign unit indices to ditto
units2[seq,cc] <- rep(units[cc], each=size)
# Not needed anymore
seq <- NULL
gc <- gc()
}
keep <- !is.na(groups)
groups <- groups[keep]
units2 <- units2[keep]
# Not needed anymore
keep <- NULL
gc <- gc()
map <- data.frame(unit=units2, group=groups, cell=cells)
verbose && exit(verbose)
verbose && exit(verbose)
class(map) <- c("UnitGroupCellMap", class(map))
if (useFileCache) {
verbose && enter(verbose, "Saving to file cache")
saveCache(map, key=key, dirs=dirs)
verbose && exit(verbose)
}
map
}, private=TRUE)
setMethodS3("getUnitGroupCellChromosomePositionMap", "ChipEffectFile", function(this, units=NULL, chromosomes=NULL, orderByPosition=TRUE, ..., force=FALSE, verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
cdf <- getCdf(this)
# Argument 'units':
ugcMap <- NULL
if (is.null(units)) {
} else if (isUnitGroupCellMap(units)) {
ugcMap <- units
units <- ugcMap[,"unit"]
}
if (!is.null(units)) {
units <- Arguments$getIndices(units, range=c(1, nbrOfUnits(cdf)))
}
units0 <- units
# Get the genome position information
gi <- getGenomeInformation(cdf)
# Argument 'chromosomes':
if (!is.null(chromosomes)) {
allChromosomes <- getChromosomes(gi)
unknown <- chromosomes[!(chromosomes %in% allChromosomes)]
if (length(unknown) > 0) {
throw("Argument 'chromosomes' contains unknown values: ",
paste(unknown, collapse=", "))
}
}
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Getting (unit, group, cell, chromosome, position) map")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Check for cached results
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Look for results in file cache
verbose && enter(verbose, "Checking cache")
chipType <- getChipType(cdf)
key <- list(method="getUnitGroupCellChromosomePositionMap",
class=class(this)[1],
chipType=chipType, units=units, ugcMap=ugcMap,
chromosomes=chromosomes, orderByPosition=orderByPosition)
dirs <- c("aroma.affymetrix", chipType)
if (!force) {
map <- loadCache(key=key, dirs=dirs)
if (!is.null(map)) {
verbose && cat(verbose, "Found cached results")
verbose && exit(verbose)
return(map)
}
}
# Select by chromosome(s)?
if (!is.null(chromosomes)) {
verbose && cat(verbose, "Units:")
verbose && str(verbose, units)
verbose && cat(verbose, "Subset by chromosomes:")
verbose && str(verbose, chromosomes)
units <- getUnitsOnChromosomes(gi, chromosomes)
verbose && cat(verbose, "Units:")
verbose && str(verbose, units)
if (!is.null(units0)) {
units <- intersect(units, units0)
}
}
verbose && cat(verbose, "Units:")
verbose && str(verbose, units)
# Get the (unit, group, cell) map?
if (!isUnitGroupCellMap(ugcMap)) {
ugcMap <- getUnitGroupCellMap(this, units=units, force=force, verbose=less(verbose, 10))
verbose && cat(verbose, "(unit, group, cell) map:")
verbose && str(verbose, ugcMap)
}
# Get the (chromosome, position) map
cpMap <- getData(gi, units=ugcMap[,"unit"], force=force, verbose=less(verbose, 10))
verbose && cat(verbose, "(chromosome, position) map:")
verbose && str(verbose, cpMap)
# Sanity check
stopifnot(nrow(ugcMap) == nrow(cpMap))
# Merge the two maps
map <- cbind(ugcMap, cpMap)
# Not needed anymore
ugcMap <- cpMap <- NULL
if (orderByPosition) {
o <- with(map, order(chromosome, physicalPosition))
map <- map[o,,drop=FALSE]
# Not needed anymore
o <- NULL
verbose && cat(verbose, "Reordered by genomic position")
}
rownames(map) <- NULL
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Save to cache
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Save only results > 50kB
if (object.size(map) > 50e3) {
saveCache(map, key=key, dirs=dirs)
verbose && cat(verbose, "Saved to file cache")
}
verbose && exit(verbose)
map
}, private=TRUE)
setMethodS3("getDataFlat", "ChipEffectFile", function(this, units=NULL, fields=c("theta", "sdTheta", "outliers"), ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Retrieving data as a flat data frame")
# Get unit-to-cell map
suppressWarnings({
map <- getUnitGroupCellMap(this, units=units, ..., verbose=less(verbose))
})
verbose && enter(verbose, "Reading data fields")
celFields <- c(theta="intensities", sdTheta="stdvs", outliers="pixels")
suppressWarnings({
data <- getData(this, indices=map[,"cell"], fields=celFields[fields])
})
rownames(data) <- seq_len(nrow(data)); # Work around?!? /HB 2006-11-28
# Decode
names <- colnames(data)
names <- gsub("intensities", "theta", names)
names <- gsub("stdvs", "sdTheta", names)
names <- gsub("pixels", "outliers", names)
colnames(data) <- names
verbose && str(verbose, data)
if ("outliers" %in% names) {
data[,"outliers"] <- as.logical(-data[,"outliers"])
}
verbose && exit(verbose)
len <- sapply(data, FUN=length)
keep <- (len == nrow(map))
data <- data[keep]
data <- as.data.frame(data)
data <- cbind(map, data)
verbose && exit(verbose)
data
}, private=TRUE)
setMethodS3("updateDataFlat", "ChipEffectFile", function(this, data, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'data':
names <- colnames(data)
namesStr <- paste(names, collapse=", ")
if (!"cell" %in% names)
throw("Argument 'data' must contain a column 'cell': ", namesStr)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose2 <- -as.integer(verbose)-2
verbose && enter(verbose, "Storing flat data to file")
# Encode?
if ("outliers" %in% names) {
data[,"outliers"] <- -as.integer(data[,"outliers"])
}
names <- gsub("theta", "intensities", names)
names <- gsub("sdTheta", "stdvs", names)
names <- gsub("outliers", "pixels", names)
colnames(data) <- names
verbose && enter(verbose, "Updating file")
indices <- data[,"cell"]
keep <- (names %in% c("intensities", "stdvs", "pixels"))
data <- data[,keep]
pathname <- getPathname(this)
pathname <- Arguments$getWritablePathname(pathname)
.updateCel(pathname, indices=indices, data, verbose=verbose2)
verbose && exit(verbose)
verbose && exit(verbose)
invisible(data)
}, private=TRUE)
setMethodS3("mergeGroups", "ChipEffectFile", function(this, fcn, fields=c("theta", "sdTheta"), ..., pathname, overwrite=FALSE, verbose=FALSE) {
# Argument 'fcn':
if (!is.function(fcn)) {
throw("Argument 'fcn' is not a function: ", class(fcn)[1])
}
# Argument 'pathname':
pathname <- Arguments$getWritablePathname(pathname, mustNotExist=!overwrite)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Merging groups")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Test 'fcn':
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && enter(verbose, "Testing merge function")
for (size in 1:10) {
y <- matrix(1000+1:(size*4), nrow=size)
yOut <- fcn(y)
if (!identical(dim(yOut), dim(y))) {
throw("Function 'fcn' must not change the dimension of the data: ",
paste(dim(yOut), collapse="x"), " != ",
paste(dim(y), collapse="x"))
}
}
verbose && exit(verbose)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Get flat (unit, group, cell) map
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
map <- getUnitGroupCellMap(this, verbose=less(verbose))
verbose && str(verbose, map)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Merge data for each unit size separately (in reverse order!)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Get all possible unit sizes (number of groups per unit)
uSizes <- sort(unique(data[,"group"]))
verbose && cat(verbose, "Different number of groups per unit identified:")
verbose && print(verbose, uSizes)
data <- getDataFlat(this, units=map, ..., verbose=less(verbose))
verbose && str(verbose, data)
for (size in rev(uSizes)) {
verbose && enter(verbose, "Unit size ", size)
# Identify the units of that size
idxs <- which(data[,"group"] == size)
unitsS <- data[idxs, "unit"]
# Get the subset of the data for such units
idxs <- which(data[,"unit"] %in% unitsS)
for (field in fields) {
# Extract signals as a matrix where each column is one unit
y <- data[idxs, field]
y <- matrix(y, nrow=size)
verbose && str(verbose, y)
y <- fcn(y)
verbose && str(verbose, y)
# Update data table
data[idxs, field] <- as.vector(y)
}
# Get the cells where to store the merged data
verbose && exit(verbose)
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Copy CEL file and update the copy
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && enter(verbose, "Storing merged data")
verbose && cat(verbose, "Pathname: ", pathname)
# Create CEL file to store results, if missing
verbose && enter(verbose, "Creating CEL file for results, if missing")
cfN <- createFrom(this, filename=pathname, path=NULL, methods="create", clear=TRUE, verbose=less(verbose))
verbose && print(verbose, cfN)
verbose && exit(verbose)
verbose && enter(verbose, "Writing merged data")
updateDataFlat(cfN, data=data, verbose=less(verbose))
verbose && exit(verbose)
verbose && exit(verbose)
verbose && exit(verbose)
cfN
}, protected=TRUE)
setMethodS3("extractMatrix", "ChipEffectFile", function(this, ..., field=c("theta", "sdTheta")) {
# Argument 'field':
field <- match.arg(field)
NextMethod("extractMatrix", field=field)
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.