DESeq2: struggle to add multiple variables to DESeqDataSet

[JonB]

Hi,

I want to add the two variables condition and treatment to my DESeqDataSet, and I thought I had succeeded:

> ddsHTSeq
class: DESeqDataSet
dim: 21680 11
exptData(0):
assays(1): counts
rownames(21680): scign000001 scign000002 ... scign021679 scign021680
rowData metadata column names(0):
colnames(11): R1D0B1 R1D2B1 ... R3DxB1 R3DzB1
colData names(2): condition treatment

I added four levels to condition and two levels to treatment and I am able to pull out DESeq results from the different levels of condition, but I can't pull out the results from treatment:

> resTreatment = results(ddsHTSeq, contrast=c("treatment", "cut", "uncut"))

Any ideas what I am doing wrong?

[Michael Love]

hi Jon,

Can you include the code you used to create ddsHTSeq?

[JonB]

Here's the code I used:

library("DESeq2")

#Set the directory where the HTSeq count files are located
directory = "/Users/jonbra/Desktop/Gene_count_files/"

#Catch the files with the pattern "GFF"
sampleFiles = grep("GFF", list.files(directory), value=TRUE)
sampleFiles

#Setup sample conditions in the same order of the files in sampleFiles
sampleCondition = c("D0", "D2", "Dx", "Dz", "D0", "D2", "Dx", "Dz", "D0", "Dx", "Dz")
sampleName = c("R1D0B1", "R1D2B1", "R1DxB1", "R1DzB1", "R2D0B1", "R2D2B1", "R2DxB1", "R2DzB1", "R3D0B1", "R3DxB1", "R3DzB1")
sampleTreatment = c("uncut", "cut", "cut", "cut", "uncut", "cut", "cut", "cut", "uncut", "cut", "cut")
sampleTable = data.frame(sampleName=sampleName, fileName=sampleFiles, condition=sampleCondition, treatment=sampleTreatment)
sampleTable

#Creates the DESeqDataSet which stores the count data for DESeq2
ddsHTSeq = DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design= ~ condition)
ddsHTSeq

colData(ddsHTSeq)$condition=factor(colData(ddsHTSeq)$condition, levels=c("D0", "D2", "Dx", "Dz"))
colData(ddsHTSeq)$treatment=factor(colData(ddsHTSeq)$treatment, levels=c("uncut", "cut"))

I tried to use this design formula instead:

ddsHTSeq = DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design= ~ treatment + condition)

But got this error message:

Error in DESeqDataSet(se, design = design, ignoreRank) :
the model matrix is not full rank, i.e. one or more variables in the design formula are linear combinations of the others

[Michael Love]

The problem you will have with modeling ~ treatment + condition is that D0 is indistinguishable from uncut, i.e. all samples which are D0 are also uncut and vice versa. That makes the model impossible to fit as specified.

If you leave out the treatment variable, you can make comparisons of each level against the other. And you can also make the comparison of D2,Dx,Dz vs D0 (so the treatment contrast you are interested) using the list argument to the contrast argument:

results(dds, 
  contrast=list(c("conditionD2","conditionDx","conditionDz"),"conditionD0"), 
  listValues=c(1/3,-1))

[JonB]

Ah, I see. That makes sense. Thanks!