I have 3 sample groups "Group": White_lion (n=25), Black_lion (n=5), White_tiger (n=28) with taxa groups.

When I run the phyloseq to DEseq2 using

diagdds = phyloseq_to_deseq2(physeq, group="letters")
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(diagdds), 1, gm_mean)
diagdds = estimateSizeFactors(diagdds, geoMeans = geoMeans)
diagdds = DESeq(diagdds, fitType="local")

res = results(diagdds)
res = res[order(res$padj, na.last=NA), ]
alpha = 0.01

sum(res$padj < alpha, na.rm=TRUE)

sigtab = res[(res$padj < alpha), ]
# Add taxonomic labels for plotting
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(physeq2)[rownames(sigtab), ], "matrix"))

dim(sigtab)

I obtain 17 differentially abundant taxons (I assume this takes the first group alphabetically, and compares to the other two independently of their group size? -- is this correct?)

But I am really interested in to compare A to B, and A to C. (White lion to black lion; and white lion to white tiger). How can I include two-group test within my group comparison?

• Should I run separately these by creating a new phyloseq only containing these features?
• Could I maybe run something like this

  diagdds = phyloseq_to_deseq2(physeq2, ~ Animal + Colour)
  

  I end up with 23 significant features,

How can I obtain differentially abundant features between paired-groups?

Thanks

You don't have to assume how the functions work, just read the documentation. The vignette discusses how results works, as does ?results. In particular see the contrast argument.

In addition, I'm skeptical of the use of DESeq2 testing for microbiome/metagenomics. So I won't comment on whether more or less DA features is indicating anything meaningful.