No result with DESeq2, but with Cuffdiff 2.1.1

[sindrle]

I have a question, I get no significantly changed genes running DESeq2. I have two groups at two time points (paired) and ran DEseq2 v1.2. DESeq2 test the difference between groups, disregarding time point, difference between time points disregarding groups & finally if one group had another effect between the two timepoints than the other.

I have also used Cuffdiff 2.1.1 testing difference between groups, disregarding time point & difference between time points disregarding groups (Cuffdiff cannot handle the last test).

Cuffdiff found over 1000 diff exp. genes and running GAGE it all made biological sense. So why cant I find anything with DEseq2?

[sindrle]

I haw tweeked abit and the results now looks ok, but differing still from cuffdiff. My study design is however one that cuffdiff cannot handle..

[Simon Anders]

If you told us what you actually did, we might be able to tell you whether it is correct. So, how have you called DESeq2?

[sindrle]

Hi!
With this run I get some genes sig. changed and I guess they are more to be trusted than Cuffdiff results due to correct design?
I will see if EdgeR also gives the same results.

library(DESeq2);
library(GenomicFeatures);
library(htSeqTools);
sampleFiles <- list.files(path="/Volumes/timemachine/HTseq_DEseq2",pattern="*.txt");
sort(sampleFiles,decreasing=TRUE);
status <- factor(c(rep("Diabetic",22), rep("Healthy",26)), levels=c("Diabetic","Healthy"));
timepoints = as.factor(c(1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2));
patients <- factor(c(88,87,86,83,82,81,75,72,71,70,13,08,01,88,87,86,83,82,81,75,72,71,70,13,08,01,79,77,76,73,67,62,61,55,53,21,04,79,77,76,73,67,62,61,55,53,21,04));
sampleTable <- data.frame(sampleName = sampleFiles, fileName = sampleFiles, status=status, timepoints=timepoints, patients=patients);
directory <- c("/Volumes/timemachine/HTseq_DEseq2/");
des <- formula(~patients + timepoints*status);
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= des);
ddsHTSeq;
dds <- estimateSizeFactors(ddsHTSeq);
dds <- estimateDispersions(dds);
dds <- nbinomWaldTest(dds, maxit = 5000)
resultsNames(dds);
statusResults <- results(dds, "status_Healthy_vs_Diabetic");
timepointsResults <- results(dds, "timepoints_2_vs_1");
statusTreatmentResults <- results(dds, "timepoints2.statusHealthy");
write.csv(as.data.frame(statusResults),file="/Volumes/timemachine/HTseq_DEseq2/1_Results/statusResults.csv");
write.csv(as.data.frame(timepointsResults),file="/Volumes/timemachine/HTseq_DEseq2/1_Results/timepointsResults.csv");
write.csv(as.data.frame(statusTreatmentResults),file="/Volumes/timemachine/HTseq_DEseq2/1_Results/statusTreatmentResults.csv");
plotDispEsts(dds)