[rozitaa]

hi rozitaa,

The issue here is that there are more animals in certain tissues. So when the R function model.matrix tries to form the model matrix internal to the DESeq() function, this results in a column of zeros for the nested animal 6, tissue A combination for example.

There is an easy solution though. In the next release of Bioconductor, which will be released 17 April 2015, and has DESeq2 v1.8, you can provide a custom model matrix to DESeq(). (You can test the development version of DESeq2 v1.7 already, but this requires downloading the development version of R, and this is not impossible but not always easy.)

So all you need to do is to remove these interaction columns which have no samples:

design <- ~ gut_microbiota + gut_microbiota:animal.nested + gut_microbiota:tissue
mm <- model.matix(design, data=colData(dds))
## remove columns of mm which have all zeros:
(all.zero <- apply(mm, 2, function(x) all(x == 0)))
mm <- mm[, !all.zero]
dds <- DESeq(dds, design=mm, betaPrior=FALSE)

We do have to turn off the log2 fold change shrinkage for this options though (betaPrior=FALSE). But the inference and results tables work the same.

Thanks a lot Micheal! I got the idea!

The only problem that I have now is I am using HTSeq-count data this my code:

design <- ~ gut_microbiota + gut_microbiota:animal.nested + gut_microbiota:tissue
mm <- model.matrix(design, data=sample_table)
(all.zero <- apply(mm, 2, function(x) all(x == 0)))
mm <- as.data.frame(mm[, !all.zero])
dds = DESeqDataSetFromHTSeqCount(sampleTable=sample_table, directory='~/htseq/', design=mm)

But I get an error:

Error in terms.default(object, data = data) : 
  no terms component nor attribute

Is it because of the way I am fitting the design matrix into the htseq-count data?

[Michael Love]

So, again, you'll only be able to do this custom model matrix in version 1.8.

packageVersion("DESeq2")

Additionally, you can only provide the matrix to the 'full' argument of DESeq(). In order to construct a DESeqDataSet object, you can use design=~1, because this will be ignored and the custom matrix will be used instead.

[rozitaa]

So, again, you'll only be able to do this custom model matrix in version 1.8.

packageVersion("DESeq2")

Additionally, you can only provide the matrix to DESeq(). In order to construct a DESeqDataSet object, you can use design=~1, because this will be ignored and the custom matrix will be used instead.

I see!! I cannot wait till April 17th since I have deadline. And it is hard for me to install the developmental version of R! Is there any other solution for my problem?

[Michael Love]

"Now I would like to see DE genes resulted from comparing of every gut_microbiota status (B, C, D, E) to the baseline (A) on each tissue separately."

Sorry, reading over your column data and question more carefully, I realized that this advice about nesting doesn't apply to your experimental design. The nesting approach works for experimental designs where the comparison of interest is nested within individuals/animal. For example if you wanted to compare tissue differences for each microbiota. However, you want to compare microbiota differences for each tissue. The problem for analysis with the nesting approach is that you cannot control for the animal-specific effects and compare across microbiota, because these variables are confounded. Whereas, animal and tissue are not confounded, which allows the nesting approach to work.

I'd recommend you subset each tissue into 3 DESeqDataSets and run a design of ~ gut_microbiota for each.

[rozitaa]

I'd recommend you subset each tissue into 3 DESeqDataSets and run a design of ~ gut_microbiota for each.

Thanks a lot Micheal!

[lbragg]

Also, I had to update that custom line of size factor code for when all rows have a zero. If you were using that line of code to avoid the issue that all rows of your count matrix have one or more zeros, you should get the edited version:

https://support.bioconductor.org/p/62246/

We're implementing a solution here which was on the list of todos for a while now, but it won't be out for 6 months. Typical RNA-Seq doesn't have this issue, but since phyloseq we are seeing more metagenomics datasets which often have no genes without a zero.

Hi there,

I have a metatranscriptomic dataset where every gene has at least one zero count. I tried the geoMeans function you've supplied via the link, but I keep getting the following error:

dds = estimateSizeFactors(dds, geoMeans=geoMeans)
Error in if (any(value <= 0)) { : missing value where TRUE/FALSE needed

This is using DESeq2 1.8.1 on R 3.2.

Happy to provide a MRE if you like.

Lauren

[Michael Love]

what does summary(geoMeans) give you?

[lbragg]

what does summary(geoMeans) give you?

summary(geoMeans)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.414   2.289   2.828   3.731   4.000  28.650

[Michael Love]

Another test (I suspect there is a column with all zeros):

table(colSums(counts(dds)) == 0)

[lbragg]

Another test (I suspect there is a column with all zeros):

table(colSums(counts(dds)) == 0)

You're right, there is one sample with all zeros. Thanks.