[BioC] Limma, decideTests

[Gordon Smyth]

 > Date: Mon, 17 Jan 2005 15:12:12 +0100
 > From: "Ingunn Berget" <ingunn.berget at umb.no>
 > Subject: [BioC] Limma, decideTests
 > To: <bioconductor at stat.math.ethz.ch>
 >
 > Hello mail-list
 >
 > Experiment
 > This microarray experiment was conducted to study gene expression in
 > bacteria at different growth conditions. There are 12 different conditions:
 >
 >> unique(conditions)
 >  [1] "gr46"     "NaCl"     "Glycerol" "HCl"      "NaOH"     "CH3COOH"
 >  [7] "EtOH"     "gr15"     "PK"       "BC7ppm"   "BC9ppm"   "EtBr"
 >
 > The one called PK is the positive control (normal conditions), whereas the
 > other conditions are "stress conditions". A common reference, used for all
 > arrays, is labelled with Cy3. On each array, one of the other conditions is
 > labelled with Cy5. There are three biologial replicates of each
 > array/condition, expect NaOH and HCl where something went wrong with one of
 > replicates. Totally 34 arrays.
 >
 > The goal is to find which genes are differentially expressed in PK and the
 > other conditions.
 >
 > ANALYSIS, copy of RCODE + some comments:
 > MAptl2 <- normalizeWithinArrays(modobj2)#modobj2: RGlist with raw data
 > D <- modelMatrix(modobj2$target,ref = "Ref")
 > GG <- MAptl2[keep,] #keep: logical index because want to use genes only, not
 > controls
 > cor<- duplicateCorrelation(GG,design=D)
 > fit <- lmFit(GG,design=D,ndups=2,correlation = cor$consensus.correlation)
 > fit <- eBayes(fit)
 > topTable(fit,n=30,adjust="fdr")
 >
 >
 > The topTable command results in a list with genes having small p-values,
 > have also tried with
 > topTable(fit,coef.=i,n=30,adjust="fdr") and different values of i
 > (i=1,2,..,12)
 >
 >> contrast.matrix <-
 >> 
makeContrasts(gr15.PK=gr15-PK,gr46.PK=gr46-PK,gr46.15=gr46-gr15,NaCl.PK=NaCl-PK,EtBr.PK=EtBr-PK,
 > +
 > 
EtOH.PK=EtOH-PK,NaOH.PK=NaOH-PK,CH2COOH.PK=CH3COOH-PK,HCl.PK=HCl-PK,BC7.PK=BC7ppm-PK,
 > + BC9.PK=BC9ppm-PK,BC9.BC7=BC9ppm-BC7ppm,Gly.PK=Glycerol-PK,levels=D)
 >>
 >> fit2 <- contrasts.fit(fit,contrast.matrix)
 >> fit2 <- eBayes(fit2)
 >>
 >> results <- decideTests(fit2,method="global")
 >> summary(results)
 >
 >    gr15.PK gr46.PK gr46.15 NaCl.PK EtBr.PK EtOH.PK NaOH.PK CH2COOH.PK HCl.PK
 > BC7.PK BC9.PK BC9.BC7 Gly.PK
 > -1
 > 0
 > 1
 >> table(results)
 > results
 >    -1     0     1
 >    13 78780    26
 >>
 >
 >
 > Questions:
 > Since the topTable command lead to a genelist with very low p-values, this
 > means that there are differentially expressed genes in the data?
 > This is differentially expressed compared to the common reference?
 > Is this contrast matrix the correct for comparing gene expression in the
 > stress conditions and PK
 > I don't understand why the summary(result) is "empty", does this means that
 > a)  there's something wrong in the code, b)  there is not enough data for
 > making this many contrasts (only 3 replicates?) c) no differential expressed
 > genes??

An entirely NULL (i.e., empty) summary table was arising when some of the 
p-values were NA (missing) and hence so were some of the TestResults. This 
bug is now fixed in limma 1.8.19 available from http://bioinf.wehi.edu.au

Another issue is why you have NA p-values. Often one can revisit the 
background correct/quality weighting scheme to avoid this.

Gordon

 > If it is of help, here is a part of the target-slot in the original RG list
 > (have not indluded all since this mail already is long)
 >
 >    Slidenumber               FileNameCy3               FileNameCy5      Cy5
 > Cy3         Name
 > 1         1279     Cy3_1279_46gr_III.txt     Cy5_1279_46gr_III.txt     gr46
 > Ref     46gr_III
 > 2         1608      Cy3_1608_NaCl_II.txt      Cy5_1608_NaCl_II.txt     NaCl
 > Ref      NaCl_II
 > 3         1609  Cy3_1609_Glycerol_II.txt  Cy5_1609_Glycerol_II.txt Glycerol
 > Ref  Glycerol_II
 > 4         1610      Cy3_1610_HCl_III.txt      Cy5_1610_HCl_III.txt      HCl
 > Ref      HCl_III
 > 5         1612     Cy3_1612_NaOH_III.txt     Cy5_1612_NaOH_III.txt     NaOH
 > Ref     NaOH_III
 > 6         1613  Cy3_1613_CH3COOH_III.txt  Cy5_1613_CH3COOH_III.txt  CH3COOH
 > Ref  CH3COOH_III
 > 7         1625      Cy3_1625_46gr_II.txt      Cy5_1625_46gr_II.txt     gr46
 > Ref      46gr_II
 >
 > all 34 rows have Ref in the Cy3 column, in the Cy5 column according to the
 > condition RNA is extracted from on each array
 >
 > Thanks for any help in advance, and sorry for long mail
 >
 > Best regards
 >
 > Ingunn