Hello All!

I am attempting to do differential abundance analysis using DESeq2 after doing the initial part of my analysis in QIIME 2. I was able to successfully export a .biom file, added sample metadata and taxonomy information, along with a tree and reference sequences using the import_biom function. I then used the code pasted with DESeq2 package version 1.18.1. However, the output always gave me the OTUID rather than the actually taxonomy... I even checked in phyloseq to see if the taxonomy was added correctly to the .biom file using tax_table(biom_file)[1000:1010, 1:7], and it returned taxonomic information... Any guidance is much appreciated!

subset <- subset_samples(biom_file, variable != "NA")
variabledds = phyloseq_to_deseq2(subset, ~ variable)

#a zero-tolerant variant of geometric mean:
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

geoMeans = apply(counts(variabledds), 1, gm_mean)
variabledds = estimateSizeFactors(variabledds, geoMeans = geoMeans)
variabledds = DESeq(variabledds, fitType="local")

alpha = 0.05

res0 = results(variabledds, contrast = c("variable", "Yes", "No"))
res0 = res0[order(res0$padj, na.last=NA), ]
sigtab0 = res0[(res0$padj < alpha), ]
sigtab0

This is a phyloseq question more than a DESeq2 question, as DESeq2 is not in control of the handoff of information. DESeq2 just receives rownames and count matrices from upstream software. And we don’t have any support for taxonomy or such structure on the rows of a DESeqDataSet, only rownames or what tabular data can go on the rowRanges.