Complex Contrast in DESeq2
Entering edit mode
jshouse ▴ 10
Last seen 19 months ago

I have an experimental design with 3 main factor effects. (data = mapped rnaseq counts (raw))

  1. Days of treatment (levels = 1 and 5)
  2. Diet (levels of BD, HFD, and MCD)
  3. Exposure (levels of Vehicle and PERC)

I created a full model where concatenation of each of these:

dds<-DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = input_dir, design=~condition.Group)

Question 1) We are not super interested in interaction effects between all 3 main effects. I went ahead and ran a full model as listed above instead of splitting day1 and day2.  My understanding from a prior post is that this is the recommendation unless one expects the variances to be significantly different from each other on days of treatment =1 vs. days of treatment =2.  Is this still accurate?

Question 2) Main effect contrasts are pretty straight-forward. For example, I used the following to compare the difference between HFD and BD when exposure = Vehicle on Day1  :

d1.hfd.veh <
  results(dds, contrast=c("condition.Group","HFD+Vehicle+1","BD+Vehicle+1"),
​          parallel = TRUE)

 My question regards a more complex contrast.  I want to compare effect of (PERC vs. Vehicle on HFD on Day1) to (PERC vs. Vehicle on BD on Day1).  This is the contrast I wrote after reading this post DESeq2: with multiple factors and interaction terms won't show all effects and I wanted to make sure my understanding is correct. <
results(dds, contrast=
  parallel = TRUE)

Question 3:) In the past, to examine PCA plots, etc..I have take normalized counts from the dds object with normcounts<-counts(dds,normalized=TRUE) and then log transformed them to use in PCA plots. Is using plotPCA on rlog(dds, blind=FALSE) better, similar, or different?

Much appreciated.



deseq2 • 1.9k views
Entering edit mode
Last seen 1 day ago
United States

The best way to set up the model is to ask yourself if it makes sense that there would be interactions between these variables. It sounds like you may want to have an interaction between Exposure and Diet, to test for differences of differences as you've stated in Question 2. And then you could have Days in the model as an added factor, but not interacting with the other two. Does this sound right?

Question 3: If you want a simple log plus pseudocount transformation, you can use normTransform() which produces a DESeqTransform object that you can run plotPCA() on. The idea behind rlog() and vst() are statistical approaches to transforming the data such that the variance is relatively stable across the mean as opposed to what can often be arbitrary decisions about cutoffs and pseudocounts. You can read the DESeq2 paper which talks about this, or look in the DESeq2 vignette section on transformations, and in particular the plots of the SD over the mean for the various transformations.

Entering edit mode

Thanks Michael for the prompt helpful response as always. I'm not sure that I want to combine the information from Day1 and Day5 when examining DE changes for Diet and Exposure. I expect quite a different set of genes to be DE from 1 day of treatment vs. 5 days of treatment, so I didn't want to include Day in the model. I am interested in effect of Diet between Day1 and Day5 for vehicle but again that is a simple contrast. Unless I'm missing something.

For question 2, was the syntax of the contrast I wrote correct to test  "effect of (PERC vs. Vehicle on HFD on Day1) to (PERC vs. Vehicle on BD on Day1)."



Entering edit mode

" I expect quite a different set of genes to be DE from 1 day of treatment vs. 5 days of treatment"

To me, this sounds like you want to use a design with interactions between all variables. You are interested in testing for the effect of Diet on Exposure, and additionally you think that this will change across Days. The interaction model makes it much easier to extract differences of differences such as "(PERC vs. Vehicle on HFD on Day1) vs (PERC vs. Vehicle on BD on Day1)", although it makes extracting the group comparisons a bit more difficult (but not so much). 

If you go ahead and use a design of ~Days*Diet*Exposure, you can report back the resultsNames, and I can explain how to use results to get out the comparisons of interest.

Entering edit mode
ddsfull<-DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = input_dir, design=~condition.Diet*condition.Exposure*condition.Days)
ddsfull <- ddsfull[ rowSums(counts(ddsfull)) > 1, ]
ddsfull$condition.Exposure<-relevel(dds$condition.Exposure, ref = "Vehicle")
ddsfull$condition.Diet<-relevel(dds$condition.Diet, ref = "BD")
ptm <- proc.time()
ddsfull<-DESeq(ddsfull,parallel = TRUE)
proc.time() - ptm

​> resultsNames(ddsfull)
 [1] "Intercept"                                                
 [3] "condition.Diet_MCD_vs_BD"                                 
 [5] "condition.Days_5_vs_1"                                    
 [7] "condition.DietMCD.condition.ExposurePERC"                 
 [9] "condition.DietMCD.condition.Days5"                        
[11] "condition.DietHFD.condition.ExposurePERC.condition.Days5" 
Entering edit mode

There are a lot of "condition." here, which takes up a lot of visual space. If you rename these columns in the colData but remove the prefix "condition.", the output is a lot easier to manage. I'm going to write the following assuming the resultsNames don't contain "condition.", as in:

​> resultsNames(ddsfull)
 [1] "Intercept" "Diet_HFD_vs_BD"                                
 [3] "Diet_MCD_vs_BD" "Exposure_PERC_vs_Vehicle"                      
 [5] "Days_5_vs_1" "DietHFD.ExposurePERC"                
 [7] "DietMCD.ExposurePERC" "DietHFD.Days5"                       
 [9] "DietMCD.Days5" "ExposurePERC.Days5"                  
[11] "DietHFD.ExposurePERC.Days5" "DietMCD.ExposurePERC.Days5"

HFD vs BD for Vehicle on Day 1 is:

results(dds, name="Diet_HFD_vs_BD")

...this is because Vehicle and Day 1 are reference levels.

The contrast (PERC vs. Vehicle on HFD on Day 1) vs (PERC vs. Vehicle on BD on Day 1) is given by the single interaction term:

results(dds, name="DietHFD.ExposurePERC")
Entering edit mode

Thanks Michael. That is much easier. What about comparing the second two levels of diet to each other? MCD vs. HFD.  I would need to set a new reference level and re-run ?  



Entering edit mode

You don't have to re-run, all the possible contrasts are there, but you have to add terms together or contrast them to extract them.

MCD vs HFD for Vehicle on Day 1 is:

results(dds, contrast=list("Diet_MCD_vs_BD","Diet_HFD_vs_BD"))

You can do the math to see how it works out:

(MCD - BD) - (HFD - BD) = MCD - HFD - BD + BD = MCD - HFD

Entering edit mode

Got it! Thanks!


Login before adding your answer.

Traffic: 404 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6