DESEQ2: Which of the following 2 contrast methods does what I want ?

[EVELE]

Hi all,

I have 2 time course experiments.

For the first experiment I have 5 treated samples (1 for the vegetative stage with no replicates,
2 replicates for the early stage of autogamy and
2 replicates for the late stage of autogamy)

For the second experiment I have 6 untreated samples (2 replicates for the vegetative stage,
2 replicates for the early stage of the autogamy cycle and
2 replicates for the late stage of the autogamy cycle)

I am interested in founding 3 sets of differentially expressed genes:

1 set: Treated vs Untreated all samples belonging to the vegetative stage,
2 set: Treated vs Untreated all samples belonging to the early stage of autogamy and
3 set: Treated vs Untreated all samples belonging to the late stage of autogamy

I have applied 2 methods in order to obtain those 3 sets and I get different results.
Apparently the methods are not equivalent but I cannot see the difference.
I would be grateful if anyone could make this difference clear to me.

This part of code is common to both methods
# 0 for vegetative,
# 1 for early autogamy,
# 2 for late autogamy
times=factor(c(rep(1,2), rep(2,2), 0,rep(1,2), rep(2,2), rep(0,2)), levels=c(0,1,2))

# factors for treated samples
# controls for untreated samples
strain = factor(c(rep("factors", 5), rep("controls",6)), levels=c("controls","factors"))

colData = data.frame( strain= strain, times=times)
rownames(colData)=colnames(duo)

1st method
dds <-DESeqDataSetFromMatrix(countData= as.matrix(duo), colData= colData, design=~ strain + times + strain:times)
dds<-DESeq(dds)

resultsNames(dds)
[1] "Intercept" "straincontrols" "strainfactors"
[4] "times0" "times1" "times2"
[7] "straincontrols.times0" "strainfactors.times0" "straincontrols.times1"
[10] "strainfactors.times1" "straincontrols.times2" "strainfactors.times2"

# In order to obtain the list of differentially expressed genes for the vegetative stage
resVEG <- results(dds, contrast=list("straincontrols.times0","strainfactors.times0"))
resVEG<- resVEG[!is.na(resVEG$padj) & resVEG$padj <= 0.05,]
list_of_differentially_expressed_genes_for_VEG<-rownames(resVEG)
length(list_of_differentially_expressed_genes_for_VEG)
822

# In order to obtain the list of differentially expressed genes for the early stage of autogamy
resEARLY <- results(dds, contrast=list("straincontrols.times1","strainfactors.times1"))
resEARLY <- resEARLY[!is.na(resEARLY$padj) & resEARLY$padj <= 0.05,]
list_of_differentially_expressed_genes_for_EARLY<-rownames(resEARLY)
length(list_of_differentially_expressed_genes_for_EARLY)
461

# In order to obtain the list of differentially expressed genes for the late stage of autogamy
resLATE <- results(dds, contrast= list("straincontrols.times2","strainfactors.times2"))
resLATE <- resLATE[!is.na(resLATE$padj) & resLATE$padj <= 0.05,]
list_of_differentially_expressed_genes_for_LATE<-rownames(resLATE)
length(list_of_differentially_expressed_genes_LATE)
924

2nd method
colData$group <- factor(paste0(colData$strain,colData$times))
colData
strain times group
ND7.T5 factors 1 factors1
ND7.T10 factors 1 factors1
ND7.T30 factors 2 factors2
ND7.T40 factors 2 factors2
ND7.V factors 0 factors0
T0.2 controls 1 controls1
T5.2 controls 1 controls1
T20.3 controls 2 controls2
T30 controls 2 controls2
V1.2 controls 0 controls0
VK1 controls 0 controls0

dds<-DESeqDataSetFromMatrix(countData= as.matrix(duo), colData= colData, design=~group)
dds<-DESeq(dds)

# In order to obtain the list of differentially expressed genes for the vegetative stage
resVEG<-results(dds, contrast=c("group", "controls0", "factors0"))
resVEG<- resVEG[!is.na(resVEG$padj) & resVEG$padj <= 0.05,]
list_of_differentially_expressed_genes_for_VEG<-rownames(resVEG)
length(list_of_differentially_expressed_genes_for_VEG)
3500

# In order to obtain the list of differentially expressed genes for the early stage of autogamy
resEARLY <- results(dds, contrast=c("group", "controls1", "factors1"))
resEARLY <- resEARLY[!is.na(resEARLY$padj) & resEARLY$padj <= 0.05,]
list_of_differentially_expressed_genes_for_EARLY<-rownames(resEARLY)
length(list_of_differentially_expressed_genes_for_EARLY)
4156

# In order to obtain the list of differentially expressed genes for the late stage of autogamy
resLATE <- results(dds, contrast= c("group", "controls2", "factors2"))
resLATE <- resLATE[!is.na(resLATE$padj) & resLATE$padj <= 0.05,]
list_of_differentially_expressed_genes_for_LATE<-rownames(resLATE)
length(list_of_differentially_expressed_genes_LATE)
3959

Thank you in advance !