[Robast-commits] r267 - in pkg/RobLoxBioC: . R inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 14 12:28:09 CET 2009
Author: stamats
Date: 2009-03-14 12:28:08 +0100 (Sat, 14 Mar 2009)
New Revision: 267
Added:
pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
Modified:
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
pkg/RobLoxBioC/man/robloxbioc.Rd
Log:
added method for BeadLevelList, package now depends on package beadarray ...
Modified: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION 2009-03-14 11:22:16 UTC (rev 266)
+++ pkg/RobLoxBioC/DESCRIPTION 2009-03-14 11:28:08 UTC (rev 267)
@@ -1,11 +1,11 @@
Package: RobLoxBioC
-Version: 0.1
-Date: 2009-03-12
+Version: 0.2
+Date: 2009-03-14
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
+Depends: R(>= 2.8.1), methods, Biobase, affy, beadarray
Author: Matthias Kohl <Matthias.Kohl at stamats.de>
Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
LazyLoad: yes
Added: pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R (rev 0)
+++ pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2009-03-14 11:28:08 UTC (rev 267)
@@ -0,0 +1,191 @@
+setMethod("robloxbioc", signature(x = "BeadLevelList"),
+ function(x, log = FALSE, imagesPerArray = 1, what = "G", probes = NULL, arrays = NULL,
+ eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 3L, mad0 = 1e-4){
+ 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)
+ pr <- pr[!nasinf]
+ finten <- finten[!nasinf]
+ binten <- rep(0, length(finten))
+ }
+ 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)
+ 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)
+ if (imagesPerArray == 1) {
+ G <- GBeadStDev <- GNoBeads <- matrix(0, nrow = noprobes, ncol = len)
+ colnames(G) <- colnames(GBeadStDev) <- colnames(GNoBeads) <- arraynms
+ if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
+ R <- RBeadStDev <- RNoBeads <- G
+ else R <- NULL
+ }
+ else if (imagesPerArray == 2) {
+ G <- GBeadStDev <- GNoBeads <- matrix(0, nrow = noprobes, ncol = (len/2))
+ colnames(G) <- colnames(GBeadStDev) <- colnames(GNoBeads) <- arraynms[seq(1, len, by = 2)]
+ if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
+ R <- RBeadStDev <- RNoBeads <- G
+ else R <- NULL
+ }
+ i <- j <- 1
+ while (j <= len) {
+ probeIDs <- as.integer(pr)
+ start <- 0
+ blah <- rmxBeadSummary(x = finten, probeIDs = probeIDs, probes = probes,
+ eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
+ steps = steps, mad0 = mad0)
+ G[, i] <- blah$fore
+ GBeadStDev[, i] <- blah$sd
+ GNoBeads[, i] <- blah$noBeads
+ 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))
+ 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)
+ R[, i] <- blah$fore
+ RBeadStDev[, i] <- blah$sd
+ RNoBeads[, i] <- blah$noBeads
+ }
+ j <- j + imagesPerArray
+ i <- i + 1
+ rm(probeIDs, blah)
+ 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]
+ binten <- rep(0, length(finten))
+ }
+ 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]
+ binten <- rep(0, length(finten))
+ ord <- order(pr)
+ pr <- pr[ord]
+ finten <- finten[ord]
+ }
+ }
+ GBeadStDev <- GBeadStDev/sqrt(GNoBeads)
+ if(!is.null(R)) RBeadStDev <- RBeadStDev/sqrt(RNoBeads)
+ if (whatelse == "R") {
+ rownames(G) <- rownames(R) <- rownames(GBeadStDev) <- rownames(RBeadStDev) <- rownames(GNoBeads) <- rownames(RNoBeads) <- probes
+ BSData <- new("NChannelSet", R = R, G = G, GBeadStDev = GBeadStDev,
+ RBeadStDev = RBeadStDev, GNoBeads = GNoBeads, RNoBeads = RNoBeads)
+ }
+ else {
+ BSData <- new("ExpressionSetIllumina")
+ assayData(BSData) <- assayDataNew(exprs = G, se.exprs = GBeadStDev,
+ NoBeads = GNoBeads, storage.mode = "list")
+ rownames(exprs(BSData)) <- rownames(se.exprs(BSData)) <- rownames(NoBeads(BSData)) <- probes
+ featureData(BSData) <- new("AnnotatedDataFrame", data = data.frame(ProbeID = probes, row.names = probes))
+ }
+ if (nrow(pData(BLData)) == len) {
+ if (imagesPerArray == 2)
+ BSData at phenoData <- new("AnnotatedDataFrame", data = pData(BLData at phenoData)[arrayord, , drop = FALSE][seq(1, len, by = 2), , drop = FALSE])
+ else BSData at phenoData <- BLData at phenoData
+ }
+ else {
+ BSData at phenoData <- new("AnnotatedDataFrame", data = data.frame(sampleName = colnames(G)))
+ }
+ if (!is.null(pData(BSData)$sampleName))
+ sampleNames(BSData) <- pData(BSData)$sampleName
+ else sampleNames(BSData) <- colnames(G)
+ if (whatelse == "R") {
+ varMetadata <- data.frame(labelDescription = colnames(BSData at phenoData@data), channel = "_ALL_")
+ BSData at phenoData <- new("AnnotatedDataFrame", data = data.frame(BSData at phenoData@data), varMetadata = varMetadata)
+ }
+ BSData at annotation <- BLData at annotation
+ if ("qcScores" %in% names(BLData at arrayInfo))
+ t <- try(BSData at BeadLevelQC <- BLData at arrayInfo$qcScores, silent = TRUE)
+ BSData
+ })
+## computation of bead summaries via robloxbioc for "matrix"
+rmxBeadSummary <- function(x, probeIDs, probes, eps, eps.lower, eps.upper, steps, mad0){
+ noBeads <- as.vector(table(probeIDs))
+ noBeads.uni <- as.integer(names(table(noBeads)))
+ len <- length(probes)
+ fore <- numeric(len)
+ SD <- numeric(len)
+ for(i in seq(along = noBeads.uni)){
+ index <- noBeads == noBeads.uni[i]
+ IDs <- probes[index]
+ temp <- matrix(x[probeIDs %in% IDs], ncol = noBeads.uni[i], byrow = TRUE)
+ rmx <- robloxbioc(temp, eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
+ steps = steps, mad0 = mad0)
+ fore[index] <- rmx[,"mean"]
+ SD[index] <- rmx[,"sd"]
+ }
+
+ return(list(fore = fore, sd = SD, noBeads = noBeads))
+}
Modified: pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2009-03-14 11:22:16 UTC (rev 266)
+++ pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2009-03-14 11:28:08 UTC (rev 267)
@@ -109,7 +109,7 @@
library(RobLoxBioC)
###########################################################
## Example 1: Analogous to "classical" MAS 5.0
-## both computations take about 38 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
+## both computations take about 55 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
###########################################################
## hgu95a
system.time(eset.hgu95a <- robloxbioc(spikein.hgu95a, normalize = TRUE, add.constant = 0))
@@ -128,7 +128,7 @@
###########################################################
## Example 2: MAS 5.0 + 32
-## both computations take about 38 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
+## both computations take about 55 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
###########################################################
## hgu95a
system.time(eset.hgu95a32 <- robloxbioc(spikein.hgu95a, normalize = TRUE, add.constant = 32))
Modified: pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
===================================================================
--- pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2009-03-14 11:22:16 UTC (rev 266)
+++ pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2009-03-14 11:28:08 UTC (rev 267)
@@ -12,8 +12,8 @@
\details{
\tabular{ll}{
Package: \tab RobLoxBioC\cr
-Version: \tab 0.1 \cr
-Date: \tab 2009-03-12 \cr
+Version: \tab 0.2 \cr
+Date: \tab 2009-03-14 \cr
Depends: \tab R (>= 2.8.1), methods, methods, Biobase, affy\cr
LazyLoad: \tab yes\cr
License: \tab LGPL-3\cr
Modified: pkg/RobLoxBioC/man/robloxbioc.Rd
===================================================================
--- pkg/RobLoxBioC/man/robloxbioc.Rd 2009-03-14 11:22:16 UTC (rev 266)
+++ pkg/RobLoxBioC/man/robloxbioc.Rd 2009-03-14 11:28:08 UTC (rev 267)
@@ -3,6 +3,7 @@
\alias{robloxbioc-methods}
\alias{robloxbioc,matrix-method}
\alias{robloxbioc,AffyBatch-method}
+\alias{robloxbioc,BeadLevelList-method}
\title{Generic Function for Preprocessing Biological Data}
\description{
@@ -18,6 +19,10 @@
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)
+
+\S4method{robloxbioc}{BeadLevelList}(x, log = FALSE, imagesPerArray = 1, what = "G", probes = NULL,
+ arrays = NULL, eps = NULL, eps.lower = 0, eps.upper = 0.1,
+ steps = 3L, mad0 = 1e-4)
}
\arguments{
\item{x}{ biological data. }
@@ -47,6 +52,16 @@
\item{delta}{ a number denoting the delta parameter; confer the MAS 5.0
PM correction algorithm. }
\item{sc}{ value at which all arrays will be scaled to. }
+ \item{log}{ if \code{TRUE}, then the log2 intensities for each bead-type are summarized}
+ \item{imagesPerArray}{ Specifies how many images (strips) there are per array.
+ Normally 1 for a SAM and 1 or 2 for a BeadChip. The images (strips) from the same array
+ will be combined so that each column in the output represents a sample. }
+ \item{what}{ character string specifying which intensities/values to summarize.
+ See \code{\link[beadarray]{getArrayData}} for a list of possibilities. }
+ \item{probes}{ Specify particular probes to summarize. If left \code{NULL} then all
+ the probes on the first array are used. }
+ \item{arrays}{ integer (scalar or vector) specifying the strips/arrays to summarize.
+ If \code{NULL}, then all strips/arrays are summarized. }
}
\details{
The optimally-robust resp. the radius-minimax (rmx) estimator for normal location
@@ -68,8 +83,12 @@
The algorithm used for Affymetrix data is similar to MAS 5.0 (cf. Affymetrix (2002)).
The main difference is the substitution of the Tukey one-step estimator by our rmx
- k-step (k >= 1) estimator. The optional scale normalization is performed as given
- in Affymetrix (2002).
+ k-step (k >= 1) estimator in the PM/MM correction step. The optional scale normalization
+ is performed as given in Affymetrix (2002).
+
+ In case of Illumina data, the rmx estimator is used to summarize the bead types.
+ The implementation for the most part was taken from function
+ \code{\link[beadarray]{createBeadSummaryData}}.
}
\value{ Return value depends on the class of \code{x}.
In case of \code{"matrix"} a matrix with columns "mean" and "sd" is returned.
@@ -98,7 +117,8 @@
\seealso{\code{\link[RobLox]{roblox}}, \code{\link[RobLox]{rowRoblox}},
\code{\link[affy]{AffyBatch-class}},
\code{\link[affy]{generateExprVal.method.mas}},
- \code{\link[Biobase]{ExpressionSet-class}} }
+ \code{\link[Biobase]{ExpressionSet-class}},
+ \code{\link[beadarray]{createBeadSummaryData}} }
\examples{
## similar to rowRoblox of package RobLox
ind <- rbinom(200, size=1, prob=0.05)
@@ -132,6 +152,10 @@
## Affymetrix scale normalization
eset1 <- robloxbioc(Dilution, normalize = TRUE)
}
+
+## using Illumina-Data
+data(BLData)
+BSData <- robloxbioc(BLData, eps.upper = 0.5)
}
\concept{normal location and scale}
\concept{infinitesimal robustness}
More information about the Robast-commits
mailing list