Hello everyone, Im just fowarding a mail I got from the MAANOVA maling list before christmas. I hope Jason and Yong dont mind. Contains code for maanova anlysis of affymetrix chips.. thought someone here may be interested.. best morten #################################################################### # affyprocessing.R - v1.2 7/14/2004 # # Template for pre-processing of Affymetrix CEL data using # R/affy (Bioconductor) and formatting for R/maanova analysis # [ref: http://www.bioconductor.org] # # Jason Affourtit (jason.affourtit at jax.org) # originally developed by Yong Woo (yhw at jax.org) # #################################################################### #read in R/affy (rma) library library(affy) #read in all CEL files in R working directory to create AffyBatch object CELData = ReadAffy() #process CEL data using rma (default options) rma.CELData = rma(CELData) #export rma data into dataframe object rma.expr = exprs(rma.CELData) #transform to put row name (Probe ID) into the first column rma.expr.df = data.frame(ProbeID=row.names(rma.expr),rma.expr) #save as a tab-delimited text file with no row names write.table(rma.expr.df,"rma.expr.dat",sep="\t",row=F,quote=F) #check order and names of samples to set up design matrix pData(CELData) #create design matrix for R/maanova analysis design.matrix=data.frame(Array=row.names(pData(rma.CELData)),Strain=c( "wt","wt","wt","mut","mut","mut"),Sample=c(1,2,3,4,5,6),Dye=c(1,1,1,1, 1,1)) #export design matrix to tab-delimited text file write.table(design.matrix,"design.dat",sep="\t",row=F,quote=F) #save R environment file save.image("ProjectName.Rdata") ############################################################# # affymaanova.R - v1.3 11/01/2004 # # Template for R/maanova analysis of data pre-processed # using rma (affyprocessing.R) # [ref: #http://www.jax.org/staff/churchill/labsite/] # # Jason Affourtit (jason.affourtit at jax.org) # originally developed by Yong Woo (yhw at jax.org) # ############################################################# #load the R/maanova library library(maanova) #read in the rma-processed experimental data and design file raw.data = read.madata("rma.expr.dat",designfile="design.dat",cloneid=1,pmt=2,spo t=F) #convert the raw data from all arrays into an madata object data = createData(raw.data,n.rep=1,log.trans=F) #make the model based on the design model.full.fix = makeModel(data=data,formula=~Strain) #fit fixed model ANOVA anova.full.fix = fitmaanova(data,model.full.fix) #verify correct arrays are being used in analysis prior to ftest model.full.fix$design #fit permutation ftest and plot histogram of Fs p-values ftest.full.fix = matest(data,model.full.fix,term="Strain",n.perm=500,shuffle.method="re sid", pval.pool=TRUE) hist(ftest.full.fix$Fs$Pvalperm, main ="Fs Pvalue histogram - ProjectName", breaks=100) #output R workspace image in a file following f-test save.image("ProjectName.Rdata") #determine ANOVA object column order to set up output below anova.full.fix$Strain.level #assign column names and calculate fold change and ratios #FoldChange and ratios default to [tester - control] or fold change of tester #be sure that you have these set correctly depending upon column order determined above #relative to control tester= 2 #column 2= mutant control=1 #column 1= wild type RelativeFoldchange=sign(anova.full.fix$Strain[,tester] - anova.full.fix$Strain[,control])*2^(abs((anova.full.fix$Strain[,tester ]-anova.full.fix$Strain[,control]))) RelativeFoldchangeName=paste("relativefoldchange",anova.full.fix$Strai n.level[tester],"relative to",anova.full.fix$Strain.level[control]) #load qvalue library and generate qvalues #[ref:http://faculty.washington.edu/~jstorey/qvalue/] library(qvalue) #calculate q-value based upon permuted p-value and summarize results ftest.full.fix.Fs.qobj=qvalue(ftest.full.fix$Fs$Pvalperm) qsummary(ftest.full.fix.Fs.qobj) #test if the order of p-value is scrambled - should yield TRUE for all rows table(ftest.full.fix.Fs.qobj$pval==ftest.full.fix$Fs$Pvalperm) #concatenate them together result.df=data.frame(data$cloneid, RelativeFoldchange, ftest.full.fix$Fs$Pvalperm, ftest.full.fix.Fs.qobj$qval) #assign names to each column names(result.df)=c("cloneid",RelativeFoldchangeName, "FsPvalue", "Qvalue.FDR") #create subset where qvalue <0.05 index.Fs=ftest.full.fix.Fs.qobj$qval<0.05 #write two files - all genes and top hits write.table(result.df,"all.genes.result.dat",sep="\t",row=F, quote=F) write.table(result.df[index.Fs,],"top.hits.result.dat",sep="\t",row=F, quote=F) #generate volcano plots jpeg("volcanoPlot.jpeg") plot(anova.full.fix$Strain[,tester]-anova.full.fix$Strain[,control], -log10(ftest.full.fix$F1$Ptab),main="volcano plot - ProjectName (q<0.05)", xlab=paste("log((",anova.full.fix$Strain.level[tester],") - ", anova.full.fix$Strain.level[control],")",sep=""),ylab="-log10(F1 Tabulated P-Value)",cex=.3,pch=3,col="blue") points(anova.full.fix$Strain[index.Fs,tester]-anova.full.fix$Strain[in dex.Fs,control], -log10(ftest.full.fix$F1$Ptab[index.Fs]),pch=1,col="red") dev.off() #output R workspace image in a file save.image("ProjectName.Rdata")