medinalab.org

  • Increase font size
  • Default font size
  • Decrease font size
Home Protocols Bioinformatics Latest MAANOVA commands

Latest MAANOVA commands

E-mail Print PDF

 

 - Mickey DeSalvo - 13 Feb 2009 
Change your working directory to folder where datafile and designfile are located.

library(maanova)

library(qvalue)

raw<-read.madata(datafile="arraydata2_filled_11feb.txt",designfile="arraydesign6_11feb.txt",arrayType=c("twoColor"),header=TRUE,spotflag=FALSE,n.rep=1,avgreps=0,log.trans=TRUE,metarow=0,metacol=0,row=0,col=0,probeid=1,intensity=2)

	reads in your data file.

model<-fitmaanova(raw,formula=~Array+Sample,method=c("REML"))

	fits maanova model to your experiment.
	above if for a fixed-effect model.
	below is for a mixed-effect model:
	model<-fitmaanova(raw,formula=~Array+Sample, random=~Array)
	if your experiment involves a dye-swap, then your formula could include "Dye" to take into account the dye effect.

output<-matest(raw,model,term="Sample",n.perm=500,test.type=c("ftest"))

	performs f-test given your data and model.
	for a t-test do the following:
	output<-matest(raw,model,term="Sample",PairContrast(6),n.perm=500,test.type=c("ttest"))

output.adj<-adjPval(output,'jsFDR')

	adjusts p-values based on storey's qvalue approach.

volcano<-volcano(output.adj,threshold=c(0.05,0.05),method=c("unadj","unadj"))

	picks DE genes and shows volcano plots.
	above is to pick DE genes at unadjusted a=0.05 for both F1 and Fs tests.
	below is to pick DE genes at fdr-adjusted a=0.01 for both F1 and Fs tests.
	volcano<-volcano(output.adj,threshold=c(0.01,0.01),method=c("fdr","fdr"))

	the following commands are necessary to generate an output file that contains CloneIDs, pvalues, and gene expression estimates for the DE genes in each sample according to the Fs test.

cluster<-macluster(model,term="Sample",idx.gene=volcano$idx.Fs,what="gene",method="kmean",kmean.ngroups=6,n.perm=100)

Datacon.kmeanFs<-consensus(cluster, 0.8, draw=FALSE)

names<-raw$probeid [volcano$idx.Fs]

Ptab<-output$Fs$Ptab [volcano$idx.Fs]

	above is choosing tabulated p-vaules.
	below is for choosing adjusted p-values:
	adjPtab<-output.adj$Fs$adjPtab [volcano$idx.Fs]

Fs<-data.frame(names,Datacon.kmeanFs$data.draw[,],Ptab)

	if using adjusted p-values change "Ptab" in the above line to "adjPtab".

write.table(Fs, "Fs.txt", row.names=FALSE, sep="\t", quote=FALSE)

	Fs.txt can be changed to generate a more descriptive table title.