[Robast-commits] r272 - in pkg/RobLoxBioC: . R inst inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 8 09:05:51 CEST 2009
Author: stamats
Date: 2009-04-08 09:05:51 +0200 (Wed, 08 Apr 2009)
New Revision: 272
Added:
pkg/RobLoxBioC/R/KolmogorovMinDist.R
pkg/RobLoxBioC/R/robloxbiocAffyBatch.R
pkg/RobLoxBioC/R/robloxbiocMatrix.R
pkg/RobLoxBioC/man/KolmogorovMinDist.Rd
Removed:
pkg/RobLoxBioC/R/robloxAffyBatch.R
pkg/RobLoxBioC/R/robloxMatrix.R
Modified:
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/NAMESPACE
pkg/RobLoxBioC/R/AllGeneric.R
pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
pkg/RobLoxBioC/inst/NEWS
pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
pkg/RobLoxBioC/inst/scripts/IlluminaExample.R
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
pkg/RobLoxBioC/man/robloxbioc.Rd
Log:
Version 0.3 of the package including methods for Kolmogorov minimum distance.
Modified: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION 2009-03-31 16:44:19 UTC (rev 271)
+++ pkg/RobLoxBioC/DESCRIPTION 2009-04-08 07:05:51 UTC (rev 272)
@@ -1,11 +1,11 @@
Package: RobLoxBioC
-Version: 0.2
-Date: 2009-03-14
+Version: 0.3
+Date: 2009-04-08
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, beadarray
+Depends: R(>= 2.8.1), methods, Biobase, affy, beadarray, distr
Author: Matthias Kohl <Matthias.Kohl at stamats.de>
Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
LazyLoad: yes
Modified: pkg/RobLoxBioC/NAMESPACE
===================================================================
--- pkg/RobLoxBioC/NAMESPACE 2009-03-31 16:44:19 UTC (rev 271)
+++ pkg/RobLoxBioC/NAMESPACE 2009-04-08 07:05:51 UTC (rev 272)
@@ -1,3 +1,3 @@
import("methods")
-exportMethods("robloxbioc")
+exportMethods("robloxbioc", "KolmogorovMinDist")
Modified: pkg/RobLoxBioC/R/AllGeneric.R
===================================================================
--- pkg/RobLoxBioC/R/AllGeneric.R 2009-03-31 16:44:19 UTC (rev 271)
+++ pkg/RobLoxBioC/R/AllGeneric.R 2009-04-08 07:05:51 UTC (rev 272)
@@ -7,3 +7,8 @@
setGeneric("robloxbioc",
function(x, ...) standardGeneric("robloxbioc"))
}
+
+if(!isGeneric("KolmogorovMinDist")){
+ setGeneric("KolmogorovMinDist",
+ function(x, D, ...) standardGeneric("KolmogorovMinDist"))
+}
Added: pkg/RobLoxBioC/R/KolmogorovMinDist.R
===================================================================
--- pkg/RobLoxBioC/R/KolmogorovMinDist.R (rev 0)
+++ pkg/RobLoxBioC/R/KolmogorovMinDist.R 2009-04-08 07:05:51 UTC (rev 272)
@@ -0,0 +1,235 @@
+###############################################################################
+## Methods for KolmogorovMinDist: computation of minimum Kolmogorov distance
+###############################################################################
+setMethod("KolmogorovMinDist", signature(x = "matrix",
+ D = "Norm"),
+ function(x, D, mad0 = 1e-4){
+ stopifnot(is.numeric(x))
+ ksdist <- function(pars, x){
+ if(any(ina <- is.na(x))) x <- x[!ina]
+ y <- .Internal(sort(x, FALSE))
+ p.y <- pnorm(y, mean = pars[1], sd = pars[2])
+ n <- length(y)
+ max(p.y - (seq_len(n)-1)/n, seq_len(n)/n - p.y)
+ }
+ Med <- rowMedians(x, na.rm = TRUE)
+ Mad <- pmax(rowMedians(abs(x-Med), na.rm = TRUE)/qnorm(0.75), mad0)
+ startPars <- cbind(Med, Mad)
+ n <- nrow(x)
+ res <- numeric(n)
+ for(i in seq_len(n)){
+ temp <- try(optim(par = startPars[i,], fn = ksdist, x = x[i,],
+ method = "L-BFGS-B", lower = c(-Inf, 1e-15),
+ upper = c(Inf, Inf))$value)
+ if(inherits(temp, "try-error")){
+ res[i] <- NA
+ }else{
+ res[i] <- temp
+ }
+ }
+ res
+ })
+
+setMethod("KolmogorovMinDist", signature(x = "AffyBatch",
+ D = "AbscontDistribution"),
+ function(x, D, bg.correct = TRUE, pmcorrect = TRUE, verbose = TRUE){
+ if(bg.correct){
+ if(verbose) cat("Background correcting ...")
+ x <- affy::bg.correct.mas(x, griddim = 16)
+ if(verbose) cat(" done.\n")
+ }
+ 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)))
+
+ intensData <- intensity(x)
+ kmd <- numeric(m*n)
+ if(pmcorrect){
+ if(verbose) cat("Calculating PM/MM ...")
+ div <- function(INDEX, x){
+ l.pm <- INDEX[,1]
+ if(ncol(INDEX) == 2)
+ l.mm <- INDEX[,2]
+ else
+ l.mm <- integer()
+ x[l.pm, , drop = FALSE]/x[l.mm, , drop = FALSE]
+ }
+ res <- lapply(INDEX, div, x = intensData)
+ if(verbose) cat(" done.\n")
+ }else{
+ if(verbose) cat("Extract PM data ...")
+ pm.only <- function(INDEX, x){
+ l.pm <- INDEX[,1]
+ x[l.pm, , drop = FALSE]
+ }
+ res <- lapply(INDEX, pm.only, x = intensData)
+ if(verbose) cat(" done.\n")
+ }
+ if(verbose) cat("Computing minimum Kolmogorov distance ...")
+ 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))
+ kmd[ind1] <- KolmogorovMinDist(log2(t(temp)), D = D)
+ }
+ if(verbose) cat(" done.\n")
+ kmd.mat <- matrix(kmd, nrow = m)
+ dimnames(kmd.mat) <- list(ids, sampleNames(x))
+ kmd.mat
+ })
+
+setMethod("KolmogorovMinDist", signature(x = "BeadLevelList",
+ D = "AbscontDistribution"),
+ function(x, D, log = FALSE, imagesPerArray = 1, what = "G", probes = NULL, arrays = NULL){
+ BLData <- x
+ arraynms <- arrayNames(BLData)
+ if(!is.null(arrays) && !is.character(arrays)) arraynms <- arraynms[arrays]
+ if(is.character(arrays)) arraynms <- which(arraynms %in% arrays)
+ len <- length(arraynms)
+ what <- match.arg(what, c("G", "R", "RG", "M", "A", "beta"))
+ whatelse <- ""
+ if(what == "RG"){
+ if(BLData at arrayInfo$channels == "two"){
+ what <- "G"
+ whatelse <- "R"
+ }else{
+ stop("Need two-channel data to calculate summary R and G values")
+ }
+ }
+ if(imagesPerArray == 1){
+ sel <- getArrayData(BLData, what = "ProbeID", array = arraynms[1]) != 0
+ pr <- getArrayData(BLData, what = "ProbeID", array = arraynms[1])[sel]
+ finten <- getArrayData(BLData, what = what, log = log, array = arraynms[1])[sel]
+ nasinf <- !is.finite(finten) | is.na(finten)
+ finten <- finten[!nasinf]
+ }
+ else if(imagesPerArray == 2){
+ if(length(arraynms)%%2 != 0)
+ stop("Need an even number of arrays when 'imagesPerArray=2'")
+ arrayord <- order(arraynms)
+ arraynms <- arraynms[arrayord]
+ tmp <- unlist(strsplit(arraynms, "_"))
+ chipnums <- tmp[seq(1, length(tmp), by = 3)]
+ pos <- tmp[seq(2, length(tmp), by = 3)]
+ stripnum <- as.numeric(tmp[seq(3, length(tmp), by = 3)])
+ check <- ((chipnums[seq(1, len, by = 2)] == chipnums[seq(2, len, by = 2)])
+ & (pos[seq(1, len, by = 2)] == pos[seq(2, len, by = 2)])
+ & (stripnum[seq(1, len, by = 2)] == stripnum[seq(2, len, by = 2)] - 1))
+ if(sum(check) != length(check)) stop("Missing arrays")
+ sel1 <- getArrayData(BLData, what = "ProbeID", array = arraynms[1]) != 0
+ sel2 <- getArrayData(BLData, what = "ProbeID", array = arraynms[2]) != 0
+ pr <- append(getArrayData(BLData, what = "ProbeID", array = arraynms[1])[sel1],
+ getArrayData(BLData, what = "ProbeID", array = arraynms[2])[sel2])
+ finten <- append(getArrayData(BLData, what = what, log = log, array = arraynms[1])[sel1],
+ getArrayData(BLData, what = what, log = log, array = arraynms[2])[sel2])
+ nasinf <- !is.finite(finten) | is.na(finten)
+ finten <- finten[!nasinf]
+ }else{
+ stop("You can only specify 1 or 2 images per array")
+ }
+ if(is.null(probes)) probes <- sort(unique(pr))
+ probes <- probes[probes > 0 & !is.na(probes)]
+ noprobes <- length(probes)
+ pr <- pr[!nasinf]
+ if (imagesPerArray == 1) {
+ G.kmd <- matrix(0, nrow = noprobes, ncol = len)
+ colnames(G.kmd) <- arraynms
+ if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
+ R.kmd <- G.kmd
+ else R.kmd <- NULL
+ }
+ else if (imagesPerArray == 2) {
+ G.kmd <- matrix(0, nrow = noprobes, ncol = (len/2))
+ colnames(G.kmd) <- arraynms[seq(1, len, by = 2)]
+ if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
+ R.kmd <- G.kmd
+ else R.kmd <- NULL
+ }
+ i <- j <- 1
+ while (j <= len) {
+ probeIDs <- as.integer(pr)
+ start <- 0
+ G.kmd[,i] <- kmdBeadLevel(x = finten, D = D, probeIDs = probeIDs, probes = probes)
+ if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[i]]]$R) && whatelse == "R") {
+ if (imagesPerArray == 1) {
+ finten <- getArrayData(BLData, what = whatelse, log = log, array = arraynms[i])[sel]
+ nasinf <- !is.finite(finten) | is.na(finten)
+ finten <- finten[!nasinf]
+ binten <- rep(0, length(finten))
+ }
+ else if (imagesPerArray == 2) {
+ finten <- append(getArrayData(BLData, what = whatelse, log = log, array = arraynms[j])[sel1],
+ getArrayData(BLData, what = whatelse, log = log, array = arraynms[j + 1])[sel2])
+ nasinf <- !is.finite(finten) | is.na(finten)
+ finten <- finten[!nasinf]
+ binten <- rep(0, length(finten))
+ }
+ start <- 0
+ R.kmd[,i] <- kmdBeadLevel(x = finten, D = D, probeIDs = probeIDs, probes = probes)
+ }
+ j <- j + imagesPerArray
+ i <- i + 1
+ rm(probeIDs)
+ gc()
+ if ((imagesPerArray == 1) && (i <= len)) {
+ sel <- getArrayData(BLData, what = "ProbeID", array = arraynms[i]) != 0
+ pr <- getArrayData(BLData, what = "ProbeID", array = arraynms[i])[sel]
+ finten <- getArrayData(BLData, what = what, log = log, array = arraynms[i])[sel]
+ nasinf <- !is.finite(finten) | is.na(finten)
+ pr <- pr[!nasinf]
+ finten <- finten[!nasinf]
+ }
+ else if ((imagesPerArray == 2) && (j < len)) {
+ sel1 <- getArrayData(BLData, what = "ProbeID", array = arraynms[j]) != 0
+ sel2 <- getArrayData(BLData, what = "ProbeID", array = arraynms[j + 1]) != 0
+ pr <- append(getArrayData(BLData, what = "ProbeID", array = arraynms[j])[sel1],
+ getArrayData(BLData, what = "ProbeID", array = arraynms[j + 1])[sel2])
+ finten <- append(getArrayData(BLData, what = what, log = log, array = arraynms[j])[sel1],
+ getArrayData(BLData, what = what, log = log, array = arraynms[j + 1])[sel2])
+ nasinf <- !is.finite(finten) | is.na(finten)
+ pr <- pr[!nasinf]
+ finten <- finten[!nasinf]
+ }
+ }
+ if (whatelse == "R") {
+ rownames(G.kmd) <- rownames(R.kmd) <- probes
+ res <- list(G = G.kmd, R = R.kmd)
+ }
+ else {
+ rownames(G.kmd) <- probes
+ res <- G.kmd
+ }
+ return(res)
+ })
+kmdBeadLevel <- function(x, D, probeIDs, probes){
+ comIDs <- intersect(probeIDs, probes)
+ x <- x[probeIDs %in% comIDs]
+ probeIDs <- probeIDs[probeIDs %in% comIDs]
+ noBeads <- as.vector(table(probeIDs))
+ noBeads.uni <- as.integer(names(table(noBeads)))
+ probes1 <- comIDs
+ len1 <- length(probes1)
+ kmd1 <- numeric(len1)
+ for(i in seq(along = noBeads.uni)){
+ index <- noBeads == noBeads.uni[i]
+ IDs <- probes1[index]
+ if(noBeads.uni[i] == 1){
+ kmd1[index] <- 0.5
+ }else{
+ temp <- matrix(x[probeIDs %in% IDs], ncol = noBeads.uni[i], byrow = TRUE)
+ kmd1[index] <- KolmogorovMinDist(temp, D = D)
+ }
+ }
+ len <- length(probes)
+ kmd <- numeric(len)
+ nas <- !(probes %in% comIDs)
+ kmd[nas] <- NA
+ kmd[!nas] <- kmd1
+
+ return(kmd)
+}
Deleted: pkg/RobLoxBioC/R/robloxAffyBatch.R
===================================================================
--- pkg/RobLoxBioC/R/robloxAffyBatch.R 2009-03-31 16:44:19 UTC (rev 271)
+++ pkg/RobLoxBioC/R/robloxAffyBatch.R 2009-04-08 07:05:51 UTC (rev 272)
@@ -1,104 +0,0 @@
-###############################################################################
-## Use robloxbioc to preprocess Affymetrix-data - comparable to MAS 5.0
-###############################################################################
-setMethod("robloxbioc", signature(x = "AffyBatch"),
- function(x, bg.correct = TRUE, pmcorrect = TRUE, normalize = FALSE,
- add.constant = 32, verbose = TRUE,
- eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 3L,
- mad0 = 1e-4, contrast.tau = 0.03, scale.tau = 10, delta = 2^(-20), sc = 500) {
- if(bg.correct){
- if(verbose) cat("Background correcting ...")
- x <- affy::bg.correct.mas(x, griddim = 16)
- if(verbose) cat(" done.\n")
- }
- 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)))
-
- intensData <- intensity(x)
- rob.est <- matrix(NA, ncol = 2, nrow = m*n)
- if(pmcorrect){
- if(verbose) cat("PM/MM correcting ...")
- 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])
- }
- 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] <- robloxbioc(t(temp), eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
- steps = 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 <- matrix(pmax.int(pps.pm - pps.im, delta), ncol = ncol(pps.pm))
- colnames(pm.corrected) <- colnames(pps.pm)
- rownames(pm.corrected) <- rownames(pps.pm)
- res[[k]] <- pm.corrected
- }
- if(verbose) cat(" done.\n")
- }else{
- if(verbose) cat("Extract PM data ...")
- pm.only <- function(INDEX, x){
- l.pm <- INDEX[,1]
- x[l.pm, , drop = FALSE]
- }
- res <- lapply(INDEX, pm.only, x = intensData)
- if(verbose) cat(" done.\n")
- }
- if(verbose) 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] <- robloxbioc(log2(t(temp)), eps = eps, eps.lower = eps.lower,
- eps.upper = eps.upper, steps = steps, mad0 = mad0)
- rob.est[ind1, 2]/sqrt(k)
- }
- if(verbose) cat(" done.\n")
- exp.mat <- 2^matrix(rob.est[,1], nrow = m) + add.constant
- se.mat <- 2^matrix(rob.est[,2], nrow = m) + add.constant
-
- 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))
- if(normalize){
- if(verbose) cat("Scale normalization ...")
- for (i in 1:ncol(exprs(eset))) {
- slg <- exprs(eset)[, i]
- sf <- sc/mean(slg, trim = 0.02)
- reported.value <- sf * slg
- exprs(eset)[, i] <- reported.value
- }
- if(verbose) cat(" done.\n")
- }
- return(eset)
- })
Deleted: pkg/RobLoxBioC/R/robloxMatrix.R
===================================================================
--- pkg/RobLoxBioC/R/robloxMatrix.R 2009-03-31 16:44:19 UTC (rev 271)
+++ pkg/RobLoxBioC/R/robloxMatrix.R 2009-04-08 07:05:51 UTC (rev 272)
@@ -1,135 +0,0 @@
-###############################################################################
-## Use robloxbioc to compute optimally robust (rmx) estimator for rows of a
-## matrix
-###############################################################################
-setMethod("robloxbioc", signature(x = "matrix"),
- function(x, eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 3L, mad0 = 1e-4){
- if(is.null(eps)){
- 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){
- warning("'eps' has to be of length 1 => only first element is used")
- eps <- eps[1]
- }
- if(!is.numeric(eps))
- stop("'eps' has to be a double in (0, 0.5]")
- 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(!is.null(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)
- })
-
-###############################################################################
-## 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)
-}
Copied: pkg/RobLoxBioC/R/robloxbiocAffyBatch.R (from rev 271, pkg/RobLoxBioC/R/robloxAffyBatch.R)
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocAffyBatch.R (rev 0)
+++ pkg/RobLoxBioC/R/robloxbiocAffyBatch.R 2009-04-08 07:05:51 UTC (rev 272)
@@ -0,0 +1,104 @@
+###############################################################################
+## Use robloxbioc to preprocess Affymetrix-data - comparable to MAS 5.0
+###############################################################################
+setMethod("robloxbioc", signature(x = "AffyBatch"),
+ function(x, bg.correct = TRUE, pmcorrect = TRUE, normalize = FALSE,
+ add.constant = 32, verbose = TRUE,
+ eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 3L,
+ mad0 = 1e-4, contrast.tau = 0.03, scale.tau = 10, delta = 2^(-20), sc = 500) {
+ if(bg.correct){
+ if(verbose) cat("Background correcting ...")
+ x <- affy::bg.correct.mas(x, griddim = 16)
+ if(verbose) cat(" done.\n")
+ }
+ 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)))
+
+ intensData <- intensity(x)
+ rob.est <- matrix(NA, ncol = 2, nrow = m*n)
+ if(pmcorrect){
+ if(verbose) cat("PM/MM correcting ...")
+ 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]/x[l.mm, , drop = FALSE])
+ }
+ 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] <- robloxbioc(t(temp), eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
+ steps = 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 <- matrix(pmax.int(pps.pm - pps.im, delta), ncol = ncol(pps.pm))
+ colnames(pm.corrected) <- colnames(pps.pm)
+ rownames(pm.corrected) <- rownames(pps.pm)
+ res[[k]] <- pm.corrected
+ }
+ if(verbose) cat(" done.\n")
+ }else{
+ if(verbose) cat("Extract PM data ...")
+ pm.only <- function(INDEX, x){
+ l.pm <- INDEX[,1]
+ x[l.pm, , drop = FALSE]
+ }
+ res <- lapply(INDEX, pm.only, x = intensData)
+ if(verbose) cat(" done.\n")
+ }
+ if(verbose) 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] <- robloxbioc(log2(t(temp)), eps = eps, eps.lower = eps.lower,
+ eps.upper = eps.upper, steps = steps, mad0 = mad0)
+ rob.est[ind1, 2] <- rob.est[ind1, 2]/sqrt(k)
+ }
+ if(verbose) cat(" done.\n")
+ exp.mat <- 2^matrix(rob.est[,1], nrow = m) + add.constant
+ 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))
+ if(normalize){
+ if(verbose) cat("Scale normalization ...")
+ for (i in 1:ncol(exprs(eset))) {
+ slg <- exprs(eset)[, i]
+ sf <- sc/mean(slg, trim = 0.02)
+ reported.value <- sf * slg
+ exprs(eset)[, i] <- reported.value
+ }
+ if(verbose) cat(" done.\n")
+ }
+ return(eset)
+ })
Property changes on: pkg/RobLoxBioC/R/robloxbiocAffyBatch.R
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
Modified: pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2009-03-31 16:44:19 UTC (rev 271)
+++ pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2009-04-08 07:05:51 UTC (rev 272)
@@ -17,13 +17,12 @@
}
}
if(imagesPerArray == 1){
- sel <- getArrayData(BLData, what = "ProbeID", array = arraynms[1]) != 0
- pr <- getArrayData(BLData, what = "ProbeID", array = arraynms[1])[sel]
+ pr <- getArrayData(BLData, what = "ProbeID", array = arraynms[1])
+ sel <- pr != 0
+ pr <- pr[sel]
finten <- getArrayData(BLData, what = what, log = log, array = arraynms[1])[sel]
nasinf <- !is.finite(finten) | is.na(finten)
- pr <- pr[!nasinf]
finten <- finten[!nasinf]
- binten <- rep(0, length(finten))
}
else if(imagesPerArray == 2){
if(length(arraynms)%%2 != 0)
@@ -42,22 +41,18 @@
sel1 <- getArrayData(BLData, what = "ProbeID", array = arraynms[1]) != 0
sel2 <- getArrayData(BLData, what = "ProbeID", array = arraynms[2]) != 0
pr <- append(getArrayData(BLData, what = "ProbeID", array = arraynms[1])[sel1],
- getArrayData(BLData, what = "ProbeID", array = arraynms[2])[sel2])
+ getArrayData(BLData, what = "ProbeID", array = arraynms[2])[sel2])
finten <- append(getArrayData(BLData, what = what, log = log, array = arraynms[1])[sel1],
- getArrayData(BLData, what = what, log = log, array = arraynms[2])[sel2])
+ getArrayData(BLData, what = what, log = log, array = arraynms[2])[sel2])
nasinf <- !is.finite(finten) | is.na(finten)
- pr <- pr[!nasinf]
finten <- finten[!nasinf]
- binten <- rep(0, length(finten))
- ord <- order(pr)
- pr <- pr[ord]
- finten <- finten[ord]
}else{
stop("You can only specify 1 or 2 images per array")
}
if(is.null(probes)) probes <- sort(unique(pr))
probes <- probes[probes > 0 & !is.na(probes)]
noprobes <- length(probes)
+ pr <- pr[!nasinf]
if (imagesPerArray == 1) {
G <- GBeadStDev <- GNoBeads <- matrix(0, nrow = noprobes, ncol = len)
colnames(G) <- colnames(GBeadStDev) <- colnames(GNoBeads) <- arraynms
@@ -87,19 +82,13 @@
finten <- getArrayData(BLData, what = whatelse, log = log, array = arraynms[i])[sel]
nasinf <- !is.finite(finten) | is.na(finten)
finten <- finten[!nasinf]
- binten <- rep(0, length(finten))
}
else if (imagesPerArray == 2) {
finten <- append(getArrayData(BLData, what = whatelse, log = log, array = arraynms[j])[sel1],
getArrayData(BLData, what = whatelse, log = log, array = arraynms[j + 1])[sel2])
nasinf <- !is.finite(finten) | is.na(finten)
finten <- finten[!nasinf]
- binten <- rep(0, length(finten))
- ord <- order(pr)
- pr <- pr[ord]
- finten <- finten[ord]
}
- start <- 0
blah <- rmxBeadSummary(x = finten, probeIDs = probeIDs, probes = probes,
eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
steps = steps, mad0 = mad0)
@@ -118,7 +107,6 @@
nasinf <- !is.finite(finten) | is.na(finten)
pr <- pr[!nasinf]
finten <- finten[!nasinf]
- binten <- rep(0, length(finten))
}
else if ((imagesPerArray == 2) && (j < len)) {
sel1 <- getArrayData(BLData, what = "ProbeID", array = arraynms[j]) != 0
@@ -130,10 +118,6 @@
nasinf <- !is.finite(finten) | is.na(finten)
pr <- pr[!nasinf]
finten <- finten[!nasinf]
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 272
More information about the Robast-commits
mailing list