Question

DESeq2 complex model matrix and result interpretations

0

Entering edit mode

sutturka • 0

@sutturka-14580

Last seen 20 months ago

United States

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 • 2.6k views

ADD COMMENT • link updated 5.6 years ago by Michael Love 43k • written 5.6 years ago by sutturka • 0

score 0 · Answer 1 · 2019-05-16

0

Entering edit mode

Michael Love 43k

@mikelove

Last seen 15 hours ago

United States

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?

ADD COMMENT • link 5.6 years ago Michael Love 43k

0

Entering edit mode

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)

head(New_sample)

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.

ADD REPLY • link 5.6 years ago sutturka • 0

1

Entering edit mode

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.

ADD REPLY • link 5.6 years ago Michael Love 43k

0

Entering edit mode

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

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

ADD REPLY • link 5.6 years ago sutturka • 0

1

Entering edit mode

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!

ADD REPLY • link 5.6 years ago Michael Love 43k