Hi,

I have a dataset with five time points and for each time point I have a set of control samples, a set A of treated samples and a set B of samples with a different treatment (these are not repetitive measures)

Now, what I want is to identify the differentially expressed genes between the control and each treatment at each time point. So at T1 - DE genes between Ctrl and Treatment A + DE genes between Ctrl and Treatment B; then at T2 - DE genes between Ctrl and Treatment A + DE genes between Ctrl and Treatment B; and so on ...

So as a full model I used Code should be placed in three backticks as shown below

ds = phyloseq_to_deseq2(Physeq_no_field, ~ Disturbance + Time.Points +  Disturbance:Time.Points )

and for my reduced model I used

ddsLRT = DESeq(ds, test="LRT", reduced=~Disturbance + Time.Points)

Is that correct? Then I looked at individual comparisons using the Wald test

# for T1 - control vs TreatmentA
resT1_treatA <- results(ddsLRT, name="TreatmentA.Time.PointsT1", test="Wald")

# for T1 - control vs TreatmentB
resT1_treatB <- results(ddsLRT, name="TreatmentB.Time.PointsT1", test="Wald")

# Continue for each individual time point
# ...

Is that correct? Or am I making a huge mistake.

Cheers!

Hi Maria,

That design should be correct! It is the same experimental design as shown in the RNA-seq workflow vignette, part 9: http://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#time-course-experiments.

Note that if you want to additionally test treatment A vs treatment B at time t, you will need to use contrast, like results(ddsLRT, contrast=list("treatA.Tt","treatB.Tt")).

Also, note that running the results function without specifying the name argument would perform a global tests, i.e, identifying genes that show a treatment-specific effect at one or more time points (after time 0).

Kind regards,

Jeroen