################################################### ### chunk number 1: ################################################### #line 206 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" options(continue=' ') ################################################### ### chunk number 2: ################################################### #line 210 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" library(marray) library(dyebias) library(dyebiasexamples) library(convert) options(stringsAsFactors = FALSE) ################################################### ### chunk number 3: ################################################### #line 223 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" datadir <- system.file("extradata", package="dyebiasexamples") datadir ################################################### ### chunk number 4: ################################################### #line 231 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" sample.file <- file.path(datadir, "Tuteja-samples.txt") targets <- read.marrayInfo(sample.file) summary(targets) ################################################### ### chunk number 5: ################################################### #line 247 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" layout <- read.marrayLayout(fname=NULL, ngr=12, ngc=4, nsr=22, nsc=19) n.spots <- maNspots(layout) summary(layout) ################################################### ### chunk number 6: ################################################### #line 264 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" first.file <- file.path(datadir, maInfo(targets)$filename[1]) first.file gnames <- read.marrayInfo(first.file, info.id=c(7,8,10), labels=10, skip=0) names(maInfo(gnames)) <- c("control.type", "reporter.group", "reporterId") maInfo(gnames)$reporterId <- as.character(maInfo(gnames)$reporterId) #### Following is not really needed (it makes the ### ``"Controls are generated from an arbitaray columns\n"'' message go way controls <- maInfo(gnames)$reporter.group controls [ controls == "Experimental"] <- "gene" maControls(layout) <- as.factor(controls) summary(gnames) ################################################### ### chunk number 7: ################################################### #line 292 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" data.raw <- read.GenePix(fnames = maInfo(targets)$filename, path = datadir, name.Gf = "GenePix:F532 Median", name.Gb ="GenePix:B532 Median", name.Rf = "GenePix:F635 Median", name.Rb = "GenePix:B635 Median", name.W = "GenePix:Flags", layout = layout, gnames = gnames, targets = targets, skip=0, sep = "\t", quote = '"', DEBUG=TRUE) ################################################### ### chunk number 8: ################################################### #line 314 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" maW(data.raw)[ maW(data.raw) == 0] <- 1 maW(data.raw)[ maW(data.raw) < 0] <- 0 maLayout(data.raw) <- layout summary(data.raw) ################################################### ### chunk number 9: ################################################### #line 329 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" normalized.data.file <- file.path(datadir, "E-MTAB-32-processed-data-1641029599.txt") data <- read.table(normalized.data.file, sep="\t", as.is=T, header=T, skip=1) names(data) <- c("reporterId", "X13685041", "X13685040" , "X13685153" ,"X13685151" ,"X13685042") ## get rid of overly long gene names: data$reporterId <- sub(pattern="ebi.ac.uk:MIAMExpress:Reporter:A-MEXP-1165\\.(P[0-9]+) .*", replacement="\\1", x=data$reporterId, perl=T, fixed=F) data$reporterId <- sub(pattern="ebi.ac.uk:MIAMExpress:Reporter:A-MEXP-1165\\.Blank.*", replacement="Blank", x=data$reporterId, perl=T, fixed=F) data.norm <- as(data.raw,"marrayNorm") ## now replace the actual data: m <- as.matrix(data[,c(2:6)]) n <- as.numeric(m) dim(n) <- dim(m) maM(data.norm) <- n ################################################### ### chunk number 10: ################################################### #line 368 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" maM(data.norm)[,c(3,4)] <- -1 *maM(data.norm)[,c(3,4)] ################################################### ### chunk number 11: ################################################### #line 377 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" maA(data.norm) <- log2( maRf(data.raw) * maGf(data.raw) ) / 2 ################################################### ### chunk number 12: ################################################### #line 404 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" iGSDBs.estimated <- dyebias.estimate.iGSDBs(data.norm, is.balanced=FALSE, reference="input", verbose=TRUE) summary(iGSDBs.estimated) ################################################### ### chunk number 13: ################################################### #line 425 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" hist(iGSDBs.estimated$dyebias, breaks=50) ################################################### ### chunk number 14: ################################################### #line 447 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" estimator.subset <- (maInfo(maGnames(data.norm))$reporter.group=="Experimental") summary(estimator.subset) ################################################### ### chunk number 15: ################################################### #line 474 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" application.subset <- ((maW(data.norm)==1) & dyebias.application.subset(data.raw=data.raw, min.SNR=1.5, maxA=15, use.background=TRUE) ) summary(as.vector(application.subset)) ################################################### ### chunk number 16: ################################################### #line 489 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" correction <- dyebias.apply.correction(data.norm=data.norm, iGSDBs=iGSDBs.estimated, estimator.subset=estimator.subset, application.subset=application.subset, verbose=TRUE) ################################################### ### chunk number 17: ################################################### #line 507 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" correction$summary[,c("slide", "avg.correction", "var.ratio", "reduction.perc", "p.value")] ################################################### ### chunk number 18: ################################################### #line 543 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" layout(matrix(1:2, nrow=1,ncol=2)) dyebias.rgplot(data=data.norm, slide=3, iGSDBs=iGSDBs.estimated, application.subset=application.subset, main="uncorrected", cex=0.2, cex.lab=0.8, cex.main=0.8, output=NULL) dyebias.rgplot(data=correction$data.corrected, slide=3, iGSDBs=iGSDBs.estimated, application.subset=application.subset, main="corrected", cex=0.2, cex.lab=0.8, cex.main=0.8, output=NULL) ################################################### ### chunk number 19: ################################################### #line 584 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" layout(matrix(1:2, nrow=1,ncol=2)) order <- dyebias.boxplot(data=data.norm, iGSDBs=iGSDBs.estimated, order=NULL, ylim=c(-1,1), application.subset=application.subset, main="uncorrected", cex=0.2, cex.lab=0.8, cex.main=0.8, output=NULL ) order <- dyebias.boxplot(data=correction$data.corrected, iGSDBs=iGSDBs.estimated, order=order, ylim=c(-1,1), application.subset=application.subset, main="corrected", cex=0.2, cex.lab=0.8, cex.main=0.8, output=NULL ) ################################################### ### chunk number 20: eval=FALSE ################################################### ## #line 660 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## library(GEOquery) ## library(marray) ## library(limma) ## library(dyebias) ## library(dyebiasexamples) ## ## options(stringsAsFactors = FALSE) ## ## ## gse.id <- "GSE9318" ## ## dir <- system.file("data", package = "dyebias") ## dir <- getwd() ## ## if (! file.exists(dir)) { ## stop(sprintf("could not find directory '%s' (to write GSE9318 data to/from)", dir)) ## } ## ## file.RData <- sprintf("%s/%s.RData", dir,gse.id) ## ## if( file.exists(file.RData) ) { ## cat(sprintf("Loading existing data from %s\n", file.RData)) ## load(file.RData) ## } else { ## if ( interactive()) { ## if (readline(prompt=sprintf("Could not find file %s.\nDo you want me to download the data from GEO?\nThis will take a while [y/N]: ", file.RData)) != "y" ) { ## stop("Exiting") ## } ## } ## cat("Downloading .soft file from GEO ...\n") ## file.soft <- getGEOfile(gse.id, destdir=dir, amount="full") ## cat("done\nParsing .soft file ...\n") ## gse <- getGEO(filename=file.soft) ## cat(sprintf("done\nSaving to %s ... ", file.RData)) ## save.image(file.RData) ## cat("done\n") ## } ## ################################################### ### chunk number 21: eval=FALSE ################################################### ## #line 727 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## ## following slides are used to obtain iGSDB estimates: ## slide.ids.est <- c( "GSM237352", "GSM237353", "GSM237354", "GSM237355" ) ## ## ## column names: ## cy3.name <- "label_ch1" ## cy5.name <- "label_ch2" ## ## ## function to recognize 'genes': ## gene.selector <- function(table){grep("^A_75_", as.character(table[["ProbeName"]]))} ## ## reporterid.name <- "ProbeName" ## M.name <- "VALUE" #normalized ## Gf.name <- "Cy3" #raw ## Gb.name <- "Cy3_Background" ## Rf.name <- "Cy5" #raw ## Rb.name <- "Cy3_Background" # yes! (error in the data: Cy5_Background == Cy5 ...) ## ## ## convert the raw data ... ## data.raw.est <- ## dyebias.geo2marray(gse=gse, ## slide.ids=slide.ids.est, ## type="raw", ## gene.selector=gene.selector, ## reporterid.name=reporterid.name, ## cy3.name=cy3.name, ## cy5.name=cy5.name, ## Rf.name=Rf.name, ## Gf.name=Gf.name, ## Rb.name=Rb.name, ## Gb.name=Gb.name, ## ) ## ## ## ... and the normalized data: ## data.norm.est <- ## dyebias.geo2marray(gse=gse, ## slide.ids=slide.ids.est, ## type="norm", ## gene.selector=gene.selector, ## reporterid.name=reporterid.name, ## cy3.name=cy3.name, ## cy5.name=cy5.name, ## M.name=M.name ## ) ## ## ## maA was not set; do that here: ## maA(data.norm.est) <- log2( maRf(data.raw.est) * maGf(data.raw.est) ) / 2 ## ## ## unswap the swapped dye-swaps: ## maM(data.norm.est)[, c(2,4)] <- - maM(data.norm.est)[, c(2,4)] ## ################################################### ### chunk number 22: eval=FALSE ################################################### ## #line 784 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## maInfo(maTargets(data.norm.est)) ## ################################################### ### chunk number 23: eval=FALSE ################################################### ## #line 797 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## info <- maInfo(maTargets(data.norm.est)) ## ## info$Cy3 <- c("wt", "wt IP", "td", "td IP" ) ## info$Cy5 <- c("wt IP", "wt", "td IP", "td") ## ## maInfo(maTargets(data.norm.est)) <- info ## ## references <- c("wt", "td") ## ## ### The estimation would then be run as: ## ### ## ### iGSDBs.estimated <- dyebias.estimate.iGSDBs(data.norm.est, is.balanced=FALSE, ## ### reference=references, verbose=TRUE) ## ################################################### ### chunk number 24: eval=FALSE ################################################### ## #line 828 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## iGSDBs.estimated <- dyebias.estimate.iGSDBs(data.norm.est, is.balanced=TRUE, verbose=TRUE) ## summary(iGSDBs.estimated) ## ################################################### ### chunk number 25: eval=FALSE ################################################### ## #line 842 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## subset <- sapply(GSMList(gse), function(x){Meta(x)$platform_id}) == "GPL5991" ## slide.ids.all <- sapply(GSMList(gse), function(x){Meta(x)$geo_accession})[subset] ## ## data.raw.all <- ## dyebias.geo2marray(gse=gse, ## slide.ids=slide.ids.all, ## type="raw", ## gene.selector=gene.selector, ## reporterid.name=reporterid.name, ## cy3.name=cy3.name, ## cy5.name=cy5.name, ## Rf.name=Rf.name, ## Gf.name=Gf.name, ## Rb.name=Rb.name, ## Gb.name=Gb.name, ## ) ## ## data.norm.all <- ## dyebias.geo2marray(gse=gse, ## slide.ids=slide.ids.all, ## type="norm", ## gene.selector=gene.selector, ## reporterid.name=reporterid.name, ## cy3.name=cy3.name, ## cy5.name=cy5.name, ## M.name=M.name ## ) ## ## maA(data.norm.all) <- log2( maRf(data.raw.all) * maGf(data.raw.all) ) / 2 ## ## swap <- c(1,2,3,5, 7,8,9) ## maM(data.norm.all)[, swap] <- - maM(data.norm.all)[, swap] ## ################################################### ### chunk number 26: eval=FALSE ################################################### ## #line 889 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## application.subset <- ( (maW(data.norm.all) == 1) & ## dyebias.application.subset(data.raw=data.raw.all, ## min.SNR=1.5, ## maxA=15, ## use.background=TRUE) ## ) ## ## correction <- dyebias.apply.correction(data.norm=data.norm.all, ## iGSDBs=iGSDBs.estimated, ## estimator.subset=T, ## application.subset=application.subset, ## verbose=FALSE) ## ## correction$summary[,c("slide", "avg.correction", "var.ratio", "reduction.perc", "p.value")] ## ################################################### ### chunk number 27: eval=FALSE ################################################### ## #line 921 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## layout(matrix(1:2, nrow=1,ncol=2)) ## ## dyebias.maplot(data=data.norm.all, ## slide=6, ## iGSDBs=iGSDBs.estimated, ## application.subset=application.subset, ## main="uncorrected", ## cex=0.2, ## cex.lab=0.8, ## cex.main=0.8, ## output=NULL) ## ## dyebias.maplot(data=correction$data.corrected, ## slide=6, ## iGSDBs=iGSDBs.estimated, ## application.subset=application.subset, ## main="corrected", ## cex=0.2, ## cex.lab=0.8, ## cex.main=0.8, ## output=NULL) ################################################### ### chunk number 28: eval=FALSE ################################################### ## #line 974 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" ## ## layout(matrix(1:2, nrow=1,ncol=2)) ## ## order <- dyebias.trendplot(data=data.norm.all, ## iGSDBs=iGSDBs.estimated, ## application.subset=application.subset, ## order=NULL, ## lwd=0.1, ## ylim=c(-0.5,0.5), ## output=NULL, ## main="uncorrected") ## ## order <- dyebias.trendplot(data=correction$data.corrected, ## iGSDBs=iGSDBs.estimated, ## application.subset=application.subset, ## order=order, ## lwd=0.1, ## ylim=c(-0.5,0.5), ## output=NULL, ## main="corrected") ## ################################################### ### chunk number 29: ################################################### #line 1044 "vignettes/dyebias/inst/doc/dyebias-vignette.Rnw" toLatex(sessionInfo())