[BioC] 2fold sigchanges from limma

[Gordon Smyth]

Dear Ivan,

Also, you can use the lfc (log-fold-change) argument to decideTests().

Best wishes
Gordon

>Date: Sat, 02 Dec 2006 09:50:31 -0500
>From: "James W. MacDonald" <jmacdon at med.umich.edu>
>Subject: Re: [BioC] 2fold sigchanges from limma
>To: Ivan Baxter <ibaxter at purdue.edu>
>Cc: bioconductor at stat.math.ethz.ch
>
>Hi Ivan,
>
>Take a look at limma2annaffy in the affycoretools package. I think this
>will do what you want in a one-line function.
>
>Best,
>
>Jim
>
>Ivan Baxter wrote:
> > Greetings-
> >
> > I am using limma to analyze a fairly simple multiple treatment
> > experiment and I would like to pull out all the genes which are
> > significantly expressed and change at least two fold. I have figured out
> > a way to do this, but it is a little complicated and I am running into
> > problems downstream when I try to format output tables using annaffy
> > (described below).  Is there a simple way to filter the output of
> > topTable to only get the genes which are significant and change more
> > than a give fold-change cutoff?
> >
> >
> > The way that I am going about this is....
> >
> > cont.matrix <- makeContrasts(
> >     EGvsC    = EG.C.C.C - C.C.C.C,
> >     ECvsC    = C.EC.C.C - C.C.C.C,
> >     EGECvsC  = EG.EC.C.C - C.C.C.C,
> >     GTvsC   = C.C.G.C - C.C.C.C,
> >     GT_APvsC   = C.C.G.A  - C.C.C.C,
> >     GT_APvsGT = C.C.G.A  - C.C.G.C,
> >     levels=design)
> > fit2 <- contrasts.fit(fit, cont.matrix)
> > fit2 <- eBayes(fit2)
> > results <- decideTests(fit2)
> > idx <- abs(fit2$coefficients) > 1
> > #make a dataframe like results where 1/-1 indicates sigchange greater
> > than 2 fold
> > comb <- data.frame(gene = rownames(results))
> > for(i in 1:length(rownames(results))){
> >     for(j in 1:length(colnames(results))){
> >
> >         if(results[i,j] != 0 & idx[i,j] == "TRUE"){
> >             comb[i,j+1] <- results[i,j]
> >            }
> >        else{comb[i,j+1] <-0}
> >     }
> > }
> >
> > #but then I want to make an output table with gene names, gene
> > annotations, fold change, pvalue and the expression values #across the
> > arrays......
> >
> >
> > cax1_genes <- unlist(as.list(comb$gene[comb$GTAPvsC != 0]))
> > syms <- unlist(mget(cax1_genes, hgu133a2GENENAME))
> > test <- match(cax1_genes, geneNames(human.eps2))
> > anncols <- aaf.handler()[c(1:4,7)]
> > anntable <- aafTableAnn(geneNames(human.eps2)[test], "hgu133a2", anncols)
> >
> > #now I need to get the fold- change and adj.p.value for each gene, and
> > here is where I run into trouble.
> > #I tried pulling out all the significant changes using topTable and then
> > pulling out the list of genes that met my criteria, but...
> > contp <- 5 # this is the contrast which is being tested
> > caxsig <- length(which(p.adjust(fit2$p.value[,contp], method = "BH") <
> > 0.05))
> > cax1sig <- topTable(fit2, coef =contp, number =caxsig, adjust.method =
> > "none", sort.by = "p", resort.by = "M")
> > testtable <- aafTable("log2 change" = format(cax1sig$M[cax1sig$ID ==
> > cax1_genes], digits = 2),
> >                                     "pval" =
> > format(cax1sig$adj.P.Val[cax1sig$ID == cax1_genes], digits = 2))
> > anntablep <- merge(anntable, testtable)
> >
> >
> > # doesn't work, I get the following error:
> >  > testtable <- aafTable("log2 change" = format(cax1sig$M[cax1sig$ID ==
> > cax1_genes], digits = 2), "pval" = format(cax1sig$adj.P.Val[cax1sig$ID
> > == cax1_genes], digits = 2))
> > Warning messages:
> > 1: longer object length
> >         is not a multiple of shorter object length in: cax1sig$ID ==
> > cax1_genes
> > 2: longer object length
> >         is not a multiple of shorter object length in: cax1sig$ID ==
> > cax1_genes
> >
> > If this is actually a good way to pull out the two fold genes, could
> > anyone tell me what I am doing wrong with this last step?
> >
> > thanks in advance.
> >
> > Ivan
> >
> >
> >
> >  > sessionInfo()
> > Version 2.3.0 (2006-04-24)
> > i386-pc-mingw32
> >
> > attached base packages:
> >  [1] "grid"      "splines"   "tools"     "methods"   "stats"
> > "graphics"  "grDevices" "utils"     "datasets"  "base"
> >
> > other attached packages:
> >      annaffy         KEGG           GO     hgu133a2 RColorBrewer
> > geneplotter     annotate       hexbin   colorspace      lattice
> > genefilter
> >      "1.4.0"     "1.12.0"     "1.12.0"     "1.12.0"      "0.2-3"
> > "1.10.0"     "1.10.0"      "1.6.0"        "0.9"     "0.13-8"     "1.10.1"
> >     survival        limma        gcrma  matchprobes         affy
> > affyio      Biobase
> >       "2.24"      "2.7.3"      "2.4.1"      "1.4.0"     "1.10.0"
> > "1.0.0"     "1.10.0"
> >
> >
> >
>
>
>--
>James W. MacDonald
>University of Michigan
>Affymetrix and cDNA Microarray Core
>1500 E Medical Center Drive
>Ann Arbor MI 48109
>734-647-5623