DESeq2 comparison failure

Question

[dec986]

Hello,

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,

-Dave

[dec986]

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.

-Dave

[Devon Ryan]

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.

[dec986]

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"))

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition),~sex+treatment)
#dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)
dds <- DESeq(dds)

#------------------------- Same Condition, different sex
write.table(file="DESeq2/female_control_vs_male_control_transcripts.tsv",sep="\t",quote=FALSE,results(dds,cooksCutoff=FALSE,
contrast=c(0,1,1,1,0,0,0,0,0,0,-1,-1,-1,0,0,0,0,0,0)/3))

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.

-Dave

[Devon Ryan]

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

Answer 1

[Devon Ryan]

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.