[BioC] Limma, decideTests

[Gordon K 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>
> Message-ID: <004b01c4fc9e$8a4d98b0$9fb12780 at ihf4651>
> Content-Type: text/plain; format=flowed; charset="iso-8859-1";
> 	reply-type=original
>
> 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??

Neither do I.  Try 'show(fit2)' to see if you actually have any data in the fitted model object.

If you want topTable() and decideTests() to agree, why not use them on the same set of contrasts?

I notice that you've chosen method="global", which is of course is different from the topTable()
method.  It is less powerful if you have a lot of related contrasts.

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