I am trying to perform differential gene analysis using limma package along with qvalue package on several Affy microarrays. Basically I am comparing two conditions, let's say treated and untreated. I found that after running:
library(limma) library(annaffy) phenoData <- read.AnnotatedDataFrame("phenodata.txt") data <- ReadAffy(celfile.path=".") combinations <- factor(paste(pData(phenoData)[,1], sep = "_")) design <- model.matrix(~combinations) fit <- lmFit(eset, design) fit <- eBayes(fit) table <- topTable(fit,coef=2, 10000, adjust="fdr") library(qvalue) pvalues <- table$P.Value hist(table$P.Value, nclass = 20) qobj <- qvalue(p = table$P.Value, lambda=0) summary(qobj) pi0 <- qobj$pi0 qvalues <- qobj$qvalues localFDR <- lfdr(p = table$P.Value, lambda = 0) write.qvalue(qobj,file="qvalue.txt",sep="\t")
q-values I obtained from "qpackage" and p.adjusted values from "fdr" adjustment for multiple testing result in the same values for hundreds of genes making only 2 genes significant. I strongly believe that is not the actual case, as you can see great fold change for many of the genes. In addition, many p values fall below 0.05 indicating the significance of the genes. I was wondering if anyone could explain to me why this is the output I get and how should I overcome it to produce an accurate list of differentially expressed genes. Thank you so much.