affymetrix in maanova, code
0
0
Entering edit mode
Morten ▴ 300
@morten-929
Last seen 9.6 years ago
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")
qvalue maanova qvalue maanova • 1.2k views
ADD COMMENT

Login before adding your answer.

Traffic: 702 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6