I've used DESeq2 previously on my 16S microbiota data, where I ran the phyloseqtodeseq function then went straight on to the deseq function, e.g.
MERGEDnonsinglet.age.dds <- phyloseqtodeseq2(MERGEDnonsinglet.age, ~ Ageclass) MERGEDnonsinglet.age.dds2 = DESeq(MERGED_nonsinglet.age.dds1, test="Wald", fitType="parametric")
...where 'MERGED_nonsinglet.age' is my phyloseq object with singletons removed (but not rarefied).
However, with recent updates to the package (V1.22.2 I think) I can no longer do this as I get the following error:
estimating size factors Error in estimateSizeFactorsForMatrix(counts(object), locfunc = locfunc, : every gene contains at least one zero, cannot compute log geometric means
I've managed to get the script running but only if I include the geomeans function:
MERGEDnonsinglet.age.dds <- phyloseqtodeseq2(MERGEDnonsinglet.age, ~ Age_class)
gmmean = function(x, na.rm=FALSE){ exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))} geoMeans = apply(counts(MERGEDnonsinglet.age.dds), 1, gm_mean)
MERGEDnonsinglet.age.dds1 = estimateSizeFactors(MERGEDnonsinglet.age.dds, geoMeans = geoMeans) MERGEDnonsinglet.age.dds2 = DESeq(MERGEDnonsinglet.age.dds1, test="Wald", fitType="parametric")
MERGEDnonsinglet.age.dds.res = results(MERGEDnonsinglet.age.dds2, cooksCutoff = F) View(MERGED_nonsinglet.age.dds.res)
...and I've read conflicting information on whether or not you should transform your data. Also the output results are a bit weird. If I look at differences between, for example, males and females, then do the same for adults vs juveniles...3 ASVs appear significant with sex and 2 ASVs with age class, but those 2 ASVs for age class are the same as those for sex - if that makes sense? Additionally, those ASVs are not what I would expect to come out - I can see some ASVs that clearly have a hugely different abundance that aren't picked up.
After I looked in to this more, I've seen that there is some debate over whether or not you should transform the data using geomeans so I'm unsure how to proceed. Is it this geomeans step that could be giving these weird results? If yes, what should I be doing differently to make it work?
Thanks for any advice you may have, Sophie