Hello,

I am working with RNA-seq data and using DESeq2 to perform differential expression analysis. My design has a single factor, diagnosis, which is a multi-level factor with six different diagnoses: 'Craniopharyngioma', 'ATRT', 'Ependymoma', 'Glioblastoma', 'Glioma', and 'Medulloblastoma'.

My design is as follows:

dds <- DESeqDataSetFromMatrix(countData = count_data, colData = meta, design = ~ diagnosis)

I am interested in comparing each diagnosis against 'Craniopharyngioma'. To do so, I tried to relevel the 'diagnosis' factor:

dds$diagnosis <- relevel(dds$diagnosis, ref = "Craniopharyngioma")

However, when I run the results() function without specifying any contrast, I get comparisons with 'ATRT' as the base level, as seen in the names of the coefficients from resultsNames(dds):

resultsNames(dds)

Output:

[1] "Intercept"                           "diagnosis_Craniopharyngioma_vs_ATRT" 
[2] "diagnosis_Ependymoma_vs_ATRT"        "diagnosis_Glioblastoma_vs_ATRT"      
[3] "diagnosis_Glioma_vs_ATRT"            "diagnosis_Medulloblastoma_vs_ATRT" 

How can I interpret this result? Does DESeq2 use the relevelled factor in its calculations, but then print the results with the old base level?

I would appreciate if someone could clarify this for me.

Thanks,
Amdom

You can make your life simply and use a contrast for each, for example:

library(DESeq2)

dds <- makeExampleDESeqDataSet()
diagnosis <- c('Craniopharyngioma', 'ATRT', 'Ependymoma', 'Glioblastoma', 'Glioma', 'Medulloblastoma')
dds$diagnosis <- factor(rep(group, each=2))
design(dds) <- ~ diagnosis
dds <- DESeq(dds)

contrasts <- lapply(setdiff(diagnosis, "Craniopharyngioma"), function(x) c("diagnosis", x, "Craniopharyngioma"))
results.list <- lapply(contrasts, function(con) results(dds, contrast = con))