Question: DESeq2 comparison failure
0
gravatar for dec986
4.1 years ago by
dec986200
United States
dec986200 wrote:

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

 

 

rna-seq R • 2.6k views
ADD COMMENTlink modified 4.1 years ago • written 4.1 years ago by dec986200

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

ADD REPLYlink written 4.1 years ago by dec986200

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 4.1 years ago by Devon Ryan92k

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

ADD REPLYlink modified 4.1 years ago • written 4.1 years ago by dec986200

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

ADD REPLYlink written 4.1 years ago by Devon Ryan92k
1
gravatar for Devon Ryan
4.1 years ago by
Devon Ryan92k
Freiburg, Germany
Devon Ryan92k 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 4.1 years ago • written 4.1 years ago by Devon Ryan92k
Please log in to add an answer.

Help
Access

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