[Robast-commits] r261 - in pkg/RobLoxBioC: . R inst inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 5 20:09:04 CET 2009
Author: stamats
Date: 2009-03-05 20:09:04 +0100 (Thu, 05 Mar 2009)
New Revision: 261
Added:
pkg/RobLoxBioC/inst/scripts/
pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
Modified:
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/R/robloxAffyBatch.R
pkg/RobLoxBioC/R/robloxMatrix.R
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
pkg/RobLoxBioC/man/robloxbioc.Rd
Log:
code seems to work correctly, affycompTable in AffymetrixExample does not work yet (because of missing dilution data).
Modified: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION 2009-03-04 08:58:34 UTC (rev 260)
+++ pkg/RobLoxBioC/DESCRIPTION 2009-03-05 19:09:04 UTC (rev 261)
@@ -1,6 +1,6 @@
Package: RobLoxBioC
Version: 0.1
-Date: 2009-03-03
+Date: 2009-03-05
Title: Infinitesimally robust estimators for preprocessing omics data
Description: Functions for the determination of optimally robust
influence curves and estimators for preprocessing omics data,
Modified: pkg/RobLoxBioC/R/robloxAffyBatch.R
===================================================================
--- pkg/RobLoxBioC/R/robloxAffyBatch.R 2009-03-04 08:58:34 UTC (rev 260)
+++ pkg/RobLoxBioC/R/robloxAffyBatch.R 2009-03-05 19:09:04 UTC (rev 261)
@@ -2,9 +2,15 @@
## Use robloxbioc to preprocess Affymetrix-data - comparable to MAS 5.0
###############################################################################
setMethod("robloxbioc", signature(x = "AffyBatch"),
- function(x, pmcorrect = "roblox", normalize = FALSE, verbose = TRUE,
- eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 1L, mad0 = 1e-4,
- contrast.tau = 0.03, scale.tau = 10, delta = 2^(-20), sc = 500) {
+ function(x, bg.correct = TRUE, pmcorrect = TRUE, normalize = FALSE,
+ add.constant = 32, verbose = TRUE,
+ eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 1L,
+ 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)
@@ -14,23 +20,18 @@
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"){
+ 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){
@@ -56,16 +57,22 @@
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 <- 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(" done.\n")
- cat("Computing expression values ...")
- }
+ if(verbose) cat("Computing expression values ...")
for(k in nr){
ind <- which(NROW == k)
temp <- matrix(do.call(rbind, res[ind]), nrow = k)
@@ -74,8 +81,8 @@
eps.upper = eps.upper, steps = 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)
+ 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))
Modified: pkg/RobLoxBioC/R/robloxMatrix.R
===================================================================
--- pkg/RobLoxBioC/R/robloxMatrix.R 2009-03-04 08:58:34 UTC (rev 260)
+++ pkg/RobLoxBioC/R/robloxMatrix.R 2009-03-05 19:09:04 UTC (rev 261)
@@ -1,39 +1,4 @@
###############################################################################
-## 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
###############################################################################
@@ -134,3 +99,37 @@
}
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)
+}
Added: pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R (rev 0)
+++ pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2009-03-05 19:09:04 UTC (rev 261)
@@ -0,0 +1,380 @@
+###############################################################################
+## Affymetrix example
+###############################################################################
+
+###############################################################################
+## References:
+## L.M. Cope1 , R.A. Irizarry, H.A. Jaffee, Z. Wu and T.P. Speed (2004):
+## "A benchmark for Affymetrix GeneChip expression measures".
+## Bioinformatics 20(3): 323-331
+##
+## R.A. Irizarry1, Z. Wu and H.A. Jaffee (2006):
+## "Comparison of Affymetrix GeneChip expression measures"
+## Bioinformatics 22(7): 789-794
+###############################################################################
+
+###############################################################################
+## Data
+## Spike-in hgu95a data:
+## http://www.biostat.jhsph.edu/~ririzarr/affycomp/spikein.tgz
+##
+## Spike-in hgu133a data:
+## http://www.biostat.jhsph.edu/~ririzarr/affycomp/hgu133spikein.tgz
+##
+## Dilution data:
+## The links:
+## http://qolotus02.genelogic.com/datasets.nsf/ (no response)
+## and
+## http://www.genelogic.com/media/studies/dilution.cfm (not found)
+## seem not to lead to the data any longer.
+## An email to the support of genelogic is still unanswered ...
+###############################################################################
+
+
+library(affy)
+library(affycomp) ## Version 1.19.4
+
+###################
+## replace with your path to hgu95a data!!!
+PATH <- "./spikein"
+###################
+
+fn <- list.celfiles(path = PATH, full.names=TRUE)
+data(spikein.phenodata)
+(spikein.hgu95a <- read.affybatch(filenames = fn,
+ phenoData = spikein.phenodata))
+
+###################
+## replace with your path to hgu133a data!!!
+PATH <- "./hgu133spikein"
+###################
+
+fn <- list.celfiles(path = PATH, full.names=TRUE)
+fn <- fn[c(seq(1,40,3), seq(2, 41, 3), seq(3, 42, 3))]
+fn <- fn[c(6:14, 1:5, 20:28, 15:19, 34:42, 29:33)]
+data(hgu133a.spikein.phenodata)
+
+## Attention:
+## Order of filenames in fn has to be identical to
+## sampleNames(hgu133a.spikein.phenodata)!!!
+
+(spikein.hgu133a <- read.affybatch(filenames = fn,
+ phenoData = hgu133a.spikein.phenodata))
+
+
+###########################################################
+## assessments for MAS 5.0 and RMA including dilution data from package affycomp
+###########################################################
+data(mas5.assessment)
+data(rma.assessment)
+data(mas5.assessment.133)
+data(rma.assessment.133)
+
+res.rma <- rma(spikein.hgu133a)
+
+library(RobLoxBioC)
+###########################################################
+## Example 1: Analogous to "classical" MAS 5.0
+## computation takes about 38 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
+###########################################################
+## hgu95a
+system.time(eset.hgu95a <- robloxbioc(spikein.hgu95a, normalize = TRUE, add.constant = 0))
+eset.hgu95a.log2 <- eset.hgu95a
+exprs(eset.hgu95a.log2) <- log2(exprs(eset.hgu95a))
+roblox.hgu95a <- assessSpikeIn(eset.hgu95a.log2, method.name = "roblox")
+
+## hgu133a
+system.time(eset.hgu133a <- robloxbioc(spikein.hgu133a, normalize = TRUE, add.constant = 0))
+eset.hgu133a.log2 <- eset.hgu133a
+exprs(eset.hgu133a.log2) <- log2(exprs(eset.hgu133a))
+roblox.hgu133a <- assessSpikeIn(eset.hgu133a.log2, method.name = "roblox")
+
+
+###########################################################
+## Example 2: MAS 5.0 + 32
+## computation takes about 38 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
+###########################################################
+## hgu95a
+system.time(eset.hgu95a32 <- robloxbioc(spikein.hgu95a, normalize = TRUE, add.constant = 32))
+eset.hgu95a.log232 <- eset.hgu95a32
+exprs(eset.hgu95a.log232) <- log2(exprs(eset.hgu95a32))
+roblox.hgu95a32 <- assessSpikeIn(eset.hgu95a.log232, method.name = "roblox + 32")
+
+## hgu133a
+system.time(eset.hgu133a32 <- robloxbioc(spikein.hgu133a, normalize = TRUE, add.constant = 32))
+eset.hgu133a.log232 <- eset.hgu133a32
+exprs(eset.hgu133a.log232) <- log2(exprs(eset.hgu133a32))
+roblox.hgu133a32 <- assessSpikeIn(eset.hgu133a.log232, method.name = "roblox + 32")
+
+
+###########################################################
+## Example 3: MAS 5.0 with PM only
+## computation takes about 19 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
+###########################################################
+## hgu95a
+system.time(eset.hgu95a.pmonly <- robloxbioc(spikein.hgu95a, bg.correct = TRUE,
+ pmcorrect = FALSE, normalize = TRUE,
+ add.constant = 0))
+eset.hgu95a.log2.pmonly <- eset.hgu95a.pmonly
+exprs(eset.hgu95a.log2.pmonly) <- log2(exprs(eset.hgu95a.pmonly))
+roblox.hgu95a.pmonly <- assessSpikeIn(eset.hgu95a.log2.pmonly, method.name = "roblox (PM)")
+
+## hgu133a
+system.time(eset.hgu133a.pmonly <- robloxbioc(spikein.hgu133a, bg.correct = TRUE,
+ pmcorrect = FALSE, normalize = TRUE,
+ add.constant = 0))
+eset.hgu133a.log2.pmonly <- eset.hgu133a.pmonly
+exprs(eset.hgu133a.log2.pmonly) <- log2(exprs(eset.hgu133a.pmonly))
+roblox.hgu133a.pmonly <- assessSpikeIn(eset.hgu133a.log2.pmonly, method.name = "roblox (PM)")
+
+
+###########################################################
+## Example 4: MAS 5.0 + 32 with PM only
+## computation takes about 19 sec on Intel P9500 (64bit Linux, 4 GByte RAM)
+###########################################################
+## hgu95a
+system.time(eset.hgu95a.pmonly32 <- robloxbioc(spikein.hgu95a, bg.correct = TRUE,
+ pmcorrect = FALSE, normalize = TRUE,
+ add.constant = 32))
+eset.hgu95a.log2.pmonly32 <- eset.hgu95a.pmonly32
+exprs(eset.hgu95a.log2.pmonly32) <- log2(exprs(eset.hgu95a.pmonly32))
+roblox.hgu95a.pmonly32 <- assessSpikeIn(eset.hgu95a.log2.pmonly32, method.name = "roblox + 32 (PM)")
+
+## hgu133a
+system.time(eset.hgu133a.pmonly32 <- robloxbioc(spikein.hgu133a, bg.correct = TRUE,
+ pmcorrect = FALSE, normalize = TRUE,
+ add.constant = 32))
+eset.hgu133a.log2.pmonly32 <- eset.hgu133a.pmonly32
+exprs(eset.hgu133a.log2.pmonly32) <- log2(exprs(eset.hgu133a.pmonly32))
+roblox.hgu133a.pmonly32 <- assessSpikeIn(eset.hgu133a.log2.pmonly32, method.name = "roblox + 32 (PM)")
+
+
+###############################################################################
+## Figure 1: The MA plot shows log fold change as a function of mean log
+## expression level. A set of 14 arrays representing a single experiment from
+## the Affymetrix spike-in data are used for this plot. A total of 13 sets of
+## fold changes are generated by comparing the first array in the set to each
+## of the others. Genes are symbolized by numbers representing the nominal
+## log2 fold change for the gene. Non-differentially expressed genes with
+## observed fold changes larger than 2 are plotted in red. All other probesets
+## are represented with black dots.
+###############################################################################
+## hgu95a
+par(mfrow = c(3, 2))
+affycompPlot(roblox.hgu95a$MA)
+affycompPlot(roblox.hgu95a32$MA)
+affycompPlot(roblox.hgu95a.pmonly$MA)
+affycompPlot(roblox.hgu95a.pmonly32$MA)
+affycompPlot(mas5.assessment$MA)
+affycompPlot(rma.assessment$MA)
+
+## hgu133a
+par(mfrow = c(3, 2))
+affycompPlot(roblox.hgu133a$MA)
+affycompPlot(roblox.hgu133a32$MA)
+affycompPlot(roblox.hgu133a.pmonly$MA)
+affycompPlot(roblox.hgu133a.pmonly32$MA)
+affycompPlot(mas5.assessment.133$MA)
+affycompPlot(rma.assessment.133$MA)
+
+
+###############################################################################
+## Figure 4a: Average observed log2 intensity plotted against nominal log2
+## concentration for each spiked-in gene for all arrays in Affymetrix spike-In
+## experiment
+###############################################################################
+## hgu95a
+par(mfrow = c(3, 2))
+affycomp.figure4a(roblox.hgu95a$Signal)
+affycomp.figure4a(roblox.hgu95a32$Signal)
+affycomp.figure4a(roblox.hgu95a.pmonly$Signal)
+affycomp.figure4a(roblox.hgu95a.pmonly32$Signal)
+affycomp.figure4a(mas5.assessment$Signal)
+affycomp.figure4a(rma.assessment$Signal)
+
+## hgu133a
+par(mfrow = c(3, 2))
+affycomp.figure4a(roblox.hgu133a$Signal)
+affycomp.figure4a(roblox.hgu133a32$Signal)
+affycomp.figure4a(roblox.hgu133a.pmonly$Signal)
+affycomp.figure4a(roblox.hgu133a.pmonly32$Signal)
+affycomp.figure4a(mas5.assessment.133$Signal)
+affycomp.figure4a(rma.assessment.133$Signal)
+
+## Comparison plot
+## hgu95a
+affycomp.compfig4a(list(roblox.hgu95a$Signal,
+ roblox.hgu95a32$Signal,
+ roblox.hgu95a.pmonly$Signal,
+ roblox.hgu95a.pmonly32$Signal,
+ mas5.assessment$Signal,
+ rma.assessment$Signal),
+ method.names = c("roblox", "roblox + 32", "roblox (PM)",
+ "roblox + 32 (PM)", "MAS 5.0", "RMA"))
+
+## hgu133a
+affycomp.compfig4a(list(roblox.hgu133a$Signal,
+ roblox.hgu133a32$Signal,
+ roblox.hgu133a.pmonly$Signal,
+ roblox.hgu133a.pmonly32$Signal,
+ mas5.assessment.133$Signal,
+ rma.assessment.133$Signal),
+ method.names = c("roblox", "roblox + 32", "roblox (PM)",
+ "roblox + 32 (PM)", "MAS 5.0", "RMA"))
+
+
+###############################################################################
+## Figure 5: A typical identification rule for differential expression filters
+## genes with fold change exceeding a given threshold. This figure shows
+## average ROC curves which offer a graphical representation of both
+## specificity and sensitivity for such a detection rule.
+## a) Average ROC curves based on comparisons with nominal fold changes
+## ranging from 2 to 4096.
+## b) As a) but with nominal fold changes equal to 2.
+###############################################################################
+## Figure 5a:
+## hgu95a
+par(mfrow = c(3, 2))
+affycomp.figure5a(roblox.hgu95a$FC)
+affycomp.figure5a(roblox.hgu95a32$FC)
+affycomp.figure5a(roblox.hgu95a.pmonly$FC)
+affycomp.figure5a(roblox.hgu95a.pmonly32$FC)
+affycomp.figure5a(mas5.assessment$FC)
+affycomp.figure5a(rma.assessment$FC)
+
+## hgu133a
+par(mfrow = c(3, 2))
+affycomp.figure5a(roblox.hgu133a$FC)
+affycomp.figure5a(roblox.hgu133a32$FC)
+affycomp.figure5a(roblox.hgu133a.pmonly$FC)
+affycomp.figure5a(roblox.hgu133a.pmonly32$FC)
+affycomp.figure5a(mas5.assessment.133$FC)
+affycomp.figure5a(rma.assessment.133$FC)
+
+## Comparison plot
+## hgu95a
+affycomp.compfig5a(list(roblox.hgu95a$FC,
+ roblox.hgu95a32$FC,
+ roblox.hgu95a.pmonly$FC,
+ roblox.hgu95a.pmonly32$FC,
+ mas5.assessment$FC,
+ rma.assessment$FC),
+ method.names = c("roblox", "roblox + 32", "roblox (PM)",
+ "roblox + 32 (PM)", "MAS 5.0", "RMA"))
+
+## hgu133a
+affycomp.compfig5a(list(roblox.hgu133a$FC,
+ roblox.hgu133a32$FC,
+ roblox.hgu133a.pmonly$FC,
+ roblox.hgu133a.pmonly32$FC,
+ mas5.assessment$FC,
+ rma.assessment$FC),
+ method.names = c("roblox", "roblox + 32", "roblox (PM)",
+ "roblox + 32 (PM)", "MAS 5.0", "RMA"))
+
+## Figure 5b:
+## hgu95a
+par(mfrow = c(3, 2))
+affycomp.figure5b(roblox.hgu95a$FC)
+affycomp.figure5b(roblox.hgu95a32$FC)
+affycomp.figure5b(roblox.hgu95a.pmonly$FC)
+affycomp.figure5b(roblox.hgu95a.pmonly32$FC)
+affycomp.figure5b(mas5.assessment$FC)
+affycomp.figure5b(rma.assessment$FC)
+
+## hgu133a
+par(mfrow = c(3, 2))
+affycomp.figure5b(roblox.hgu133a$FC)
+affycomp.figure5b(roblox.hgu133a32$FC)
+affycomp.figure5b(roblox.hgu133a.pmonly$FC)
+affycomp.figure5b(roblox.hgu133a.pmonly32$FC)
+affycomp.figure5b(mas5.assessment.133$FC)
+affycomp.figure5b(rma.assessment.133$FC)
+
+## Comparison plot
+## hgu95a
+affycomp.compfig5b(list(roblox.hgu95a$FC,
+ roblox.hgu95a32$FC,
+ roblox.hgu95a.pmonly$FC,
+ roblox.hgu95a.pmonly32$FC,
+ mas5.assessment$FC,
+ rma.assessment$FC),
+ method.names = c("roblox", "roblox + 32", "roblox (PM)",
+ "roblox + 32 (PM)", "MAS 5.0", "RMA"))
+
+## hgu133a
+affycomp.compfig5b(list(roblox.hgu133a$FC,
+ roblox.hgu133a32$FC,
+ roblox.hgu133a.pmonly$FC,
+ roblox.hgu133a.pmonly32$FC,
+ mas5.assessment$FC,
+ rma.assessment$FC),
+ method.names = c("roblox", "roblox + 32", "roblox (PM)",
+ "roblox + 32 (PM)", "MAS 5.0", "RMA"))
+
+
+###############################################################################
+## Figure 6:
+## a) Observed log fold changes plotted against nominal log fold changes. The
+## dashed lines represent highest, 25th highest, 100th highest, 25th
+## percentile, 75th percentile, smallest 100th, smallest 25th, and smallest
+## log fold change for the genes that were not differentially expressed.
+## b) Like a) but the observed fold changes were calculated for spiked in
+## genes with nominal concentrations no higher than 2pM.
+###############################################################################
+## Figure 6a:
+## hgu95a
+par(mfrow = c(3, 2))
+affycomp.figure6a(roblox.hgu95a$FC)
+affycomp.figure6a(roblox.hgu95a32$FC)
+affycomp.figure6a(roblox.hgu95a.pmonly$FC)
+affycomp.figure6a(roblox.hgu95a.pmonly32$FC)
+affycomp.figure6a(mas5.assessment$FC)
+affycomp.figure6a(rma.assessment$FC)
+
+## hgu133a
+par(mfrow = c(3, 2))
+affycomp.figure6a(roblox.hgu133a$FC)
+affycomp.figure6a(roblox.hgu133a32$FC)
+affycomp.figure6a(roblox.hgu133a.pmonly$FC)
+affycomp.figure6a(roblox.hgu133a.pmonly32$FC)
+affycomp.figure6a(mas5.assessment.133$FC)
+affycomp.figure6a(rma.assessment.133$FC)
+
+## Figure 6b:
+## hgu95a
+par(mfrow = c(3, 2))
+affycomp.figure6b(roblox.hgu95a$FC)
+affycomp.figure6b(roblox.hgu95a32$FC)
+affycomp.figure6b(roblox.hgu95a.pmonly$FC)
+affycomp.figure6b(roblox.hgu95a.pmonly32$FC)
+affycomp.figure6b(mas5.assessment$FC)
+affycomp.figure6b(rma.assessment$FC)
+
+## hgu133a
+par(mfrow = c(3, 2))
+affycomp.figure6b(roblox.hgu133a$FC)
+affycomp.figure6b(roblox.hgu133a32$FC)
+affycomp.figure6b(roblox.hgu133a.pmonly$FC)
+affycomp.figure6b(roblox.hgu133a.pmonly32$FC)
+affycomp.figure6b(mas5.assessment.133$FC)
+affycomp.figure6b(rma.assessment.133$FC)
+
+
+###############################################################################
+## Table
+###############################################################################
+## hgu95a
+round(tableAll(roblox.hgu95a, roblox.hgu95a32, roblox.hgu95a.pmonly,
+ roblox.hgu95a.pmonly32, mas5.assessment, rma.assessment), 4)
+
+## hgu133a
+round(tableAll(roblox.hgu133a, roblox.hgu133a32, roblox.hgu133a.pmonly,
+ roblox.hgu133a.pmonly32, mas5.assessment.133, rma.assessment.133), 4)
+
+## smaller table, more informative ...
+## hgu95a
+affycompTable(roblox.hgu95a, mas5.assessment, rma.assessment)
+
+## hgu133a
+affycompTable(roblox.hgu133a, mas5.assessment.133, rma.assessment.133)
Property changes on: pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
===================================================================
--- pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2009-03-04 08:58:34 UTC (rev 260)
+++ pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2009-03-05 19:09:04 UTC (rev 261)
@@ -13,7 +13,7 @@
\tabular{ll}{
Package: \tab RobLoxBioC\cr
Version: \tab 0.1 \cr
-Date: \tab 2009-03-03 \cr
+Date: \tab 2009-03-05 \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-04 08:58:34 UTC (rev 260)
+++ pkg/RobLoxBioC/man/robloxbioc.Rd 2009-03-05 19:09:04 UTC (rev 261)
@@ -14,7 +14,8 @@
\S4method{robloxbioc}{matrix}(x, eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 1L, mad0 = 1e-4)
-\S4method{robloxbioc}{AffyBatch}(x, pmcorrect = "roblox", normalize = "FALSE", verbose = TRUE,
+\S4method{robloxbioc}{AffyBatch}(x, bg.correct = TRUE, pmcorrect = TRUE, normalize = FALSE,
+ add.constant = 32, verbose = TRUE,
eps = NULL, eps.lower = 0, eps.upper = 0.1, steps = 1L, mad0 = 1e-4,
contrast.tau = 0.03, scale.tau = 10, delta = 2^(-20), sc = 500)
}
@@ -29,10 +30,15 @@
upper bound for the amount of gross errors. See details below. }
\item{steps}{ positive integer. k-step is used to compute the optimally robust estimator. }
\item{mad0}{ scale estimate used if computed MAD is equal to zero}
- \item{pmcorrect}{ method used for PM correction; so far only "roblox" and \code{pm.only}
- are implemented. The algorithm is comparable to the algorithm of MAS 5.0; confer
- \code{\link[affy]{pmcorrect.mas}}. }
+ \item{bg.correct}{ if \code{TRUE} MAS 5.0 background correction is performed;
+ confer \code{\link[affy]{bg.correct.mas}}. }
+ \item{pmcorrect}{ method used for PM correction; \code{TRUE} calls an algorithm which is
+ comparable to the algorithm of MAS 5.0; confer \code{\link[affy]{pmcorrect.mas}}.
+ If \code{FALSE} only the PM intensities are used. }
\item{normalize}{ logical: if \code{TRUE}, Affymetrix scale normalization is performed. }
+ \item{add.constant}{ constant added to the MAS 5.0 expression values before the normalization
+ step. Improves the variance of the measure one no longer devides by numbers close to 0
+ when computing fold-changes. }
\item{verbose}{ logical: if \code{TRUE}, some messages are printed. }
\item{contrast.tau}{ a number denoting the contrast tau parameter; confer the MAS 5.0
PM correction algorithm. }
@@ -57,12 +63,13 @@
the case, try to find a rough estimate for the amount of gross errors, such that
it lies between \code{eps.lower} and \code{eps.upper}.
- If \code{eps} is missing, the radius-minimax (rmx) estimator in sense of
+ If \code{eps} is \code{NULL}, the radius-minimax (rmx) estimator in sense of
Rieder et al. (2001, 2008), respectively Section 2.2 of Kohl (2005) is used.
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.
+ k-step (k >= 1) estimator. The optional scale normalization is performed as given
+ in Affymetrix (2002).
}
\value{ Return value depends on the class of \code{x}.
In case of \code{"matrix"} a matrix with columns "mean" and "sd" is returned.
@@ -107,6 +114,8 @@
## using Affymetrix-Data
## confer example to generateExprVal.method.mas
+## A more worked out example can be found in the scripts folder
+## of the package.
data(SpikeIn)
probes <- pm(SpikeIn)
mas <- generateExprVal.method.mas(probes)
@@ -120,6 +129,7 @@
if(require(affydata)){
data(Dilution)
eset <- robloxbioc(Dilution)
+ ## Affymetrix scale normalization
eset1 <- robloxbioc(Dilution, normalize = TRUE)
}
}
More information about the Robast-commits
mailing list