###########################################################################/**
# @RdocClass MatNormalization
#
# @title "The MatNormalization class"
#
# \description{
# @classhierarchy
#
# This class represents a normalization method that corrects for systematic
# effects in the probe intensities due to differences in the number of
# A, C, G, and T:s and the match scores according to MAT [1].
# }
#
# @synopsis
#
# \arguments{
# \item{...}{Arguments passed to the constructor of
# @see "AbstractProbeSequenceNormalization".}
# \item{unitsToFit}{The units from which the normalization curve should
# be estimated. If @NULL, all are considered.}
# \item{model}{A @character string specifying the model used to fit
# the base-count effects.}
# \item{nbrOfBins}{The number of bins to use for the variance smoothing step.}
# }
#
# \section{Fields and Methods}{
# @allmethods "public"
# }
#
# \section{Requirements}{
# This class requires that an aroma probe sequence file and aroma
# match scores file is available for the chip type.
# }
#
# @author "MR"
#
# \references{
# [1] Johnson WE, Li W, Meyer CA, Gottardo R, Carroll JS, Brown M, Liu XS.
# \emph{Model-based analysis of tiling-arrays for ChIP-chip}, PNAS, 2006.
# }
#*/###########################################################################
setConstructorS3("MatNormalization", function(..., unitsToFit=NULL, model=c("lm"), nbrOfBins=200) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'model':
model <- match.arg(model)
# Argument 'nbrOfBins':
nbrOfBins <- Arguments$getInteger(nbrOfBins, range=c(1,Inf))
extend(AbstractProbeSequenceNormalization(..., unitsToFit=unitsToFit), "MatNormalization",
.model = model,
.scaleResiduals = TRUE,
.nbrOfBins = as.integer(nbrOfBins)
)
})
setMethodS3("getAromaCellMatchScoreFile", "MatNormalization", function(this, ..., force=FALSE) {
apm <- this$.apm
if (force || is.null(apm)) {
dataSet <- getInputDataSet(this)
cdf <- getCdf(dataSet)
chipType <- getChipType(cdf, fullname=FALSE)
nbrOfCells <- nbrOfCells(cdf)
apm <- AromaCellMatchScoreFile$byChipType(chipType, nbrOfCells=nbrOfCells, ...)
this$.apm <- apm
}
apm
}, protected=TRUE)
setMethodS3("getAsteriskTags", "MatNormalization", function(this, collapse=NULL, ...) {
tags <- NextMethod("getAsteriskTags", collapse=NULL)
# Add model tag?
model <- this$.model
tags <- c(tags, model)
# Collapse?
tags <- paste(tags, collapse=collapse)
tags
}, protected=TRUE)
setMethodS3("getParameters", "MatNormalization", function(this, ...) {
# Get parameters from super class
params <- NextMethod("getParameters")
params <- c(params, list(
model = this$.model,
nbrOfBins = this$.nbrOfBins
))
params
}, protected=TRUE)
setMethodS3("getDesignMatrix", "MatNormalization", function(this, cells=NULL, model=NULL, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'cells':
if (is.null(cells)) {
} else {
# Validated below...
}
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Retrieving design matrix")
verbose && cat(verbose, "Cells:")
verbose && str(verbose, cells)
verbose && enter(verbose, "Locating probe-sequence annotation data")
# Locate AromaCellSequenceFile holding probe sequences
aps <- getAromaCellSequenceFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Locating match scores annotation data")
# Locate AromaCellMatchScoreFile holding match scores
apm <- getAromaCellMatchScoreFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Reading sequence matrix")
sm <- readSequenceMatrix(aps, cells=cells, verbose=verbose)
verbose && exit(verbose)
verbose && enter(verbose, "Reading match scores")
ms <- readColumns(apm, rows=cells, verbose=verbose)
verbose && exit(verbose)
verbose && enter(verbose, "Constructing design matrix")
nT <- rowSums(sm == "T")
G <- (sm == "G")+0
A <- (sm == "A")+0
C <- (sm == "C")+0
designMatrix <- cbind(nT, A, C, G, rowSums(A)^2, rowSums(C)^2, rowSums(G)^2, nT^2, log(as.integer(ms[,1])))
# Garbage collect
# Not needed anymore
nT <- G <- A <- C <- ms <- sm <- NULL
gc <- gc()
verbose && print(verbose, gc)
verbose && str(verbose, designMatrix)
verbose && cat(verbose, "object.size(designMatrix): ",
object.size(designMatrix))
verbose && exit(verbose)
verbose && exit(verbose)
designMatrix
}, private=TRUE)
setMethodS3("fitOne", "MatNormalization", function(this, df, ram=NULL, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'ram':
ram <- getRam(aromaSettings, ram)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Locating probe-sequence annotation data")
# Locate AromaCellSequenceFile holding probe sequences
aps <- getAromaCellSequenceFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Locating match scores annotation data")
# Locate AromaCellMatchScoreFile holding match scores
apm <- getAromaCellMatchScoreFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Reading 'non-missing' cells to fit")
cellsToFit <- which( !(isMissing(aps) | isMissing(apm)) )
nbrOfCells <- length(cellsToFit)
verbose && cat(verbose, "Cells to fit:")
verbose && str(verbose, cellsToFit)
verbose && exit(verbose)
# this code adopted from Dave Fourniers 17/08/2007 post to r-help
# mailing list entitled "[R] Linear models over large datasets"
# Calculate the number of cells to process per chunk. It should:
# (1) increase with the 'ram' option.
# (2) there is only on array, so it should be independent of everything.
cellsPerChunk <- ceiling(ram * 1e6 + 1L)
cellsPerChunk <- Arguments$getInteger(cellsPerChunk, range=c(1,Inf))
## In aroma.affymetrix v1.9.4 and before, the number of cells processed
## per chunk increase with the total number of cells. /HB 2011-02-15
## cellsPerChunk <- ceiling(nbrOfCells/numChunks) + 1
nbrOfChunks <- ceiling(nbrOfCells / cellsPerChunk)
verbose && cat(verbose, "Number cells per chunk: ", cellsPerChunk)
verbose && enter(verbose, "Reading signals to fit")
y <- extractMatrix(df, cells=cellsToFit, verbose=less(verbose, 10))
verbose && exit(verbose)
verbose && enter(verbose, "Log2 transforming signals")
y <- log2(y)
verbose && cat(verbose, "Target log2 probe signals:")
verbose && str(verbose, y)
verbose && exit(verbose)
start <- 0L
xtx <- xty <- 0
chunk <- 1L
while (start < nbrOfCells) {
verbose && enter(verbose, sprintf("Chunk #%d of %d", chunk, nbrOfChunks))
verbose && enter(verbose, "Working on indices over range")
from <- start + 1L
to <- min(start+cellsPerChunk, nbrOfCells)
indSubset <- (from:to)
rng <- c(from, to)
rng <- rng / nbrOfCells
verbose && cat(verbose, sprintf("[%g,%g]", rng[1], rng[2]))
verbose && exit(verbose)
verbose && enter(verbose, "Reading design matrix for this iteration")
X <- getDesignMatrix(this, cells=cellsToFit[indSubset], verbose=verbose)
verbose && exit(verbose)
verbose && enter(verbose, "Calculating cross products")
xtx <- xtx + crossprod(X)
xty <- xty + crossprod(X, y[indSubset])
verbose && exit(verbose)
start <- start + cellsPerChunk
chunk <- chunk + 1L
verbose && exit(verbose)
} # while (...)
verbose && enter(verbose, "Solving normal equations")
#fit <- list(xtx=xtx, xty=xty) #,beta=solve(xtx, xty))
verbose && exit(verbose)
fit <- list(beta=solve(xtx, xty), scaleResiduals=this$.scaleResiduals)
fit
}, protected=TRUE)
setMethodS3("predictOne", "MatNormalization", function(this, fit, ram=NULL, ..., verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'ram':
ram <- getRam(aromaSettings, ram)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Locating probe-sequence annotation data")
# Locate AromaCellSequenceFile holding probe sequences
aps <- getAromaCellSequenceFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Locating match scores annotation data")
# Locate AromaCellMatchScoreFile holding match scores
apm <- getAromaCellMatchScoreFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Allocating mu vector")
mu <- double(nbrOfCells(aps))
verbose && str(verbose, mu)
verbose && exit(verbose)
verbose && enter(verbose, "Reading 'non-missing' cells to predict")
cellsToPredict <- which( !(isMissing(aps) | isMissing(apm)) )
nbrOfCells <- length(cellsToPredict)
verbose && cat(verbose, "Cells to predict:")
verbose && str(verbose, cellsToPredict)
verbose && exit(verbose)
# Calculate the number of cells to process per chunk. It should:
# (1) increase with the 'ram' option.
# (2) there is only on array, so it should be independent of everything.
cellsPerChunk <- ceiling(ram * 1e6 + 1L)
cellsPerChunk <- Arguments$getInteger(cellsPerChunk, range=c(1,Inf))
## In aroma.affymetrix v1.9.4 and before, the number of cells processed
## per chunk increase with the total number of cells. /HB 2011-02-15
## cellsPerChunk <- ceiling(nbrOfCells/numChunks) + 1
nbrOfChunks <- ceiling(nbrOfCells / cellsPerChunk)
verbose && cat(verbose, "Number cells per chunk: ", cellsPerChunk)
start <- 0L
chunk <- 1L
while (start < nbrOfCells) {
verbose && enter(verbose, sprintf("Chunk #%d of %d", chunk, nbrOfChunks))
verbose && enter(verbose, "Working on indices over range")
from <- start + 1L
to <- min(start+cellsPerChunk, nbrOfCells)
indSubset <- (from:to)
rng <- c(from, to)
rng <- rng / nbrOfCells
verbose && cat(verbose, sprintf("[%g,%g]", rng[1], rng[2]))
verbose && exit(verbose)
verbose && enter(verbose, "Reading design matrix for this iteration")
X <- getDesignMatrix(this, cells=cellsToPredict[indSubset], verbose=verbose)
verbose && exit(verbose)
verbose && enter(verbose, "Calculating mu")
mu[ cellsToPredict[indSubset] ] <- X %*% fit$beta
verbose && exit(verbose)
start <- start + cellsPerChunk
chunk <- chunk + 1L
verbose && exit(verbose)
} # while (...)
# Not needed anymore
X <- indSubset <- NULL
gc <- gc()
#nbrOfBins <- this$.nbrOfBins
#q <- quantile(mu[cellsToPredict],prob=(0:nbrOfBins)/nbrOfBins)
#cuts<-cut(mu[cellsToPredict],breaks=q,labels=1:(length(q)-1)) # define
#ss<-split(data.frame(resid),cuts)
#ssvar<-sapply(ss,var)
#v<-ssvar[as.character(cuts)]
#for(j in 1:length(b))
#rr<-resid/sqrt(v)
# Return results
mu
}, protected=TRUE)
###########################################################################/**
# @RdocMethod process
#
# @title "Normalizes the data set"
#
# \description{
# @get "title".
# }
#
# @synopsis
#
# \arguments{
# \item{...}{Not used.}
# \item{force}{If @TRUE, data already normalized is re-normalized,
# otherwise not.}
# \item{verbose}{See @see "R.utils::Verbose".}
# }
#
# \value{
# Returns a @double @vector.
# }
#
# \seealso{
# @seeclass
# }
#*/###########################################################################
setMethodS3("process", "MatNormalization", function(this, ..., ram=NULL, force=FALSE, verbose=FALSE) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'ram':
ram <- getRam(aromaSettings, ram)
# Argument 'verbose':
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
verbose && enter(verbose, "Normalizing data set for probe-sequence effects")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Already done?
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if (!force && isDone(this)) {
verbose && cat(verbose, "Already normalized")
verbose && exit(verbose)
outputDataSet <- getOutputDataSet(this)
return(invisible(outputDataSet))
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Setup
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Get input data set
ds <- getInputDataSet(this)
# Get (and create) the output path
outputPath <- getPath(this)
verbose && enter(verbose, "Locating probe-sequence annotation data")
# Locate AromaCellSequenceFile holding probe sequences
aps <- getAromaCellSequenceFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Locating match scores annotation data")
# Locate AromaCellMatchScoreFile holding match scores
apm <- getAromaCellMatchScoreFile(this, verbose=less(verbose, 5))
verbose && exit(verbose)
verbose && enter(verbose, "Reading 'non-missing' cells to fit")
cellsToFit <- which( !(isMissing(aps) | isMissing(apm)) )
verbose && cat(verbose, "Cells to fit:")
verbose && str(verbose, cellsToFit)
verbose && exit(verbose)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Normalize all arrays simultaneously
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
nbrOfArrays <- length(ds)
verbose && enter(verbose, "Normalizing ", nbrOfArrays, " arrays")
verbose && enter(verbose, "Path: ", outputPath)
df <- getOneFile(ds)
nbrOfCells <- nbrOfCells(df)
xtx <- 0
xtyList <- as.list(double(nbrOfArrays))
nbrOfCells <- length(cellsToFit)
# Calculate the number of cells to process per chunk. It should:
# (1) increase with the 'ram' option.
# (2) the arrays are processed sequentially, i.e. it is constant
# in number of arrays.
cellsPerChunk <- ceiling(ram * 1e6 + 1L)
cellsPerChunk <- Arguments$getInteger(cellsPerChunk, range=c(1,Inf))
## In aroma.affymetrix v1.9.4 and before, the number of cells processed
## per chunk increase with the total number of cells. /HB 2011-02-15
## cellsPerChunk <- ceiling(nbrOfCells/numChunks) + 1
nbrOfChunks <- ceiling(nbrOfCells / cellsPerChunk)
verbose && cat(verbose, "Number cells per chunk: ", cellsPerChunk)
idxs <- 1:nbrOfCells
head <- 1:cellsPerChunk
count <- 1L
while (length(idxs) > 0) {
verbose && enter(verbose, sprintf("Fitting chunk #%d of %d", count, nbrOfChunks))
if (length(idxs) < cellsPerChunk) {
head <- 1:length(idxs)
}
cc <- idxs[head]
verbose && cat(verbose, "Cells: ")
verbose && str(verbose, cellsToFit[cc])
verbose && enter(verbose, "Reading design matrix")
X <- getDesignMatrix(this, cells=cellsToFit[cc], verbose=verbose)
verbose && exit(verbose)
verbose && enter(verbose, "Calculating cross product X'X")
xtx <- xtx + crossprod(X)
verbose && exit(verbose)
verbose && enter(verbose, "Calculating cross product X'y for each array")
for (ii in seq_len(nbrOfArrays)) {
df <- ds[[ii]]
verbose && enter(verbose, sprintf("Array #%d ('%s') of %d",
ii, getName(df), nbrOfArrays))
y <- extractMatrix(df, cells=cellsToFit[cc], verbose=verbose)
y <- log2(y)
verbose && cat(verbose, "Target log2 probe signals:")
verbose && str(verbose, y)
xtyList[[ii]] <- xtyList[[ii]] + crossprod(X, y)
# Not needed anymore
y <- NULL
verbose && exit(verbose)
} # for (ii ...)
verbose && exit(verbose)
# Not needed anymore
X <- NULL
# Next chunk
idxs <- idxs[-head]
count <- count + 1L
verbose && exit(verbose)
} # while (...)
verbose && enter(verbose, "Solving for each array")
fits <- lapply(xtyList, FUN=function(xty) {
list(beta=solve(xtx, xty))
})
verbose && exit(verbose)
# Not needed anmore
# Not needed anymore
xtx <- xtyList <- NULL
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Save model fits
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
for (ii in seq_len(nbrOfArrays)) {
df <- ds[[ii]]
verbose && enter(verbose, sprintf("Array #%d ('%s') of %d",
ii, getName(df), nbrOfArrays))
fullname <- getFullName(df)
# Store model fit
verbose && enter(verbose, "Saving model fit")
# Store fit and parameters (in case someone are interested in looking
# at them later; no promises of backward compatibility though).
filename <- sprintf("%s,fit.RData", fullname)
fitPathname <- Arguments$getWritablePathname(filename,
path=outputPath, ...)
saveObject(fits[[ii]], file=fitPathname)
verbose && str(verbose, fits[[ii]], level=-50)
verbose && exit(verbose)
# Garbage collect
gc <- gc()
verbose && print(verbose, gc)
verbose && exit(verbose)
} # for (ii ...)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Create output CEL files
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && enter(verbose, "Creating (temporary) CEL files")
pathnamesT <- listenv()
for (ii in seq_len(nbrOfArrays)) {
df <- ds[[ii]]
verbose && enter(verbose, sprintf("Array #%d ('%s') of %d",
ii, getName(df), nbrOfArrays))
fullname <- getFullName(df)
filename <- sprintf("%s.CEL", fullname)
pathname <- Arguments$getWritablePathname(filename, path=outputPath, ...)
# Write to a temporary file (allow rename of existing one if forced)
isFile <- isFile(pathname)
pathnameT <- pushTemporaryFile(pathname, isFile=isFile, verbose=verbose)
# Create CEL file to store results, if missing
verbose && enter(verbose, "Creating CEL file for results, if missing")
createFrom(df, filename=pathnameT, path=NULL, verbose=less(verbose))
verbose && exit(verbose)
# Record temporary filename and reuse below
pathnamesT[[ii]] <- pathnameT
verbose && exit(verbose)
} ## for (ii ...)
## Assert proper temporary files
pathnamesT <- unlist(pathnamesT)
stopifnot(length(pathnamesT) == nbrOfArrays)
verbose && exit(verbose)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Predict
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
nbrOfCells <- length(cellsToFit)
# Calculate the number of cells to process per chunk. It should:
# (1) increase with the 'ram' option.
# (2) the arrays are processed sequentially, i.e. it is constant
# in number of arrays.
cellsPerChunk <- ceiling(ram * 1e6 + 1L)
cellsPerChunk <- Arguments$getInteger(cellsPerChunk, range=c(1,Inf))
## In aroma.affymetrix v1.9.4 and before, the number of cells processed
## per chunk increase with the total number of cells. /HB 2011-02-15
## cellsPerChunk <- ceiling(nbrOfCells/numChunks) + 1
nbrOfChunks <- ceiling(nbrOfCells / cellsPerChunk)
verbose && cat(verbose, "Number cells per chunk: ", cellsPerChunk)
start <- 0L
xtx <- 0
mu <- vector("list", length=nbrOfArrays)
count <- 1L
while (start < nbrOfCells) {
verbose && enter(verbose, sprintf("Fitting chunk #%d of %d", count, nbrOfChunks))
verbose && enter(verbose, "Working on indices over range")
from <- start + 1L
to <- min(start+cellsPerChunk, nbrOfCells)
indSubset <- (from:to)
rng <- c(from, to)
rng <- rng / nbrOfCells
verbose && cat(verbose, sprintf("[%g,%g]", rng[1], rng[2]))
verbose && exit(verbose)
verbose && enter(verbose, "Set of cells")
cellsChunk <- cellsToFit[indSubset]
verbose && str(verbose, cellsChunk)
verbose && exit(verbose)
verbose && enter(verbose, "Reading design matrix for this iteration")
X <- getDesignMatrix(this, cells=cellsChunk, verbose=verbose)
verbose && exit(verbose)
verbose && enter(verbose, "Calculating model fits")
xtx <- xtx + crossprod(X)
verbose && enter(verbose, "Processing ", nbrOfArrays, " arrays")
for (ii in seq_len(nbrOfArrays)) {
df <- ds[[ii]]
verbose && enter(verbose, sprintf("Array #%d ('%s') of %d",
ii, getName(df), nbrOfArrays))
mu <- X %*% fits[[ii]]$beta
mu <- as.double(mu)
verbose && str(verbose, mu)
verbose && enter(verbose, "Updating temporary CEL file")
pathnameT <- pathnamesT[ii]
verbose2 <- as.logical(verbose)
.updateCel(pathnameT, indices=cellsChunk, intensities=2^mu, verbose=verbose2)
verbose && exit(verbose)
verbose && exit(verbose)
} # for (ii ...)
verbose && exit(verbose)
# Not needed anymore
mu <- NULL
# Garbage collect
gc <- gc()
verbose && print(verbose, gc)
verbose && exit(verbose)
start <- start + cellsPerChunk
#start <- nbrOfCells + 1
count <- count + 1L
verbose && exit(verbose)
} # while (...)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Scale residuals
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && enter(verbose, "Scaling residuals")
nbrOfBins <- this$.nbrOfBins
verbose && cat(verbose, "Number of bins: ", nbrOfBins)
probs <- (0:nbrOfBins) / nbrOfBins
verbose && cat(verbose, "Quantile probabilities:")
verbose && str(verbose, probs)
cutLabels <- seq_len(nbrOfBins)
res <- listenv()
for (ii in seq_len(nbrOfArrays)) {
df <- ds[[ii]]
verbose && enter(verbose, "Binning predicted values, calculating and scaling residuals")
pathnameT <- pathnamesT[ii]
res[[ii]] %<-% {
## Original probe intensities
y <- extractMatrix(df, cells=cellsToFit, verbose=verbose)
y <- log2(y)
## Prediction
mu <- .readCel(pathnameT, indices=cellsToFit, readOutliers=FALSE, readHeader=FALSE, readMasked=FALSE, verbose=less(verbose,10))$intensities
mu <- log2(mu)
## Residuals
r <- y - mu
## Normalize residuals
q <- quantile(mu, probs=probs)
cuts <- cut(mu, breaks=q, labels=cutLabels) # define
ss <- split(r, cuts, drop=FALSE)
ssvar <- sapply(ss, FUN=var)
v <- ssvar[as.character(cuts)]
r <- r / sqrt(v)
r <- as.double(r)
#return(list(y=y,mu=mu,r=r))
verbose && enter(verbose, "Updating temporary CEL file")
verbose2 <- as.logical(verbose)
.updateCel(pathnameT, indices=cellsToFit, intensities=2^r, verbose=verbose2)
verbose && exit(verbose)
# Not needed anymore
q <- ss <- ssvar <- v <- r <- y <- NULL
gc <- gc()
verbose && print(verbose, gc)
# Rename temporary file
pathname <- popTemporaryFile(pathnameT, verbose=verbose)
## Create checksum file
dfZ <- getChecksumFile(pathname)
pathname
} ## %<-%
verbose && exit(verbose)
} # for (ii ...)
## Resolve futures
res <- as.list(res)
res <- NULL
verbose && exit(verbose)
outputDataSet <- getOutputDataSet(this, force=TRUE)
verbose && exit(verbose)
invisible(outputDataSet)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.