[Robast-commits] r438 - branches/robast-0.9/pkg/ROptEstOld/R branches/robast-0.9/pkg/RobAStBase/R branches/robast-0.9/pkg/RobLoxBioC/R branches/robast-0.9/pkg/RobLoxBioC/inst/scripts pkg/ROptEstOld/R pkg/RobAStBase/R pkg/RobLoxBioC/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 2 16:19:30 CET 2010
Author: ruckdeschel
Date: 2010-12-02 16:19:30 +0100 (Thu, 02 Dec 2010)
New Revision: 438
Added:
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R
Modified:
branches/robast-0.9/pkg/ROptEstOld/R/AllPlot.R
branches/robast-0.9/pkg/ROptEstOld/R/infoPlot.R
branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R
branches/robast-0.9/pkg/RobAStBase/R/comparePlot.R
branches/robast-0.9/pkg/RobAStBase/R/infoPlot.R
branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R
branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
pkg/ROptEstOld/R/AllPlot.R
pkg/ROptEstOld/R/infoPlot.R
pkg/RobAStBase/R/AllPlot.R
pkg/RobAStBase/R/comparePlot.R
pkg/RobAStBase/R/infoPlot.R
pkg/RobLoxBioC/R/AffySimStudyFunction.R
pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
Log:
added argument no.readonly = TRUE in all assignments of type opar <- par();
[RobLoxBioC] had forgotten to also check new items from trunk when creating branch robast-0.9
now added AffymetrixReproducibility.R, AffymetrixSimStudy.R, IlluminaReproducibility.R
from [trunc] scripts folder; @Matthias pls check whether I missed something else...
Modified: branches/robast-0.9/pkg/ROptEstOld/R/AllPlot.R
===================================================================
--- branches/robast-0.9/pkg/ROptEstOld/R/AllPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/ROptEstOld/R/AllPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -37,7 +37,7 @@
w0 <- options("warn")
options(warn = -1)
- opar <- par()
+ opar <- par(no.readonly = TRUE)
devNew()
nrows <- trunc(sqrt(dims))
ncols <- ceiling(dims/nrows)
Modified: branches/robast-0.9/pkg/ROptEstOld/R/infoPlot.R
===================================================================
--- branches/robast-0.9/pkg/ROptEstOld/R/infoPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/ROptEstOld/R/infoPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -49,7 +49,7 @@
ncols <- ceiling(dims/nrows)
w0 <- options("warn")
options(warn = -1)
- opar <- par()
+ opar <- par(no.readonly = TRUE)
devNew()
par(mfrow = c(nrows, ncols))
for(i in 1:dims){
Modified: branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -173,7 +173,7 @@
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
- opar <- par()
+ opar <- par(no.readonly = TRUE)
opar$cin <- opar$cra <- opar$csi <- opar$cxy <- opar$din <- NULL
on.exit(par(opar))
if (!withSweave)
Modified: branches/robast-0.9/pkg/RobAStBase/R/comparePlot.R
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/R/comparePlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/RobAStBase/R/comparePlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -226,7 +226,7 @@
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
- opar <- par()
+ opar <- par(no.readonly = TRUE)
opar$cin <- opar$cra <- opar$csi <- opar$cxy <- opar$din <- NULL
if(mfColRow) on.exit(par(opar))
Modified: branches/robast-0.9/pkg/RobAStBase/R/infoPlot.R
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/R/infoPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/RobAStBase/R/infoPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -235,7 +235,7 @@
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
- opar <- par()
+ opar <- par(no.readonly = TRUE)
opar$cin <- opar$cra <- opar$csi <- opar$cxy <- opar$din <- NULL
if(mfColRow) on.exit(par(opar))
# if (!withSweave)
@@ -375,7 +375,7 @@
do.call(lines, args=c(list(x.vec, absInfo, type = plty,
lty = lty, lwd = lwd, col = col), dotsL))
if(with.legend)
- legend(legend.location[[1]],
+ legend(legend.location[[1]],
legend = c("class. opt. IC", objectc),
bg = legend.bg,
lty = c(ltyI, lty), col = c(colI, col),
@@ -420,7 +420,7 @@
do.call(lines, args = c(list(x.vec, yc.vec, type = plty,
lty = ltyI, col = colI, lwd = lwdI), dotsL))
if(with.legend)
- legend(legend.location[[i+in1to.draw]],
+ legend(legend.location[[i+in1to.draw]],
bg = legend.bg,
legend = c("class. opt. IC", objectc),
col = c(colI, col), lwd = c(lwdI, lwd),
Modified: branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -72,7 +72,7 @@
abline(h = 0)
boxplot(Ergebnis2, col = myCol[c(1,2,4)], pch = 20, main = "Scale")
abline(h = 1)
- op <- par(mar = rep(2, 4))
+ op <- par(mar = rep(2, 4), no.readonly = TRUE)
plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
legend("center", c("ML", "Med/MAD", "biweight", "rmx"),
fill = myCol, ncol = 4, cex = 1.5)
Modified: branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2010-12-02 14:24:42 UTC (rev 437)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -79,7 +79,7 @@
abline(h = 0)
boxplot(Ergebnis2, col = myCol, pch = 20, main = "Scale")
abline(h = 1)
- op <- par(mar = rep(2, 4))
+ op <- par(mar = rep(2, 4), no.readonly = TRUE)
plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
legend("center", c("ML", "Med/MAD", "Illumina", "rmx"),
fill = myCol, ncol = 4, cex = 1.5)
Added: branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R (rev 0)
+++ branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -0,0 +1,177 @@
+###############################################################################
+## MAS 5.0 vs. robloxbioc for Uni-RNA samples
+###############################################################################
+
+## load MAQC-I data
+library(MAQCsubsetAFX)
+data(refA)
+data(refB)
+data(refC)
+data(refD)
+
+## Minimum Kolmogorov distance
+## takes about 50 min for each reference set (Core i5 520M with 8 GByte RAM)
+library(RobLoxBioC)
+system.time(minKD.A <- KolmogorovMinDist(refA, Norm()))
+system.time(minKD.B <- KolmogorovMinDist(refB, Norm()))
+system.time(minKD.C <- KolmogorovMinDist(refC, Norm()))
+system.time(minKD.D <- KolmogorovMinDist(refD, Norm()))
+
+## load the results for random normal samples from R-forge
+con <- url("http://robast.r-forge.r-project.org/data/minKD_norm.RData")
+load(file = con)
+close(con)
+
+uni.n <- rep(c(11, 16), 4)
+
+#######################################
+## Figure 4 in Kohl and Deigner (2010)
+#######################################
+resA <- split(as.vector(minKD.A$dist), as.vector(minKD.A$n))[c(4,8)]
+resB <- split(as.vector(minKD.B$dist), as.vector(minKD.B$n))[c(4,8)]
+resC <- split(as.vector(minKD.C$dist), as.vector(minKD.C$n))[c(4,8)]
+resD <- split(as.vector(minKD.D$dist), as.vector(minKD.D$n))[c(4,8)]
+#setEPS(height = 6, width = 9)
+#postscript(file = "Figure4.eps")
+par(mar = c(4, 4, 3, 1))
+boxplot(c(resA, resB, resC, resD), main = "Minimum Kolmogorov distance",
+ ylim = c(0, 0.49), at = c(1, 2, 4, 5, 7, 8, 10, 11), xlim = c(0.5, 14.5),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size",
+ pch = 20)
+boxplot(minKD.norm[,c(7,12)], at = c(13, 14), add = TRUE, pch = 20)
+lines(c(1,2), 1/(2*c(11,16)), lwd = 2)
+lines(c(4,5), 1/(2*c(11,16)), lwd = 2)
+lines(c(7,8), 1/(2*c(11,16)), lwd = 2)
+lines(c(10,11), 1/(2*c(11,16)), lwd = 2)
+lines(c(13,14), 1/(2*c(11,16)), lwd = 2)
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+abline(v = c(3, 6, 9, 12), lwd = 1.5)
+text(c(1.5, 4.5, 7.5, 10.5, 13.5), rep(0.48, 5), labels = c("refA", "refB", "refC", "refD", "normal"), font = 2)
+legend("bottomleft", legend = "minimal possible distance", lty = 1,
+ bg = "white", cex = 0.8)
+dev.off()
+
+## Reference set B
+pdf(file = "minKD_refB.pdf")
+par(mfrow = c(1, 2))
+boxplot(res, main = "Reference set B", ylim = c(0, 0.45),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size")
+lines(1:length(uni.n), 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+
+boxplot(minKD.norm[,c(7,12)], main = "Normal samples", ylim = c(0, 0.45),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size")
+lines(1:length(uni.n), 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+dev.off()
+
+## Reference set C
+pdf(file = "minKD_refC.pdf")
+par(mfrow = c(1, 2))
+boxplot(res, main = "Reference set C", ylim = c(0, 0.45),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size")
+lines(1:length(uni.n), 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+
+boxplot(minKD.norm[,c(7,12)], main = "Normal samples", ylim = c(0, 0.45),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size")
+lines(1:length(uni.n), 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+dev.off()
+
+## Reference set D
+pdf(file = "minKD_refD.pdf")
+par(mfrow = c(1, 2))
+boxplot(res, main = "Reference set D", ylim = c(0, 0.45),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size")
+lines(1:length(uni.n), 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+
+boxplot(minKD.norm[,c(7,12)], main = "Normal samples", ylim = c(0, 0.45),
+ ylab = "minimum Kolmogorov distance", xlab = "sample size")
+lines(1:length(uni.n), 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
+abline(h = c(0.1, 0.15), lty = 2, lwd = 1.5)
+dev.off()
+
+
+## MAS 5.0
+## takes about 4.5 minutes for each reference set (Core i5 520M with 8 GByte RAM)
+system.time(res.mas5.A <- mas5(refA))
+system.time(res.mas5.B <- mas5(refB))
+system.time(res.mas5.C <- mas5(refC))
+system.time(res.mas5.D <- mas5(refD))
+
+## MAS rmx
+## takes about 45 seconds for each reference set (Core i5 520M with 8 GByte RAM)
+library(RobLoxBioC)
+system.time(res.rmx.A <- robloxbioc(refA, normalize = TRUE, add.constant = 0))
+system.time(res.rmx.B <- robloxbioc(refB, normalize = TRUE, add.constant = 0))
+system.time(res.rmx.C <- robloxbioc(refC, normalize = TRUE, add.constant = 0))
+system.time(res.rmx.D <- robloxbioc(refD, normalize = TRUE, add.constant = 0))
+
+## Spearman correlations
+cor.mas5.A <- cor(exprs(res.mas5.A), method = "spearman")
+cor.rmx.A <- cor(exprs(res.rmx.A), method = "spearman")
+(diff.A <- cor.rmx.A-cor.mas5.A)
+(rel.A <- cor.rmx.A/cor.mas5.A)
+
+cor.mas5.B <- cor(exprs(res.mas5.B), method = "spearman")
+cor.rmx.B <- cor(exprs(res.rmx.B), method = "spearman")
+(diff.B <- cor.rmx.B-cor.mas5.B)
+(rel.B <- cor.rmx.B/cor.mas5.B)
+
+cor.mas5.C <- cor(exprs(res.mas5.C), method = "spearman")
+cor.rmx.C <- cor(exprs(res.rmx.C), method = "spearman")
+(diff.C <- cor.rmx.C-cor.mas5.C)
+(rel.C <- cor.rmx.C/cor.mas5.C)
+
+cor.mas5.D <- cor(exprs(res.mas5.D), method = "spearman")
+cor.rmx.D <- cor(exprs(res.rmx.D), method = "spearman")
+(diff.D <- cor.rmx.D-cor.mas5.D)
+(rel.D <- cor.rmx.D/cor.mas5.D)
+
+100*(range(rel.A[col(rel.A) > row(rel.A)])-1)
+100*(range(rel.B[col(rel.B) > row(rel.B)])-1)
+100*(range(rel.C[col(rel.C) > row(rel.C)])-1)
+100*(range(rel.D[col(rel.D) > row(rel.D)])-1)
+range(diff.A[col(diff.A) > row(diff.A)])
+range(diff.B[col(diff.B) > row(diff.B)])
+range(diff.C[col(diff.C) > row(diff.C)])
+range(diff.D[col(diff.D) > row(diff.D)])
+
+## Pearson correlations of log2-transformed data
+cor.mas5.A1 <- cor(log2(exprs(res.mas5.A)))
+cor.rmx.A1 <- cor(log2(exprs(res.rmx.A)))
+(diff.A1 <- cor.rmx.A1-cor.mas5.A1)
+(rel.A1 <- cor.rmx.A1/cor.mas5.A1)
+
+cor.mas5.B1 <- cor(log2(exprs(res.mas5.B)))
+cor.rmx.B1 <- cor(log2(exprs(res.rmx.B)))
+(diff.B1 <- cor.rmx.B1-cor.mas5.B1)
+(rel.B1 <- cor.rmx.B1/cor.mas5.B1)
+
+cor.mas5.C1 <- cor(log2(exprs(res.mas5.C)))
+cor.rmx.C1 <- cor(log2(exprs(res.rmx.C)))
+(diff.C1 <- cor.rmx.C1-cor.mas5.C1)
+(rel.C1 <- cor.rmx.C1/cor.mas5.C1)
+
+cor.mas5.D1 <- cor(log2(exprs(res.mas5.D)))
+cor.rmx.D1 <- cor(log2(exprs(res.rmx.D)))
+(diff.D1 <- cor.rmx.D1-cor.mas5.D1)
+(rel.D1 <- cor.rmx.D1/cor.mas5.D1)
+
+100*(range(rel.A1[col(rel.A1) > row(rel.A1)])-1)
+100*(range(rel.B1[col(rel.B1) > row(rel.B1)])-1)
+100*(range(rel.C1[col(rel.C1) > row(rel.C1)])-1)
+100*(range(rel.D1[col(rel.D1) > row(rel.D1)])-1)
+range(diff.A1[col(diff.A1) > row(diff.A1)])
+range(diff.B1[col(diff.B1) > row(diff.B1)])
+range(diff.C1[col(diff.C1) > row(diff.C1)])
+range(diff.D1[col(diff.D1) > row(diff.D1)])
+
Added: branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R (rev 0)
+++ branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -0,0 +1,100 @@
+###############################################################################
+## Simulation study comparing Tukey's biweight with the rmx estimator
+###############################################################################
+
+## fixed variables
+n <- 11
+M <- 1e5
+eps.lower <- 0
+eps.upper <- 0.05
+seed <- 123 ## due to 1e5 replications the influence of the seed is neglectable
+
+eps <- 0.01
+contD <- Norm(0, 9)
+(res1 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Td(df = 3)
+(res2 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Cauchy()
+(res3 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(3, 1)
+(res4 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(10, 1)
+(res5 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Dirac(1.51)
+(res6 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Dirac(1000)
+(res7 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+
+eps <- 0.02
+contD <- Norm(0, 9)
+(res11 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Td(df = 3)
+(res12 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Cauchy()
+(res13 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(3, 1)
+(res14 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(10, 1)
+(res15 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Dirac(1.51)
+(res16 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Dirac(1000)
+(res17 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+
+eps <- 0.04
+contD <- Norm(0, 9)
+(res21 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Td(df = 3)
+(res22 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Cauchy()
+(res23 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(3, 1)
+(res24 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(10, 1)
+(res25 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Norm(1.51)
+(res26 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
+contD <- Dirac(1000)
+(res27 <- AffySimStudy(n = n, M = M, eps = eps, seed = seed,
+ eps.lower = eps.lower, eps.upper = eps.upper,
+ contD = contD))
Added: branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R (rev 0)
+++ branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -0,0 +1,268 @@
+###############################################################################
+## Illumina's default method vs. robloxbioc for technical replicates
+###############################################################################
+
+###############################################################################
+## References:
+## Dunning, M.J., Smith, M.L., Ritchie, M.E., Tavare, S.:
+## beadarray: R classes and methods for Illumina bead-based data.
+## Bioinformatics 2007, 23(16):2183-4.
+##
+## M.J. Dunning, N.L. Barbosa-Morais, A.G. Lynch, S. Tavaré and M.E. Ritchie:
+## Statistical issues in the analysis of Illumina data.
+## BMC Bioinformatics 2008, 9:85.
+###############################################################################
+
+###############################################################################
+## Data:
+## Can be obtained via
+## http://www.compbio.group.cam.ac.uk/Resources/spike/index.html
+###############################################################################
+
+## Load the required packages
+library(beadarray)
+library(RobLoxBioC)
+
+###############################################################################
+## Extract all *.zip file to directory "SpikeInData".
+## Copy spike_targets.txt to directory "SpikeInData".
+##
+## Code to read the bead level data from the directory "SpikeInData"
+##
+## NOTE: reading in the raw data for the entire experiment requires at
+## least 4Gb of RAM for each processing method.
+###############################################################################
+
+###########################################################
+## Read targets information
+targets <- read.table("./SpikeInData/spike_targets.txt",header=TRUE)
+arraynms <- as.character(targets$ArrayNo)
+
+## Use sharpened, subtracted data from text files
+spikeInData <- readIllumina(path = "./SpikeInData", arrayNames=arraynms[1:2],
+ useImages=FALSE, textType=".csv")
+#save(spikeInData, compress = TRUE, file = "spikeInData.RData")
+#load(file = "spikeInData.RData")
+
+## takes about 80 sec (Core i5 520M with 8 GByte RAM)
+system.time(res.ill <- createBeadSummaryData(spikeInData, log = TRUE, imagesPerArray = 2))
+#save(res.ill, file = "SpikeInData_Ill.RData")
+#load(file = "SpikeInData_Ill.RData")
+
+## takes about 490 sec (Core i5 520M with 8 GByte RAM)
+system.time(res.rmx <- robloxbioc(spikeInData, imagesPerArray = 2))
+#save(res.rmx, file = "SpikeInData_rmx.RData")
+#load(file = "SpikeInData_rmx.RData")
+
+###########################################################
+## From Dunning et al. (2008):
+## "The spikes were added at concentrations of 1000, 300, 100, 30, 10 and 3 pM
+## on the six arrays from the first four BeadChips. A further four chips were
+## hybridised with spikes at concentrations of 1, 0.3, 0.1, 0.03, 0.01 and 0 pM.
+## The spikes on a given array were all added at the same concentration. Each
+## concentration was allocated to the same position on all replicate BeadChips.
+## For example, 1000 pM was always array 1 on a chip and 300 pM was array 2 and
+## so on."
+###########################################################
+
+A <- c(1,7,13,19)
+cor.ill.A <- cor(exprs(res.ill[,A]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.A <- cor(exprs(res.rmx[,A]), method = "spearman", use = "pairwise.complete.obs")
+(diff.A <- cor.rmx.A-cor.ill.A)
+(rel.A <- cor.rmx.A/cor.ill.A)
+range(diff.A[col(diff.A) > row(diff.A)])
+100*(range(rel.A[col(rel.A) > row(rel.A)])-1)
+
+cor.ill.A1 <- cor(exprs(res.ill[,A]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.A1 <- cor(exprs(res.rmx[,A]), method = "pearson", use = "pairwise.complete.obs")
+(diff.A1 <- cor.rmx.A1-cor.ill.A1)
+(rel.A1 <- cor.rmx.A1/cor.ill.A1)
+range(diff.A1[col(diff.A1) > row(diff.A1)])
+100*(range(rel.A1[col(rel.A1) > row(rel.A1)])-1)
+
+
+B <- c(2,8,14,20)
+cor.ill.B <- cor(exprs(res.ill[,B]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.B <- cor(exprs(res.rmx[,B]), method = "spearman", use = "pairwise.complete.obs")
+(diff.B <- cor.rmx.B-cor.ill.B)
+(rel.B <- cor.rmx.B/cor.ill.B)
+range(diff.B[col(diff.B) > row(diff.B)])
+100*(range(rel.B[col(rel.B) > row(rel.B)])-1)
+
+cor.ill.B1 <- cor(exprs(res.ill[,B]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.B1 <- cor(exprs(res.rmx[,B]), method = "pearson", use = "pairwise.complete.obs")
+(diff.B1 <- cor.rmx.B1-cor.ill.B1)
+(rel.B1 <- cor.rmx.B1/cor.ill.B1)
+range(diff.B1[col(diff.B1) > row(diff.B1)])
+100*(range(rel.B1[col(rel.B1) > row(rel.B1)])-1)
+
+
+C <- c(3,9,15,21)
+cor.ill.C <- cor(exprs(res.ill[,C]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.C <- cor(exprs(res.rmx[,C]), method = "spearman", use = "pairwise.complete.obs")
+(diff.C <- cor.rmx.C-cor.ill.C)
+(rel.C <- cor.rmx.C/cor.ill.C)
+range(diff.C[col(diff.C) > row(diff.C)])
+100*(range(rel.C[col(rel.C) > row(rel.C)])-1)
+
+cor.ill.C1 <- cor(exprs(res.ill[,C]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.C1 <- cor(exprs(res.rmx[,C]), method = "pearson", use = "pairwise.complete.obs")
+(diff.C1 <- cor.rmx.C1-cor.ill.C1)
+(rel.C1 <- cor.rmx.C1/cor.ill.C1)
+range(diff.C1[col(diff.C1) > row(diff.C1)])
+100*(range(rel.C1[col(rel.C1) > row(rel.C1)])-1)
+
+
+D <- c(4,10,16,22)
+cor.ill.D <- cor(exprs(res.ill[,D]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.D <- cor(exprs(res.rmx[,D]), method = "spearman", use = "pairwise.complete.obs")
+(diff.D <- cor.rmx.D-cor.ill.D)
+(rel.D <- cor.rmx.D/cor.ill.D)
+range(diff.D[col(diff.D) > row(diff.D)])
+100*(range(rel.D[col(rel.D) > row(rel.D)])-1)
+
+cor.ill.D1 <- cor(exprs(res.ill[,D]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.D1 <- cor(exprs(res.rmx[,D]), method = "pearson", use = "pairwise.complete.obs")
+(diff.D1 <- cor.rmx.D1-cor.ill.D1)
+(rel.D1 <- cor.rmx.D1/cor.ill.D1)
+range(diff.D1[col(diff.D1) > row(diff.D1)])
+100*(range(rel.D1[col(rel.D1) > row(rel.D1)])-1)
+
+
+E <- c(5,11,17,23)
+cor.ill.E <- cor(exprs(res.ill[,E]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.E <- cor(exprs(res.rmx[,E]), method = "spearman", use = "pairwise.complete.obs")
+(diff.E <- cor.rmx.E-cor.ill.E)
+(rel.E <- cor.rmx.E/cor.ill.E)
+range(diff.E[col(diff.E) > row(diff.E)])
+100*(range(rel.E[col(rel.E) > row(rel.E)])-1)
+
+cor.ill.E1 <- cor(exprs(res.ill[,E]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.E1 <- cor(exprs(res.rmx[,E]), method = "pearson", use = "pairwise.complete.obs")
+(diff.E1 <- cor.rmx.E1-cor.ill.E1)
+(rel.E1 <- cor.rmx.E1/cor.ill.E1)
+range(diff.E1[col(diff.E1) > row(diff.E1)])
+100*(range(rel.E1[col(rel.E1) > row(rel.E1)])-1)
+
+
+F <- c(6,12,18,24)
+cor.ill.F <- cor(exprs(res.ill[,F]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.F <- cor(exprs(res.rmx[,F]), method = "spearman", use = "pairwise.complete.obs")
+(diff.F <- cor.rmx.F-cor.ill.F)
+(rel.F <- cor.rmx.F/cor.ill.F)
+range(diff.F[col(diff.F) > row(diff.F)])
+100*(range(rel.F[col(rel.F) > row(rel.F)])-1)
+
+cor.ill.F1 <- cor(exprs(res.ill[,F]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.F1 <- cor(exprs(res.rmx[,F]), method = "pearson", use = "pairwise.complete.obs")
+(diff.F1 <- cor.rmx.F1-cor.ill.F1)
+(rel.F1 <- cor.rmx.F1/cor.ill.F1)
+range(diff.F1[col(diff.F1) > row(diff.F1)])
+100*(range(rel.F1[col(rel.F1) > row(rel.F1)])-1)
+
+
+#######################################
+## second series of dilutions
+#######################################
+A <- c(25,31,37,43)
+cor.ill.A <- cor(exprs(res.ill[,A]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.A <- cor(exprs(res.rmx[,A]), method = "spearman", use = "pairwise.complete.obs")
+(diff.A <- cor.rmx.A-cor.ill.A)
+(rel.A <- cor.rmx.A/cor.ill.A)
+range(diff.A[col(diff.A) > row(diff.A)])
+100*(range(rel.A[col(rel.A) > row(rel.A)])-1)
+
+cor.ill.A1 <- cor(exprs(res.ill[,A]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.A1 <- cor(exprs(res.rmx[,A]), method = "pearson", use = "pairwise.complete.obs")
+(diff.A1 <- cor.rmx.A1-cor.ill.A1)
+(rel.A1 <- cor.rmx.A1/cor.ill.A1)
+range(diff.A1[col(diff.A1) > row(diff.A1)])
+100*(range(rel.A1[col(rel.A1) > row(rel.A1)])-1)
+
+
+B <- c(26,32,38,44)
+cor.ill.B <- cor(exprs(res.ill[,B]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.B <- cor(exprs(res.rmx[,B]), method = "spearman", use = "pairwise.complete.obs")
+(diff.B <- cor.rmx.B-cor.ill.B)
+(rel.B <- cor.rmx.B/cor.ill.B)
+range(diff.B[col(diff.B) > row(diff.B)])
+100*(range(rel.B[col(rel.B) > row(rel.B)])-1)
+
+cor.ill.B1 <- cor(exprs(res.ill[,B]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.B1 <- cor(exprs(res.rmx[,B]), method = "pearson", use = "pairwise.complete.obs")
+(diff.B1 <- cor.rmx.B1-cor.ill.B1)
+(rel.B1 <- cor.rmx.B1/cor.ill.B1)
+range(diff.B1[col(diff.B1) > row(diff.B1)])
+100*(range(rel.B1[col(rel.B1) > row(rel.B1)])-1)
+
+
+C <- c(27,33,39,45)
+cor.ill.C <- cor(exprs(res.ill[,C]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.C <- cor(exprs(res.rmx[,C]), method = "spearman", use = "pairwise.complete.obs")
+(diff.C <- cor.rmx.C-cor.ill.C)
+(rel.C <- cor.rmx.C/cor.ill.C)
+range(diff.C[col(diff.C) > row(diff.C)])
+100*(range(rel.C[col(rel.C) > row(rel.C)])-1)
+100*(1/rel.C[col(rel.C) > row(rel.C)]-1)
+
+
+cor.ill.C1 <- cor(exprs(res.ill[,C]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.C1 <- cor(exprs(res.rmx[,C]), method = "pearson", use = "pairwise.complete.obs")
+(diff.C1 <- cor.rmx.C1-cor.ill.C1)
+(rel.C1 <- cor.rmx.C1/cor.ill.C1)
+range(diff.C1[col(diff.C1) > row(diff.C1)])
+100*(range(rel.C1[col(rel.C1) > row(rel.C1)])-1)
+100*(1/rel.C1[col(rel.C1) > row(rel.C1)]-1)
+
+
+D <- c(28,34,40,46)
+cor.ill.D <- cor(exprs(res.ill[,D]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.D <- cor(exprs(res.rmx[,D]), method = "spearman", use = "pairwise.complete.obs")
+(diff.D <- cor.rmx.D-cor.ill.D)
+(rel.D <- cor.rmx.D/cor.ill.D)
+range(diff.D[col(diff.D) > row(diff.D)])
+100*(range(rel.D[col(rel.D) > row(rel.D)])-1)
+100*(1/rel.D[col(rel.D) > row(rel.D)]-1)
+
+cor.ill.D1 <- cor(exprs(res.ill[,D]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.D1 <- cor(exprs(res.rmx[,D]), method = "pearson", use = "pairwise.complete.obs")
+(diff.D1 <- cor.rmx.D1-cor.ill.D1)
+(rel.D1 <- cor.rmx.D1/cor.ill.D1)
+range(diff.D1[col(diff.D1) > row(diff.D1)])
+100*(range(rel.D1[col(rel.D1) > row(rel.D1)])-1)
+100*(1/rel.D1[col(rel.D1) > row(rel.D1)]-1)
+
+
+E <- c(29,35,41,47)
+cor.ill.E <- cor(exprs(res.ill[,E]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.E <- cor(exprs(res.rmx[,E]), method = "spearman", use = "pairwise.complete.obs")
+(diff.E <- cor.rmx.E-cor.ill.E)
+(rel.E <- cor.rmx.E/cor.ill.E)
+range(diff.E[col(diff.E) > row(diff.E)])
+100*(range(rel.E[col(rel.E) > row(rel.E)])-1)
+
+cor.ill.E1 <- cor(exprs(res.ill[,E]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.E1 <- cor(exprs(res.rmx[,E]), method = "pearson", use = "pairwise.complete.obs")
+(diff.E1 <- cor.rmx.E1-cor.ill.E1)
+(rel.E1 <- cor.rmx.E1/cor.ill.E1)
+range(diff.E1[col(diff.E1) > row(diff.E1)])
+100*(range(rel.E1[col(rel.E1) > row(rel.E1)])-1)
+
+
+F <- c(30,36,42,48)
+cor.ill.F <- cor(exprs(res.ill[,F]), method = "spearman", use = "pairwise.complete.obs")
+cor.rmx.F <- cor(exprs(res.rmx[,F]), method = "spearman", use = "pairwise.complete.obs")
+(diff.F <- cor.rmx.F-cor.ill.F)
+(rel.F <- cor.rmx.F/cor.ill.F)
+range(diff.F[col(diff.F) > row(diff.F)])
+100*(range(rel.F[col(rel.F) > row(rel.F)])-1)
+
+cor.ill.F1 <- cor(exprs(res.ill[,F]), method = "pearson", use = "pairwise.complete.obs")
+cor.rmx.F1 <- cor(exprs(res.rmx[,F]), method = "pearson", use = "pairwise.complete.obs")
+(diff.F1 <- cor.rmx.F1-cor.ill.F1)
+(rel.F1 <- cor.rmx.F1/cor.ill.F1)
+range(diff.F1[col(diff.F1) > row(diff.F1)])
+100*(range(rel.F1[col(rel.F1) > row(rel.F1)])-1)
+
+
+
Modified: pkg/ROptEstOld/R/AllPlot.R
===================================================================
--- pkg/ROptEstOld/R/AllPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ pkg/ROptEstOld/R/AllPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -37,7 +37,7 @@
w0 <- options("warn")
options(warn = -1)
- opar <- par()
+ opar <- par(no.readonly = TRUE)
devNew()
nrows <- trunc(sqrt(dims))
ncols <- ceiling(dims/nrows)
Modified: pkg/ROptEstOld/R/infoPlot.R
===================================================================
--- pkg/ROptEstOld/R/infoPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ pkg/ROptEstOld/R/infoPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -49,7 +49,7 @@
ncols <- ceiling(dims/nrows)
w0 <- options("warn")
options(warn = -1)
- opar <- par()
+ opar <- par(no.readonly = TRUE)
devNew()
par(mfrow = c(nrows, ncols))
for(i in 1:dims){
Modified: pkg/RobAStBase/R/AllPlot.R
===================================================================
--- pkg/RobAStBase/R/AllPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ pkg/RobAStBase/R/AllPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -173,7 +173,7 @@
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
- opar <- par()
+ opar <- par(no.readonly = TRUE)
opar$cin <- opar$cra <- opar$csi <- opar$cxy <- opar$din <- NULL
on.exit(par(opar))
if (!withSweave)
Modified: pkg/RobAStBase/R/comparePlot.R
===================================================================
--- pkg/RobAStBase/R/comparePlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ pkg/RobAStBase/R/comparePlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -226,7 +226,7 @@
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
- opar <- par()
+ opar <- par(no.readonly = TRUE)
opar$cin <- opar$cra <- opar$csi <- opar$cxy <- opar$din <- NULL
if(mfColRow) on.exit(par(opar))
Modified: pkg/RobAStBase/R/infoPlot.R
===================================================================
--- pkg/RobAStBase/R/infoPlot.R 2010-12-02 14:24:42 UTC (rev 437)
+++ pkg/RobAStBase/R/infoPlot.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -235,7 +235,7 @@
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
- opar <- par()
+ opar <- par(no.readonly = TRUE)
opar$cin <- opar$cra <- opar$csi <- opar$cxy <- opar$din <- NULL
if(mfColRow) on.exit(par(opar))
# if (!withSweave)
Modified: pkg/RobLoxBioC/R/AffySimStudyFunction.R
===================================================================
--- pkg/RobLoxBioC/R/AffySimStudyFunction.R 2010-12-02 14:24:42 UTC (rev 437)
+++ pkg/RobLoxBioC/R/AffySimStudyFunction.R 2010-12-02 15:19:30 UTC (rev 438)
@@ -41,7 +41,7 @@
ind <- sample(1:M, min(M, 20))
if(plot1) dev.new()
print(
- stripplot(rep(1:min(M, 20), each = n) ~ as.vector(Mre[ind,]),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 438
More information about the Robast-commits
mailing list