TAXDB <- 'sbdi-gtdb_R22-RS232-1'
GTDB_TAXRANKS <- c('domain', 'phylum', 'class', 'order', 'family', 'genus', 'species')
ORDERED_TAXRANKS <- c('domain', 'kingdom', 'phylum', 'class', 'order', 'family', 'genus', 'species', 'strain', 'subspecies')
BLAST6_COLNAMES <- c('qseqid', 'sseqid', 'pident', 'length', 'mismatch', 'gapopen', 'qstart', 'qend', 'sstart', 'send', 'evalue', 'bitscore')
BLAST6_COLTYPES <- 'ccdiiiiiiidd'rank <- function(ds) {
ds %>%
mutate(
rank = case_when(
!is.na(species) ~ 'species',
!is.na(genus) ~ 'genus',
!is.na(family) ~ 'family',
!is.na(order) ~ 'order',
!is.na(class) ~ 'class',
!is.na(phylum) ~ 'phylum',
!is.na(domain) ~ 'domain'
)
)
}taxonomy <- read_tsv(
'data/lmo16s.sbdi-gtdb.taxonomy.tsv',
col_types = cols(.default = col_character(), confidence = col_double())
) %>%
anti_join(
read_tsv(
'data/lmo16s.silva.taxonomy.tsv',
col_types = cols(.default = col_character(), confidence = col_double())
) %>%
filter(Class == 'Chloroplast'),
by = join_by(ASV_ID)
) %>%
rename_with(str_to_lower) %>%
rank()
samples <- read_sheet('https://docs.google.com/spreadsheets/d/1Pxw-pXlSCCEORKOEwm9GrDxHVLca_ipwgpDdpphgARA/edit?gid=288998864#gid=288998864') %>%
filter(project == 'LMO', ! replaced) %>%
mutate(
filter_fraction = as.character(filter_fraction),
replicate = as.character(replicate)
) %>%
unnest(date)
counts <- read_tsv('data/lmo16s.tsv.gz', col_types = 'cci') %>%
semi_join(taxonomy, by = join_by(asv_id)) %>%
group_by(sample) %>%
mutate(relab = count/sum(count)) %>%
ungroup()
pplace <- read_tsv(Sys.glob('data/lmo16s.*.sbdi-gtdb.pplace.tsv'), col_types = 'cddddc') %>%
rename_with(str_to_lower) %>%
rename(asv_id = name) %>%
mutate(taxopath = str_remove_all(taxopath, '[a-z]__')) %>%
separate(taxopath, into = GTDB_TAXRANKS, sep = ';', fill = 'right') %>%
semi_join(taxonomy, by = join_by(asv_id)) %>%
rank()I compared the DADA2 method assignTaxonomy with phylogenetic placement performed by nf-core/phyloplace using the same GTDB-SBDI database and the full set of LMO 16S ASV sequences. The latter turned out to assign substantially more ASVs at the more specific ranks species and genus. At the level of individual ASVs, I could observe movements between ranks in both directions in terms of specificity, i.e. from more to less specific in the DADA2 assignment and phylogenetic placement respectively. The addSpecies method of DADA2, based on exact matches to reference sequences, assigned a very small proportion of ASVs to species. Overall, phylogenetic placement appears to outperform the sequence composition-based DADA2 algorithm.
I assigned taxonomy to the ASVs from the LMO time series using the SBDI-GTDB database (release R11-RS232; Lundin and Andersson (2021)) using the standard DADA2 method assignTaxonomy (Callahan et al. (2016)) and phylogenetic placement with nf-core/phyloplace (Lundin (2025)), both called as part of nf-core/ampliseq (Straub et al. (2020)). After removal of ASVs assigned to the “Chloroplast” class by SILVA (Quast et al. (2013)), 199619 remained.
I also counted the number of ASVs assigned to species by DADA2’s addSpecies method that is based on exact matches to reference sequences with unambiguous species assignment.
This analysis was run with R version 4.6.0 (2026-04-24). Versions of packages used can be found in Table 1.
cite_packages(output = "table", pkgs = "Session", out.dir = getwd()) %>%
kable()DADA2 assigned 198879 ASVs at least at the domain level, while nf-core/phyloplace placed 199234. (The alternative addSpecies method of DADA2 assigned only 3279 ASVs to species based on exact matches to unambiguous reference sequences.) Phylogenetic placement proved to be substantially more specific than the DADA2 assignTaxonomy (Figure 1 A) with 111678 ASVs placed at species level with nf-core/phyloplace compared with only 75323 assigned to species by DADA2. Movement between ranks occurred both from and to more specific, with both algorithms (Figure 1 B). Placement species were recruited from all ranks in the assignTaxonomy, while species in the latter were mostly assigned either to species or genus in the former (with a few exceptions moving to order).
(Surprisingly few ASVs were placed at the family rank.)
p0 <- taxonomy %>%
filter(!is.na(rank)) %>%
mutate(source = 'DADA2-assignTaxonomy') %>%
count(source, rank) %>%
union(
pplace %>%
semi_join(taxonomy, by = join_by(asv_id)) %>%
mutate(source = 'nf-core/phyloplace') %>%
count(source,rank)
) %>%
mutate(rank = factor(rank, levels = GTDB_TAXRANKS, ordered = TRUE)) %>%
ggplot(aes(x = rank, y = n, fill = source)) +
geom_col(position = 'dodge') +
geom_label(aes(y = n + 1e4, label = n)) +
labs(x = '', y = '', fill = '')
p1 <- taxonomy %>% transmute(asv_id, DADA2 = rank) %>%
full_join(
pplace %>% transmute(asv_id, phyloplace = rank),
by = join_by(asv_id)
) %>%
make_long(DADA2, phyloplace) %>%
ggplot(aes(x = x, next_x = next_x, node = node, next_node = next_node, fill = factor(node, levels = GTDB_TAXRANKS, ordered = TRUE))) +
geom_sankey(colour = 'black') +
labs(x = '', y = '', fill = '') +
scale_fill_brewer(palette = 'Paired')
p0 / p1 +
plot_annotation(tag_levels = 'A') +
plot_layout(heights = c(2, 3))
# vim: sw=2