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)
#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)
#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")