# ------------------------------------------
# Example analyzing the 6 arrays in the
# AffySpikeU95Subset data set
# Loading the data
data(AffySpikeU95Subset)
# Defining design and contrast matrix
group<-factor(rep(1:2,each=3))
design<-model.matrix(~group-1)
contrast<-matrix(c(1,-1),1,2)
# Analyzing with an AffyBatch object as input
model1<-plw(AffySpikeU95Subset,design=design,contrast=contrast,
epsilon=0.01)
## Look at fitted vs observed density for log(s2)
varHistPlot(model1)
## Look at fitted curve for scale parameter
scaleParameterPlot(model1)
## Selecting top genes
topRankSummary(model1,nGenes=10)
## Plotting t-statistics and log2FC for top genes
par(mfrow=c(1,2))
plotSummaryT(model1,nGenes=20)
plotSummaryLog2FC(model1,nGenes=20)
###---------------------------------------
# Analyzing with BG-adjusted and normalized PM data
pm1<-pm(bg.correct.rma(AffySpikeU95Subset, bgtype = 2))
pm2<-matrix(.C("qnorm_c", as.double(as.vector(pm1)),
as.integer(nrow(pm1)),
as.integer(ncol(pm1)))[[1]],
nrow(pm1),ncol(pm1))
data<-log2(pm2)
probenames<-probeNames(AffySpikeU95Subset)
model2<-plw(data,design=design,contrast=contrast,
probenames=probenames,epsilon=0.01)
###---------------------------------------
# Model1 and model2 should give identical result
# For example identical top ranking:
range(topRankSummary(model1)$t-
topRankSummary(model2)$t,na.rm=TRUE)
Run the code above in your browser using DataLab