[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