I have a dataset with 3-factor levels (A, B, C). If I use the data of A, B, C to perform a pairwise comparison between A and B, I got a result very different from that when I just use the data of A and B.

Is there a bug in DESeq() of DESeq2?

In the example,

deg_output.tsv has 1857 non-NA padj deg_output2.tsv has 16209 non-NA padj

this is a quite big difference.

The data is here.

https://filebin.net/izu2idggjyorzuof

The metadata is here.

$ cat meta_data.tsv
sample_id   treatment
A0116300    A
A0216455    A
A0316316    A
A0416309    A
B0516310    B
B0616470    B
B0716468    B
B0816450    B
C0916315    C
C1016460    C
C1116453    C
C1216301    C

The code is here.

==> main.R <==
# vim: set noexpandtab tabstop=2:

count_file='count_data.tsv'
sample_file='meta_data.tsv'
condA='A'
condB='B'

read.tsv = function(file, header = T, sep = "\t", quote = "", check.names = F, comment.char = "", stringsAsFactors = F, ...) {
  read.table(file, header = header, sep = sep, quote = quote, check.names = check.names, comment.char = comment.char, stringsAsFactors = stringsAsFactors, ...)
}

sample_table = read.table(sample_file, header=T)
count_table = read.table(count_file, header=T)[, sample_table[[1]], drop=F]
count_table = count_table[rowSums(count_table) > 0, ]
design = as.formula(sprintf('~%s', names(sample_table)[[2]]))

suppressWarnings(suppressPackageStartupMessages(library(DESeq2)))
cds = suppressWarnings(
  DESeqDataSetFromMatrix(
    countData = count_table
    , colData = sample_table
    , design = design
    )
  )

cds = DESeq(cds)
res=data.frame(
  results(cds
    , contrast=c(names(sample_table)[[2]], condA, condB)
    )
  )

write.table(
  cbind(symbol=rownames(res), res)
    , quote = F, sep = "\t", row.names = F
  , file='deg_output.tsv'
  )

==> main2.R <==
# vim: set noexpandtab tabstop=2:

count_file='count_data.tsv'
sample_file='meta_data.tsv'
condA='A'
condB='B'

read.tsv = function(file, header = T, sep = "\t", quote = "", check.names = F, comment.char = "", stringsAsFactors = F, ...) {
  read.table(file, header = header, sep = sep, quote = quote, check.names = check.names, comment.char = comment.char, stringsAsFactors = stringsAsFactors, ...)
}

sample_table = read.table(sample_file, header=T)[1:8, ]
count_table = read.table(count_file, header=T)[, sample_table[[1]], drop=F]
count_table = count_table[rowSums(count_table) > 0, ]
design = as.formula(sprintf('~%s', names(sample_table)[[2]]))

suppressWarnings(suppressPackageStartupMessages(library(DESeq2)))
cds = suppressWarnings(
  DESeqDataSetFromMatrix(
    countData = count_table
    , colData = sample_table
    , design = design
    )
  )

cds = DESeq(cds)
res=data.frame(
  results(cds
    , contrast=c(names(sample_table)[[2]], condA, condB)
    )
  )

write.table(
  cbind(symbol=rownames(res), res)
    , quote = F, sep = "\t", row.names = F
  , file='deg_output2.tsv'
  )

Nope, not a bug. This is a property of linear models in general. This is also a FAQ in the vignette.