Entering edit mode
4.4 years ago
Lepomis_8
▴
30
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")