Question: DESeq2 complex model matrix and result interpretations
0
5 months ago by
sutturka0
sutturka0 wrote:

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 countData and 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)

[1] "Intercept"                      "ResponseSensitive"
[3] "ResponseResistant.pt.n2"        "ResponseSensitive.pt.n2"
[5] "ResponseResistant.pt.n3"        "ResponseSensitive.pt.n3"
[7] "ResponseResistant.pt.n4"        "ResponseSensitive.pt.n4"
[9] "ResponseResistant.pt.n5"        "ResponseSensitive.pt.n5"
[11] "ResponseResistant.pt.n6"        "ResponseSensitive.pt.n6"
[13] "ResponseSensitive.pt.n7"        "ResponseSensitive.pt.n8"
[15] "ResponseSensitive.pt.n9"        "ResponseResistant.Treatmentpre"
[17] "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.

deseq2 • 144 views
modified 5 months ago by Michael Love25k • written 5 months ago by sutturka0
Answer: DESeq2 complex model matrix and result interpretations
0
5 months ago by
Michael Love25k
United States
Michael Love25k wrote:

A few notes:

The re-factoring has to happen before you run model.matrix. Right now you are getting pre vs post, but if you fix this, it will be post vs pre.

The second and third results tables are correct.

The first results table can't really be formed with this design, because you want to control for patient. You can't control for patient and at the same time compare across the response groups, due to the nesting of individuals in groups.

I don't understand your last comparison of interest? Can you restate what you want to compare?

Thank you for the response. Much appreciated.

(1) I did the relevel before building the model.matrix and I get the post_vs_pre now

resultsNames(dds)

[1] "Intercept"                       "ResponseResistant"
[3] "ResponseSensitive.pt.n2"         "ResponseResistant.pt.n2"
[5] "ResponseSensitive.pt.n3"         "ResponseResistant.pt.n3"
[7] "ResponseSensitive.pt.n4"         "ResponseResistant.pt.n4"
[9] "ResponseSensitive.pt.n5"         "ResponseResistant.pt.n5"
[11] "ResponseSensitive.pt.n6"         "ResponseResistant.pt.n6"
[13] "ResponseSensitive.pt.n7"         "ResponseSensitive.pt.n8"
[15] "ResponseSensitive.pt.n9"         "ResponseSensitive.Treatmentpost"
[17] "ResponseResistant.Treatmentpost"


(2) If I understand you correctly, the results table Resistant vs Sensitive can't be formed with this design.

res1 = results(dds, contrast=list("ResponseSensitive"))


The result table above (after relevel) correspond to Post vs Pre? But then, why the resultName is ResponseSensitive.

(3) I must perform Resistant vs Sensitive comparison. I removed the Treatment variable from the experiment design, and the build the model.matrix as below. I think this will give me the desired comparison with patient effect?

New_sample = select(samples, Response, Patient)

Response    Patient pt.n
Resistant     1      1
Resistant     1      1
Resistant     2      2
Resistant     2      2
Resistant     3      3
Resistant     3      3
Resistant     4      4
Resistant     4      4
Resistant     5      5
Resistant     5      5
Resistant     6      6
Resistant     6      6
Sensitive     7      1
Sensitive     7      1
Sensitive     8      2
Sensitive     8      2
Sensitive     9      3
Sensitive     9      3
Sensitive    10      4
Sensitive    10      4
Sensitive    11      5
Sensitive    11      5
Sensitive    12      6
Sensitive    12      6
Sensitive    13      7
Sensitive    13      7
Sensitive    14      8
Sensitive    14      8
Sensitive    15      9
Sensitive    15      9

m1 = model.matrix(~ Response + Response:pt.n, New_sample)


Thanks.

1

You can't control for the patient and compare across the response, because these are perfectly confounded. If you want to control across the response, remove the patient term, and use ~treatment + response.

Thanks. Can you please answer the last question - what is being compared with:

results(dds, contrast=list("ResponseResistant"))

1

That depends on the design. When you are adding the variables per patient and the patient response in the design, I believe that the resistant vs sensitive coefficient is just for a single pair of samples (whichever happen to be the first alphabetically by their ID). This is not what you want!