[Robast-commits] r257 - in pkg: . RobLoxBioC RobLoxBioC/R RobLoxBioC/inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 2 15:19:35 CET 2009
Author: stamats
Date: 2009-03-02 15:19:35 +0100 (Mon, 02 Mar 2009)
New Revision: 257
Added:
pkg/RobLoxBioC/
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/NAMESPACE
pkg/RobLoxBioC/R/
pkg/RobLoxBioC/R/robloxAffyBatch.R
pkg/RobLoxBioC/R/robloxMatrix.R
pkg/RobLoxBioC/R/sysdata.rda
pkg/RobLoxBioC/inst/
pkg/RobLoxBioC/inst/CITATION
pkg/RobLoxBioC/inst/NEWS
pkg/RobLoxBioC/man/
Log:
added first version of new package RobLoxBioC
Added: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION (rev 0)
+++ pkg/RobLoxBioC/DESCRIPTION 2009-03-02 14:19:35 UTC (rev 257)
@@ -0,0 +1,13 @@
+Package: RobLoxBioC
+Version: 0.1
+Date: 2009-03-02
+Title: Infinitesimally robust estimators for preprocessing omics data
+Description: Functions for the determination of optimally robust
+ influence curves and estimators for preprocessing omics data,
+ in particular gene expression data.
+Depends: R(>= 2.8.1), methods, Biobase, affy
+Author: Matthias Kohl <Matthias.Kohl at stamats.de>
+Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
+LazyLoad: yes
+License: LGPL-3
+URL: http://robast.r-forge.r-project.org/
Property changes on: pkg/RobLoxBioC/DESCRIPTION
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/RobLoxBioC/NAMESPACE
===================================================================
--- pkg/RobLoxBioC/NAMESPACE (rev 0)
+++ pkg/RobLoxBioC/NAMESPACE 2009-03-02 14:19:35 UTC (rev 257)
@@ -0,0 +1,3 @@
+import("methods")
+
+exportMethods("robloxbioc")
Property changes on: pkg/RobLoxBioC/NAMESPACE
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/RobLoxBioC/R/robloxAffyBatch.R
===================================================================
--- pkg/RobLoxBioC/R/robloxAffyBatch.R (rev 0)
+++ pkg/RobLoxBioC/R/robloxAffyBatch.R 2009-03-02 14:19:35 UTC (rev 257)
@@ -0,0 +1,86 @@
+###############################################################################
+## Use robloxbioc to preprocess Affymetrix-data - comparable to MAS 5.0
+###############################################################################
+setMethod("robloxbioc", signature(x = "AffyBatch"),
+ function(x, pmcorrect = "roblox", verbose = TRUE,
+ eps, eps.lower = 0, eps.upper = 0.1, steps = 1L, mad0 = 1e-4,
+ contrast.tau = 0.03, scale.tau = 10, delta = 2^(-20)) {
+ n <- length(x)
+ ids <- featureNames(x)
+ m <- length(ids)
+
+ CDFINFO <- getCdfInfo(x)
+ INDEX <- sapply(ids, get, envir = CDFINFO)
+ NROW <- unlist(lapply(INDEX, nrow))
+ nr <- as.integer(names(table(NROW)))
+
+ diff.log2 <- function(INDEX, x){
+ l.pm <- INDEX[,1]
+ if(ncol(INDEX) == 2)
+ l.mm <- INDEX[,2]
+ else
+ l.mm <- integer()
+ log2(x[l.pm, , drop = FALSE]) - log2(x[l.mm, , drop = FALSE])
+ }
+ pm.only <- function(INDEX, x){
+ l.pm <- INDEX[,1]
+ x[l.pm, , drop = FALSE]
+ }
+
+ intensData <- intensity(x)
+ rob.est <- matrix(NA, ncol = 2, nrow = m*n)
+ if(verbose) cat("PM/MM correcting ...")
+ if(pmcorrect == "roblox"){
+ res <- lapply(INDEX, diff.log2, x = intensData)
+ rob.est1 <- matrix(NA, ncol = 2, nrow = m*n)
+ for(k in nr){
+ ind <- which(NROW == k)
+ temp <- matrix(do.call(rbind, res[ind]), nrow = k)
+ ind1 <- as.vector(sapply(seq_len(n)-1, function(x, ind, m){ ind + x*m }, ind = ind, m = m))
+ rob.est1[ind1, 1:2] <- robloxEM(t(temp), eps = eps, eps.lower = eps.lower, k = steps,
+ mad0 = mad0)
+ }
+ sb <- matrix(rob.est1[,1], nrow = m)
+ for(k in seq_len(m)){
+ IDX <- INDEX[[k]]
+ l.pm <- IDX[,1]
+ if(ncol(IDX) == 2){
+ l.mm <- IDX[,2]
+ }else{
+ l.mm <- integer()
+ }
+ pps.pm <- intensData[l.pm, , drop = FALSE]
+ pps.mm <- intensData[l.mm, , drop = FALSE]
+ pps.im <- pps.mm
+ l <- t(t(pps.mm >= pps.pm) & (sb[k,] > contrast.tau))
+ pps.im[l] <- t(t(pps.pm)/2^sb[k,])[l]
+ l <- t(t(pps.mm >= pps.pm) & (sb[k,] <= contrast.tau))
+ pps.im[l] <- t(t(pps.pm)/2^(contrast.tau/(1 + (contrast.tau - sb[k,])/scale.tau)))[l]
+ pm.corrected <- pmax.int(pps.pm - pps.im, delta)
+ # pm.corrected[pm.corrected < delta] <- delta
+ res[[k]] <- pm.corrected
+ }
+ }else{
+ res <- lapply(INDEX, pm.only, x = intensData)
+ }
+ if(verbose){
+ cat(" done.\n")
+ cat("Computing expression values ...")
+ }
+ for(k in nr){
+ ind <- which(NROW == k)
+ temp <- matrix(do.call(rbind, res[ind]), nrow = k)
+ ind1 <- as.vector(sapply(seq_len(n)-1, function(x, ind, m){ ind + x*m }, ind = ind, m = m))
+ rob.est[ind1, 1:2] <- robloxEM(log2(t(temp)), eps = eps, eps.lower = eps.lower, k = steps, mad0 = mad0)
+ }
+ if(verbose) cat(" done.\n")
+ exp.mat <- 2^matrix(rob.est[,1], nrow = m)
+ se.mat <- 2^matrix(rob.est[,2], nrow = m)
+
+ dimnames(exp.mat) <- list(ids, sampleNames(x))
+ dimnames(se.mat) <- list(ids, sampleNames(x))
+ eset <- new("ExpressionSet", phenoData = phenoData(x),
+ experimentData = experimentData(x), exprs = exp.mat,
+ se.exprs = se.mat, annotation = annotation(x))
+ return(eset)
+ })
Property changes on: pkg/RobLoxBioC/R/robloxAffyBatch.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/RobLoxBioC/R/robloxMatrix.R
===================================================================
--- pkg/RobLoxBioC/R/robloxMatrix.R (rev 0)
+++ pkg/RobLoxBioC/R/robloxMatrix.R 2009-03-02 14:19:35 UTC (rev 257)
@@ -0,0 +1,150 @@
+###############################################################################
+## computation of k-step construction in case x is a matrix
+###############################################################################
+.onestep.locsc.matrix <- function(x, initial.est, A1, A2, a, b){
+ mean <- initial.est[,1]
+ sd <- initial.est[,2]
+ u <- A1*(x-mean)/sd^2
+ v <- A2*(((x-mean)/sd)^2-1)/sd - a
+ ind <- b/sqrt(u^2 + v^2) <= 1
+ IC1 <- rowMeans(u*(ind*b/sqrt(u^2 + v^2) + !ind), na.rm = TRUE)
+ IC2 <- rowMeans(v*(ind*b/sqrt(u^2 + v^2) + !ind), na.rm = TRUE)
+ IC <- cbind(IC1, IC2)
+ return(initial.est + IC)
+}
+.kstep.locsc.matrix <- function(x, initial.est, A1, A2, a, b, mean, k){
+ est <- .onestep.locsc.matrix(x = x, initial.est = initial.est, A1 = A1, A2 = A2, a = a, b = b)
+ if(k > 1){
+ for(i in 2:k){
+ A1 <- est[,2]^2*A1/initial.est[,2]^2
+ A2 <- est[,2]^2*A2/initial.est[,2]^2
+ a <- est[,2]*a/initial.est[,2]
+ b <- est[,2]*b/initial.est[,2]
+ initial.est <- est
+ est <- .onestep.locsc.matrix(x = x, initial.est = est, A1 = A1, A2 = A2, a = a, b = b)
+ }
+ }
+ A1 <- est[,2]^2*A1/initial.est[,2]^2
+ A2 <- est[,2]^2*A2/initial.est[,2]^2
+ a <- est[,2]*a/initial.est[,2]
+ b <- est[,2]*b/initial.est[,2]
+
+ return(est)
+}
+
+
+###############################################################################
+## Use robloxbioc to compute optimally robust (rmx) estimator for rows of a
+## matrix
+###############################################################################
+setMethod("robloxbioc", signature(x = "matrix"),
+ function(x, eps, eps.lower = 0, eps.upper = 0.1, steps = 1L, mad0 = 1e-4){
+ if(missing(x))
+ stop("'x' is missing with no default")
+ if(is.data.frame(x))
+ x <- data.matrix(x)
+ else
+ x <- as.matrix(x)
+ if(!is.matrix(x))
+ stop("'x' has to be a matrix resp. convertable to a matrix by 'as.matrix'
+ or 'data.matrix'")
+
+ if(missing(eps) && missing(eps.lower) && missing(eps.upper)){
+ eps.lower <- 0
+ eps.upper <- 0.5
+ }
+ if(missing(eps)){
+ if(!missing(eps.lower) && missing(eps.upper))
+ eps.upper <- 0.5
+ if(missing(eps.lower) && !missing(eps.upper))
+ eps.lower <- 0
+ if(length(eps.lower) != 1 || length(eps.upper) != 1)
+ stop("'eps.lower' and 'eps.upper' have to be of length 1")
+ if(!is.numeric(eps.lower) || !is.numeric(eps.upper) || eps.lower >= eps.upper)
+ stop("'eps.lower' < 'eps.upper' is not fulfilled")
+ if((eps.lower < 0) || (eps.upper > 0.5))
+ stop("'eps.lower' and 'eps.upper' have to be in [0, 0.5]")
+ }else{
+ if(length(eps) != 1)
+ stop("'eps' has to be of length 1")
+ if(eps == 0)
+ stop("'eps = 0'! => use functions 'mean' and 'sd' for estimation")
+ if((eps < 0) || (eps > 0.5))
+ stop("'eps' has to be in (0, 0.5]")
+ }
+ if(length(steps) != 1){
+ warning("'steps' has to be of length 1 => only first element is used!")
+ steps <- steps[1]
+ }
+ if(steps < 1){
+ stop("'steps' has to be some positive integer value")
+ }
+ steps <- as.integer(steps)
+
+ mean <- rowMedians(x, na.rm = TRUE)
+ sd <- rowMedians(abs(x-mean), na.rm = TRUE)/qnorm(0.75)
+ if(any(sd == 0)){
+ warning("Some of the initial scale estimates were 0 => set to 'mad0'")
+ sd[sd == 0] <- mad0
+ }
+
+ if(!missing(eps)){
+ r <- sqrt(ncol(x))*eps
+ if(r > 10){
+ b <- sd*1.618128043
+ const <- 1.263094656
+ A2 <- b^2*(1+r^2)/(1+const)
+ A1 <- const*A2
+ a <- -0.6277527697*A2/sd
+ mse <- A1 + A2
+ }else{
+ A1 <- sd^2*.getA1.locsc(r)
+ A2 <- sd^2*.getA2.locsc(r)
+ a <- sd*.geta.locsc(r)
+ b <- sd*.getb.locsc(r)
+ mse <- A1 + A2
+ }
+ robEst <- .kstep.locsc.matrix(x = x, initial.est = cbind(mean, sd),
+ A1 = A1, A2 = A2, a = a, b = b, k = steps)
+ colnames(robEst) <- c("mean", "sd")
+ }else{
+ sqrtn <- sqrt(ncol(x))
+ rlo <- sqrtn*eps.lower
+ rup <- sqrtn*eps.upper
+ if(rlo > 10){
+ r <- (rlo + rup)/2
+ }else{
+ r <- uniroot(.getlsInterval, lower = rlo+1e-8, upper = rup,
+ tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
+ }
+ if(r > 10){
+ b <- sd*1.618128043
+ const <- 1.263094656
+ A2 <- b^2*(1+r^2)/(1+const)
+ A1 <- const*A2
+ a <- -0.6277527697*A2/sd
+ mse <- A1 + A2
+ }else{
+ A1 <- sd^2*.getA1.locsc(r)
+ A2 <- sd^2*.getA2.locsc(r)
+ a <- sd*.geta.locsc(r)
+ b <- sd*.getb.locsc(r)
+ mse <- A1 + A2
+ }
+ if(rlo == 0){
+ ineff <- (A1 + A2 - b^2*r^2)/(1.5*sd^2)
+ }else{
+ if(rlo > 10){
+ ineff <- 1
+ }else{
+ A1lo <- sd^2*.getA1.locsc(rlo)
+ A2lo <- sd^2*.getA2.locsc(rlo)
+ ineff <- (A1 + A2 - b^2*(r^2 - rlo^2))/(A1lo + A2lo)
+ }
+ }
+ robEst <- .kstep.locsc.matrix(x = x, initial.est = cbind(mean, sd),
+ A1 = A1, A2 = A2, a = a, b = b, k = steps)
+ colnames(robEst) <- c("mean", "sd")
+ }
+ return(robEst)
+ })
Added: pkg/RobLoxBioC/R/sysdata.rda
===================================================================
(Binary files differ)
Property changes on: pkg/RobLoxBioC/R/sysdata.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/RobLoxBioC/inst/CITATION
===================================================================
--- pkg/RobLoxBioC/inst/CITATION (rev 0)
+++ pkg/RobLoxBioC/inst/CITATION 2009-03-02 14:19:35 UTC (rev 257)
@@ -0,0 +1,19 @@
+if(!exists("meta") || is.null(meta)) meta <- packageDescription("RobLoxBioC")
+year <- sub("-.*", "", meta$Date)
+note <- sprintf("R package version %s", meta$Version)
+
+citHeader("To cite package RobLoxBioC in publications use:")
+
+citEntry(entry="Manual",
+ title = "RobLoxBioC: Infinitesimally robust estimators for preprocessing omics data",
+ author = personList(as.person("M. Kohl")),
+ language = "English",
+ year = year,
+ note = note,
+ type = "R package",
+ url = "http://robast.r-forge.r-project.org/",
+ textVersion = paste("Kohl, M.",
+ sprintf("(%s).", year),
+ "RobLoxBioC: Infinitesimally robust estimators for preprocessing omics data.",
+ paste(note, ".", sep = ""),
+ "URL http://robast.r-forge.r-project.org/"))
Property changes on: pkg/RobLoxBioC/inst/CITATION
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/RobLoxBioC/inst/NEWS
===================================================================
--- pkg/RobLoxBioC/inst/NEWS (rev 0)
+++ pkg/RobLoxBioC/inst/NEWS 2009-03-02 14:19:35 UTC (rev 257)
@@ -0,0 +1,8 @@
+###############################################################################
+## News: to package RobLoxBioC
+###############################################################################
+
+#######################################
+version 0.1
+#######################################
++ start of development
Property changes on: pkg/RobLoxBioC/inst/NEWS
___________________________________________________________________
Name: svn:executable
+ *
More information about the Robast-commits
mailing list