Hi all,
I'd like to use DESeq2 for treatment-specific differences over time.
Here is my colData with a common time 0 for Mock, treatment1 and treatment 2:
> as.data.frame(colData(cds))
times treatments
Control_a 00 Mock
Control_b 00 Mock
Control_c 00 Mock
Mock_03hr_a 03 Mock
Mock_03hr_b 03 Mock
Mock_03hr_c 03 Mock
Mock_06hr_a 06 Mock
Mock_06hr_b 06 Mock
Mock_06hr_c 06 Mock
Mock_12hr_a 12 Mock
Mock_12hr_b 12 Mock
Mock_12hr_c 12 Mock
treatment1_03hr_a 03 treatment1
treatment1_03hr_b 03 treatment1
treatment1_03hr_c 03 treatment1
treatment1_06hr_a 06 treatment1
treatment1_06hr_b 06 treatment1
treatment1_06hr_c 06 treatment1
treatment1_12hr_a 12 treatment1
treatment1_12hr_b 12 treatment1
treatment1_12hr_c 12 treatment1
treatment2_03hr_a 03 treatment2
treatment2_03hr_b 03 treatment2
treatment2_03hr_c 03 treatment2
treatment2_06hr_a 06 treatment2
treatment2_06hr_b 06 treatment2
treatment2_06hr_c 06 treatment2
treatment2_12hr_a 12 treatment2
treatment2_12hr_b 12 treatment2
treatment2_12hr_c 12 treatment2
When I tried this:
cds <- DESeqDataSetFromMatrix(countData = data, colData = colData, design = ~ treatments+times+treatments:times)
I got the error message:
Error in checkFullRank(modelMatrix) : the model matrix is not full rank, so the model cannot be fit as specified. One or more variables or interaction terms in the design formula are linear combinations of the others and must be removed.
I've been searching for the reasons and solutions like this post: With DESeq2 "Not full rank" Error with design ~ line + time + condition, and went through its comments. However, since my data is unpaired, using:
colData$times.nested = factor(c(rep(1,3),rep(rep(2:5,each=3),3))) cds <- DESeqDataSetFromMatrix(countData = data, colData = colData, design = ~ treatments + times + treatments:times.nested + treatments:times)
or
mm = model.matrix(~ treatments + times + treatments:times, colData(cds))
cds = DESeq(cds, full=mm, betaPrior=FALSE)
gave me the same error.
How do I capture differentially expressed genes of each treatment vs. Mock over time with common time 0 incorporated?
Best regards,
Ting
Hi Michael,
With your solution, it works well now! Thank you so much for the prompt reply.
Hi Michael,
As I checked the comparisons of the analysis, it showed:
Does that shows there's no statistical comparisons between treatments as you said in A: deseq2 multifactor timeseries help ? Can genes with p.adj<0.05 in "times03.treatmentstreatment2" be interpreted as significant differences between treatment2 and Mock at times03 or just significant differences in treatment2 at times03 over time? If it's the latter case, is there other way to do the differential analysis between one of the treatments and Mock over time including common time 0? I once thought copy-paste the data of common time 0 to each treatment, thus the data set will be paired. But I guess this will change size factors?
Or sacrifice time 0 by removing the common time 0 and then just compare treatments over times03, times06, and times12.
times03.treatmentstreatment2 => treatment 2 vs mock at time 3
To compare treatment 2 vs treatment 1 at time 3, you would use the 'contrast' argument: