Entering edit mode
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}}
