Dear list, I performed a 9 array Affymetrix microarray experiment with three genotypes and three replicates of each. I am trying to do GCRMA normalization, filtering based on MAS5 calls, and then use limma with contrasts (and BH to adjust the p-values) to identify differentially expressed genes. I have no problems with the normalization and filtering, but I am having trouble with the limma aspect, specifically generating the results of the contrasts with BH and writing the results to a table. In the past, I have only used limma with experiments involving two genotypes. I have pasted my R console session below. The table I get from the below commands has four columns, one with the affy IDs and one for each of the three contrasts, but the only entries are 1,0, and -1. Any advice is much appreciated. Thank you. Matthew ----------------------------------------------------- Matthew R. Willmann, Ph.D. Research Associate, Poethig Lab University of Pennsylvania Department of Biology 433 S. University Avenue Philadelphia, PA 19104 Lab phone: 215-898-8916 Cell: 508-243-2495 Fax: 215-898-8780 ---------------------------------------------------------------------- ---------------------------------------------------------------------- ---------------------------------------- > dir() [1] "Grandcentral (3864) WT1 ATH1.CEL" [2] "Grandcentral (3865) WT2 ATH1.CEL" [3] "Grandcentral (3866) WT3 ATH1.CEL" [4] "Grandcentral (3867) cct1 ATH1.CEL" [5] "Grandcentral (3868) cct2 ATH1.CEL" [6] "Grandcentral (3869) cct3 ATH1.CEL" [7] "Grandcentral (3870) gct1 ATH1.CEL" [8] "Grandcentral (3871) gct2 ATH1.CEL" [9] "Grandcentral (3872) gct3 ATH1.CEL" [10] "R Console.txt" [11] "R Console2.txt" [12] "R ConsoleFIRST.txt" [13] "StewartarraysLIMMAresults_01072010.txt" [14] "StewartarraysLIMMAresults_01072010.txt.fmt" [15] "StewartarraysLIMMAresults_01072010_A.txt" [16] "StewartarraysLIMMAresults_01072010_B.txt" [17] "StewartarraysLIMMAresults_01072010_C.txt" [18] "StewartarraysLIMMAresults_01072010_D.txt" [19] "StewartarraysLIMMAresults_01072010_E.txt" [20] "StewartarraysLIMMAresults_01072010_E.txt.fmt" [21] "csgtest.txt" [22] "dataStewartgcrma01072010.txt" > library(affy) Loading required package: Biobase Loading required package: tools Welcome to Bioconductor Vignettes contain introductory material. To view, type 'openVignette()'. To cite Bioconductor, see 'citation("Biobase")' and for packages 'citation(pkgname)'. Loading required package: affyio Loading required package: preprocessCore > library(gcrma) Loading required package: matchprobes Loading required package: splines > library(limma) > data.Stewart=RedAffy(filenames=Cel.files) Error: could not find function "RedAffy" > data.Stewart=ReadAffy(filenames=Cel.files) Error in AllButCelsForReadAffy(..., filenames = filenames, widget = widget, : object "Cel.files" not found > Cel.files=list.files(pattern=".CEL") > Cel.files [1] "Grandcentral (3864) WT1 ATH1.CEL" "Grandcentral (3865) WT2 ATH1.CEL" [3] "Grandcentral (3866) WT3 ATH1.CEL" "Grandcentral (3867) cct1 ATH1.CEL" [5] "Grandcentral (3868) cct2 ATH1.CEL" "Grandcentral (3869) cct3 ATH1.CEL" [7] "Grandcentral (3870) gct1 ATH1.CEL" "Grandcentral (3871) gct2 ATH1.CEL" [9] "Grandcentral (3872) gct3 ATH1.CEL" > data.Stewart=ReadAffy(filenames=Cel.files) > eset<-gcrma(data,fast=F) Error in function (classes, fdef, mtable) : unable to find an inherited method for function "indexProbes", for signature "function", "character" > eset<-gcrma(data.Stewart,fast=F) Adjusting for optical effect.........Done. Computing affinities.Done. Adjusting for non-specific binding.........Done. Normalizing Calculating Expression > data.pma<-mas5calls(data.Stewart) Getting probe level data... Computing p-values Making P/M/A Calls > data.pma.exprs=exprs(data.pma) > index.9arrays=grep(".CEL",colnames(data.pma.exprs)) > numP=apply(data.pma.exprs[,index.9arrays]=="P",1,sum) > gene.select=which(numP!=0) > length(gene.select) [1] 17920 > data.wk=eset[gene.select,] > write.table(data.wk,"csgtest08.txt",sep="\t") > design <- model.matrix(~ 0+factor(c(1,1,1,2,2,2,3,3,3))) > colnames(design) <- c("group1", "group2", "group3") > fit <- lmFit(eset, design) > contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design) > fit2 <- contrasts.fit(fit, contrast.matrix) > fit2 <- eBayes(fit2) > results <- decideTests(fit2, method="separate") > objects() [1] "Cel.files" "affinity.spline.coefs" "contrast.matrix" [4] "data.Stewart" "data.pma" "data.pma.exprs" [7] "data.wk" "design" "eset" [10] "fit" "fit2" "gene.select" [13] "index.9arrays" "numP" "results" > write.table(results,"StewartarraysLIMMAresults_01082010.txt",sep="\t ",col.names=NA) >