Hi all,

I have RNA-seq data from three groups (rapid, slow, and control) at three timepoints (T1, T2, T3). Each group has three replicates except for the control group. After performing the time course experiments in DESeq2, I got over 1000 DEGs.

Now, my simple question is:

Are those DEGs (1) behaving in a condition-specific way (aka rapid is expressed differentially from slow since T1 + rapid and slow are both differ from control) or (2) those DEGs just behaving differently from the control but there's no different between rapid and slow groups.

Here is my script:

countdata=read.table("all.genes.txt", sep=" ", header=TRUE, row.names=1)
condition <- factor(c(rep("control",2), rep("rapid",3), rep("slow",3),
                      rep("control",2), rep("rapid",3), rep("slow",3),
                      rep("control",2), rep("rapid",3), rep("slow",3)))
timepoints <- factor(c(rep("t1",8), rep("t2",8), rep("t3",8)))
sampleTable <- data.frame(condition = as.factor(condition),
                          timepoints = as.factor(timepoints))
d.deseq<-DESeqDataSetFromMatrix(countData = countdata,
                                colData = sampleTable,
                                design = ~condition+timepoints+condition:timepoints)
ddsTC<-DESeq(d.deseq, test="LRT", reduced = ~ condition+timepoints)
ddsTC$condition <- relevel(ddsTC$condition, ref="control")
     resTC<-results(ddsTC, alpha = 0.05)

"control" is already the reference level of condition thanks to alphabetical sorting, but if it were not you would need to relevel before calling DESeq().

These DEG should exhibit changes across time which are condition specific: either rapid or slow or both exhibit a different expression profile over time compared to control.