[BioC] Limma :: decideTest question

Gordon K Smyth smyth at wehi.EDU.AU
Sun Oct 17 00:50:28 CEST 2010


Dear Lorenzo,

limma uses all the arrays in your experiment to estimate the residual 
standard errors, not just the arrays that you specify in your contrast. 
In your second analysis, it appears that you are analysing only a subset 
of the arrays.  If you do this, all the results will change, because the 
number of residual degrees of freedom uses to estimate the standard errors 
is reduced.  Normally I recommend that you keep all the data together when 
using limma, as in your first analysis, not analyse subsets of arrays 
piece-meal, as in your second analysis.  Generally, analysing all the 
arrays together in one lmFit give you more power to detect differential 
expression.

Best wishes
Gordon

> Date: Thu, 14 Oct 2010 15:08:12 +0100
> From: Lorenzo Bomba <lory.bomb at gmail.com>
> To: bioconductor at stat.math.ethz.ch
> Subject: [BioC] Limma :: decideTest question
> Message-ID: <E1905B89-3ED3-484D-88AD-DA891066C790 at gmail.com>
> Content-Type: text/plain; charset="us-ascii"
>
> Hi all,
>
> I would like to ask a question about the function decideTest:
>
> I'm analysing microarray data and I'm using a the following commands to get the differentially express genes:
>
> #RDE is a  MAlist object
>
> target<-read.table(file="TargetALL.txt", header=TRUE)
> TS <- factor(target$Diet, levels=c("CTR0","CTR","TRT"))
> design <- model.matrix(~0+TS)
> colnames(design) <- c("CTR0","CTR","TRT")
> fit3_1 <- lmFit(RDE,design,weights=RDE$weights)
> cont.matrix <- makeContrasts(CTR0vsCTR = CTR0-CTR, CTR0vsTRT= CTR0-TRT, CTRvsTRT= CTR-TRT,levels=design)
> fit3_1c <- contrasts.fit(fit3_1, cont.matrix)
> fit3 <- eBayes(fit3_1c)
> TableDEGfit3T <- topTable(fit3, adjust="BH", number= 29767)
>
> # follow the "TargetALL.txt" file that I use to create the design matrix
> FileName	Diet
> S1M	CTR0
> S2	CTR0
> S3	CTR0
> S4	CTR0
> S6	CTR0
> S10M	CTR
> S11	CTR
> S12	CTR
> S13	CTR
> S14	CTR
> S15	CTR
> S17rt	TRT
> S18	TRT
> S19	TRT
> S20	TRT
> S21	TRT
> S22	TRT
>
> #after I make a SEPARATE decide test I got the number of gene in common between CTR0vsCTR and CTR0vsTRT contrast like in the venn #diagramm following
>
> results <- decideTests(fit3,p.value=0.01)
> #Venn Diagrams
> jpeg("Venndigramsnocov.jpg")
> vennDiagram(results)
> dev.off()
>
> -------------- next part --------------
>
> after that I tried to do separate analysis of the two control using this commands
>
> #CONTRAST:  CTR0vsCTR
> target<-read.table(file="TargetALL.txt", header=TRUE)
> TS <-factor(target$Diet,levels=c("CTR0","CTR"))
> design <- model.matrix(~-1+TS)
> colnames(design) <- c("CTR0","CTR")
> fit4_1 <- lmFit(RDEctr0ctr,design,weights=RDEctr0ctr$weights)
> cont.matrix <- makeContrasts(CTR0vsCTR= CTR0-CTR,levels=design)
> fit4_1c <- contrasts.fit(fit4_1, cont.matrix)
> fit4 <- eBayes(fit4_1c)
> TableDEGfit4T <- topTable(fit4, adjust="BH", number= 29767)
>
>
> #CONTRAST:  CTR0vsTRT
> target<-read.table(file="TargetALL.txt", header=TRUE)
> TS <-factor(target$Diet,levels=c("CTR0","TRT"))
> design <- model.matrix(~-1+TS)
> colnames(design) <- c("CTR0","TRT")
> fit5_1 <- lmFit(RDEctr0ctr,design,weights=RDEctr0ctr$weights)
> cont.matrix <- makeContrasts(CTR0vsTRT= CTR0-TRT,levels=design)
> fit5_1c <- contrasts.fit(fit5_1, cont.matrix)
> fit5 <- eBayes(fit5_1c)
> TableDEGfit5T <- topTable(fit5, adjust="BH", number= 29767)
>
> and than I tried to figured out what are the gene in common that I've suppose to be the same number of the vennDiagram: 6472 genes in common
>
> p001_4<-which(TableDEGfit4T$adj.P.Val <= 0.01)
> p001_5<-which(TableDEGfit5T$adj.P.Val <= 0.01)
> Descr4_001<-TableDEGfit4T$Description[p001_4]
> Descr5_001<-TableDEGfit5T$Description[p001_5]
> com <- match(Descr4_001, Descr5_001)
> length(which(com != "NA"))
>
> after This I got this number 4176 instead 6472
> I cannot understand because in the decideTest the contrast are take separately .........
> Thanks in advance for the help!
>
> Lorenzo Bomba

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}



More information about the Bioconductor mailing list