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 ?