Hi Mike and community,
I’m using theDESeq2 package for my RNA-seq analysis. I would like to get feedback on the choice of my models and applied contrasts. If there is a better model to answer our questions, I would appreciate any suggestions.
I use a 3-level factor analysis:
Cultivar factor with 2 levels (Jonathan, Prima)
Time factor with 3 levels (T1, T2, T3)
Treatment factor with 2 levels (Infected, Control) with 3 biological replicates each.
Experiment: two cultivars have been inoculated witha pathogenic fungus, and sampleswere taken for sequencing at day 5 (time point1), 15 (time point 2), and 30 (time point3). Control samples have been treated in the same manner, but inoculation has been done with water. One of the purpose of control samples was to remove genes that were expressed due to physical cut of the plant.
1) Effect of treatment on cultivar
2) Effect of time on cultivar
Complication: In order to compare the effect of treatment within and across cultivars, we have to reduce infected samples with control samples before making any comparison for each time point and cultivar. Therefore, I have nested cultivar with time point for my first question, and run full model using Wald test where I have used only interaction terms to compare effect of treatment within and between time points for each cultivar.
I have been wondering if using main effect + interaction term could give me the same control over treatment within each cultivar and time point without nesting?
1) Effect of treatment on cultivar:
My reference levels are in alphabetical order. Jonathan_t1, Control.
full~Cultivar + Cultivar:Treatment
# "Intercept" "Cultivar_Jonathan_t2_vs_Jonathan_t1" "Cultivar_Jonathan_t3_vs_Jonathan_t1"
# "Cultivar_Prima_t1_vs_Jonathan_t1" "Cultivar_Prima_t2_vs_Jonathan_t1" "Cultivar_Prima_t3_vs_Jonathan_t1"
# "CultivarJonathan_t1.TreatmentInf" "CultivarJonathan_t2.TreatmentInf" "CultivarJonathan_t3.TreatmentInf"
# "CultivarPrima_t1.TreatmentInf" "CultivarPrima_t2.TreatmentInf" "CultivarPrima_t3.TreatmentInf"
In order to compare across cultivars while controlling for specific time points, I applied this contrast:
contrast = list(c("CultivarJonathan_t1.TreatmentInf","CultivarPrima_t1.TreatmentInf")), test = 'Wald', alpha = 0.01)
name = "CultivarJonathan_t1.TreatmentInf", test = 'Wald', alpha = 0.01
2) Effect of time on cultivar:
I want to detect genes that are differentially expressed due to time within cultivar and across cultivars, I didn’t nest any variables in the model, I have followed this post: DESeq2 time series contrasts between timepoints.
full ~Cultivar + Treatment + Time + Cultivar:Time
reduced ~Cultivar + Treatment + Time, test="LRT"
 "Intercept" "Cultivar_Prima_vs_Jonathan" "Treatment_Inf_vs_Cont" "Time_T2_vs_T1" "Time_T3_vs_T1"
 "CultivarPrima.TimeT2" "CultivarPrima.TimeT3"
SampleName Treatment Time Cultivar
1 Jonathan Inf T1 Jonathan_t1
2 Jonathan Inf T1 Jonathan_t1
3 Jonathan Inf T1 Jonathan_t1
4 Jonathan Cont T1 Jonathan_t1
5 Jonathan Cont T1 Jonathan_t1
6 Jonathan Cont T1 Jonathan_t1
7 Prima Inf T1 Prima_t1
8 Prima Inf T1 Prima_t1
9 Prima Inf T1 Prima_t1
10 Prima Cont T1 Prima_t1
11 Prima Cont T1 Prima_t1
12 Prima Cont T1 Prima_t1
13 Jonathan Inf T2 Jonathan_t2
14 Jonathan Inf T2 Jonathan_t2
15 Jonathan Inf T2 Jonathan_t2
16 Jonathan Cont T2 Jonathan_t2
17 Jonathan Cont T2 Jonathan_t2
18 Jonathan Cont T2 Jonathan_t2
19 Prima Inf T2 Prima_t2
20 Prima Inf T2 Prima_t2
21 Prima Inf T2 Prima_t2
22 Prima Cont T2 Prima_t2
23 Prima Cont T2 Prima_t2
24 Prima Cont T2 Prima_t2
25 Jonathan Inf T3 Jonathan_t3
26 Jonathan Inf T3 Jonathan_t3
27 Jonathan Inf T3 Jonathan_t3
28 Jonathan Cont T3 Jonathan_t3
29 Jonathan Cont T3 Jonathan_t3
30 Jonathan Cont T3 Jonathan_t3
31 Prima Inf T3 Prima_t3
32 Prima Inf T3 Prima_t3
33 Prima Inf T3 Prima_t3
34 Prima Cont T3 Prima_t3
35 Prima Cont T3 Prima_t3
36 Prima Cont T3 Prima_t3