[Robast-commits] r268 - in pkg/RobLoxBioC: R inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 15 14:13:42 CET 2009
Author: stamats
Date: 2009-03-15 14:13:42 +0100 (Sun, 15 Mar 2009)
New Revision: 268
Added:
pkg/RobLoxBioC/inst/scripts/IlluminaExample.R
Modified:
pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
Log:
added Illumina example based on file by Dunning and Ritchie ...
Modified: pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2009-03-14 11:28:08 UTC (rev 267)
+++ pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2009-03-15 13:13:42 UTC (rev 268)
@@ -172,20 +172,40 @@
})
## computation of bead summaries via robloxbioc for "matrix"
rmxBeadSummary <- function(x, probeIDs, probes, eps, eps.lower, eps.upper, steps, mad0){
+ comIDs <- intersect(probeIDs, probes)
+ x <- x[probeIDs %in% comIDs]
+ probeIDs <- probeIDs[probeIDs %in% comIDs]
noBeads <- as.vector(table(probeIDs))
noBeads.uni <- as.integer(names(table(noBeads)))
+ probes1 <- comIDs
+ len1 <- length(probes1)
+ fore1 <- numeric(len1)
+ SD1 <- numeric(len1)
+ for(i in seq(along = noBeads.uni)){
+ index <- noBeads == noBeads.uni[i]
+ IDs <- probes1[index]
+ if(noBeads.uni[i] == 1){
+ fore1[index] <- x[probeIDs %in% IDs]
+ SD1[index] <- mad0
+ }else{
+ 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)
+ fore1[index] <- rmx[,"mean"]
+ SD1[index] <- rmx[,"sd"]
+ }
+ }
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"]
- }
+ noBeads1 <- numeric(len)
+ nas <- !(probes %in% comIDs)
+ fore[nas] <- NA
+ fore[!nas] <- fore1
+ SD[nas] <- NA
+ SD[!nas] <- SD1
+ noBeads1[nas] <- 0
+ noBeads1[!nas] <- noBeads
- return(list(fore = fore, sd = SD, noBeads = noBeads))
+ return(list(fore = fore, sd = SD, noBeads = noBeads1))
}
Added: pkg/RobLoxBioC/inst/scripts/IlluminaExample.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/IlluminaExample.R (rev 0)
+++ pkg/RobLoxBioC/inst/scripts/IlluminaExample.R 2009-03-15 13:13:42 UTC (rev 268)
@@ -0,0 +1,1583 @@
+###############################################################################
+## Illumina Example
+###############################################################################
+
+###############################################################################
+## 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
+###############################################################################
+
+###############################################################################
+## This example is based on the R code of Mark Dunning and Matt Ritchie
+## available under
+## http://www.compbio.group.cam.ac.uk/Resources/spike/scripts/Analysis.R
+##
+## This file was slightly adapted and code for the computation of
+## rmx-estimators was added.
+###############################################################################
+
+########################################
+## Analysis of Illumina in Spike data set
+##
+## November 2007
+##
+## Mark Dunning and Matt Ritchie
+########################################
+
+## Load the required packages
+library(beadarray)
+library(gplots)
+library(RColorBrewer)
+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.
+## For this reason the data is read in sequentially saved, then removed
+## before moving on to the next method
+###############################################################################
+
+
+###########################################################
+## Read targets information
+targets <- read.table("./SpikeInData/spike_targets.txt",header=TRUE)
+
+
+###########################################################
+## Read text and tif files from all BeadChips from spike experiment
+## sharpen images, no local background subtraction
+bld.sharpen.nobgc <- readIllumina(path = "./SpikeInData", textType=".csv",
+ arrayNames=targets$ArrayNo,
+ imageManipulation="sharpen",
+ backgroundMethod="none")
+save(bld.sharpen.nobgc, file="bld.sharpen.nobgc.rda")
+rm(bld.sharpen.nobgc); gc()
+
+
+
+## Take the values found in the bead level text files as sharpened, subtracted intensities
+bld.sharpen.subtract <- readIllumina(path = "./SpikeInData", textType=".csv",
+ arrayNames=targets$ArrayNo,
+ useImages=FALSE,
+ backgroundMethod="none")
+save(bld.sharpen.subtract, file="bld.sharpen.subtract.rda")
+rm(bld.sharpen.subtract); gc()
+
+## sharpen images, with local background subtraction and
+## normal-exponential convolution model (avoids negative background
+## corrected intensities)
+bld.sharpen.normexp <- readIllumina(path = "./SpikeInData", textType=".csv",
+ arrayNames=targets$ArrayNo,
+ useImages=FALSE,
+ backgroundMethod="normexp")
+save(bld.sharpen.normexp, file="bld.sharpen.normexp.rda")
+rm(bld.sharpen.normexp); gc()
+
+
+
+
+
+## no sharpening, no local background subtraction
+bld.nosharpen.nobgc <- readIllumina(path = "./SpikeInData", textType=".csv",
+ arrayNames=targets$ArrayNo,
+ imageManipulation="none",
+ backgroundMethod="none")
+save(bld.nosharpen.nobgc, file="bld.nosharpen.nobgc.rda")
+rm(bld.nosharpen.nobgc); gc()
+
+## no sharpening with local background subtraction
+bld.nosharpen.subtract <- readIllumina(path = "./SpikeInData", textType=".csv",
+ arrayNames=targets$ArrayNo,
+ imageManipulation="none",
+ backgroundMethod="subtract")
+save(bld.nosharpen.subtract, file="bld.nosharpen.subtract.rda")
+rm(bld.nosharpen.subtract); gc()
+
+
+###########################################################
+## Figure 1 - plots of foreground and background intensities
+###########################################################
+## Read text and tif files from first array only (to save time and memory)
+## sharpen images, no local background subtraction
+B <- readIllumina(path = "./SpikeInData", textType=".csv",
+ arrayNames=targets$ArrayNo[1:12],
+ imageManipulation="sharpen",
+ backgroundMethod="none")
+
+## sharpen images, with local background subtraction
+B.bc <- backgroundCorrect(B, method="subtract")
+
+## Setup colours
+pal <- brewer.pal(6, "Set1")
+
+ylim <- c(8.5,10.5)
+names <- c("Array1 Strip1", "Array1 Strip2", "Array2 Strip1", "Array2 Strip2",
+ "Array3 Strip1", "Array3 Strip2", "Array4 Strip1", "Array4 Strip2",
+ "Array5 Strip1", "Array5 Strip2", "Array6 Strip1", "Array6 Strip2")
+
+pdf("Figure1.pdf", width=12, height=8)
+par(mfrow=c(1,2), mar=c(6,3,1.5,0.15), oma=c(0,1,0,0))
+boxplotBeads(B, names=names, las=2, outline=FALSE, col=rep(pal, each=2),
+ main="Raw Foreground", what="G", ylim=ylim, ylab="", medlwd=1)
+boxplotBeads(B, names=names, las=2, outline=FALSE, col=rep(pal, each=2),
+ main="Raw Background", what="Gb", ylim=ylim, ylab="", medlwd=1)
+mtext("log2 intensity", side=2, outer=TRUE)
+dev.off()
+
+
+###########################################################
+## Derive low-level properties of Illumina data quoted in the text
+## Number of negative corrected intensities from first BeadChip
+negs <- numeric(12)
+for(i in 1:12)
+ negs[i] <- sum(B.bc[[i]]$G<0)
+
+## Percentage negative from first BeadChip - sharpened
+round(negs/numBeads(B.bc)*100,2)
+# [1] 1.04 0.76 0.93 0.68 0.92 0.68 0.94 0.82 0.84 0.74 0.90 0.82
+# old figures [1] 7.87 7.53 7.32 6.95 6.94 6.29 6.18 5.91 5.54 5.25 4.80 4.92
+
+## quartiles
+quantile(round(negs/numBeads(B.bc)*100,2), c(0,0.25,0.5,0.75,1))
+# 0% 25% 50% 75% 100%
+# 0.6800 0.7550 0.8300 0.9225 1.0400
+
+## How many missing beads (contaminated left-hand edge - 0 intensities)
+artefact = NULL
+for(i in 1:12)
+ artefact[i] = sum(B[[i]]$G==0)
+
+## Percentage negative from first BeadChip - sharpened
+round(artefact/numBeads(B)*100,2)
+# [1] 7.32 7.21 6.88 6.71 6.49 6.08 5.70 5.56 5.16 4.96 4.42 4.67
+
+## quartiles
+quantile(round(artefact/numBeads(B)*100,2), c(0,0.25,0.5,0.75,1))
+# 0% 25% 50% 75% 100%
+# 4.4200 5.1100 5.8900 6.7525 7.3200
+
+## Calculate median background level on both original and log2-scale
+bgorig = NULL
+bglog = NULL
+for(i in 1:12) {
+ bgorig[[i]] = getArrayData(B, wh="Gb", log=FALSE, array=i)
+ bglog[[i]] = getArrayData(B, wh="Gb", log=TRUE, array=i)
+}
+
+## log2-scale
+sapply(bglog, FUN="median", na.rm=TRUE)
+# [1] 9.308339 9.308339 9.308339 9.306062 9.306062 9.306062 9.308339 9.306062 9.308339 9.306062 9.308339 9.308339
+
+## original scale
+sapply(bgorig, FUN="median")
+# [1] 634 634 634 633 633 633 634 633 634 633 634 634
+
+
+###########################################################
+## Number of negative corrected intensities from full data set - sharpened
+load("bld.sharpen.subtract.rda")
+narrays <- length(arrayNames(bld.sharpen.subtract))
+
+negs <- numeric(12)
+for(i in 1:narrays)
+ negs[i] = sum(bld.sharpen.subtract[[i]]$G<0)
+
+## Percentage negative from full data set
+round(negs/numBeads(bld.sharpen.subtract)*100,2)
+# [1] 1.04 0.76 0.93 0.68 0.92 0.68 0.94 0.82 0.84 0.74 0.90 0.82 0.29 0.13 0.54
+#[16] 0.15 0.41 0.22 0.29 0.15 0.38 0.14 0.41 0.13 0.61 0.19 0.46 0.13 0.52 0.18
+#[31] 0.38 0.21 0.38 0.17 0.40 0.16 0.29 0.10 0.48 0.11 0.28 0.12 0.32 0.15 0.39
+#[46] 0.15 0.32 0.14 0.50 0.17 0.46 0.26 0.64 0.22 0.57 0.16 0.52 0.34 0.49 0.22
+#[61] 0.61 0.14 0.63 0.29 0.61 0.20 0.54 0.14 0.39 0.20 0.48 0.19 0.50 0.28 0.79
+#[76] 0.32 0.75 0.30 0.87 0.21 0.63 0.32 0.81 0.30 0.22 0.18 0.41 0.22 0.57 0.14
+#[91] 0.46 0.23 0.22 0.12 0.55 0.13
+
+# quartiles
+quantile(round(negs/numBeads(bld.sharpen.subtract)*100,2), c(0,0.25,0.5,0.75,1))
+# 0% 25% 50% 75% 100%
+#0.100 0.190 0.320 0.555 1.040
+
+
+###########################################################
+## Create bead summary data for each pre-processing option using
+## the default Illumina method of removing outliers > 3 MADs
+## from the median before averaging the intensities
+
+## sharpen images, no local background subtraction
+load("bld.sharpen.nobgc.rda")
+bsd.sharpen.nobgc <- createBeadSummaryData(bld.sharpen.nobgc, log=TRUE,
+ imagesPerArray=2, what="G",
+ method="illumina", n=3)
+rm(bld.sharpen.nobgc)
+
+## sharpen images, with local background subtraction
+load("bld.sharpen.subtract.rda")
+bsd.sharpen.subtract <- createBeadSummaryData(bld.sharpen.subtract, log=TRUE,
+ imagesPerArray=2, what="G",
+ method="illumina", n=3)
+bsd.sharpen.subtract.raw <- createBeadSummaryData(bld.sharpen.subtract, log=FALSE,
+ imagesPerArray=2, what="G",
+ method="illumina", n=3)
+rm(bld.sharpen.subtract)
+
+
+## sharpen images, with local background subtraction and
+## normal-exponential convolution model (avoids negative background
+## corrected intensities)
+load("bld.sharpen.normexp.rda")
+bsd.sharpen.normexp <- createBeadSummaryData(bld.sharpen.normexp, log=TRUE,
+ imagesPerArray=2, what="G",
+ method="illumina", n=3)
+rm(bld.sharpen.normexp)
+
+## no sharpening, no local background subtraction
+load("bld.nosharpen.nobgc.rda")
+bsd.nosharpen.nobgc <- createBeadSummaryData(bld.nosharpen.nobgc, log=TRUE, imagesPerArray=2, what="G", method="illumina", n=3)
+rm(bld.nosharpen.nobgc)
+
+## no sharpening, with local background subtraction
+load("bld.nosharpen.subtract.rda")
+bsd.nosharpen.subtract = createBeadSummaryData(bld.nosharpen.subtract, log=TRUE, imagesPerArray=2, what="G", method="illumina", n=3)
+rm(bld.nosharpen.subtract)
+
+## Save summary data
+save(bsd.sharpen.nobgc, bsd.sharpen.subtract, bsd.sharpen.normexp, bsd.nosharpen.nobgc, bsd.nosharpen.subtract, file="bsd.log.ill.rda")
+
+
+###########################################################
+## Figure 2 - plot results from simulation study (varying number of outliers
+## were introduced at the saturation level to assess how many outliers can be
+## tolerated before serious bias is introduced in to the ## expression measures)
+
+targets <- read.table("./SpikeInData/spike_targets.txt", header=TRUE, sep=" ")
+arraynms <- as.character(targets$ArrayNo)
+narrays <- length(arraynms)
+
+## introduce outliers in 2nd BeadChip, at varying %
+## Use sharpened, subtracted data from text files
+bld.sharpen.array2 <- readIllumina(path = "./SpikeInData", arrayNames=arraynms[13:24],
+ useImages=FALSE, textType=".csv")
+
+## percentage of beads to be saturated
+per.out <- c(0, 0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4)
+
+nchips <- 12
+arraynms <- arrayNames(bld.sharpen.array2)
+
+## Saturated
+ill.summ.array2 <- list(NoBeads=list(), exprs=list(), se.exprs=list())
+ill.summ2.array2 <- med.summ.array2 <- mean.summ.array2 <- trim.summ.array2 <- win.summ.array2 <- ill.summ.array2
+rmx.summ.array2 <- ill.summ.array2
+
+set.seed(06092007)
+for(l in 1:length(per.out)) {
+ cat(l, "of 9\n")
+ simdataarray2 <- copyBeadLevelList(bld.sharpen.array2)
+ if(per.out[l]!=0) {
+ for(i in 1:nchips) {
+ cat(i, "of 12\n")
+ ## no gross error model!
+# nsamp <- round(per.out[l]*numBeads(simdataarray2, array=i),0)
+# ind <- sample(1:numBeads(simdataarray2, array=i), nsamp)
+ ## gross error model!
+# ind <- as.logical(rbinom(numBeads(simdataarray2, array=i), prob = per.out[l], size = 1))
+# ind <- (1:numBeads(simdataarray2, array=i))[ind]
+# simdataarray2 at beadData[[arraynms[i]]]$G[ind] <- 2^16
+ ## problem: there are bead types where more than 50% of the values are contaminated
+ ## => there is no meaningful estimator which can handle this!
+ sel <- which(simdataarray2 at beadData[[arraynms[i]]]$ProbeID != 0)
+ pr <- simdataarray2 at beadData[[arraynms[i]]]$ProbeID[sel]
+ probes <- sort(unique(pr))
+ indices <- NULL
+ for(j in seq(along = probes)){
+ ind <- pr == probes[j]
+ ## < 50% of values contaminated
+ repeat{
+ cont <- as.logical(rbinom(sum(ind), prob = per.out[l], size = 1))
+ if(sum(cont) < sum(ind)/2) break
+ }
+ indices <- c(indices, sel[ind][cont])
+ }
+ simdataarray2 at beadData[[arraynms[i]]]$G[indices] <- 2^16
+ }
+ }
+ tmp.ill <- createBeadSummaryData(simdataarray2, method="illumina", log=FALSE, n=3, imagesPerArray=2)
+ tmp.ill2 <- createBeadSummaryData(simdataarray2, method="illumina", log=FALSE, n=2, imagesPerArray=2)
+ tmp.med <- createBeadSummaryData(simdataarray2, method="median", log=FALSE, imagesPerArray=2)
+ tmp.mean <- createBeadSummaryData(simdataarray2, method="mean", log=FALSE, imagesPerArray=2)
+ tmp.trim <- createBeadSummaryData(simdataarray2, method="trim", trim=0.1, log=FALSE, imagesPerArray=2)
+ tmp.win <- createBeadSummaryData(simdataarray2, method="winsorize", trim=0.1, log=FALSE, imagesPerArray=2)
+ tmp.rmx <- robloxbioc(simdataarray2, imagesPerArray=2, eps = max(per.out[l], 0.05))
+ ill.summ.array2$NoBeads[[l]] = NoBeads(tmp.ill)
+ ill.summ.array2$exprs[[l]] = exprs(tmp.ill)
+ ill.summ.array2$se.exprs[[l]] = se.exprs(tmp.ill)
+ ill.summ2.array2$NoBeads[[l]] = NoBeads(tmp.ill2)
+ ill.summ2.array2$exprs[[l]] = exprs(tmp.ill2)
+ ill.summ2.array2$se.exprs[[l]] = se.exprs(tmp.ill2)
+ med.summ.array2$NoBeads[[l]] = NoBeads(tmp.med)
+ med.summ.array2$exprs[[l]] = exprs(tmp.med)
+ med.summ.array2$se.exprs[[l]] = se.exprs(tmp.med)
+ mean.summ.array2$NoBeads[[l]] = NoBeads(tmp.mean)
+ mean.summ.array2$exprs[[l]] = exprs(tmp.mean)
+ mean.summ.array2$se.exprs[[l]] = se.exprs(tmp.mean)
+ win.summ.array2$NoBeads[[l]] = NoBeads(tmp.win)
+ win.summ.array2$exprs[[l]] = exprs(tmp.win)
+ win.summ.array2$se.exprs[[l]] = se.exprs(tmp.win)
+ trim.summ.array2$NoBeads[[l]] = NoBeads(tmp.trim)
+ trim.summ.array2$exprs[[l]] = exprs(tmp.trim)
+ trim.summ.array2$se.exprs[[l]] = se.exprs(tmp.trim)
+ rmx.summ.array2$NoBeads[[l]] = NoBeads(tmp.rmx)
+ rmx.summ.array2$exprs[[l]] = exprs(tmp.rmx)
+ rmx.summ.array2$se.exprs[[l]] = se.exprs(tmp.rmx)
+ rm(simdataarray2, tmp.ill, tmp.med, tmp.mean, tmp.trim, tmp.win, tmp.rmx)
+}
+save(ill.summ.array2, ill.summ2.array2, med.summ.array2, mean.summ.array2,
+ win.summ.array2, trim.summ.array2, rmx.summ.array2, file="sim.summary.bias.raw.rda")
+
+## calculate bias - assume original, complete data set gave true values
+bias.ill.array2 <- list(NoBeads=list(), exprs=list(), se.exprs=list())
+bias.ill2.array2 <- bias.med.array2 <- bias.mean.array2 <- bias.trim.array2 <- bias.win.array2 <- bias.ill.array2
+bias.rmx.array2 <- bias.ill.array2
+for(l in 1:length(per.out)) {
+ cat(l, "\n")
+ bias.ill.array2$exprs[[l]] <- ill.summ.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.ill.array2$se.exprs[[l]] <- ill.summ.array2$se.exprs[[l]]^2*ill.summ.array2$NoBeads[[l]]
+ bias.ill.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-ill.summ.array2$NoBeads[[l]]
+
+ bias.ill2.array2$exprs[[l]] <- ill.summ2.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.ill2.array2$se.exprs[[l]] <- ill.summ2.array2$se.exprs[[l]]^2*ill.summ2.array2$NoBeads[[l]]
+ bias.ill2.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-ill.summ2.array2$NoBeads[[l]]
+
+ bias.med.array2$exprs[[l]] <- med.summ.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.med.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-med.summ.array2$NoBeads[[l]]
+
+ bias.mean.array2$exprs[[l]] <- mean.summ.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.mean.array2$se.exprs[[l]] <- mean.summ.array2$se.exprs[[l]]^2*mean.summ.array2$NoBeads[[l]]
+ bias.mean.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-mean.summ.array2$NoBeads[[l]]
+
+ bias.trim.array2$exprs[[l]] <- trim.summ.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.trim.array2$se.exprs[[l]] <- trim.summ.array2$se.exprs[[l]]^2*trim.summ.array2$NoBeads[[l]]
+ bias.trim.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-trim.summ.array2$NoBeads[[l]]
+
+ bias.win.array2$exprs[[l]] <- win.summ.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.win.array2$se.exprs[[l]] <- win.summ.array2$se.exprs[[l]]^2*win.summ.array2$NoBeads[[l]]
+ bias.win.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-win.summ.array2$NoBeads[[l]]
+
+ bias.rmx.array2$exprs[[l]] <- rmx.summ.array2$exprs[[l]]-mean.summ.array2$exprs[[1]]
+ bias.rmx.array2$se.exprs[[l]] <- rmx.summ.array2$se.exprs[[l]]^2*rmx.summ.array2$NoBeads[[l]]
+ bias.rmx.array2$NoBeads[[l]] <- mean.summ.array2$NoBeads[[1]]-rmx.summ.array2$NoBeads[[l]]
+}
+
+## calculate average bias from `replicate arrays
+avebias.ill.array2 <- avebias.ill2.array2 <- avebias.mean.array2 <- avebias.med.array2 <- avebias.trim.array2 <- avebias.win.array2 <- NULL
+avebias.rmx.array2 <- NULL
+for(l in 1:length(per.out)) {
+ avebias.ill.array2[l] <- mean(as.vector(bias.ill.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.ill2.array2[l] <- mean(as.vector(bias.ill2.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.med.array2[l] <- mean(as.vector(bias.med.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.mean.array2[l] <- mean(as.vector(bias.mean.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.win.array2[l] <- mean(as.vector(bias.win.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.trim.array2[l] <- mean(as.vector(bias.trim.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.rmx.array2[l] <- mean(as.vector(bias.rmx.array2$exprs[[l]]), na.rm=TRUE)
+}
+
+numout.ill.array2 <- numout.ill2.array2 <- numout.mean.array2 <- numout.med.array2 <- numout.trim.array2 <- numout.win.array2 <- NULL
+numout.rmx.array2 <- NULL
+for(l in 1:length(per.out)) {
+ numout.ill.array2[l] <- mean(as.vector(bias.ill.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.ill2.array2[l] <- mean(as.vector(bias.ill2.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.med.array2[l] <- mean(as.vector(bias.med.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.mean.array2[l] <- mean(as.vector(bias.mean.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.win.array2[l] <- mean(as.vector(bias.win.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.trim.array2[l] <- mean(as.vector(bias.trim.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.rmx.array2[l] <- mean(as.vector(bias.rmx.array2$NoBeads[[l]]/mean.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+}
+
+## calculate variance
+avevar.ill.array2 <- avevar.ill2.array2 <- avevar.mean.array2 <- avevar.med.array2 <- avevar.trim.array2 <- avevar.win.array2 <- NULL
+avevar.rmx.array2 <- NULL
+for(l in 1:length(per.out)) {
+ avevar.ill.array2[l] <- mean(as.vector(bias.ill.array2$se.exprs[[l]]), na.rm=TRUE)
+ avevar.ill2.array2[l] <- mean(as.vector(bias.ill2.array2$se.exprs[[l]]), na.rm=TRUE)
+# avevar.med.array2[l] <- mean(as.vector(bias.med.array2$se.exprs[[l]]), na.rm=TRUE)
+ avevar.mean.array2[l] <- mean(as.vector(bias.mean.array2$se.exprs[[l]]), na.rm=TRUE)
+ avevar.win.array2[l] <- mean(as.vector(bias.win.array2$se.exprs[[l]]), na.rm=TRUE)
+ avevar.trim.array2[l] <- mean(as.vector(bias.trim.array2$se.exprs[[l]]), na.rm=TRUE)
+ avevar.rmx.array2[l] <- mean(as.vector(bias.rmx.array2$se.exprs[[l]]), na.rm=TRUE)
+}
+
+lwd <- 2
+sel <- 1:9
+x <- per.out[sel]*100
+myCol <- brewer.pal(5, "Set1")
+
+pdf("Figure2.pdf", width=8, height=5)
+par(mfrow=c(1,2))
+par(mar=c(3,4,1.5,0.15), oma=c(1,0,0,0))
+plot(x, avebias.ill.array2[sel], lwd=lwd, type="l", xlab="", ylab="average bias", main="A", ylim=c(0,25000), col = myCol[1])
+points(x, avebias.med.array2[sel], type="l", col=myCol[2], lwd=lwd)
+points(x, avebias.trim.array2[sel], type="l", col=myCol[3], lwd=lwd)
+points(x, avebias.mean.array2[sel], type="l", col=myCol[4], lwd=lwd)
+points(x, avebias.rmx.array2[sel], type="l", col=myCol[5], lty = 2, lwd=lwd)
+legend("topleft",legend=c("Illumina", "median", "trimmed mean", "mean", "rmx"), col=myCol, lwd=2)
+#plot(x, numout.ill.array2[sel], lwd=lwd, type="l", xlab="", ylab="% observations removed by summary method", main="(b)", ylim=c(0,40))
+#points(x, numout.med.array2[sel], type="l", col=2, lwd=lwd)
+#points(x, numout.trim.array2[sel], type="l", col=3, lwd=lwd)
+#points(x, numout.mean.array2[sel], type="l", col=4, lwd=lwd)
+#abline(0,1,col="gray", lty=2)
+plot(x, log2(avevar.ill.array2)[sel], lwd=lwd, type="l", xlab="", ylab=expression(log[2]("average variance")), main="B", ylim=c(0,30),col = myCol[1])
+points(x, log2(avevar.trim.array2)[sel], type="l", col=myCol[3], lwd=lwd)
+points(x, log2(avevar.mean.array2)[sel], type="l", col=myCol[4], lwd=lwd)
+points(x, log2(avevar.rmx.array2)[sel], type="l", col=myCol[5], lwd=lwd)
+mtext("% outliers simulated", side=1, outer=TRUE)
+dev.off()
+
+
+###############################################################################################################################################################
+###############################################################################################################################################################
+## below this point not yet adapted to use of robloxbioc
+###############################################################################################################################################################
+###############################################################################################################################################################
+
+
+## log-scale
+## Saturated
+ill.log.summ.array2 <- list(NoBeads=list(), exprs=list(), se.exprs=list())
+ill.log.summ2.array2 <- med.log.summ.array2 <- mean.log.summ.array2 <- trim.log.summ.array2 <- win.log.summ.array2 <- ill.log.summ.array2
+ill.log.rmx.array2 <- ill.log.summ.array2
+
+set.seed(06092007)
+for(l in 1:length(per.out)) {
+ cat(l, "\n")
+ simdataarray2 <- copyBeadLevelList(bld.sharpen.array2)
+ if(per.out[l]!=0) {
+ for(i in 1:nchips) {
+ nsamp <- round(per.out[l]*numBeads(simdataarray2, array=i),0)
+ ind <- sample(1:numBeads(simdataarray2, array=i), nsamp)
+ simdataarray2 at beadData[[arraynms[i]]]$G[ind] <- 2^16
+ }
+ }
+ tmp.ill.log <- createBeadSummaryData(simdataarray2, method="illumina", log=TRUE, n=3, imagesPerArray=2)
+ tmp.ill2.log <- createBeadSummaryData(simdataarray2, method="illumina", log=TRUE, n=2, imagesPerArray=2)
+ tmp.med.log <- createBeadSummaryData(simdataarray2, method="median", log=TRUE, imagesPerArray=2)
+ tmp.mean.log <- createBeadSummaryData(simdataarray2, method="mean", log=TRUE, imagesPerArray=2)
+ tmp.trim.log <- createBeadSummaryData(simdataarray2, method="trim", trim=0.1, log=TRUE, imagesPerArray=2)
+ tmp.win.log <- createBeadSummaryData(simdataarray2, method="winsorize", trim=0.1, log=TRUE, imagesPerArray=2)
+ tmp.rmx.log <- robloxbioc(simdataarray2, log = TRUE, imagesPerArray=2, eps = max(per.out[l], 0.05))
+ ill.log.summ.array2$NoBeads[[l]] <- NoBeads(tmp.ill.log)
+ ill.log.summ.array2$exprs[[l]] <- exprs(tmp.ill.log)
+ ill.log.summ.array2$se.exprs[[l]] <- se.exprs(tmp.ill.log)
+ ill.log.summ2.array2$NoBeads[[l]] <- NoBeads(tmp.ill2.log)
+ ill.log.summ2.array2$exprs[[l]] <- exprs(tmp.ill2.log)
+ ill.log.summ2.array2$se.exprs[[l]] <- se.exprs(tmp.ill2.log)
+ med.log.summ.array2$NoBeads[[l]] <- NoBeads(tmp.med.log)
+ med.log.summ.array2$exprs[[l]] <- exprs(tmp.med.log)
+ med.log.summ.array2$se.exprs[[l]] <- se.exprs(tmp.med.log)
+ mean.log.summ.array2$NoBeads[[l]] <- NoBeads(tmp.mean.log)
+ mean.log.summ.array2$exprs[[l]] <- exprs(tmp.mean.log)
+ mean.log.summ.array2$se.exprs[[l]] <- se.exprs(tmp.mean.log)
+ win.log.summ.array2$NoBeads[[l]] <- NoBeads(tmp.win.log)
+ win.log.summ.array2$exprs[[l]] <- exprs(tmp.win.log)
+ win.log.summ.array2$se.exprs[[l]] <- se.exprs(tmp.win.log)
+ trim.log.summ.array2$NoBeads[[l]] <- NoBeads(tmp.trim.log)
+ trim.log.summ.array2$exprs[[l]] <- exprs(tmp.trim.log)
+ trim.log.summ.array2$se.exprs[[l]] <- se.exprs(tmp.trim.log)
+ rmx.log.summ.array2$NoBeads[[l]] <- NoBeads(tmp.rmx.log)
+ rmx.log.summ.array2$exprs[[l]] <- exprs(tmp.rmx.log)
+ rmx.log.summ.array2$se.exprs[[l]] <- se.exprs(tmp.rmx.log)
+ rm(simdataarray2, tmp.ill.log, tmp.med.log, tmp.mean.log, tmp.trim.log, tmp.win.log, tmp.rmx.log)
+}
+save(ill.log.summ.array2, ill.log.summ2.array2, med.log.summ.array2, mean.log.summ.array2,
+ win.log.summ.array2, trim.log.summ.array2, rmx.log.summ.array2, file="sim.summary.bias.log.rda")
+
+bias.ill.log.array2 <- list(NoBeads=list(), exprs=list(), se.exprs=list())
+bias.ill2.log.array2 <- bias.med.log.array2 <- bias.mean.log.array2 <- bias.trim.log.array2 <- bias.win.log.array2 <- bias.ill.log.array2
+bias.rmx.log.array2 <- bias.ill.log.array2
+for(l in 1:length(per.out)) {
+ cat(l, "\n")
+ bias.ill.log.array2$exprs[[l]] <- ill.log.summ.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.ill.log.array2$se.exprs[[l]] <- ill.log.summ.array2$se.exprs[[l]]^2*ill.log.summ.array2$NoBeads[[l]]
+ bias.ill.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-ill.log.summ.array2$NoBeads[[l]]
+
+ bias.ill2.log.array2$exprs[[l]] <- ill.log.summ2.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.ill2.log.array2$se.exprs[[l]] <- ill.log.summ2.array2$se.exprs[[l]]^2*ill.log.summ2.array2$NoBeads[[l]]
+ bias.ill2.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-ill.log.summ2.array2$NoBeads[[l]]
+
+ bias.med.log.array2$exprs[[l]] <- med.log.summ.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.med.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-med.log.summ.array2$NoBeads[[l]]
+
+ bias.mean.log.array2$exprs[[l]] <- mean.log.summ.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.mean.log.array2$se.exprs[[l]] <- mean.log.summ.array2$se.exprs[[l]]^2*mean.log.summ.array2$NoBeads[[l]]
+ bias.mean.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-mean.log.summ.array2$NoBeads[[l]]
+
+ bias.trim.log.array2$exprs[[l]] <- trim.log.summ.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.trim.log.array2$se.exprs[[l]] <- trim.log.summ.array2$se.exprs[[l]]^2*trim.log.summ.array2$NoBeads[[l]]
+ bias.trim.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-trim.log.summ.array2$NoBeads[[l]]
+
+ bias.win.log.array2$exprs[[l]] <- win.log.summ.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.win.log.array2$se.exprs[[l]] <- win.log.summ.array2$se.exprs[[l]]^2*win.log.summ.array2$NoBeads[[l]]
+ bias.win.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-win.log.summ.array2$NoBeads[[l]]
+
+ bias.rmx.log.array2$exprs[[l]] <- rmx.log.summ.array2$exprs[[l]]-mean.log.summ.array2$exprs[[1]]
+ bias.rmx.log.array2$se.exprs[[l]] <- rmx.log.summ.array2$se.exprs[[l]]^2*rmx.log.summ.array2$NoBeads[[l]]
+ bias.rmx.log.array2$NoBeads[[l]] <- mean.log.summ.array2$NoBeads[[1]]-rmx.log.summ.array2$NoBeads[[l]]
+}
+
+avebias.ill.log.array2 <- avebias.ill2.log.array2 <- avebias.mean.log.array2 <- avebias.med.log.array2 <- avebias.trim.log.array2 <- avebias.win.log.array2 <- NULL
+avebias.rmx.log.array2 <- NULL
+for(l in 1:length(per.out)) {
+ avebias.ill.log.array2[l] <- mean(as.vector(bias.ill.log.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.ill2.log.array2[l] <- mean(as.vector(bias.ill2.log.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.med.log.array2[l] <- mean(as.vector(bias.med.log.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.mean.log.array2[l] <- mean(as.vector(bias.mean.log.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.win.log.array2[l] <- mean(as.vector(bias.win.log.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.trim.log.array2[l] <- mean(as.vector(bias.trim.log.array2$exprs[[l]]), na.rm=TRUE)
+ avebias.rmx.log.array2[l] <- mean(as.vector(bias.rmx.log.array2$exprs[[l]]), na.rm=TRUE)
+}
+
+numout.ill.log.array2 <- numout.ill2.log.array2 <- numout.mean.log.array2 <- numout.med.log.array2 <- numout.trim.log.array2 <- numout.win.log.array2 <- NULL
+numout.rmx.log.array2 <- NULL
+for(l in 1:length(per.out)) {
+ numout.ill.log.array2[l] <- mean(as.vector(bias.ill.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.ill2.log.array2[l] <- mean(as.vector(bias.ill2.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.med.log.array2[l] <- mean(as.vector(bias.med.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.mean.log.array2[l] <- mean(as.vector(bias.mean.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.win.log.array2[l] <- mean(as.vector(bias.win.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.trim.log.array2[l] <- mean(as.vector(bias.trim.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+ numout.rmx.log.array2[l] <- mean(as.vector(bias.rmx.log.array2$NoBeads[[l]]/mean.log.summ.array2$NoBeads[[1]]), na.rm=TRUE)*100
+}
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 268
More information about the Robast-commits
mailing list