DESeq2 compare within a condition

Question

This might be a really stupid question, but I can't figure out how to do this. I've read through the DESeq2 vignette and manual pages but couldn't find an answer.

I have a bunch of samples split up into different conditions (eg. celltypes and disease state). I would like to make a comparison between the two possible values within one condition, but only for the samples with a specific value of the other condition.

For example, with the following samplesheet:

patient | phenotype | type
--------+-----------+-----
1       | healthy   | A
1       | healthy   | B
1       | sick      | A
1       | sick      | B
2       | healthy   | A
2       | sick      | A
2       | sick      | B

I would like to compare "healthy" vs "sick", but only for type "A".

Currently, I have the following code, but this will also include the type "B" sample:

dds <- DESeqDataSetFromMatrix(countData=counts, colData=design, design = ~ patient + phenotype + type)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("phenotype", "healthy", "sick"))

Any idea how to accomplish this?

Answer 1

Welcome to the world of all possible pairwise comparisons in DESeq2 !

So the easiest way is to create another group in your data frame, first I simulate some sensible counts:

counts = counts(makeExampleDESeqDataSet(m=48))
design = expand.grid(id=1:3,
phenotype=rep(c("sick","healthy"),each=2),
type = rep(c("A","B"),2))

And you make a group that is a combination of phenotype and type:

design$group = paste(design$phenotype,design$type,sep="_")

Run DESeq2, including id:

dds <- DESeqDataSetFromMatrix(countData=counts, colData=design, 
design = ~ id + group)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("group", "healthy_A", "sick_A"))

So this compares only healthy_A with sick_A.. hopefully this is what you need... The other option (normal statistical) is to fit a model with interaction, but I think this is easier to explain and interpret.

Another question you might have is why not do two separate analysis for A and B, I think it's sometimes better if you model the variance with all available data, given the limited number of replicates.