I am working on a 18-patients dataset. I am trying to calculate DE genes with both EdgeR and DESeq2 for further analyses.
Its a 2-factor design (Status, Fraction). I want to calculate DE genes of different Fractions using the Status as covariate. After DE analysis, logFC are comparable between edgeR and DESeq2, but significant genes (FDR<0.05) are different. EdgeR identifies 500 genes, DESeq2 found 1500. Among those 500 DE genes found with edgeR are ALSO DE with DESeq2. This is a huge problem, since edgeR/DESeq2 are supposed to be superimposable. Any suggestion?
y <- DGEList(counts=counts,group=Fraction) keep <- rowSums(y$counts) >= 10 y$counts <- y$counts[keep,] y <- calcNormFactors(y, method="TMM") design <- model.matrix( ~ Status + Fraction) y <- estimateDisp(y,design) fit <- glmQLFit(y,design) qlf <- glmQLFTest(fit,coef=3)
dds <- DESeqDataSetFromMatrix(countData = counts,colData=summary,design= ~ Status + Fraction) keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep,] dds <- DESeq(dds) res <- results(dds)