1
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.