I'm trying to analyze RNA-Seq data for an experiment which tests 4 groups (same genotype) against 4 different conditions, and am interested in identifying differentially expressed genes among the conditions. I've gone through many different example scripts but I'm not sure that I'm doing it right! Each script is different and gives me different results (which mostly give ~12000 DEGs). I just want to get DEG's among treatments at an FDR of 0.05 and fold change of 2 or -2. Any help is appreciated!! I don't know if any or all steps are necessary or correct.

# Assemble expression set
eset=ExpressionSet(assayData = exprs,phenoData = pData,experimentData = experimentdata)

edgeR dataset:

    #create dataset
    dge=DGEList(counts=exprs(eset), group=pData(eset)$Condition

    #normalize by total count
    dge=calcNormFactors(dge)

    #estimate dispersion parameter for GLM
    dge=estimateGLMTrendedDisp(dge, method="power")

    #specify design matrix 
    design.mat=model.matrix(~ 0 + dge$samples$group)
    colnames(design.mat)=c("Green", "Yellow", "Blue","Control")

    # model fit
    fit.edgeR=glmFit(dge, design.mat)

    # Differential expression
    contrasts.edgeR=makeContrasts(Green-Yellow,Green-Blue,Green-Control,Yellow-Blue,Yellow-Control,Blue-Control, levels=design.mat)

    #Likelihood ratio test
    lrt.edgeR=glmLRT(fit.edgeR, contrast=contrasts.edgeR)

    #access results table
    edgeR_results=lrt.edgeR$table
    sig.edgeR=decideTestsDGE(lrt.edgeR, adjust.method="BH", p.value = p.threshold, lfc=2)

DESeq dataset:

#create dataset
dds=DESeqDataSetFromMatrix(countData = exprs(eset), colData = pData(eset),design = ~ Condition ) 

#run DESeq
dds=DESeq(dds)

#view results
res=results(dds, name="Condition_Blue_vs_Control", alpha = 0.05, lfcThreshold = 2 ,altHypothesis="greaterAbs")
res=results(dds, name="Condition_Yellow_vs_Control", alpha = 0.05, lfcThreshold = 2 ,altHypothesis="greaterAbs")
res=results(dds, name="Condition_Green_vs_Control", alpha = 0.05, lfcThreshold = 2 ,altHypothesis="greaterAbs")