Hi Mali, On 6/5/2011 2:36 AM, mali salmon wrote: > OK, I have noticed that when I formulate the model in a different way, I get > the expected results, and the expected columns are printed. > The model that works fine is: > >> TS<- paste(targets$Strain, targets$Treatment, sep=".") >> TS<- factor(TS, levels=c("WT.U","WT.S","Mu.U","Mu.S")) >> design<- model.matrix(~0+TS) >> colnames(design)<- levels(TS) >> fit<- lmFit(filtered, design) >> cont.matrix<- makeContrasts(WT.SvsU=WT.S-WT.U, Mu.SvsU=Mu.S-Mu.U, > Diff=(Mu.S-Mu.U)-(WT.S-WT.U),levels=design) >> fit2<- contrasts.fit(fit, cont.matrix) >> fit2<- eBayes(fit2) > > Then when I run "limma2annaffy" I get the correct columns. > > However, when I use a different approach (but equivalent according to the > user guide) below, the wrong columns are printed when I use > "limma2annaffy": >> Strain<- factor(targets$Strain, levels=c("WT","Mu")) >> Treatment<- factor(targets$Treatment, levels=c("U","S")) >> design<- model.matrix(~Strain*Treatment) >> fit<- lmFit(filtered, design) >> cont.matrix<- cbind(WT.SvsU=c(0,0,1,0),Mu.SvsU=c(0,0,1,1),Diff=c(0,0,0,1)) >> fit2<- contrasts.fit(fit, cont.matrix) >> fit2<- eBayes(fit2) > > Any idea what I did wrong? You didn't do anything wrong per se, if anything it is a shortcoming of limma2annaffy(). In order to make the code simpler, I require that the model be specified using a cell means model (as you did above where it worked for you), rather than a factor effects model (the parameterization you tried that didn't work). Maybe this isn't stated as explicitly as it should be. I'll take a look at the documentation and see if I can make some changes. Best, Jim > Thank > Mali > > > On Sun, Jun 5, 2011 at 8:54 AM, mali salmon<shalmom1 at="" gmail.com=""> wrote: > >> Dear Jim >> Thanks to your help I can generate HTML lists now, but the expression is >> missing for some of the samples. >> I also tried just to output lists of genes that are significant in each of >> the contrasts, and out of 10 samples, the expression of only 2-5 is printed >> (in both the txt and html files). >> Below is the function I used >>> limma2annaffy(filtered, fit2, design,cont.matrix, annotation(eset),adjust >> = "fdr",anncols = aaf.handler(chip="yeast2.db")[-c(5,6)], >> fldfilt=1,tstat=TRUE,pfilt = 0.05,text=TRUE,expression = TRUE) >> >> For "Diff" contrast genes only the expression of "Mu.S.1", "Mu.S.2" is >> printed >> For "Mu.SvsU" and "WT.SvsU" contrasts only the expression of "WT.S.1" >> "WT.S.2" "WT.S.3" "Mu.S.1" "Mu.S.2" is printed >> >> Below are the targets, contrast matrix and design objects >> >>> target >> FileName Strain Treatment >> 1 WT.U.1 WT U >> 2 WT.U.2 WT U >> 3 WT.U.3 WT U >> 4 WT.S.1 WT S >> 5 WT.S.2 WT S >> 6 WT.S.3 WT S >> 7 Mu.U.1 Mu U >> 8 Mu.U.2 Mu U >> 9 Mu.S.1 Mu S >> 10 Mu.S.2 Mu S >> >> >>> cont.matrix >> WT.SvsU Mu.SvsU Diff >> [1,] 0 0 0 >> [2,] 0 0 0 >> [3,] 1 1 0 >> [4,] 0 1 1 >> >>> design >> (Intercept) StrainMu TreatmentS StrainMu:TreatmentS >> 1 1 0 0 0 >> 2 1 0 0 0 >> 3 1 0 0 0 >> 4 1 0 1 0 >> 5 1 0 1 0 >> 6 1 0 1 0 >> 7 1 1 0 0 >> 8 1 1 0 0 >> 9 1 1 1 1 >> 10 1 1 1 1 >> >> Thanks >> Mali >> >> >> >> >> >> On Fri, Jun 3, 2011 at 6:51 AM, mali salmon<shalmom1 at="" gmail.com=""> wrote: >> >>> Thanks a lot Jim for your prompt support >>> Mali >>> >>> >>> On Thu, Jun 2, 2011 at 8:36 PM, James W. MacDonald<jmacdon at="" med.umich.edu="">>>> wrote: >>> >>>> Hi Mali, >>>> >>>> >>>> On 6/2/2011 3:00 PM, mali salmon wrote: >>>> >>>>> Thanks Jim for your reply. >>>>> I repeated the analysis without resetting the annotation slot but I >>>>> still >>>>> get the same error. >>>>> Below is the code I used: >>>>> >>>>> library(limma) >>>>>> library(affy) >>>>>> library(affycoretools) >>>>>> data<- ReadAffy() >>>>>> eset<- rma(data) >>>>>> #remove low variance genes >>>>>> library(genefilter) >>>>>> eset2<-varFilter(eset,var.func=IQR,var.cutoff=0.1) >>>>>> #remove control probes (AFFX) >>>>>> filtered<-featureFilter(eset2, require.entrez=F, remove.dupEntrez=F) >>>>>> targets<-readTargets("targets.txt") >>>>>> TS<- paste(targets$Strain, targets$Treatment, sep=".") >>>>>> Strain<- factor(targets$Strain, levels=c("WT","Mu")) >>>>>> Treatment<- factor(targets$Treatment, levels=c("U","S")) >>>>>> design<- model.matrix(~Strain*Treatment) >>>>>> fit<- lmFit(filtered, design) >>>>>> cont.matrix<- >>>>>> cbind(WT.SvsU=c(0,0,1,0),Mu.SvsU=c(0,0,1,1),Diff=c(0,0,0,1)) >>>>>> fit2<- contrasts.fit(fit, cont.matrix) >>>>>> fit2<- eBayes(fit2) >>>>>> results<- decideTests(fit2,lfc=1) >>>>>> vennSelect(filtered, design, results, cont.matrix, fit2) >>>>>> >>>>> >>>> The problem is that you are by default trying to get annotation for >>>> >>>>> aaf.handler()[c(1:3,6:7, 9:12)] >>>> [1] "Probe" "Symbol" "Description" "GenBank" >>>> [5] "Gene" "UniGene" "PubMed" "Gene Ontology" >>>> [9] "Pathway" >>>> >>>> But since this is a yeast chip, not all of these are available. For this >>>> sort of issue I have included a '...' argument to vennSelect() that allows >>>> you to pass arguments to lower level functions, so you can just do something >>>> like >>>> >>>> library(annaffy) ## we need this to expose the aaf.handler() function >>>> vennSelect(filtered, design, results, cont.matrix, fit2, anncols = >>>> aaf.handler(chip = "yeast2.db")) >>>> >>>> Which will then give you all possible annotations for the yeast2.db chip, >>>> which are >>>> >>>>> aaf.handler(chip = "yeast2.db") >>>> [1] "Probe" "Description" "Chromosome" >>>> [4] "Chromosome Location" "PubMed" "Gene Ontology" >>>> [7] "Pathway" >>>> >>>> Unfortunately, you will still get the error you see, which comes from the >>>> call for Gene Ontology data, so if you amend to >>>> >>>> vennSelect(filtered, design, results, cont.matrix, fit2, anncols = >>>> aaf.handler(chip = "yeast2.db")[-6]) >>>> >>>> Then it will work. I will contact Colin Smith, the maintainer of annaffy, >>>> separately to see if he can track down the error you see with GO terms. >>>> >>>> Best, >>>> >>>> Jim >>>> >>>> >>>> >>>>> Loading required package: yeast2.db >>>>> Loading required package: org.Sc.sgd.db >>>>> Error in sqliteExecStatement(con, statement, bind.data) : >>>>> RS-DBI driver: (error in statement: cannot join using column gene_id - >>>>> column not present in both tables) >>>>> >>>>> Thanks >>>>> Mali >>>>> >>>>> >>>>> >>>>> On Thu, Jun 2, 2011 at 2:45 PM, James W. MacDonald< >>>>> jmacdon at med.umich.edu>wrote: >>>>> >>>>> Hi Mali, >>>>>> >>>>>> >>>>>> On 6/2/2011 7:40 AM, mali salmon wrote: >>>>>> >>>>>> Dear list >>>>>>> I am trying to use affycoretools for the analysis of affy yeast >>>>>>> arrays. >>>>>>> When I tried to use vennSelect in order generate HTML files containing >>>>>>> the >>>>>>> results of decideTests I got the following error: >>>>>>> >>>>>>> annotation(eset)<-"yeast2.db" >>>>>>> >>>>>>>> >>>>>>>> >>>>>>> Good thing you didn't snip the above line from your code. In the >>>>>> future, it >>>>>> is better to give us more than what you think we need, rather than >>>>>> less. >>>>>> >>>>>> The problem here is that you are resetting the annotation slot for your >>>>>> ExpressionSet (unnecessarily, and incorrectly). Not sure why you are >>>>>> doing >>>>>> this; did you see some code somewhere that lead you to believe this is >>>>>> either desired or necessary? If so, please tell us so we can hopefully >>>>>> get >>>>>> it corrected. >>>>>> >>>>>> Note that the annotation slot for an ExpressionSet doesn't imply >>>>>> annotation >>>>>> of (in this case) probeset IDs to genes, etc. The annotation slot >>>>>> contains >>>>>> name of the package that has the mappings of probes to probesets that >>>>>> were >>>>>> used to created the summarized values in that ExpressionSet. This is >>>>>> nice to >>>>>> have if you save the object, and wonder later how exactly the data were >>>>>> summarized. >>>>>> >>>>>> So start back at the line preceding the line above where you reset the >>>>>> annotation slot, and try again, without resetting the annotation slot. >>>>>> >>>>>> Best, >>>>>> >>>>>> Jim >>>>>> >>>>>> >>>>>> ... >>>>>> >>>>>>> >>>>>>> fit2<- eBayes(fit2) >>>>>>>> results<- decideTests(fit2,lfc=1) >>>>>>>> vennSelect(filtered,design,results,cont.matrix,fit2) >>>>>>>> >>>>>>>> Error in sqliteExecStatement(con, statement, bind.data) : >>>>>>> RS-DBI driver: (error in statement: cannot join using column gene_id >>>>>>> - >>>>>>> column not present in both tables) >>>>>>> >>>>>>> Any idea why I get this error? >>>>>>> Thanks >>>>>>> Mali >>>>>>> >>>>>>> sessionInfo() >>>>>>> >>>>>>>> >>>>>>>> R version 2.12.1 (2010-12-16) >>>>>>> Platform: i386-redhat-linux-gnu (32-bit) >>>>>>> >>>>>>> locale: >>>>>>> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C >>>>>>> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 >>>>>>> [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 >>>>>>> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C >>>>>>> [9] LC_ADDRESS=C LC_TELEPHONE=C >>>>>>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C >>>>>>> >>>>>>> attached base packages: >>>>>>> [1] grid stats graphics grDevices utils datasets >>>>>>> methods >>>>>>> [8] base >>>>>>> >>>>>>> other attached packages: >>>>>>> [1] hgu133a2.db_2.4.5 org.Hs.eg.db_2.4.6 annaffy_1.22.0 >>>>>>> [4] affycoretools_1.22.0 KEGG.db_2.4.5 GO.db_2.4.5 >>>>>>> [7] gplots_2.8.0 caTools_1.11 bitops_1.0-4.1 >>>>>>> [10] gdata_2.8.1 gtools_2.6.2 yeast2.db_2.4.5 >>>>>>> [13] org.Sc.sgd.db_2.4.6 RSQLite_0.9-4 DBI_0.2-5 >>>>>>> [16] annotate_1.28.1 AnnotationDbi_1.12.0 genefilter_1.32.0 >>>>>>> [19] yeast2cdf_2.7.0 affy_1.28.0 Biobase_2.10.0 >>>>>>> [22] limma_3.6.9 >>>>>>> >>>>>>> loaded via a namespace (and not attached): >>>>>>> [1] affyio_1.18.0 biomaRt_2.6.0 Biostrings_2.18.2 >>>>>>> [4] Category_2.16.1 gcrma_2.22.0 GOstats_2.16.0 >>>>>>> [7] graph_1.28.0 GSEABase_1.12.2 IRanges_1.8.9 >>>>>>> [10] preprocessCore_1.12.0 RBGL_1.26.0 RCurl_1.5-0 >>>>>>> [13] splines_2.12.1 survival_2.36-2 tcltk_2.12.1 >>>>>>> [16] tools_2.12.1 XML_3.2-0 xtable_1.5-6 >>>>>>> >>>>>>> [[alternative HTML version deleted]] >>>>>>> >>>>>>> _______________________________________________ >>>>>>> Bioconductor mailing list >>>>>>> Bioconductor at r-project.org >>>>>>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>>>>> Search the archives: >>>>>>> http://news.gmane.org/gmane.science.biology.informatics.conductor >>>>>>> >>>>>>> >>>>>> -- >>>>>> James W. MacDonald, M.S. >>>>>> Biostatistician >>>>>> Douglas Lab >>>>>> University of Michigan >>>>>> Department of Human Genetics >>>>>> 5912 Buhl >>>>>> 1241 E. Catherine St. >>>>>> Ann Arbor MI 48109-5618 >>>>>> 734-615-7826 >>>>>> ********************************************************** >>>>>> Electronic Mail is not secure, may not be read every day, and should >>>>>> not be >>>>>> used for urgent or sensitive issues >>>>>> >>>>>> >>>>> >>>> -- >>>> James W. MacDonald, M.S. >>>> Biostatistician >>>> Douglas Lab >>>> University of Michigan >>>> Department of Human Genetics >>>> 5912 Buhl >>>> 1241 E. Catherine St. >>>> Ann Arbor MI 48109-5618 >>>> 734-615-7826 >>>> ********************************************************** >>>> Electronic Mail is not secure, may not be read every day, and should not >>>> be used for urgent or sensitive issues >>>> >>> >>> >> > -- James W. MacDonald, M.S. Biostatistician Douglas Lab University of Michigan Department of Human Genetics 5912 Buhl 1241 E. Catherine St. Ann Arbor MI 48109-5618 734-615-7826 ********************************************************** Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues