Hello,
I have a question regarding DESeq2 and contrasts, and I may simply be over-thinking the problem. I have a dataset with 3 groups, and two treatments
Sample Condition Group
sample1 Control Bact1
sample2 Control Bact1
sample3 Control Bact1
sample4 Treated Bact1
sample5 Treated Bact1
sample6 Treated Bact1
sample7 Control Bact2
sample8 Control Bact2
sample9 Control Bact2
sample10 Treated Bact2
sample11 Treated Bact2
sample12 Treated Bact2
sample13 Control Bact3
sample14 Control Bact3
sample15 Control Bact3
sample16 Treated Bact3
sample17 Treated Bact3
sample18 Treated Bact3
The main goal in the experiment is to determine if strain Bact1 has differentially expressed genes after treatment compared to Bact2. Previously, I did an analysis of Bact1 vs Bact2, (using the table above, minus Bact3) and was set up as:
dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design = ~ group + cond + group:cond)
res<- results(name="groupBact2.condTreated")
# shows interaction term, or how the treatment condition is different across Bact2 v. Bact1
res_Bact2 <- results(contrast=list(c("condition_Treated_vs_Control", "groupBact2.condTreated")))
# shows treatment condition in Bact2
# alternatively I can use dds$group <- relevel(dds$group, ref="Bact2") to look at Bact1
Now, there is a third strain in the analysis, and here is where I'm stuck. I tried creating the table above, using the design
design = ~ group + cond + group:cond
# called DESeq2() with betaprior=FALSE
and it ran successfully without errors, but the resultsNames() table was
[1] "Intercept" "group_BACT2_vs_BACT1" "group_BACT3_vs_BACT1"
[4] "cond_Trt_vs_Ctrl" "groupBACT2.condTrt" "groupBACT3.condTrt"
which doesn't seem to be what I want, because checking the attributes with
attr(dds, "ModelMatrix")
shows that the variables "group_..." doesn't include the treatment variables, and "groupBACT2.condTrt" is only the treated groups (I think?). Contrasting "groupBACT2.condTrt" and "group_BACT2_vs_BACT1" shows an output, but I'm not sure if it's doing what I'd want, i.e. contrasting the difference between Bact3 v. Bact2 in the treatment condition.
If I combine the factors into groups with
dds$grp <- factor(paste0(dds$group, dds$cond))
design(dds) <- ~ grp
it outputs
[1] "Intercept" "grpall_Bact1Trt_vs_Bact1Ctrl"
[3] "grpall_Bact2Ctrl_vs_Bact2Ctrl" "grpall_Bact2Trt_vs_CBact1Ctrl"
[5] "grpall_Bact3Ctrl_vs_Bact1Ctrl" "grpall_Bact3Trt_vs_Bact1Ctrl"
Which is still not what I'm looking for. Am I on the wrong track entirely, or should I simply relevel() the group factors to assign the Treatment as the reference, or something else?
Thanks!
