############################################################################## ### ### ### file: twoColorAnalysisCode.R ### ### Author: B. M. Bolstad ### ### Aim: Implement pre-processing analysis for two color microarray dataset ### ############################################################################# ### ### ### Use code is designed to work with the BioC 1.8 release set of packages ### on R-2.3.x. it may or may not work on later versions. ### ### library(limma) library(marray) library(GEOquery) library(arrayQuality) library(convert) ############################################################################## ### ### ### Initial Code for retrieving data from GEO and setting up the basic ### structures holding the raw data. ### ### ### ############################################################################# #### #### This code should be used for retrieving the data automatically from #### the GEO. Depending on the speed of your connection this may take quite awhile #### ### download information about the chip design gpl <- getGEO("GPL1089") ### download the sample data. Each is about 4 megabytes big gsm.1 <- getGEO("GSM24215") gsm.2 <- getGEO("GSM24216") gsm.3 <- getGEO("GSM24217") gsm.4 <- getGEO("GSM24218") gsm.5 <- getGEO("GSM24219") gsm.6 <- getGEO("GSM24220") gsm.7 <- getGEO("GSM24221") gsm.8 <- getGEO("GSM24224") gsm.9 <- getGEO("GSM24225") gsm.10 <- getGEO("GSM24229") gsm.11 <- getGEO("GSM24230") gsm.12 <- getGEO("GSM24231") gsm.13 <- getGEO("GSM24232") gsm.14 <- getGEO("GSM24233") gsm.15 <- getGEO("GSM24234") ### Save the raw sample data to tab delimited text files in case we need them later. ### Once you have saved them you will not need to re-download them if you intend to ### repeat the analysis write.table(Table(gsm.1),quote=FALSE,sep="\t",file=paste(Meta(gsm.1)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.1)[,1]) write.table(Table(gsm.2),quote=FALSE,sep="\t",file=paste(Meta(gsm.2)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.2)[,1]) write.table(Table(gsm.3),quote=FALSE,sep="\t",file=paste(Meta(gsm.3)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.3)[,1]) write.table(Table(gsm.4),quote=FALSE,sep="\t",file=paste(Meta(gsm.4)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.4)[,1]) write.table(Table(gsm.5),quote=FALSE,sep="\t",file=paste(Meta(gsm.5)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.5)[,1]) write.table(Table(gsm.6),quote=FALSE,sep="\t",file=paste(Meta(gsm.6)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.6)[,1]) write.table(Table(gsm.7),quote=FALSE,sep="\t",file=paste(Meta(gsm.7)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.7)[,1]) write.table(Table(gsm.8),quote=FALSE,sep="\t",file=paste(Meta(gsm.8)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.8)[,1]) write.table(Table(gsm.9),quote=FALSE,sep="\t",file=paste(Meta(gsm.9)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.9)[,1]) write.table(Table(gsm.10),quote=FALSE,sep="\t",file=paste(Meta(gsm.10)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.10)[,1]) write.table(Table(gsm.11),quote=FALSE,sep="\t",file=paste(Meta(gsm.11)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.11)[,1]) write.table(Table(gsm.12),quote=FALSE,sep="\t",file=paste(Meta(gsm.12)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.12)[,1]) write.table(Table(gsm.13),quote=FALSE,sep="\t",file=paste(Meta(gsm.13)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.13)[,1]) write.table(Table(gsm.14),quote=FALSE,sep="\t",file=paste(Meta(gsm.14)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.14)[,1]) write.table(Table(gsm.15),quote=FALSE,sep="\t",file=paste(Meta(gsm.15)$geo_accession,".raw",sep=""),row.names=FALSE,col.names=Columns(gsm.15)[,1]) ### ### Writing out information about the chip layout. ### write.table(Table(gpl),quote=FALSE,sep="\t",file="gpl1089chiplayout.txt",row.names=FALSE) ### Meta data giving information about each slide sample.info <- rbind(do.call("c",Meta(gsm.1)), do.call("c",Meta(gsm.2)), do.call("c",Meta(gsm.3)), do.call("c",Meta(gsm.4)), do.call("c",Meta(gsm.5)), do.call("c",Meta(gsm.6)), do.call("c",Meta(gsm.7)), do.call("c",Meta(gsm.8)), do.call("c",Meta(gsm.9)), do.call("c",Meta(gsm.10)), do.call("c",Meta(gsm.11)), do.call("c",Meta(gsm.12)), do.call("c",Meta(gsm.13)), do.call("c",Meta(gsm.14)), do.call("c",Meta(gsm.15))) ### ### Write it out in case we need it later. ### write.table(sample.info,quote=FALSE,sep="\t",file="sampleinfo.txt",row.names=FALSE) ### ### This next section is setting up the initial data structures we will use. ### ### Although it is a little confusing there are two frameworks for ### storing two color microarray data. One for limma and one for marray ### ### From limma: ### RGList - for storing raw R, G, Rb, Gb and associated information ### MAList - for storing processed M and A values and associated information ### ### From marray: ### marrayRaw - for storing raw R, G, Rb, Gb and associated information ### marrayNorm - for storing processed M and A values and associated information ### ### Because each package provides slightly different functionality ### we will have to use both frameworks. Luckily the convert package ### allows us to convert an object from one framework to the other. ### ### to prevent confusion in the following code we will pre-prend all ### variable names with the object type. eg RGList.*, MAList.*, marrayRaw.* ### and marrayNorm.* ### ### set up the spot types ### This depends on the particular slide design. ### this particular array has some empty spots as well as negative ### and positive controls ### spot.types <- rbind( c("probes" , "*" , "*" , "black"), c("Empty" , "Empty" , "*" , "yellow"), c("Empty" , "EMPTY" , "*" , "yellow"), c("Negative", "AT*" , "*" , "orange"), c("Positive", "M200009348", "*" , "brown"), c("Positive", "M200003425", "*" , "brown") ) spot.types <- as.data.frame(spot.types) colnames(spot.types) <- c("SpotType","Operon_ID","ID","Color") ### read the data into a structure that can hold the raw R,G, Rb, Gb signal intensities RGList.rawdata <- read.maimages(dir(pattern=("GSM.*raw")),columns=list(R="F635 Mean", G = "F532 Mean",Rb="B635 Median",Gb="B532 Median")) ### append additional information to the data structure which describes the array design ### that was used and the samples. RGList.rawdata$genes <- Table(gpl) RGList.rawdata$printer <- getLayout(RGList.rawdata$genes) RGList.rawdata$targets <- as.data.frame(sample.info) RGList.rawdata$genes$Status <- controlStatus(spot.types,RGList.rawdata) ### convert from an RGList to marrayRaw marrayRaw.raw <- as(RGList.rawdata,"marrayRaw") ############################################################################## ### ### ### Basic quality analysis and exploration of the raw data ### ### ### ############################################################################# ### ### ### Our first examination is to look at the spatial distribution of the ### Red and Green background intensities ### ### Note that the logs are to get a logarithmic color-space ### the zlim are to make the plots strictly comparable ### ### par(mfrow=c(2,1)) imageplot(log2(RGList.rawdata$Rb[,1]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.1)$geo_accession),zlim=c(6.5,10)) imageplot(log2(RGList.rawdata$Gb[,1]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.1)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,2]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.2)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,2]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.2)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,3]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.3)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,3]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.3)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,4]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.4)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,4]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.4)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,5]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.5)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,5]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.5)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,6]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.6)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,6]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.6)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,7]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.7)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,7]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.7)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,8]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.8)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,8]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.8)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,9]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.9)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,9]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.9)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,10]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.10)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,10]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.10)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,11]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.11)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,11]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.11)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,12]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.12)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,12]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.12)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,13]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.13)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,13]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.13)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,14]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.14)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,14]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.14)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Rb[,15]+1), RGList.rawdata$printer, low="white", high="red",main=paste("Red channel background for",Meta(gsm.15)$geo_accession),zlim=c(6.5,11)) imageplot(log2(RGList.rawdata$Gb[,15]+1), RGList.rawdata$printer, low="white", high="darkgreen",main=paste("Green channel background for",Meta(gsm.15)$geo_accession),zlim=c(6.5,11)) ### ### ### lets look at signal noise ratios in the red and green channels ### ### These plots also have smoothed histograms for the control spots. ### Clearly we would like to see the positive and negative controls ### well separated. and the probes should have higher signal to noise ### ratio than the empty spots. ### ### colcode <- setCtlCol(marrayRaw.raw[,1]) par(mfrow=c(1,2)) qpS2N(marrayRaw.raw[,1],channel="red",colcode) qpS2N(marrayRaw.raw[,1],channel="green",colcode) qpS2N(marrayRaw.raw[,2],channel="red",colcode) qpS2N(marrayRaw.raw[,2],channel="green",colcode) qpS2N(marrayRaw.raw[,3],channel="red",colcode) qpS2N(marrayRaw.raw[,3],channel="green",colcode) qpS2N(marrayRaw.raw[,4],channel="red",colcode) qpS2N(marrayRaw.raw[,4],channel="green",colcode) qpS2N(marrayRaw.raw[,5],channel="red",colcode) qpS2N(marrayRaw.raw[,5],channel="green",colcode) qpS2N(marrayRaw.raw[,6],channel="red",colcode) qpS2N(marrayRaw.raw[,6],channel="green",colcode) qpS2N(marrayRaw.raw[,7],channel="red",colcode) qpS2N(marrayRaw.raw[,7],channel="green",colcode) qpS2N(marrayRaw.raw[,8],channel="red",colcode) qpS2N(marrayRaw.raw[,8],channel="green",colcode) qpS2N(marrayRaw.raw[,9],channel="red",colcode) qpS2N(marrayRaw.raw[,9],channel="green",colcode) qpS2N(marrayRaw.raw[,10],channel="red",colcode) qpS2N(marrayRaw.raw[,10],channel="green",colcode) qpS2N(marrayRaw.raw[,11],channel="red",colcode) qpS2N(marrayRaw.raw[,11],channel="green",colcode) qpS2N(marrayRaw.raw[,12],channel="red",colcode) qpS2N(marrayRaw.raw[,12],channel="green",colcode) qpS2N(marrayRaw.raw[,13],channel="red",colcode) qpS2N(marrayRaw.raw[,13],channel="green",colcode) qpS2N(marrayRaw.raw[,14],channel="red",colcode) qpS2N(marrayRaw.raw[,14],channel="green",colcode) qpS2N(marrayRaw.raw[,15],channel="red",colcode) qpS2N(marrayRaw.raw[,15],channel="green",colcode) ### this function will produce numerical values for mean and var of S2N ### ratios in the Red and Green channels Calc.SN <- function(object){ S2N.R <- log(maRf(object), 2) - log(maRb(object),2) S2N.G <- log(maGf(object), 2) - log(maGb(object),2) SN <- cbind(apply(S2N.R,2,mean.na), apply(S2N.R,2,var.na), apply(S2N.G,2,mean.na), apply(S2N.G,2,var.na)) SN <- data.frame(SN) colnames(SN) <- c("mean SN.Red", "var SN.Red","mean SN.Green", "var SN.Green") SN } Calc.SN(marrayRaw.raw) ### ### MAplots for each slide with-out and with local-background correction ### ### The MAplots where local-background correction has taken place are ### very noisy in the low intensities. For the analysis of this dataset ### we will choose to do no background correction ### ### RGList.bgc <- backgroundCorrect(RGList.rawdata,"none") marrayRaw.bgc <-as(RGList.bgc,"marrayRaw") ### ### Do the plots side by side. With no background on the left ### and local-background on the right. ### ### We use the same ylim to make the plots comparable. ### ### Note that by default plotMA does local background correction ### which is why plotMA() on an RGList.rawdata is noisy ### we created the RGList.bgc object by carrying out no background correction par(mfrow=c(1,2)) plotMA(RGList.bgc,array=1,ylim=c(-5,5)) plotMA(RGList.rawdata,array=1,ylim=c(-5,5)) plotMA(RGList.bgc,array=2,ylim=c(-5,5)) plotMA(RGList.rawdata,array=2,ylim=c(-5,5)) plotMA(RGList.bgc,array=3,ylim=c(-5,5)) plotMA(RGList.rawdata,array=3,ylim=c(-5,5)) plotMA(RGList.bgc,array=4,ylim=c(-5,5)) plotMA(RGList.rawdata,array=4,ylim=c(-5,5)) plotMA(RGList.bgc,array=5,ylim=c(-5,5)) plotMA(RGList.rawdata,array=5,ylim=c(-5,5)) plotMA(RGList.bgc,array=6,ylim=c(-5,5)) plotMA(RGList.rawdata,array=6,ylim=c(-5,5)) plotMA(RGList.bgc,array=7,ylim=c(-5,5)) plotMA(RGList.rawdata,array=7,ylim=c(-5,5)) plotMA(RGList.bgc,array=8,ylim=c(-5,5)) plotMA(RGList.rawdata,array=8,ylim=c(-5,5)) plotMA(RGList.bgc,array=9,ylim=c(-5,5)) plotMA(RGList.rawdata,array=9,ylim=c(-5,5)) plotMA(RGList.bgc,array=10,ylim=c(-5,5)) plotMA(RGList.rawdata,array=10,ylim=c(-5,5)) plotMA(RGList.bgc,array=11,ylim=c(-5,5)) plotMA(RGList.rawdata,array=11,ylim=c(-5,5)) plotMA(RGList.bgc,array=12,ylim=c(-5,5)) plotMA(RGList.rawdata,array=12,ylim=c(-5,5)) plotMA(RGList.bgc,array=13,ylim=c(-5,5)) plotMA(RGList.rawdata,array=13,ylim=c(-5,5)) plotMA(RGList.bgc,array=14,ylim=c(-5,5)) plotMA(RGList.rawdata,array=14,ylim=c(-5,5)) plotMA(RGList.bgc,array=15,ylim=c(-5,5)) plotMA(RGList.rawdata,array=15,ylim=c(-5,5)) ### ### MAplots with print-tip (grid) specific loess smoothers) ### for the non-background corrected data ### par(mfrow=c(1,1)) plot(marrayRaw.bgc[,1],main=Meta(gsm.1)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,2],main=Meta(gsm.2)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,3],main=Meta(gsm.3)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,4],main=Meta(gsm.4)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,5],main=Meta(gsm.5)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,6],main=Meta(gsm.6)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,7],main=Meta(gsm.7)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,8],main=Meta(gsm.8)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,9],main=Meta(gsm.9)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,10],main=Meta(gsm.10)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,11],main=Meta(gsm.11)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,12],main=Meta(gsm.12)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,13],main=Meta(gsm.13)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,14],main=Meta(gsm.14)$geo_accession,ylim=c(-5,5.5)) plot(marrayRaw.bgc[,15],main=Meta(gsm.15)$geo_accession,ylim=c(-5,5.5)) ### ### Boxplot M values for each grid group. ### for the non-background corrected data ### boxplot(marrayRaw.bgc[,1],main=Meta(gsm.1)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,2],main=Meta(gsm.2)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,3],main=Meta(gsm.3)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,4],main=Meta(gsm.4)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,5],main=Meta(gsm.5)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,6],main=Meta(gsm.6)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,7],main=Meta(gsm.7)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,8],main=Meta(gsm.8)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,9],main=Meta(gsm.9)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,10],main=Meta(gsm.10)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,11],main=Meta(gsm.11)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,12],main=Meta(gsm.12)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,13],main=Meta(gsm.13)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,14],main=Meta(gsm.14)$geo_accession,ylim=c(-2,2)) boxplot(marrayRaw.bgc[,15],main=Meta(gsm.15)$geo_accession,ylim=c(-2,2)) ### ### Comparing M values across arrays ### par(las=2) boxplot(marrayRaw.bgc,main="M across arrays") ### ### Now we will create spatial plots of the M (and A) values ### to make the plots nicer we will image the rank of ### the M value. Uneven spatial distributions indicate ### potential problems. ### MAList.nonorm.bgc <- normalizeWithinArrays(RGList.bgc,method="none") tmpNorm <- as(MAList.nonorm.bgc, "marrayNorm") tmpNorm@maM <- apply(maM(tmpNorm),2,rank) tmpNorm@maA <- apply(maA(tmpNorm),2,rank) col.forA <- maPalette(low="white",high="black",k=100) par(mfrow=c(1,2)) image(tmpNorm[,1], xvar = "maM",main=paste(Meta(gsm.1)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,1], xvar = "maA",main=paste(Meta(gsm.1)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,2], xvar = "maM",main=paste(Meta(gsm.2)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,2], xvar = "maA",main=paste(Meta(gsm.2)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,3], xvar = "maM",main=paste(Meta(gsm.3)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,3], xvar = "maA",main=paste(Meta(gsm.3)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,4], xvar = "maM",main=paste(Meta(gsm.4)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,4], xvar = "maA",main=paste(Meta(gsm.4)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,5], xvar = "maM",main=paste(Meta(gsm.5)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,5], xvar = "maA",main=paste(Meta(gsm.5)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,6], xvar = "maM",main=paste(Meta(gsm.6)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,6], xvar = "maA",main=paste(Meta(gsm.6)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,7], xvar = "maM",main=paste(Meta(gsm.7)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,7], xvar = "maA",main=paste(Meta(gsm.7)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,8], xvar = "maM",main=paste(Meta(gsm.8)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,8], xvar = "maA",main=paste(Meta(gsm.8)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,9], xvar = "maM",main=paste(Meta(gsm.9)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,9], xvar = "maA",main=paste(Meta(gsm.9)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,10], xvar = "maM",main=paste(Meta(gsm.10)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,10], xvar = "maA",main=paste(Meta(gsm.10)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,11], xvar = "maM",main=paste(Meta(gsm.11)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,11], xvar = "maA",main=paste(Meta(gsm.11)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,12], xvar = "maM",main=paste(Meta(gsm.12)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,12], xvar = "maA",main=paste(Meta(gsm.12)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,13], xvar = "maM",main=paste(Meta(gsm.13)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,13], xvar = "maA",main=paste(Meta(gsm.13)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,14], xvar = "maM",main=paste(Meta(gsm.14)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,14], xvar = "maA",main=paste(Meta(gsm.14)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) image(tmpNorm[,15], xvar = "maM",main=paste(Meta(gsm.15)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,15], xvar = "maA",main=paste(Meta(gsm.15)$geo_accession,"Ranked A"),bar=FALSE,col=col.forA) ############################################################################## ### ### ### Within Slide normalization and investigation of its effects ### ### ### ############################################################################# ### ### Normalize the R and G channel intensities using the printtiploess method ### Returns an object with ### MAList.norm.bgc <- normalizeWithinArrays(RGList.bgc,method = "printtiploess") marrayNorm.norm.bgc <-as(MAList.norm.bgc,"marrayNorm") ### ### Examine the MA plots after within slide normalization. ### par(mfrow=c(1,1)) plot(marrayNorm.norm.bgc[,1],main=Meta(gsm.1)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,2],main=Meta(gsm.2)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,3],main=Meta(gsm.3)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,4],main=Meta(gsm.4)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,5],main=Meta(gsm.5)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,6],main=Meta(gsm.6)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,7],main=Meta(gsm.7)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,8],main=Meta(gsm.8)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,9],main=Meta(gsm.9)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,10],main=Meta(gsm.10)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,11],main=Meta(gsm.11)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,12],main=Meta(gsm.12)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,13],main=Meta(gsm.13)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,14],main=Meta(gsm.14)$geo_accession,ylim=c(-5,7)) plot(marrayNorm.norm.bgc[,15],main=Meta(gsm.15)$geo_accession,ylim=c(-5,7)) ### ### Boxplot M values for each grid group after normalization ### boxplot(marrayNorm.norm.bgc[,1],main=Meta(gsm.1)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,2],main=Meta(gsm.2)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,3],main=Meta(gsm.3)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,4],main=Meta(gsm.4)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,5],main=Meta(gsm.5)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,6],main=Meta(gsm.6)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,7],main=Meta(gsm.7)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,8],main=Meta(gsm.8)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,9],main=Meta(gsm.9)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,10],main=Meta(gsm.10)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,11],main=Meta(gsm.11)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,12],main=Meta(gsm.12)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,13],main=Meta(gsm.13)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,14],main=Meta(gsm.14)$geo_accession,ylim=c(-2,2)) boxplot(marrayNorm.norm.bgc[,15],main=Meta(gsm.15)$geo_accession,ylim=c(-2,2)) ### ### Comparing M values across arrays after normalization ### boxplot(marrayNorm.norm.bgc,ylim=c(-1,1),main="M across arrays") ### ### Now we will create spatial plots of the M (and A) values ### to make the plots nicer we will image the rank of ### the M value. Uneven spatial distributions indicate ### potential problems. ### tmpNorm <- marrayNorm.norm.bgc tmpNorm@maM <- apply(maM(tmpNorm),2,rank) par(mfrow=c(1,1)) image(tmpNorm[,1], xvar = "maM",main=paste(Meta(gsm.1)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,2], xvar = "maM",main=paste(Meta(gsm.2)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,3], xvar = "maM",main=paste(Meta(gsm.3)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,4], xvar = "maM",main=paste(Meta(gsm.4)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,5], xvar = "maM",main=paste(Meta(gsm.5)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,6], xvar = "maM",main=paste(Meta(gsm.6)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,7], xvar = "maM",main=paste(Meta(gsm.7)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,8], xvar = "maM",main=paste(Meta(gsm.8)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,9], xvar = "maM",main=paste(Meta(gsm.9)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,10], xvar = "maM",main=paste(Meta(gsm.10)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,11], xvar = "maM",main=paste(Meta(gsm.11)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,12], xvar = "maM",main=paste(Meta(gsm.12)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,13], xvar = "maM",main=paste(Meta(gsm.13)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,14], xvar = "maM",main=paste(Meta(gsm.14)$geo_accession,"Ranked M"),bar=FALSE) image(tmpNorm[,15], xvar = "maM",main=paste(Meta(gsm.15)$geo_accession,"Ranked M"),bar=FALSE) ############################################################################## ### ### ### Between Slide normalization and investigation of its effects ### ### ### ############################################################################# MAList.between.norm.bgc <- normalizeBetweenArrays(MAList.norm.bgc,method = "quantile") marrayNorm.between.norm.bgc <-as(MAList.between.norm.bgc,"marrayNorm") boxplot(marrayNorm.between.norm.bgc,ylim=c(-1,1)) boxplot(apply(maM(marrayRaw.bgc),1,var.na),apply(maM(marrayNorm.norm.bgc),1,var.na),apply(maM(marrayNorm.between.norm.bgc),1,var.na)) ### ### Density plots of the single channel intensities ### ### Before any adjustment ### After within slide adjustment ### After within slide adjustment followed by between slide adjustment. ### par(mfrow=c(1,3)) plotDensities(RGList.bgc) plotDensities(MAList.norm.bgc) plotDensities(MAList.between.norm.bgc)