Lab3b.R
################################################### ### chunk number 1: loadpacks ################################################### library(affy) library(affydata)
################################################### ### chunk number 2: loaddata ################################################### data(Dilution)
################################################### ### chunk number 3: dil ################################################### class(Dilution) slotNames(Dilution) Dilution annotation(Dilution)
################################################### ### chunk number 4: pheno ################################################### phenoData(Dilution) pData(Dilution)
################################################### ### chunk number 5: expr ################################################### e<-exprs(Dilution) nrow(Dilution)*ncol(Dilution) dim(e)
################################################### ### chunk number 6: pmmm ################################################### PM<-pm(Dilution) dim(PM) PM[1:5,]
################################################### ### chunk number 7: AffyID ################################################### gnames<-geneNames(Dilution) length(gnames) gnames[1:5] nrow(e)/length(gnames)
################################################### ### chunk number 8: sub ################################################### dil1<-Dilution[1] class(dil1) dil1 cel1<-Dilution[[1]] class(cel1) cel1
################################################### ### chunk number 9: image1 ################################################### # Log transformation image(Dilution[1])
################################################### ### chunk number 10: image2 ################################################### # No transformation image(cel1)
################################################### ### chunk number 11: boxplot ################################################### boxplot(Dilution,col=c(2,2,3,3))
################################################### ### chunk number 12: hist ################################################### hist(Dilution, type="l", col=c(2,2,3,3), lty=rep(1:2,2), lwd=3)
################################################### ### chunk number 13: normalization ################################################### Dil10 <- normalize(Dilution[1:2]) ##conc. group 10 micrograms Dil20 <- normalize(Dilution[3:4]) ##conc. groun 20 mictograms normDil <- merge(Dil20,Dil10)
################################################### ### chunk number 14: hist2 ################################################### boxplot(normDil,col=c(2,2,3,3))
################################################### ### chunk number 15: methods ################################################### bgcorrect.methods normalize.AffyBatch.methods pmcorrect.methods express.summary.stat.methods
################################################### ### chunk number 16: rma ################################################### rmaDil<-rma(normDil,normalize=FALSE) class(rmaDil)
################################################### ### chunk number 17: cdf ################################################### annotation(Dilution) data(hgu95av2cdf) pnames<-ls(env=hgu95av2cdf) length(gnames) gnames[1:5] get(gnames[1],env=hgu95av2cdf)
################################################### ### chunk number 18: loc ################################################### plocs<-indexProbes(Dilution,which="both") plocs[[1]] pmindex(Dilution,genenames=gnames[1], xy=TRUE) pmindex(Dilution,genenames=gnames[1])
################################################### ### chunk number 19: PMvMM ################################################### plot(mm(Dilution[1]),pm(Dilution[1]),pch=".",log="xy") abline(0,1,col="red")
################################################### ### chunk number 20: ################################################### data(SpikeIn) #loads objects: concentrations and SpikeIn sampleNames(SpikeIn)
################################################### ### chunk number 21: ################################################### pms <- pm(SpikeIn) mms <- mm(SpikeIn)
##pms follow concentration par(mfrow=c(1,2)) concentrations <- matrix(as.numeric(sampleNames(SpikeIn)),20,12,byrow=TRUE) matplot(concentrations,pms,log="xy",main="PM",ylim=c(30,20000)) lines(concentrations[1,],apply(pms,2,mean),lwd=3) ##so do mms, but no as much matplot(concentrations,mms,log="xy",main="MM",ylim=c(30,20000)) lines(concentrations[1,],apply(mms,2,mean),lwd=3)
################################################### ### chunk number 22: ################################################### if(.Platform$OS.type == "unix") { library(mAffyu95A) data(affyprobeids) data(affylocs) }
################################################### ### chunk number 23: ################################################### if( .Platform$OS.type=="unix") { Index <- affylocs[,1]+ affylocs[,2]*ncol(Dilution) + 1 probenames <- probeNames(Dilution) myindex <- unlist(pmindex(Dilution)) probenames <- probenames[match(Index,myindex)] seq <- getprobes(seq(along=Index)) tmp <- sapply(seq,cgcontent) tmp <- matrix(unlist(tmp),nrow=4) atc <- tmp[1,]+tmp[2,] gcc <- tmp[3,]+tmp[4,] }
################################################### ### chunk number 24: ################################################### if( .Platform$OS.type=="unix") { y <- intensity(Dilution)[Index,1] boxplot(split(y, gcc),xlab="GC content", ylab="Iintensity", log="y",las=1, range=0,col="grey") }
################################################### ### chunk number 25: ################################################### if(.Platform$OS.type=="unix") { avggcc <- tapply(gcc, probenames, mean) barplot(table(round(avggcc)),xlab="Average GC content in probeset",las=1) }
################################################### ### chunk number 26: ################################################### if(.Platform$OS.type=="unix") { ##pm - mm pmi <- unlist(pmindex(Dilution)) mmi <- unlist(mmindex(Dilution)) pmIndex <- match(pmi,Index)
y <- intensity(Dilution)[pmi,1]-intensity(Dilution)[mmi,1] x <- gcc[pmIndex] boxplot(split(y,x),xlab="GC content",ylab="PM-MM",las=1, range=0,col="grey",ylim=c(-300,500)) }
################################################### ### chunk number 27: ################################################### if( .Platform$OS.type=="unix") { ##the following command will take a few minutes mas5Dil <- mas5(Dilution,normalize=FALSE)
boxplot(split(exprs(mas5Dil)[names(avggcc),1],round(avggcc)), log="y",range=0,las=1,yaxt="n", xlab="Average GC content in probe set", ylab="MAS 5.0 expression",col="grey")
axis(2,c(.1,10,1000),c("0.01","10","1000"),las=1) }
################################################### ### chunk number 28: ################################################### if( .Platform$OS.type=="unix") {
RMA <- rma(Dilution)
boxplot(split(exprs(rmaDil)[names(avggcc),1],round(avggcc)), log="y",range=0,las=1,yaxt="n", xlab="Average GC content in probe set", ylab="RMA (pm-only) expression",col="grey")
axis(2,c(.1,10,1000),c("0.01","10","1000"),las=1) }
################################################### ### chunk number 29: ################################################### if( .Platform$OS.type=="unix") {
rnadeg <- AffyRNAdeg(Dilution) summaryAffyRNAdeg(rnadeg) plotAffyRNAdeg(rnadeg) }
################################################### ### chunk number 30: ################################################### if( .Platform$OS.type=="unix") {
scale01 <- function(x){r <- range(x); x-r[1]/(r[2]-r[1]) } locs <- tapply(as.double(affylocs[,3]), probenames, scale01) locs <- unlist(locs)
plot(locs, locs, type="n",xlab="5' <-----> 3'\nRelative position", ylab="Relative log intesity", yaxt="n", xaxt="n", ylim=c(-.02,.35)) legend(.1,.3,c("normal RNA","degraded RNA"), col=c("blue","red"), lty=c(2,1)) ##lets just show array 1 i <- 1 o <- sample(length(gcc), 5000) pms <- log2(intensity(Dilution)[Index, i]) res <- tapply(pms, probenames, function(x) x-median(x)) tmp <- loess(unlist(res)~gcc, subset=o, span=4/5, degree=1) oo <- seq(1,length(tmp$x), len=100) lines(sort(tmp$x)[oo], tmp$fitted[order(tmp$x)][oo], col="red") }