Dear Bioconductor community,

I have timecourse RNA-seq data.
Animals received 4 differents treatments (group = GpA, GpB, GpCD, GpE).
There are 6 animals in each group.
For each animal,blood were collacted before any treatment (D0) and at 3 time points after treatment: D3, D7 and D28.
So, the experiment could be described as follows :

Animal    Time    Group Rep
M657    0    GpA 1
M657    3    GpA 1
M657    7    GpA 1
M657    28    GpA 1
M658    0    GpA 2
M658    3    GpA 2
M658    7    GpA 2
M658    28    GpA 2
.
.
.
.
D627    0    GpB 1
D627    3    GpB 1
D627    7    GpB 1
D627    28    GpB 1
D628    0    GpB 2
D628    3    GpB 2
D628    7    GpB 2
D628    28    GpB 2
.
.
.
.

With rep being a factor value for Animal-Group which distinguishes the individual nested within a group.

I want to determine :
1/ Which genes respond at either the Day3, Day7 or day28 in the different group
for exemple :
GpA : Day3vsDay0 Day7vsDay0, Day28vsDay0
do I have to analyse each Group separatly ? Or I can analyse them all together ?
Because each subject yields a value for each time point, can I take "Animal" effect into my model ? For this, I use block option in duplicateCorrelation().

Here is what I did :

TS <- paste(Design_clean$Group, Design_clean$Time, sep=".")
Design_clean_simple <-data.frame(Group=TS,Replicate=factor(Design_clean$Rep),row.names=row.names(Design_clean))
Blocked_Design_model<-model.matrix(~0+Group,data=Design_clean_simple)
dge <- DGEList(counts=Count)
dge.norm <- calcNormFactors(dge)

v.wts.bl<- voomWithQualityWeights(dge.norm, Blocked_Design_model, plot=TRUE)
corfit.wts.bl <- duplicateCorrelation(v.wts.bl,Blocked_Design_model,block=Design_clean_simple$Replicate)
v.wts.bl.corfit <- voomWithQualityWeights(dge.norm, Blocked_Design_model, plot=TRUE, correlation=corfit.wts.bl$consensus)
fit.wts.bl.corfit <- lmFit(v.wts.bl.corfit, Blocked_Design_model,block=Design_clean_simple$Replicate,correlation=corfit.wts.bl$consensus)

To identify DEG at Day3, Day7 and Day28 when compared to baseline (Day0) :

contrast.matrix.blocked<-makeContrasts(GpA.D3=GroupGpA.3-GroupGpA.0, GpA.D7=GroupGpA.7-GroupGpA.0, GpA.D28=GroupGpA.28-GroupGpA.0
 levels=colnames(Blocked_Design_model))

fit2.wts.bl.corfit <- contrasts.fit(fit.wts.bl.corfit, contrast.matrix.blocked)
fit2.wts.bl.corfit <- eBayes(fit2.wts.bl.corfit)

2/ Which genes respond differently over time in the different groups ie :
Which genes respond differently over time in the GpA relative to the GpB?

So, to answer this question, here is the contrast matrix, as an exemple
Which genes respond differently at Day3 in the GpA relative to the GpB :

contrast.matrix.blocked<-makeContrasts(GpAvsGpB.D3 = (GroupA.3-GroupA.0)-(GroupB.3-GroupB.0),levels=colnames(Blocked_Design_model))
it2.wts.bl.corfit <- contrasts.fit(fit.wts.bl.corfit, contrast.matrix.blocked)
fit2.wts.bl.corfit <- eBayes(fit2.wts.bl.corfit)

Am I doing right ?

This looks fine to me. A few suggestions:

• Make sure you've filtered out low-abundance before running calcNormFactors, etc.
• I usually like running eBayes with robust=TRUE to protect against outliers.
• I would put Animal in the design matrix (and resolve any linear dependencies) rather than blocking on it with duplicateCorrelation. The latter is slightly anticonservative, and you don't need it as you are not comparing directly between different animals. Your comparisons are either within animals (GpA.D3) or between log-fold changes (GpAvsGpB.D3).

Your new code (posted 9 March 2018) looks correct, except that I would insert the line

keep <- filterByExpr(dge, Blocked_Design_model)
dge <- dge[keep,,keep.lib.sizes=FALSE]

before calcNormFactors. That was also Aaron's no. 1 suggestion.

The function filterByExpr was added to edgeR recently, on 26 Dec 2017, so you may need to re-install edgeR to get it.

With the approach taken, your experiment is just a one-way layout, so forming any contrasts between the groups or times is very straightforward, as I think you already appreciate.

Of course, as always, it is good practice to examine plots at each stage of your analysis, and check that the consensus correlation is positive, rather than just running all the code through as a block.