#### The support.bioconductor.org 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
0
13 months ago by
Mike10
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
modified 13 months ago by Michael Love21k • written 13 months ago by Mike10
Answer: DESeq2 experimental design, 2 controls and 2 treated
1
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

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"?

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

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"))

summary(results.A)
summary(results.B)
summary(results.contrast)
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.

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.

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.