function or argument to subset taxa for plot_cladogram

[ethanbass]

Hi,

I was wondering if it would be possible to include an option for subsetting taxa after running a differential expression analysis. In the example for plot_cladogram you subset out the Firmicutes before running the analysis, but what if we want to plot only the Firmicutes from an analysis that was run on the whole dataset? Or, to give a more tractable example, is there a way to extract only the Clostridiales from the mm_lefse object created in the example?

Example code:

data(kostic_crc)

kostic_crc_small <- phyloseq::subset_taxa(
  kostic_crc,
  Phylum %in% c("Firmicutes")
)

mm_lefse <- run_lefse(
  kostic_crc_small,
  wilcoxon_cutoff = 0.01,
  group = "DIAGNOSIS",
  kw_cutoff = 0.01,
  multigrp_strat = TRUE,
  lda_cutoff = 4
)

plot_cladogram(mm_lefse, color = c("darkgreen", "red"))

Here's my attempt to subset the phyloseq object and plot the resulting tree (without the annotations):

ps <- microbiomeMarker:::create_ps_from_mm(mm_lefse, only_marker = FALSE)
tax <- as.data.frame(stringr::str_split_fixed(ps@tax_table@.Data,"\\|",n=7))
colnames(tax) <- microbiomeMarker:::available_ranks
ps <- phyloseq::subset_taxa(ps,  tax$Order == "o__Clostridiales")

microbiomeMarker:::get_treedata_phyloseq(ps) %>%
  microbiomeMarker:::generate_taxa_tree(size = branch_size)

but I get the error

Error in child_list[[i]] : 
  attempt to select less than one element in integerOneIndex
In addition: Warning message:
In y[tip.idx] <- 1:Ntip * step :
  number of items to replace is not a multiple of replacement length"

I'd be grateful for any help!

[ethanbass]

I figured out that part of the problem is that I'm wiping out all of the nodes further up in the tree.

This yields better results by retaining the nodes with no order (but the results still look a little wonky):

ps <- microbiomeMarker:::create_ps_from_mm(mm_lefse, only_marker = FALSE)
tax <- as.data.frame(stringr::str_split_fixed(ps@tax_table@.Data,"\\|",n=7))
colnames(tax) <- microbiomeMarker:::available_ranks
ps <- phyloseq::subset_taxa(ps,  tax$Order %in% c("o__Clostridiales",""))

microbiomeMarker:::get_treedata_phyloseq(ps) %>%
  microbiomeMarker:::generate_taxa_tree(size = branch_size)