I will try to correct my code again, and I will post it afterwards.

Ok thank you. So I should have the same filtering strategy for all the methods? The fact that DESEQ2 does independent filtering does not affect the filtering strategy?

Thank you for your help! I am posting the scripts used, as well. DESEQ2ff<-list.files(path="/readspergene/", pattern = "*ReadsPerGene.out.tab$", full.names = TRUE ) counts.files <-lapply(ff, read.table, skip=4) counts <- as.data.frame( sapply( counts.files, function(x) x[ , 3 ] ) ) samples<-c("1","2","3","4","5","6","7","8","9","10","11","12") colnames(counts) <- samples library(DESeq2) group<-factor(x=c("A","A","A","B","B","B","B","A","A","B","A", "B"), levels=c("A","B")) colData<-as.data.frame(group) rownames(colData)<-samples countData<-counts colnames(countData)<-samples rownames(countData)<-counts.files[[1]][["V1"]] countData<-as.matrix(countData) dds <- DESeqDataSetFromMatrix(countData=countData, colData=colData, design=~group) keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep,] plotMA(res, ylim=c(-2,2)) resLFC <- lfcShrink(dds, coef="group_B_vs_A", type="apeglm") plotMA(resLFC, ylim=c(-2,2)) dds <- DESeq(dds) res2 <- results( dds, contrast = c("group", "A", "B") ) summary(res2) res2Sig<- subset(res2, padj < 0.05) res2Sig_fold<-subset(res2Sig, abs(log2FoldChange)>1 LIMMAVOOM ff<-list.files(path=="/readspergene/", pattern = "*ReadsPerGene.out.tab", full.names = TRUE ) files<-as.character(list.files(path="/readspergene/")) dge<-readDGE(files, path="/readspergene/", columns=c(1,3), skip=4) colnames(dge)<-c("1","2","3","4","5","6","7","8","9","10","11","12") group<-as.factor(c("A","A","A","B","B","B","B","A","A","B","A", "B")) geneid<-rownames(dge$counts) library(AnnotationDbi) library(org.Hs.eg.db) library(Homo.sapiens) genes<-select(Homo.sapiens,keys=geneid, keytype="ENSEMBL",columns=c("SYMBOL","GENENAME")) genes <- genes[!duplicated(genes$GENENAME),] dge$genes<-genes cpm <- cpm(dge) lcpm <- cpm(dge, log=TRUE) summary(lcpm) group<-as.factor(group) group <- relevel(group, "B") design<-model.matrix(~group) keep.exprs <- filterByExpr(dge,design ) rownames(design)<-colnames(dge) colnames(design)<-gsub("group","",colnames(design)) dge <- dge[keep.exprs,, keep.lib.sizes=FALSE] dge <- calcNormFactors(dge, method = "TMM") plotMDS(lcpm, col= as.numeric(group)) v<-voom(dge, design, plot=TRUE) vfit <- lmFit(v, design) efit <- eBayes(vfit, robust=TRUE) plotSA(efit, main="Final model: Mean-variance trend") summary(decideTests(efit)) DGE.results_limma <- topTable (efit, number = Inf , adjust.method = "BH", sort.by = "P") EDGER ff<-list.files(path="/home/readspergene/", pattern = "*ReadsPerGene.out.tab", full.names = TRUE ) counts.files <-lapply(ff, read.table, skip=4) counts <- as.data.frame( sapply( counts.files, function(x) x[ , 3 ] ) ) samples<-c("1","2","3","4","5","6","7","8","9","10","11","12") colnames(counts) <- samples group<-factor(x=c("A","A","A","B","B","B","B","A","A","B","A","B"), levels=c("A","B")) colData<-as.data.frame(group) rownames(colData)<-samples countData<-counts colnames(countData)<-samples rownames(countData)<-counts.files[[1]][["V1"]] countData<-as.matrix(countData) edgeR.DGElist <- DGEList (counts = countData , group = group) design <- model.matrix(~group) keep <- filterByExpr(edgeR.DGElist, design) edgeR.DGElist <- edgeR.DGElist [ keep ,,keep.lib.sizes=FALSE] edgeR.DGElist <- calcNormFactors ( edgeR.DGElist , method = "TMM" ) plotMDS(edgeR.DGElist) edgeR.DGElist <- estimateDisp ( edgeR.DGElist,design ) edger_fit <- glmFit ( edgeR.DGElist , design ) edger_lrt <- glmLRT (edger_fit) DGE.results_edgeR <- topTags ( edger_lrt , n = Inf, sort.by = "PValue" , adjust.method = "BH") DGE.results_edgeR<-as.data.frame(DGE.results_edgeR) signgenes<-subset(DGE.results_edgeR,FDR<0.05)