DESeq2 group comparison
Entering edit mode
rbronste ▴ 60
Last seen 4.7 years ago

Quick question about running some 4-way comparisons, have the following set up:

table(sampleSex, sampleTreatment)

dds<-DESeqDataSet(se, design= ~batch + sex + treatment)

dds$group <- factor(paste0(dds$sex, dds$treatment))

design(dds) <- ~0 + group

dds <- dds[ rowSums(counts(dds)) > 1, ]

dds <- DESeq(dds, betaPrior = FALSE)

res1<-results(dds, format = c("GRanges"), contrast=c(1, -1/3, -1/3, -1/3)

res2 <- subset(res1, padj<.001)

res2_sorted = res2[order(res2$padj), ]

In doing it this way I get results that show the first member of the group contrast (1) as having the smallest number of counts in the matrix I am importing as the basis for DESeq2 - however I would like to have the opposite with the first member of the group having the most counts and the remaining members having the fewest. Not sure how to change the design to reflect that, thanks!

deseq2 deseq differential binding analysis • 1.1k views
Entering edit mode
Last seen 2 hours ago
United States

Sorting by (adjusted) pvalue has no preference for the sign of the estimated coefficient.

The fact that you observe more or one sign or the other at the top of the list indicates that these just have smaller pvalue. 

If you want a one sided test, see lfcThreshold in the vignette or in ?results.

Entering edit mode

Thanks for the info. I am fine sorting by the (adjusted) pvalue, Im moreso interested in matching the way I do the 4-way analysis here to the following which is what I am doing for straightforward pairwise comparisons (where I don't do the onesided lfcThreshold):


dds<-DESeqDataSet(se, design= ~batch + treatment)

dds$group <- factor(paste0(dds$treatment))

dds$treatment<-relevel(dds$treatment, ref="vehicle")

design(dds) <- ~batch + group

dds <- dds[ rowSums(counts(dds)) > "1", ]

dds <- DESeq(dds, betaPrior = FALSE)

Me_pw<-results(dds, format = c("GRanges"), contrast=c("group", "vehicle", "EB"))

Me_pw_padj.001 <- subset(Me_pw, padj<.001)

Me_pw_padj.001_sorted = Me_pw_padj.001[order(Me_pw_padj.001$padj), ]

Entering edit mode

Is there a remaining question? I think I answered why you see one or the other sign of LFC at the top of the table sorting by pvalue. If you are only interested in one sign of coefficient, use a one sided test. If you want to look at the genes where the LFC is positive or negative, subset the table.


Login before adding your answer.

Traffic: 659 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6