[Robast-commits] r410 - in pkg/RobLoxBioC: . inst inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 27 13:11:33 CEST 2010
Author: stamats
Date: 2010-08-27 13:11:33 +0200 (Fri, 27 Aug 2010)
New Revision: 410
Added:
pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R
pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R
pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R
Removed:
pkg/RobLoxBioC/inst/scripts/AffySimStudy.R
Modified:
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/inst/NEWS
pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
pkg/RobLoxBioC/inst/scripts/IlluminaExample.R
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
Log:
some modifications/additions needed for the revised version of the corresponding manuscript.
Modified: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION 2010-08-25 03:02:36 UTC (rev 409)
+++ pkg/RobLoxBioC/DESCRIPTION 2010-08-27 11:11:33 UTC (rev 410)
@@ -1,6 +1,6 @@
Package: RobLoxBioC
-Version: 0.7.1
-Date: 2010-04-01
+Version: 0.7.2
+Date: 2010-08-27
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.
Modified: pkg/RobLoxBioC/inst/NEWS
===================================================================
--- pkg/RobLoxBioC/inst/NEWS 2010-08-25 03:02:36 UTC (rev 409)
+++ pkg/RobLoxBioC/inst/NEWS 2010-08-27 11:11:33 UTC (rev 410)
@@ -8,6 +8,15 @@
information)
#######################################
+## version 0.7.2
+#######################################
++ renamed AffySimStudy.R to AffymetrixSimStudy.R
++ update of Figure numbers in IlluminaExample.R
++ added scripts AffymetrixTechnicalReplicates.R and
+ IlluminaTechnicalReplicates.R
+
+
+#######################################
## version 0.7.1
#######################################
+ corrected bug in AffySimStudy and IlluminaSimStudy
Deleted: pkg/RobLoxBioC/inst/scripts/AffySimStudy.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffySimStudy.R 2010-08-25 03:02:36 UTC (rev 409)
+++ pkg/RobLoxBioC/inst/scripts/AffySimStudy.R 2010-08-27 11:11:33 UTC (rev 410)
@@ -1,100 +0,0 @@
-###############################################################################
-## 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))
Modified: pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2010-08-25 03:02:36 UTC (rev 409)
+++ pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2010-08-27 11:11:33 UTC (rev 410)
@@ -156,7 +156,7 @@
col = c(rep(NA, 3), grey(0.6), NA, grey(0.4), NA, NA, grey(0.4), rep(NA, 3),
grey(0.6), NA))
abline(v = c(6, 10), lwd = 1.5)
-text(c(3, 8, 14), rep(0.48, 3), labels = c("HGU95A", "HGU133A", "Normal Samples"),
+text(c(3, 8, 13.5), rep(0.48, 3), labels = c("HGU95A", "HGU133A", "Normal Samples"),
font = 2)
lines(1:5, 1/(2*uni.n[8:12]), lwd = 2)
lines(7:9, 1/(2*uni.n[c(6,8,12)]), lwd = 2)
Added: pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R (rev 0)
+++ pkg/RobLoxBioC/inst/scripts/AffymetrixReproducibility.R 2010-08-27 11:11:33 UTC (rev 410)
@@ -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)])
+
Copied: pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R (from rev 409, pkg/RobLoxBioC/inst/scripts/AffySimStudy.R)
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R (rev 0)
+++ pkg/RobLoxBioC/inst/scripts/AffymetrixSimStudy.R 2010-08-27 11:11:33 UTC (rev 410)
@@ -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))
Modified: pkg/RobLoxBioC/inst/scripts/IlluminaExample.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/IlluminaExample.R 2010-08-25 03:02:36 UTC (rev 409)
+++ pkg/RobLoxBioC/inst/scripts/IlluminaExample.R 2010-08-27 11:11:33 UTC (rev 410)
@@ -83,7 +83,7 @@
close(con)
#######################################
-## Figure 4 in Kohl and Deigner (2010)
+## Figure 5 in Kohl and Deigner (2010)
#######################################
res1 <- split(as.vector(minKD.Illumina$dist), as.vector(minKD.Illumina$n))[30:50]
res2 <- split(as.vector(minKD.Illumina.log$dist), as.vector(minKD.Illumina.log$n))[30:50]
@@ -91,7 +91,7 @@
uni.n <- rep(30:50, 3)
#setEPS(height = 6, width = 9)
-#postscript(file = "Figure4.eps")
+#postscript(file = "Figure5.eps")
par(mar = c(4, 4, 3, 1))
plot(0, 0, type = "n", ylim = c(-0.01, 0.4), xlim = c(0.5, 65.5),
panel.first = abline(h = seq(0, 0.35, by = 0.05), lty = 2, col = "grey"),
@@ -116,13 +116,13 @@
#dev.off()
## Comparison of quantiles
-## Figure 5 in Kohl and Deigner (2010)
+## Figure 6 in Kohl and Deigner (2010)
res1 <- split(as.vector(minKD.Illumina$dist), as.vector(minKD.Illumina$n))[15:65]
res2 <- split(as.vector(minKD.Illumina.log$dist), as.vector(minKD.Illumina.log$n))[15:65]
res3 <- lapply(as.data.frame(minKD.Illumina.norm), function(x) x)[6:56]
#setEPS(height = 6, width = 9)
-#postscript(file = "Figure5.eps")
+#postscript(file = "Figure6.eps")
par(mar = c(4, 4, 3, 1))
plot(15:65, sapply(res3, quantile, prob = 0.99), type = "l", lwd = 2, xlab = "sample size",
ylab = "quantile of mimimum Kolmogorov distances",
@@ -166,9 +166,9 @@
ill.SD <- assessSpikeInSD(res.ill, genenames = genenames, method.name = "Illumina")
rmx.SD <- assessSpikeInSD(res.rmx, genenames = genenames, method.name = "rmx estimator")
-## Figure 6 in Kohl and Deigner (2010)
+## Figure 7 in Kohl and Deigner (2010)
setEPS(height = 6, width = 9)
-postscript(file = "Figure6.eps")
+postscript(file = "Figure7.eps")
plot(ill.SD$xsmooth, ill.SD$ysmooth, type = "l", xlab = "mean log expression",
ylab = "mean SD", main = "Spike-in data of Dunning et al. (2008)", lwd = 2,
panel.first = abline(h = c(0.1, 0.12, 0.14, 0.16), v = seq(6, 16, 2), lty = 2, col = "grey"))
Added: pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R (rev 0)
+++ pkg/RobLoxBioC/inst/scripts/IlluminaReproducibility.R 2010-08-27 11:11:33 UTC (rev 410)
@@ -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")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 410
More information about the Robast-commits
mailing list