Question

DESeq2 experimental design, 2 controls and 2 treated

0

Entering edit mode

Mike ▴ 10

@mike-12142

Last seen 3.2 years ago

Canada

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

ADD COMMENT • link updated 6.9 years ago by Michael Love 43k • written 6.9 years ago by Mike ▴ 10

score 1 · Answer 1 · 2018-01-20

1

Entering edit mode

Michael Love 43k

@mikelove

Last seen 13 hours ago

United States

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 COMMENT • link 6.9 years ago • updated 6.8 years ago Michael Love 43k

0

Entering edit mode

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 REPLY • link 6.8 years ago Mike ▴ 10

0

Entering edit mode

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

ADD REPLY • link 6.8 years ago Michael Love 43k

0

Entering edit mode

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.

ADD REPLY • link 6.8 years ago Mike ▴ 10

0

Entering edit mode

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 REPLY • link 6.8 years ago Michael Love 43k

0

Entering edit mode

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:

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 REPLY • link 6.8 years ago Mike ▴ 10

0

Entering edit mode

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 REPLY • link 6.8 years ago Michael Love 43k