I am working with a complex experiment design shown here. I followed the DESeq2 manual and suggestions in another post to remove the unpaired samples from my analysis and got the DESeq2 running without any error. I need further suggestions to verify that my understanding/interpretations about the comparisons are correct.
Below is the code I used to build the model matrix and perform analysis.
FYI - Variables
samples in the code below correspond to raw counts and the experiment design as shown above.
#Add patients IDs specific to Response r_vec = rep(1:6, each = 2) s_vec = rep(1:9, each = 2) samples$pt.n <- factor(c(r_vec, s_vec)) #build Model Matrix m1 = model.matrix(~ Response + Response:pt.n + Response:Treatment, samples) #Remove All zero columns all.zero <- apply(m1, 2, function(x) all(x==0)) all.zero idx <- which(all.zero) m1 <- m1[,-idx] colnames(m1) dds<- DESeqDataSetFromMatrix(countData = countData, colData = samples, design = m1) dds$Response = relevel(dds$Response, "Sensitive") dds$Treatment = relevel(dds$Treatment, "pre") dds <- DESeq(dds) resultsNames(dds)  "Intercept" "ResponseSensitive"  "ResponseResistant.pt.n2" "ResponseSensitive.pt.n2"  "ResponseResistant.pt.n3" "ResponseSensitive.pt.n3"  "ResponseResistant.pt.n4" "ResponseSensitive.pt.n4"  "ResponseResistant.pt.n5" "ResponseSensitive.pt.n5"  "ResponseResistant.pt.n6" "ResponseSensitive.pt.n6"  "ResponseSensitive.pt.n7" "ResponseSensitive.pt.n8"  "ResponseSensitive.pt.n9" "ResponseResistant.Treatmentpre"  "ResponseSensitive.Treatmentpre"
Can you please confirm that my comparison interpretations below are correct as deduced from the resultsNames?
# Resistant (Treatment) vs Sensitive (Control) res1 = results(dds, contrast=list("ResponseSensitive")) # Post (Treatment) vs Pre (Control) within the group "Sensitive" res1 = results(dds, contrast=list("ResponseSensitive.Treatmentpre")) # Post (Treatment) vs Pre (Control) within the group "Resistant" res1 = results(dds, contrast=list("ResponseResistant.Treatmentpre"))
Next, I also want to compare Post vs Pre considering patients effect. I think, I need to build another model matrix with adding new pt.n column numbered as per the Treatment. Is this correct?
m1 = model.matrix(~ Treatment + Treatment:pt.n + Treatment:Response, samples)
Thank you for the help.