mikemc/speedyseq: speedyseq v0.2.0

Breaking changes

• The default ordering of new taxa output by tax_glom() is different from previous versions and from phyloseq::tax_glom() in phyloseq objects that do not have phylogenetic trees. See "Minor improvements and fixes" for more information.

New features New general-purpose vectorized taxa-merging function

The new merge_taxa_vec() function provides a vectorized version of phyloseq::merge_taxa() that can quickly merge arbitrary groups of taxa and now forms the basis of all other merging functions. phyloseq::merge_taxa() takes a phyloseq object or component object x and a set of taxa eqtaxa and merges them into a single taxon. In place of the eqtaxa argument, speedyseq's merge_taxa_vec() takes a vector group of length ntaxa(physeq) that defines how all the taxa in x should be merged into multiple new groups. Its syntax and behavior is patterned after that of base::rowsum(), which it uses to do the merging in the OTU table. When aiming to merge a large number of taxa into a smaller but still large number of groups, it is much faster to do all the merging with one call to merge_taxa_vec() than to loop through many calls to merge_taxa().

A practical example is clustering amplicon sequence variants (ASVs) into OTUs defined by a given similarity threshold. Suppose we have a phyloseq object ps that has the ASV sequences stored in its refseq slot. We can cluster the ASV sequences into 97% OTUs using the DECIPHER package with

dna <- refseq(ps)
nproc <- 1 # Increase to use multiple processors
aln <- DECIPHER::AlignSeqs(dna, processors = nproc)
d <- DECIPHER::DistanceMatrix(aln, processors = nproc)
clusters <- DECIPHER::IdClusters(
  d, 
  method = "complete",
  cutoff = 0.03, # corresponds to 97% OTUs
  processors = nproc
)

Next we merge_taxa_vec() to get the merged phyloseq object,

ps0 <- merge_taxa_vec(
  ps, 
  group = clusters$cluster,
  tax_adjust = 2
)

The names of the new taxa are set to the name of the most abundant taxon within each group (the same behavior as the tax_glom() and tip_glom() functions). Future versions will likely have a names argument to control the naming behavior.

The tax_adjust argument controls how NAs and within-group disagreements in taxonomy are handled to determine the taxonomy of the new taxa (see help(merge_taxa_vec) for details). An example of the difference between tax_adjust = 1 (phyloseq::merge_taxa() behavior) and tax_adjust = 2 can be seen in the following example from the new tip_glom() documentation,

data(GlobalPatterns)

set.seed(20190421)
ps <- prune_taxa(sample(taxa_names(GlobalPatterns), 2e2), GlobalPatterns)

ps1 <- tip_glom(ps, 0.1, tax_adjust = 1)
ps2 <- tip_glom(ps, 0.1, tax_adjust = 2)
tax_table(ps1)[c(108, 136, 45),]
#> Taxonomy Table:     [3 taxa by 7 taxonomic ranks]:
#>        Kingdom    Phylum           Class                 Order               
#> 578831 "Bacteria" "Bacteroidetes"  "Sphingobacteria"     "Sphingobacteriales"
#> 2801   "Bacteria" "Planctomycetes" "Planctomycea"        "Pirellulales"      
#> 185581 "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Oceanospirillales" 
#>        Family Genus            Species                             
#> 578831 NA     "Niabella"       NA                                  
#> 2801   NA     "Rhodopirellula" NA                                  
#> 185581 "OM60" NA               "marinegammaproteobacteriumHTCC2080"
tax_table(ps2)[c(108, 136, 45),]
#> Taxonomy Table:     [3 taxa by 7 taxonomic ranks]:
#>        Kingdom    Phylum           Class                 Order               
#> 578831 "Bacteria" "Bacteroidetes"  "Sphingobacteria"     "Sphingobacteriales"
#> 2801   "Bacteria" "Planctomycetes" "Planctomycea"        "Pirellulales"      
#> 185581 "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Oceanospirillales" 
#>        Family Genus Species
#> 578831 NA     NA    NA     
#> 2801   NA     NA    NA     
#> 185581 "OM60" NA    NA

Faster and lower-memory implementation of phyloseq::tip_glom()

The new tip_glom() function provides a speedy version of phyloseq::tip_glom(). This function performs a form of indirect phylogenetic merging of taxa using the phylogenetic tree in phy_tree(physeq) by 1) using the tree to create a distance matrix, 2) performing hierarchical clustering on the distance matrix, and 3) defining new taxonomic groups by cutting the dendrogram at the height specified by the h parameter. Speedyseq's tip_glom() provides a faster and less memory-intensive alternative to phyloseq::tip_glom() through the use of vectorized merging (via merge_taxa_vec()) and faster and lower-memory phylogenetic-distance computation (via get_all_pairwise_distances() from the castor package).

Speedyseq's tip_glom() also has the new tax_adjust argument, which is passed on to merge_taxa_vec(). It is set to 1 by default for phyloseq compatibility and should give identical results to phyloseq in this case.

For phyloseq compatibility, the default clustering function is left as cluster::agnes(). However, equivalent but faster results can be obtained by using the hclust function from base R with the method == "average" option.

Direct phylogenetic merging with tree_glom()

The new tree_glom() function performs direct phylogenetic merging of taxa. This function is much faster and arguably more intuitive than tip_glom(). Advantages of direct merging over the indirect merging of tip_glom() are

1. A merged group of taxa correspond to a clade in the original tree being collapsed to a single taxon.
2. The resolution parameter that controls the degree of merging has units in terms of the tree's branch lengths, making it potentially more biologically meaningful than the h parameter in tip_glom().
3. The distance-matrix computation and hierarchical clustering in tip_glom() can be skipped, making tree_glom() much faster and less memory intensive than tip_glom() when the number of taxa is large.

tree_glom() uses functions from the castor package to determine which clades are to be merged using one of three criteria. The default behavior is to merge a clade if the maximum distance from a node to each of its tips is less than the distance resolution.

data(GlobalPatterns)
ps1 <- subset_taxa(GlobalPatterns, Phylum == "Chlamydiae")
ntaxa(ps1)
ps2 <- tree_glom(ps1, 0.05)
ntaxa(ps2)

library(dplyr)
library(ggtree)
library(cowplot)

plot1 <- phy_tree(ps1) %>%
  ggtree +
  geom_tiplab() +
  geom_label(aes(x = branch, label = round(branch.length, 4)))
plot2 <- phy_tree(ps2) %>%
  ggtree +
  geom_tiplab() +
  geom_label(aes(x = branch, label = round(branch.length, 4)))
plot_grid(plot1, plot2)

Minor improvements and fixes

• Fixed bug that applied to taxonomic merge functions when an object named new_tax_mat exists outside the function environment; described in Issue #31

• Merging functions now maintain the original order of new taxa by default, except in phyloseq objects with phylogenetic trees (for which order is and has always been determined by how archetypes are ordered in phy_tree(ps)$tip.label). This behavior can lead to different taxa orders from past speedyseq versions and from phyloseq::tax_glom() function. However, it makes the resulting taxa order more predictable. New taxa can be be reordered according to group or taxonomy in merge_taxa_vec() and tax_glom() by setting reorder = TRUE.

• Merging/glom functions now work on relevant phyloseq components as well as phyloseq objects