Question: top table or anova?
0
gravatar for Konika Chawla
4.8 years ago by
Konika Chawla20 wrote:
Hi I have a microarray data comprised of 3 treatment conditions (N,M,I) and 2 time points (4dpi and 5dpi). I used limma - ebayes fit and toptable with coefficient for a specific contrast to get the differentially expressed genes and adjusted P value. I also did a ANOVA test, and I want to know which one is better to identify differentially expressed genes. CODE: for eBAYES and Toptable design <- model.matrix(~ 0+factor(c(1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6))) colnames(design) <- c("N_4dpi", "M_4dpi", "I_4dpi","N_5dpi", "M_5dpi", "I_5dpi") fit <- lmFit(eset, design) contrast.matrix <- makeContrasts("I_4dpi-M_4dpi","I_5dpi-M_5dpi","I_4dpi-N_4dpi","I_5dpi- N_5dpi","M_4dpi-N_4dpi","M_5dpi-N_5dpi", levels=design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) 2) #If we do one comparison at a time I/M 4dpi tab_all<-topTable(fit2, coef=1,adjust="BH",number=100,p.value=1) head(tab_all) write.table(tab_all,"I-M-4dpi.txt",sep="\t",quote=FALSE) This gives zero significant genes. ############################################################ CODE for ANOVA targets <- readTargets("target_data.txt") data <- ReadAffy(filenames=targets$filename) experiment <- targets$experiment time <- targets$time treatment <- targets$treatment ##normalization eset <- rma(data) exprs.eset <- exprs(eset) exprs.eset.df <- data.frame(exprs.eset) aof <- function(x) { m<-data.frame(time, treatment, x); anova(aov(x ~ time + treatment + time * treatment, m)) } anovaresults <- apply(exprs.eset, 1, aof) This gives over a thousand genes with treatment_Pr.F <0.05 , Does this need to be corrected for multiple testing, why is the huge difference? Appreciate your help. -- With Regards, Konika Chawla NTNU, Norway Phone +4772821344 [[alternative HTML version deleted]]
microarray limma • 783 views
ADD COMMENTlink modified 4.8 years ago by James W. MacDonald50k • written 4.8 years ago by Konika Chawla20
Answer: top table or anova?
0
gravatar for James W. MacDonald
4.8 years ago by
United States
James W. MacDonald50k wrote:
Hi Konika Chawla, Gordon Smyth has published multiple papers that outline why limma is superior to individual tests for low-replication microarray studies. You can find those references in the limma User's Guide (you can type limmaUsersGuide() at an R prompt after loading the limma package). And as you note, this isn't a fair comparison. You don't say how many measurements you are making, but assuming somewhere around 25K genes, you expect 1250 significant results (unadjusted p-value < 0.05) under the null hypothesis. So having that many 'significant' genes isn't saying much. You could use p.adjust() on the vector of p-values to get BH adjusted p-values if you want to make the comparisons more reasonable. Best, Jim On 7/16/2014 1:30 PM, Konika Chawla wrote: > Hi > I have a microarray data comprised of 3 treatment conditions (N,M,I) and > 2 time points (4dpi and 5dpi). > I used limma - ebayes fit and toptable with coefficient for a specific > contrast to get the differentially expressed genes and adjusted P value. > I also did a ANOVA test, and I want to know which one is better to > identify differentially expressed genes. > > CODE: for eBAYES and Toptable > > design <- model.matrix(~ 0+factor(c(1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6))) > colnames(design) <- c("N_4dpi", "M_4dpi", "I_4dpi","N_5dpi", "M_5dpi", > "I_5dpi") > fit <- lmFit(eset, design) > contrast.matrix <- > makeContrasts("I_4dpi-M_4dpi","I_5dpi-M_5dpi","I_4dpi-N_4dpi ","I_5dpi-N_5dpi","M_4dpi-N_4dpi","M_5dpi-N_5dpi", > levels=design) > fit2 <- contrasts.fit(fit, contrast.matrix) > fit2 <- eBayes(fit2) > 2) #If we do one comparison at a time I/M 4dpi > tab_all<-topTable(fit2, coef=1,adjust="BH",number=100,p.value=1) > head(tab_all) > write.table(tab_all,"I-M-4dpi.txt",sep="\t",quote=FALSE) > > This gives zero significant genes. > ############################################################ > > CODE for ANOVA > > > targets <- readTargets("target_data.txt") > data <- ReadAffy(filenames=targets$filename) > experiment <- targets$experiment > time <- targets$time > treatment <- targets$treatment > > > ##normalization > > eset <- rma(data) > exprs.eset <- exprs(eset) > exprs.eset.df <- data.frame(exprs.eset) > aof <- function(x) { > m<-data.frame(time, treatment, x); > anova(aov(x ~ time + treatment + time * treatment, m)) > } > > anovaresults <- apply(exprs.eset, 1, aof) > > This gives over a thousand genes with treatment_Pr.F <0.05 , Does this > need to be corrected for multiple testing, why is the huge difference? > Appreciate your help. > > -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099
ADD COMMENTlink written 4.8 years ago by James W. MacDonald50k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 300 users visited in the last hour