Microarray

• http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12613

exp.txt

Filename        Sample  Group   Description
CEL/GSM315919.CEL       ICD1    ICD     ICD-1
CEL/GSM315920.CEL       ICD2    ICD     ICD-2
CEL/GSM315921.CEL       ICD3    ICD     ICD-3
CEL/GSM315922.CEL       FoxJ1_1 FoxJ1   ICD-FoxJ1-1
CEL/GSM315923.CEL       FoxJ1_2 FoxJ1   ICD-FoxJ1-2

affy-RMA.R

library(affy)
 
dataset_name <- 'GSE12613_foxj1'
 
exp_table <- read.table(file="exp.txt",header=T,stringsAsFactors=FALSE,sep="\t")
files_vector <- as.vector(exp_table$Filename,mode='character')
samples_vector <- as.vector(exp_table$Sample,mode='character')
 
affybatch_raw <- ReadAffy(filenames=files_vector,sampleNames=samples_vector)
affybatch_corrected <- bg.correct(affybatch_raw, method='rma')
affybatch_normalized <- normalize(affybatch_corrected, method='quantiles')
eset_rma <- rma(affybatch_raw)
write.exprs(eset_rma, file=paste(dataset_name,'.eset_rma.txt',sep=''))

affy-plot.R

dataset_name <- 'GSE12613_foxj1'
 
library(affy)
 
load( paste(dataset_name,'.affybatch_raw',sep='') )
load( paste(dataset_name,'.affybatch_normalized',sep='') )
load( paste(dataset_name,'.eset_rma',sep='') )
 
## Boxplot
png(filename=paste("plot/",dataset_name,".box_raw.png",sep=''),
      width=1024,height=640);
par(mar=c(10,4,4,2))
boxplot( as.data.frame(log2(intensity(affybatch_raw))), las=2, mex=0.1,
          main=paste(dataset_name,'(raw)',sep=''),
          ylab="log2 intensity")
dev.off();
 
png(filename=paste("plot/",dataset_name,".box_rma.png", sep=''),
      width=1024,height=640);
par(mar=c(10,4,4,2))
boxplot( as.data.frame(exprs(eset_rma)), las=2, mex=0.1,
      main=paste(dataset_name,"(normalized - RMA)") )
dev.off();
 
## RMA, log2 intensity
rma_intensity <- as.data.frame(eset_rma);
rma_dist <- dist( as.matrix(rma_intensity) )
rma_clust <- hclust( rma_dist, method="average")
png(filename=paste("plot/",dataset_name,".rma_log_clust.png",sep=''),
        width=1024,height=640);
 
## RMA, log2 Rank
rma_rank <- matrix(data=NA, nrow=nrow(eset_rma), ncol=ncol(eset_rma) );
nsample <- ncol(eset_rma)
for(i in 1:nsample) {
  rma_rank[,i] <- rank(exprs(eset_rma[,i]) );
}
colnames(rma_rank) <- sampleNames(eset_rma)
rma_rank_dist <- dist( t(rma_rank) )
rma_rank_clust <- hclust(rma_rank_dist, method="average")
png(filename=paste("plot/",dataset_name,'.rma_rank_clust.png',sep=''),
      width=1024, height=640);
plot(rma_rank_clust,
      main=paste(dataset_name,":hclust_average(rma,log2_rank)",sep=''))
dev.off()
 
dataset_name <- 'GSE12613_foxj1'
 
library(affy)
library(affyPLM)
 
load( paste(dataset_name,'.affybatch_raw',sep='') )
 
affybatch_raw_PLM <- fitPLM(affybatch_raw)
 
sample_names <- sampleNames(affybatch_raw)
for(i in 1:length(sample_names)) {
  filename_img <- paste(sample_names[i],".PLM_resids.png", sep="")
  png(filename=filename_img, width=1024, height=1024);
  image( affybatch_raw_PLM, which=i, type="resids", add.legend=TRUE )
  dev.off()
}
 
png(filename=paste(dataset_name,'.RLE.png',sep=''),
              width=1024,height=800)
par(mar=c(10,4,4,2))
RLE(affybatch_raw_PLM, fin=c(1024, 600), las=2, mex=0.1,
    main=paste("Relative Log Expression(RLE) for ",dataset_name, sep=''))
dev.off()