**60**wrote:

Hello, I'm having a hard time understanding the time course experimental work flow in DESeq2. I am using the F1000 paper as a template (Love MI, Anders S, Kim V and Huber W. RNA-Seq workflow: gene-level exploratory analysis and differential expression [version 1; referees: 2 approved]. F1000Research 2015, 4:1070 (doi: 10.12688/f1000research.7035.1))

Experiment: I have 2 genotypes, 3 timepoints per genotype. 4 replicates in each genotype-condition pair (3 in one). I have a counts table generated in HTSeq - 23 samples.

```
coldata <- data.frame(sample=colnames(data),
genotype=as.factor(c(rep("het",4), rep("wt",4), rep("het",4), rep("wt",4),rep("het",3), rep("wt",4))),
time=as.factor(c(rep("es",8), rep("np",8),rep("nsc",7))))
```

sample genotype time
1 es1 het es

2 es2 het es

3 es3 het es

4 es4 het es

5 es5 wt es

6 es6 wt es

7 es7 wt es

8 es8 wt es

9 np1 het np

10 np2 het np

11 np3 het np

12 np4 het np

13 np5 wt np

14 np6 wt np

15 np7 wt np

16 np8 wt np

17 nsc1 het nsc

18 nsc2 het nsc

19 nsc4 het nsc

20 nsc5 wt nsc

21 nsc6 wt nsc

22 nsc7 wt nsc

23 nsc8 wt nsc

I have tried it 2 different ways, with slightly different results. My questions are: 1. What is the difference between the two? I have used both elements as factors. 2. It's not clear how the relevel is worning under these conditions 3. Using the same results table, can I compare the intermediate timepoints, i.e, 1st and 2nd then 2nd and 3rd, generating a esults table for each, or would I have to separate the analysis entirely?

Thanks!

Here goes what I have done:

```
dds <- DESeqDataSetFromMatrix(countData=data, colData=coldata, design=~ time + genotype:time + genotype)
colData(dds)$genotype <- relevel(colData(dds)$genotype, "wt")
ddsTC <- DESeqDataSetFromMatrix(countData=data, colData=coldata,
design=~ time + genotype:time + genotype)
ddsTC <- DESeq(ddsTC, test="LRT", reduced = ~ time + genotype)
resTC <- results(ddsTC)
resTC$symbol <- mcols(ddsTC)$symbol
head(resTC[order(resTC$padj),],4)
```

Results: log2 fold change (MLE): timensc.genotypewt LRT p-value: '~ time + genotype:time + genotype' vs '~ time + genotype' DataFrame with 4 rows and 6 columns

```
summary(resTC)
```

out of 29372 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 9, 0.031% LFC < 0 (down) : 20, 0.068% outliers [1] : 6, 0.02% low counts [2] : 9434, 32% (mean count < 3) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results

```
resultsNames(ddsTC)
```

[1] "Intercept" "genotype_wt_vs_het" "condition_np_vs_es"

[4] "condition_nsc_vs_es" "genotypewt.conditionnp" "genotypewt.conditionnsc"

Second iteration:

```
resTC2 <- results(ddsTC, contrast=c("genotype","het","wt"))
head(resTC2[order(resTC2$padj),],4)
```

log2 fold change (MLE): genotype het vs wt LRT p-value: '~ time + genotype:time + genotype' vs '~ time + genotype' DataFrame with 4 rows and 6 columns

```
summary(resTC2)
```

out of 29372 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 16, 0.054% LFC < 0 (down) : 13, 0.044% outliers [1] : 6, 0.02% low counts [2] : 9434, 32% (mean count < 3)

```
> sessionInfo()
```

R version 3.2.2 (2015-08-14) Platform: x86_64-apple-darwin13.4.0 (64-bit) Running under: OS X 10.10.5 (Yosemite)

attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets [8] methods base

other attached packages:
[1] sva_3.18.0 genefilter_1.52.1

[3] mgcv_1.8-12 nlme_3.1-126

[5] PoiClaClu_1.0.2 hexbin_1.27.1

[7] RColorBrewer_1.1-2 vsn_3.38.0

[9] ggplot2_2.1.0 pheatmap_1.0.8

[11] DESeq2_1.10.1 RcppArmadillo_0.6.100.0.0
[13] Rcpp_0.12.4 SummarizedExperiment_1.0.2
[15] Biobase_2.30.0 GenomicRanges_1.22.4

[17] GenomeInfoDb_1.6.3 IRanges_2.4.8

[19] S4Vectors_0.8.11 BiocGenerics_0.16.1

[21] maptools_0.8-39 sp_1.2-2