[BioC] Limma with contrasts

[Matthew Willmann]

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)