When running DESeq2 with the same expression data I get different results depending on the order of the columns in my count data matrix. Is this expected behaviour or am I messing something up?

My code looks like this:

countData <- read.table("matrix")
colData <- read.table("samples", header=F)
colnames(colData) <- c("", "genotype", "treatment")
colData$treatment <- relevel(colData$treatment, ref="control")
dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ treatment)

I have two replicates per condition. Depending if my matrix file has the columns in the order control,control,treatment,treatment or treatment,control,treatment,control I get completely different results. I would think that the information from the colData data frame is actually there to make the column order unimportant in countData.

Thanks in advance for any help.

This should not happen, so I'd go over the code carefully. I would make sure that the columns of counts are correctly lined up with the rows of colData. Usually this is taken care of, because we do the counting or read in count files using a column from colData. If you are providing a matrix and a colData which gives the filenames of the BAM files, it's up to you to ensure that these two objects match.

With only four files, you can try to keep track of which is which by total number of reads aligned to genes:

colSums(countData)