Hi:
I performed a microarray (5 biological replicates) and analyzed using
Limma.  The software gives me the M value as a combination from the 5
biological replicates microarray.  In order to submit the data to GEO
database we need the M-value for each array ( 5 M-values).  How can I do
that? I'll appreciate your help.
 Below there is a copy of the sequence for the analysis:

library(limma)
library(geneplotter)

###Lee targets y image files
targets=readTargets("maTargetFKS12.txt")
RG=read.maimages(files=cbind(targets$FileNameCy3,targets$FileNameCy5),
  source="imagene")

##lear gal
gal=readGAL("013384_D_20050601.gal")
gal=gal[1:nrow(RG$R),]
##asignar el gal
RG$genes=gal

###definir los Spot Types y colores
spottypes <- readSpotTypes() ##usamos el default que es SpotTypes
RG$genes$Status <- controlStatus(spottypes, RG)

###definir el layour
printer= list(ngrid.r=1,ngrid.c=1,nspot.r=107,nspot.c=101)
RG$printer=structure(printer, class = "PrintLayout")

###normalize: background y normlize within
RGoriginal = RG
RG <- backgroundCorrect(RGoriginal, method="normexp", offset=50)
MA = normalizeWithinArrays(RG,method="loess",bc.method="none")

##escoger los genes nadamas
geneIndex = which(MA$genes$Status=="gene")

##preparar el disenyo para fit
design = modelMatrix(targets,ref="WT")

###ver todos los MA plots como QC: todos OK!
for(i in 1:6){
  plotMA(MA[geneIndex,i],main=targets$Name[i],ylim=range(MA$M[geneIndex,]))
  scan()
}

###estimar los paramteros usando eBayes (moderated t-stat) y calcular FDR
fit1 = lmFit(MA,design)
fit1 = eBayes(fit1)
tt=topTable(fit1,adjust="fdr",number=3000)
tt=tt[tt$P.Value<=0.01,]

###Escoger genes para incluir en las graficas
Index=as.numeric(rownames(tt))
interestGenes <-
c("CHS3","SLT2","PIR3","ECM4","ECM13","SPI1","ORF:YHR097C","MYO1","PIR3","HOG1","SWE1")
interestIndex <- lapply(interestGenes,function(x)
which(MA$genes$GeneName==x))
names(interestIndex) <- interestGenes

##MA-PLOT
plot(fit1$Amean[geneIndex],fit1$coef[geneIndex],pch=".",col="grey",ylab="M",xlab="A")
points(fit1$Amean[Index],fit1$coef[Index],pch=16,col="blue",cex=.5)
for(i in seq(along=interestGenes)){

text(fit1$Amean[interestIndex[[i]]],fit1$coef[interestIndex[[i]]],interestGenes[i],cex=.5)
}

##VOLCANO
plot(fit1$coef[geneIndex],-log10(fit1$p.value[geneIndex]),pch=".",col="grey",xlab="M",ylab="-log_10
p-value")
points(fit1$coef[Index],-log10(fit1$p.value[Index]),pch=16,col="blue",cex=.5)
for(i in seq(along=interestGenes)){

text(fit1$coef[interestIndex[[i]]],-log10(fit1$p.value[interestIndex[[i]]]),interestGenes[i],cex=.5)
}


###hacer xls file
todos = topTable(fit1,adjust="fdr",number=10807)
write.csv(todos,file="myo1.csv")
-- 
José F. Rodríguez Quiñones
Graduate student
School of Medicine
Medical Sciences Campus
University of Puerto Rico

	[[alternative HTML version deleted]]

