The editor has been updated to markdown! Please see more info at: Tutorial: Updated Support Site Editor

Question: DESeq2 experimental design, 2 controls and 2 treated
gravatar for Mike
13 months ago by
Mike10 wrote:

I have cell type A (untreated vs treated) and cell type B (untreated vs treated), 3 replicates for each:

sample cell_type treatment combined
1 A untreated A_untreated
2 A untreated A_untreated
3 A untreated A_untreated
4 A treated A_treated
5 A treated A_treated
6 A treated A_treated
7 B untreated B_untreated
8 B untreated B_untreated
9 B untreated B_untreated
10 B treated B_treated
11 B treated B_treated
12 B treated B_treated

Within each cell type, to compare the treated samples to their corresponding untreated controls I'm using:

dds <- DESeqDataSet(se, design =~ combined)
de <- DESeq(dds)
results.A <- results(de, contrast=c("combined","A_treated","A_untreated"))
results.B <- results(de, contrast=c("combined","B_treated","B_untreated"))

I think I'm doing this part correctly.

Now I want to know how/if the two cell types respond differently to treatment. I don't think this is right, because it ignores the untreated controls:

results.AB <- results(de, contrast=c("combined","B_treated","A_treated")

I think I need to first normalize (if that's the right term) each treatment to its corresponding control and then compare between cell types. What is the proper way to do this (what design formula and how to extract the results), is it simply this?:

dds <- DESeqDataSet(se, design =~ cell_type + treatment)
de <- DESeq(dds)
results <- results(de)

Thank you.


deseq2 • 327 views
ADD COMMENTlink modified 13 months ago by Michael Love21k • written 13 months ago by Mike10
Answer: DESeq2 experimental design, 2 controls and 2 treated
gravatar for Michael Love
13 months ago by
Michael Love21k
United States
Michael Love21k wrote:

You can use a design of ~celltype + celltype:treatment, then use the following to pull out the comparisons you want:

results(dds, name="celltypeB.treatmenttreated") # individual effects

results(dds, constrast=list("celltypeB.treatmenttreated","celltypeA.treatmenttreated")) # contrast
ADD COMMENTlink modified 13 months ago • written 13 months ago by Michael Love21k

Thank you for the reply. When I use design of ~celltype:treatment

dds <- DESeqDataSet(se, design = ~cell_type:treatment)

I get a matrix is not full rank error. I looked at the "Model matrix not full rank" in the vignette, should I follow the instructions under "Group-specific condition effects, individuals nested within groups"?

ADD REPLYlink written 13 months ago by Mike10

Oops I wrote too quickly. I revised the above answer now.

ADD REPLYlink written 13 months ago by Michael Love21k

Thanks I used the modified command and it worked. I get very few genes changed in the contrast which may be correct (PCA shows the samples are similar) but just want to ensure I'm doing this right:

I first had to relevel treatment:

se$treatment <- relevel(se$treatment, "untreated")

Then analyzed as follows:

dds <- DESeqDataSet(se, design = ~celltype + celltype:treatment)
results.A <- results(dds, name="celltypeA.treatmenttreated")
results.B <- results(dds, name="celltypeB.treatmenttreated")
results.contrast <- results(dds, contrast=list("celltypeB.treatmenttreated","celltypeA.treatmenttreated"))

out of 41972 with nonzero total read count
adjusted p-value < 0.1
  up down outliers lowcounts
results.A 1548, 3.7% 2116, 5% 24, 0.057% 24080, 57%
results.B 2037, 4.9% 2079, 5% 24, 0.057% 22476, 54%
results.contrast 3, 0.0071% 0, 0% 24, 0.057% 0, 0%

There are very few changes in results.contrast.  There are also no lowcounts, is that because of the reason explained here? A: DESeq2 independent filtering fails to filter low counts

Here's an example of one gene I would have expected to be significantly changed (decreased in B.treated) in results.contrast, does it make sense (statistically) that it's not significantly changed?

Sample Normalized counts (rounded)
A.untreated.1 17711
A.untreated.2 20089
A.untreated.3 21664
A.treated.1 14773
A.treated.2 18725
A.treated.3 17930
B.untreated.1 16903
B.untreated.2 19346
B.untreated.3 19460
B.treated.1 9807
B.treated.2 11073
B.treated.3 10019


  log2fc pvalue padj
results.A -0.21 0.07206 0.23238
results.B -0.85 2.948E-13 1.980E-10
results.contrast -0.64 0.00010 0.999854

 Thanks for the help.

ADD REPLYlink written 12 months ago by Mike10

Code all seems correct. P value adjustment depends on the distribution of pvalues and on the filtering, so you can get different behavior per comparison. Some users choose to just pick the same count threshold for all results tables, and turn off independent filtering. You can do this by filtering the dds before you call results() and specifying independentFiltering=FALSE.

ADD REPLYlink written 12 months ago by Michael Love21k

I checked the distribution of pvalues and the results for 'contrast' don't look right.

hist(results.A$pvalue[results.A$baseMean > 1], breaks = 0:20/20, col = "grey50", border = "white")
hist(results.B$pvalue[results.B$baseMean > 1], breaks = 0:20/20, col = "grey50", border = "white")
hist(results.contrast$pvalue[results.contrast$baseMean > 1], breaks = 0:20/20, col = "grey50", border = "white")

results.A, results.B, results.contrast:


I see this is referred to as a hill-shaped histogram, I'll have to read up on this but if you have any immediate advice on what this means or what to try next that'd be appreciated.

I also tried pre-filtering (tried values from 1 to 100), eg. for 10:

dds.filtered <- dds[ rowSums(counts(dds)) > 10, ]

results.A.filtered <- results(dds.filtered, name="celltypeA.treatmenttreated", independentFiltering = FALSE)

results.B.filtered <- results(dds.filtered, name="celltypeB.treatmenttreated", independentFiltering = FALSE)

results.contrast.filtered <- results(dds.filtered, contrast=list("celltypeB.treatmenttreated","celltypeA.treatmenttreated"), independentFiltering = FALSE)

The results for A and B changed somewhat but contrast still showed only the same 3 significantly changed.

ADD REPLYlink modified 12 months ago • written 12 months ago by Mike10

You can play with the pvalues to try to pull out differences but also let the data and current results speak to you: as you said, the groups are similar in a PCA plot. 

ADD REPLYlink written 12 months ago by Michael Love21k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 271 users visited in the last hour