RNA-Seq Time series analysis in EdgeR

Question

[josmantorres]

Hello Biostars members, I have an RNA-Seq experiment from insect antennae across different data points (0,2,4,6 and 8 days) with 7 replicates per data point (except t6 with 6 replicates). I am analysing using EdgeR to identify differentially expressed genes. I have used the following EdgeR script in my analysis: 1) I created the DGEList object, data normalization and created the design matrix of my experiment

samples_matrix
group <- factor(c('t0','t0','t0','t0','t0','t0','t0','t2','t2','t2','t2','t2','t2','t2','t4','t4','t4','t4','t4','t4','t4','t6','t6','t6','t6','t6','t6','t8','t8','t8','t8','t8','t8','t8'))
y <- DGEList(counts=samples_matrix, group=group)
y <- calcNormFactors(y)
design <- model.matrix(~0+group, data=y$samples)
colnames(design) <- levels(y$samples$group)

The design matrix is composed by 1 and 0 and it has a total of 5 columns (with the different time points) and 34 rows with the samples.

2) I estimated the dispersion and fitted to Quasilikehood ratio F model

y <- estimateDisp(y,design)
fit <- glmQLFit(y,design, robust=TRUE)

3) I defined my contrast. We performed those contrasts that are consecutives (my.constrats1) and all possible combinations (my.contrast2).

my.contrasts1 <- makeContrasts(t2vst0=t2-t0, t4vst2=t4-t2, t6vst4=t6-t4, t8vst6=t8-t6, levels=design)
my.contrasts2 <- makeContrasts(t2vst0=t2-t0, t4vst0=t4-t0, t6vst0=t6-t0, t8vst0=t8-t0, t4vst2=t4-t2, t6vst2=t6-t2, t8vst2=t8-t2, t6vst4=t6-t4, t8vst4=t8-t4, t8vst6=t8-t6, levels=design)

4) I performed the QLMFtest, filtered using HTSfilter package (script not included in this post) and corrections for both approaches. After that, I performed pairwise comparisons between the different contrasts defined in the point 3 (I did not include this part of the analysis in the post).

res1 <- glmQLFTest(fit, contrast=my.contrasts1)
res_est1 <- topTags(res1, n=Inf,  adjust.method="BH", sort.by="PValue")
res2 <- glmQLFTest(fit, contrast=my.contrasts2)
res_est2 <- topTags(res2, n=Inf,  adjust.method="BH", sort.by="PValue")

Both res_est objetcs contains logFc comparisons for the specified contrasts and the corresponding logCPM, F, PValue and FDR values. However, I have some doubts:

1. The information in both "rest_est" objects was "Coefficient: LR test on 4 degrees of freedom", however, the number of comparisons were 4 and 10 and then the degrees of freedom has to be different.
2. When I check rest_est tables, obsereved that logCPM, F, PValue and FDR values are the same in both approaches and then the differentially expressed genes obtained using decideTestsDGE was also the same (5066 not sig and 6120 were sig).

I do not know if I have something wrong in my script or I missed some step in my analysis. Besides, I would like to know which approach would be more correct, because including all the possible time comparisons increases the number of comparisons and it will affect the statistical analysis (more corrections will be necessary) and maybe I lost some genes with small changes.

Thanks a lot for your help

Jose