Question: DESeq2 comparison failure
gravatar for dec986
5.0 years ago by
United States
dec986230 wrote:


I have data that looks like this:

transcript - Gene Symbol - Female Control 1 - Female Control 2 - Female Control 3  - Female low 1 - Female low 2 - Female low 3 - Female High 1 - Female High 2 - Female High 3 - Male Control 1 - Male Control 2    Male Control 3 - Male low 1 - Male low 2 - Male low 3  - Male high 1 - Male high  2 - Male high 3

condition <- factor(c(rep("Female Control",3),rep("Female Low",3),rep("Female High",3),rep("Male Control",3),rep("Male Low",3),rep("Male High",3)))

I make the comparisons in DESeq2 by this command, which works

write.table(file="DESeq2/male_control_vs_low_transcripts.tsv",sep="\t",quote=FALSE,results(dds, contrast=c("condition","Male Control","Male Low"),cooksCutoff=FALSE))

however, this fails:

write.table(file="DESeq2/female_vs_male_transcripts.tsv",sep="\t",quote=FALSE,results(dds,cooksCutoff=FALSE, contrast=c("condition","Female","Male")))

Why does the 1st command work but the 2nd command fail?  How can I fix this?

the error given:

Error in rowSums(cts.sub == 0) :
  'x' must be an array of at least two dimensions
Calls: write.table ... cleanContrast -> contrastAllZeroCharacter -> rowSums
Execution halted

doesn't help me at all.

thanks very much,




rna-seq R • 3.0k views
ADD COMMENTlink modified 5.0 years ago • written 5.0 years ago by dec986230

Hi Devon,

thanks very much for your help!  Unfortunately I am having a very hard time setting up the comparisons in DESeq2.  I'm missing some important piece of information somewhere.

How can I do a 2-level comparison, like "male control" vs. "female control"? I don't see how to do this from examples and the manual.


ADD REPLYlink written 5.0 years ago by dec986230

If you want to do comparisons like that then it's typically easiest to make groups like you already did and then use contrasts in a vector form for things like "male" vs. "female". For example, in your original design, the contrast for "male" vs. female is something like c(0,-1,-1,-1,1,1,1)/3.

ADD REPLYlink written 5.0 years ago by Devon Ryan96k

thanks very much for your help, unfortunately I am still getting errors.

My script is so far:

CountTable = read.table("gene.count",sep="\t",row.names=1,header=T)
CountTable <- CountTable[,2:ncol(CountTable)]#remove first column
#need to make CountTable into a DESeq2 object
countData <- as.matrix(CountTable)
condition <- data.frame(sex=rep(c("Female","Male"), each=9), treatment=c("Control","Control","Control","Low","Low","Low","High","High","High","Control","Control","Control","Low","Low","Low","High","High","High"))
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition),~sex+treatment)
#dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)
dds <- DESeq(dds)

#------------------------- Same Condition, different sex

but this gives an error:

In fitNbinomGLMs(objectNZ[fitidx, , drop = FALSE], alpha_hat = alpha_hat[fitidx],  :
  1rows had non-positive estimates of variance for coefficients
Error in cleanContrast(object, contrast, expanded = isExpanded, listValues = listValues,  :
  numeric contrast vector should have one element for every element of 'resultsNames(object)'

and this object (dds) shows:

> resultsNames(dds)
[1] "Intercept"         "sexFemale"         "sexMale"          
[4] "treatmentControl"  "treatmentHigh" "treatmentLow"

But I can't figure out how I can set up the male vs. female comparison based on this, doesn't the comparison have to sum to 0?

My best guess is

contrast = c(0,1,-1,1,0,0)

but I don't see how this focuses on control.

Maybe the error is in

dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition),~sex+treatment)


thanks for any help.  I am totally lost on this.


ADD REPLYlink modified 5.0 years ago • written 5.0 years ago by dec986230

Yeah, like I said, that'll only work in your original paramaterization, not the more traditional factorial design.

ADD REPLYlink written 5.0 years ago by Devon Ryan96k
gravatar for Devon Ryan
5.0 years ago by
Devon Ryan96k
Freiburg, Germany
Devon Ryan96k wrote:

Male and Female aren't factors in your experiment. If you had instead done this:

df <- data.frame(gender=rep(c("Female","Male"), each=9), treatment=rep(c("control","low","high","control","low","high"), each=3))

and then used a model like ~gender+treatment then you'd be able to do that.

ADD COMMENTlink modified 5.0 years ago • written 5.0 years ago by Devon Ryan96k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1585 users visited in the last hour