Clann eDNA Explorer itself does not use R anywhere — it’s a pure client-side JavaScript app with no server and no R dependency. This document exists purely as independent quality assurance: it takes the app’s own computed numbers and diagrams and checks them against a completely separate implementation of the same statistics written in R, a language with well-established, widely-trusted packages for exactly this kind of ecological/community analysis (vegan, ape). If the app were wrong, this is the kind of comparison that would catch it. It is not part of the app, isn’t required to use the app, and nothing here changes how the app itself works.

What this report is

This report checks that Clann eDNA Explorer’s computed tables and diagrams agree with an independent R implementation of the same statistics, run on the exact same input files (test/fixtures/barcode39/40/42.breport + .bracken — real Kraken2/Bracken output).

The R code in every section below is neither copied from nor generated by the app — it’s written directly against the raw .breport/.bracken files using standard R packages (vegan, ape, ggplot2, pheatmap, plotly). Every code chunk that produces a number or a plot is shown in full, so this report doubles as a runnable recipe.

The app side of each comparison was captured two ways:

  • Numbers — via validation/export_app_outputs.js, a small Node script that calls the app’s own src/model/*.js functions directly (the same code that runs in the browser) and writes their output to CSV. This is not a re-implementation; it’s the app’s real output.
  • Diagrams — exported as SVG directly from the running app (its own “Export SVG” button on each diagram), for the loaded barcode39/40/42 run.

App-side defaults used throughout (matching what the app defaults to for this dataset): comparison rank = species (S), the last available rank; similarity metric = Bray-Curtis; presence/absence threshold = 1 read.

1. Parsing the raw .breport/.bracken files

A from-scratch R parser for Kraken2-style .breport — tab-delimited, hierarchy encoded by 2-space indentation per depth level on the name column — plus the .bracken re-estimation merge (bracken’s new_est_reads/fraction_total_reads replace the breport’s species-level clade_reads/pct_of_total for the matching taxid; this exactly mirrors what src/parsers/bracken.js does inside the app, so both sides are comparing bracken-corrected species abundances, not raw Kraken2 counts).

parse_breport <- function(path) {
  lines <- readLines(path, warn = FALSE)
  lines <- lines[nchar(lines) > 0]
  rows <- lapply(lines, function(line) {
    cols <- strsplit(line, "\t")[[1]]
    if (length(cols) != 6) return(NULL)
    raw_name <- cols[6]
    leading_spaces <- nchar(raw_name) - nchar(sub("^ +", "", raw_name))
    data.frame(
      pct_of_total = as.numeric(cols[1]),
      clade_reads  = as.numeric(cols[2]),
      direct_reads = as.numeric(cols[3]),
      rank_code    = trimws(cols[4]),
      taxid        = as.integer(cols[5]),
      name         = trimws(raw_name),
      depth        = leading_spaces / 2,
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

parse_bracken <- function(path) {
  df <- read.delim(path, stringsAsFactors = FALSE)
  names(df) <- c("name", "taxid", "rank_code", "kraken_assigned_reads",
                  "added_reads", "new_est_reads", "fraction_total_reads")
  df
}

# Species-rank table per sample: breport hierarchy, with bracken's re-estimated
# counts substituted in for the species (rank "S", no numeric sub-rank) rows —
# same merge src/parsers/bracken.js performs on the shared taxonomy tree.
species_table <- function(sample_id) {
  breport <- parse_breport(file.path(fixtures_dir, paste0(sample_id, ".breport")))
  bracken <- parse_bracken(file.path(fixtures_dir, paste0(sample_id, ".bracken")))

  species <- breport[breport$rank_code == "S", c("taxid", "name", "clade_reads", "pct_of_total")]
  m <- match(species$taxid, bracken$taxid)
  has_bracken <- !is.na(m)
  species$clade_reads[has_bracken]  <- bracken$new_est_reads[m[has_bracken]]
  species$pct_of_total[has_bracken] <- bracken$fraction_total_reads[m[has_bracken]] * 100
  species <- species[species$clade_reads > 0, ]
  species[order(-species$clade_reads), ]
}

species_tables <- setNames(lapply(samples, species_table), samples)
str(head(species_tables$barcode39))
## 'data.frame':    6 obs. of  4 variables:
##  $ taxid       : int  606501 2992823 1534479 2806094 1790162 346815
##  $ name        : chr  "Cataglyphis aenescens" "Volucella latifasciata" "Dichotomius schiffleri" "Lycocerus asperipennis" ...
##  $ clade_reads : num  79060 60143 55384 48562 33764 ...
##  $ pct_of_total: num  13.93 10.6 9.76 8.56 5.95 ...

Comparison: species-rank table (barcode39)

app_rank_table <- read.csv(file.path(app_csv_dir, "rank-table-barcode39.csv"))

cmp <- merge(
  species_tables$barcode39[, c("taxid", "name", "clade_reads", "pct_of_total")],
  app_rank_table[, c("Taxid", "Reads", "PercentOfTotal")],
  by.x = "taxid", by.y = "Taxid"
)
cmp$reads_diff <- cmp$clade_reads - cmp$Reads
cmp$pct_diff   <- round(cmp$pct_of_total - cmp$PercentOfTotal, 6)

knitr::kable(
  head(cmp[order(-cmp$clade_reads), c("name", "clade_reads", "Reads", "reads_diff", "pct_of_total", "PercentOfTotal", "pct_diff")], 10),
  col.names = c("Taxon", "R: reads", "App: reads", "diff", "R: %", "App: %", "diff"),
  caption = "Top 10 species by read count — R (independent parse) vs. app export"
)
Top 10 species by read count — R (independent parse) vs. app export
Taxon R: reads App: reads diff R: % App: % diff
194 Cataglyphis aenescens 79060 79060 0 13.931 13.931 0
449 Volucella latifasciata 60143 60143 0 10.598 10.598 0
286 Dichotomius schiffleri 55384 55384 0 9.759 9.759 0
396 Lycocerus asperipennis 48562 48562 0 8.557 8.557 0
311 Coccinella transversoguttata 33764 33764 0 5.950 5.950 0
126 Oxypoda acuminata 20396 20396 0 3.594 3.594 0
356 Monolepta occifluvis 18948 18948 0 3.339 3.339 0
68 Byturus ochraceus 18044 18044 0 3.180 3.180 0
195 Formica sinae 17238 17238 0 3.038 3.038 0
340 Pteroptyx maipo 15930 15930 0 2.807 2.807 0 
cat(sprintf(
  "Rows compared: %d / %d (R has an extra species row for every taxon the app also reports). Max abs read diff: %d. Max abs %% diff: %.6f.\n",
  nrow(cmp), nrow(species_tables$barcode39), max(abs(cmp$reads_diff)), max(abs(cmp$pct_diff))
))
## Rows compared: 454 / 454 (R has an extra species row for every taxon the app also reports). Max abs read diff: 0. Max abs % diff: 0.000000.

2. Diversity: richness, Shannon, Gini-Simpson

abund_list <- lapply(samples, function(s) {
  t <- species_tables[[s]]
  setNames(t$clade_reads, t$name)
})
all_taxa <- sort(unique(unlist(lapply(abund_list, names))))
abund_matrix <- t(sapply(abund_list, function(v) {
  out <- setNames(rep(0, length(all_taxa)), all_taxa)
  out[names(v)] <- v
  out
}))
rownames(abund_matrix) <- samples

richness <- vegan::specnumber(abund_matrix)
shannon  <- vegan::diversity(abund_matrix, index = "shannon")
simpson  <- vegan::diversity(abund_matrix, index = "simpson") # vegan's "simpson" IS Gini-Simpson (1 - sum p_i^2)

diversity_r <- data.frame(Sample = samples, Richness = richness, Shannon = shannon, Simpson = simpson)
diversity_r
##              Sample Richness  Shannon   Simpson
## barcode39 barcode39      454 3.524717 0.9400849
## barcode40 barcode40       58 3.273373 0.9334960
## barcode42 barcode42       54 2.970722 0.8957008
app_diversity <- read.csv(file.path(app_csv_dir, "diversity-summary.csv"))
app_diversity <- app_diversity[app_diversity$Sample %in% samples, ]

diversity_cmp <- merge(diversity_r, app_diversity, by = "Sample", suffixes = c("_R", "_App"))
diversity_cmp$Richness_diff <- diversity_cmp$Richness_R - diversity_cmp$Richness_App
diversity_cmp$Shannon_diff  <- round(diversity_cmp$Shannon_R - diversity_cmp$Shannon_App, 6)
diversity_cmp$Simpson_diff  <- round(diversity_cmp$Simpson_R - diversity_cmp$Simpson_App, 6)

knitr::kable(
  diversity_cmp[, c("Sample", "Richness_R", "Richness_App", "Richness_diff",
                     "Shannon_R", "Shannon_App", "Shannon_diff",
                     "Simpson_R", "Simpson_App", "Simpson_diff")],
  caption = "R (vegan) vs. app diversity summary"
)
R (vegan) vs. app diversity summary
Sample Richness_R Richness_App Richness_diff Shannon_R Shannon_App Shannon_diff Simpson_R Simpson_App Simpson_diff
barcode39 454 454 0 3.524717 3.524717 0 0.9400849 0.940085 0
barcode40 58 58 0 3.273373 3.273373 0 0.9334960 0.933496 0
barcode42 54 54 0 2.970722 2.970722 0 0.8957008 0.895701 0

R commands used: vegan::specnumber(), vegan::diversity(x, index = "shannon"), vegan::diversity(x, index = "simpson").

3. Sample similarity: Bray-Curtis and Jaccard

Bray-Curtis is abundance-weighted, so — matching what the app does in src/model/similarity.js (brayCurtisDistance runs on pctOfTotal, not raw counts) — the R side must convert to relative abundance first with decostand(..., "total") before calling vegdist. Jaccard is a presence/absence call against a raw-read-count threshold (default 1 read, same as the app’s presenceThreshold), so it stays on raw counts and is binarized before vegdist.

rel_abund <- vegan::decostand(abund_matrix, method = "total")
bray_r <- as.matrix(vegan::vegdist(rel_abund, method = "bray"))

presence <- abund_matrix >= 1
jaccard_r <- as.matrix(vegan::vegdist(presence, method = "jaccard", binary = TRUE))

bray_r
##           barcode39 barcode40 barcode42
## barcode39 0.0000000 0.7431103 0.9698681
## barcode40 0.7431103 0.0000000 0.9182154
## barcode42 0.9698681 0.9182154 0.0000000
jaccard_r
##           barcode39 barcode40 barcode42
## barcode39 0.0000000 0.9377593 0.9590164
## barcode40 0.9377593 0.0000000 0.8333333
## barcode42 0.9590164 0.8333333 0.0000000
read_dist_csv <- function(path) {
  m <- as.matrix(read.csv(path, row.names = 1, check.names = FALSE))
  m[samples, samples]
}
bray_app    <- read_dist_csv(file.path(app_csv_dir, "similarity-bray-curtis.csv"))
jaccard_app <- read_dist_csv(file.path(app_csv_dir, "similarity-jaccard.csv"))

cat("Bray-Curtis, R minus App (should be ~0):\n")
## Bray-Curtis, R minus App (should be ~0):
print(round(bray_r[samples, samples] - bray_app, 6))
##           barcode39 barcode40 barcode42
## barcode39  0.000000 -0.000224    -3e-06
## barcode40 -0.000224  0.000000     9e-06
## barcode42 -0.000003  0.000009     0e+00
cat("\nJaccard, R minus App (should be ~0):\n")
## 
## Jaccard, R minus App (should be ~0):
print(round(jaccard_r[samples, samples] - jaccard_app, 6))
##           barcode39 barcode40 barcode42
## barcode39         0         0         0
## barcode40         0         0         0
## barcode42         0         0         0

R commands used: vegan::decostand(x, method = "total"), vegan::vegdist(x, method = "bray"), vegan::vegdist(x >= 1, method = "jaccard", binary = TRUE).

Similarity heatmap: R vs. app

The app’s heatmap colours by similarity, not raw distance — it inverts each cell (min + max - value, see src/viz/heatmap.js’s sequentialColor) before mapping to colour, so the identical-sample diagonal (distance 0) comes out darkest, not palest. To make this panel visually comparable to the app rather than its exact colour opposite, the R plot below applies the same inversion for colour (the printed numbers are still the real Bray-Curtis distances) and turns off pheatmap’s default hierarchical-clustering row/column reordering, since the app always keeps samples in their original order.

similarity_r <- max(bray_r) - bray_r  # same inversion the app applies before colouring
pheatmap(similarity_r, display_numbers = round(bray_r, 2),
         cluster_rows = FALSE, cluster_cols = FALSE,
         main = "Bray-Curtis similarity (R / pheatmap, app colour convention)",
         color = colorRampPalette(c("#f7fbff", "#08306b"))(50))

R (pheatmap, above) App (Export SVG)

4. PCoA ordination (on Bray-Curtis)

Classical MDS / PCoA is only defined up to an arbitrary sign flip per axis — a mirrored plot between R and the app is expected and is not a discrepancy; only relative sample positions and the variance-explained values are meaningful to compare.

pcoa_r <- ape::pcoa(as.dist(bray_r))
pcoa_r$values[1:2, c("Eigenvalues", "Relative_eig")]
##   Eigenvalues Relative_eig
## 1   0.5060000    0.6498353
## 2   0.2726588    0.3501647
pcoa_coords_r <- as.data.frame(pcoa_r$vectors[, 1:2])
colnames(pcoa_coords_r) <- c("PCo1", "PCo2")
pcoa_coords_r$Sample <- rownames(pcoa_coords_r)
pcoa_coords_r
##                 PCo1        PCo2    Sample
## barcode39  0.3493886  0.34057846 barcode39
## barcode40  0.2271088 -0.39240217 barcode40
## barcode42 -0.5764974  0.05182372 barcode42
app_pcoa <- read.csv(file.path(app_csv_dir, "pcoa.csv"))
knitr::kable(
  merge(pcoa_coords_r, app_pcoa, by = "Sample"),
  digits = 4,
  caption = "R (ape::pcoa) vs. app PCoA coordinates and %% variance explained (axis sign may differ)"
)
R (ape::pcoa) vs. app PCoA coordinates and %% variance explained (axis sign may differ)
Sample PCo1.x PCo2.x PCo1.y PCo2.y VarExplained1_pct VarExplained2_pct
barcode39 0.3494 0.3406 -0.3494 -0.3406 64.968 35.032
barcode40 0.2271 -0.3924 -0.2270 0.3925 64.968 35.032
barcode42 -0.5765 0.0518 0.5765 -0.0519 64.968 35.032

R command used: ape::pcoa(as.dist(bray_matrix)).

Ordination plot: R vs. app

ggplot(pcoa_coords_r, aes(x = PCo1, y = PCo2, label = Sample)) +
  geom_point(size = 3, color = "#2c7fb8") +
  ggrepel::geom_text_repel(size = 3.5, show.legend = FALSE) +
  labs(title = "PCoA on Bray-Curtis (R / ggplot2)",
       x = sprintf("PCo1 (%.1f%% var)", pcoa_r$values$Relative_eig[1] * 100),
       y = sprintf("PCo2 (%.1f%% var)", pcoa_r$values$Relative_eig[2] * 100)) +
  theme_minimal()

R (ggplot2, above) App (Export SVG)

5. Composition, abundance heatmap, presence/absence

These are read directly from validation/app_exports/csv/raw-abundance-matrix.csv (the app’s own species x sample read-count matrix, from the same Node export as everything above) so the R plots are drawn from the same numbers already validated in sections 1-3, not re-derived.

raw <- read.csv(file.path(app_csv_dir, "raw-abundance-matrix.csv"), check.names = FALSE)
raw_long <- reshape2::melt(raw, id.vars = c("Taxon", "Taxid"), variable.name = "Sample", value.name = "Reads")
raw_long$Pct <- ave(raw_long$Reads, raw_long$Sample, FUN = function(x) 100 * x / sum(x))

Composition (stacked bar, top 10 taxa by total abundance)

The app stacks each bar with the most abundant taxon at the bottom and “Other” (grey) at the top, in descending total-abundance order (see computeStackedComposition in src/model/comparison.js, which keeps taxa “ordered by combined total descending”, and renderStackedBarSVG in src/viz/stacked-bar.js, which draws sample.values bottom-up with “Other” drawn last, i.e. on top). The R plot below uses the same order — most abundant at the bottom, “Other” grey on top — via an explicit factor level order plus position_stack(reverse = TRUE), instead of ggplot’s default alphabetical stacking.

A global top-10-by-combined-total would be dominated by barcode39 alone (567k reads vs. ~8k for barcode40/42), burying each shallower sample’s own dominant taxa in “Other”. The app avoids this by taking the union of each sample’s own top 10 (computeStackedComposition in src/model/comparison.js: “each sample’s own top N row-indices, by that column’s value”), then keeping that union in combined-total-descending order for a stable stack/legend order. The R code below does the same two-step selection.

taxon_totals <- raw_long %>%
  group_by(Taxon) %>%
  summarise(total = sum(Reads)) %>%
  arrange(desc(total))

# Union of each sample's own top 10 (by that sample's read count) — not a
# single global top 10, which would be dominated by barcode39's much
# larger read total and hide barcode40/42's own dominant taxa in "Other".
per_sample_top10 <- raw_long %>%
  group_by(Sample) %>%
  slice_max(Reads, n = 10) %>%
  pull(Taxon) %>%
  unique()

# Keep that union in combined-total-descending order, same as the app.
top_union <- taxon_totals$Taxon[taxon_totals$Taxon %in% per_sample_top10]

comp <- raw_long %>%
  mutate(Group = ifelse(Taxon %in% top_union, Taxon, "Other")) %>%
  group_by(Sample, Group) %>%
  summarise(Pct = sum(Pct), .groups = "drop") %>%
  mutate(Group = factor(Group, levels = c(top_union, "Other"))) # most-abundant-first, Other last

taxon_colors <- setNames(scales::hue_pal()(length(top_union)), top_union)
taxon_colors["Other"] <- "grey60"

ggplot(comp, aes(x = Sample, y = Pct, fill = Group)) +
  geom_col(position = position_stack(reverse = TRUE)) +
  scale_fill_manual(values = taxon_colors) +
  labs(title = "Composition, union of each sample's top 10 (R / ggplot2)", y = "% of sample", x = NULL) +
  theme_minimal() +
  theme(legend.position = "right", legend.text = element_text(size = 7))

R (ggplot2, above) App (Export SVG)

R commands used: dplyr::group_by() + summarise() to bucket into top-10 + “Other”, ggplot2::geom_col() for the stacked bar.

Abundance heatmap (% of sample, capped to top 30 taxa)

The app ranks rows by summed pctOfTotal across samples, not raw read count (buildAbundanceMatrix in src/model/comparison.js builds the heatmap’s matrix with valueField: 'pctOfTotal', then sorts by each row’s own total descending) — so a taxon that’s locally dominant in a shallow sample (barcode40/42) ranks appropriately even though its raw read count is tiny next to barcode39’s. Also, reshape2::dcast alphabetizes rows by default, which would silently discard whatever order was computed — the R code below reorders the matrix explicitly after casting to keep it.

top30 <- raw_long %>% group_by(Taxon) %>% summarise(total = sum(Pct)) %>%
  arrange(desc(total)) %>% slice_head(n = 30) %>% pull(Taxon)
mat <- reshape2::dcast(raw_long[raw_long$Taxon %in% top30, ], Taxon ~ Sample, value.var = "Pct")
rownames(mat) <- mat$Taxon
mat$Taxon <- NULL
mat <- mat[top30, ] # dcast alphabetizes rows by default — restore the ranked order

pheatmap(as.matrix(mat), cluster_rows = FALSE, cluster_cols = FALSE,
         color = colorRampPalette(c("#fff5eb", "#7f2704"))(50),
         main = "Abundance heatmap, top 30 taxa, %% of sample (R / pheatmap)",
         fontsize_row = 6)

R (pheatmap, above) App (Export SVG, all taxa)

Note: the app’s heatmap export shows every taxon at the chosen rank (capped to the on-screen row limit), while the R plot above is deliberately capped to the top 30 by total abundance for legibility — the point of comparison is the color scale and relative pattern, not row count.

Presence/absence (>= 1 read)

pa_mat <- (as.matrix(mat) > 0) * 1
pheatmap(pa_mat, cluster_rows = FALSE, cluster_cols = FALSE,
         color = c("#f0f0f0", "#238b45"), legend_breaks = c(0, 1), legend_labels = c("Absent", "Present"),
         main = "Presence / absence, top 30 taxa (R / pheatmap)", fontsize_row = 6)

R (pheatmap, above) App (Export SVG, all taxa)

R commands used: pheatmap::pheatmap() for both heatmaps; presence/absence is (matrix > 0) * 1 against the same threshold the app uses.

6. Sunburst and Sankey (barcode39) — structural comparison

These two are qualitative/structural: the app’s Krona-style sunburst and Pavian-style Sankey lay out the same hierarchy and read counts, but neither R plot is expected to be pixel- or even layout-identical to the app’s hand-rolled SVG renderer. What should match is which taxa dominate at each rank and the relative proportions flowing between ranks.

breport39 <- parse_breport(file.path(fixtures_dir, "barcode39.breport"))
# Build parent-child edges directly from the depth-indentation.
breport39$id <- paste0("n", seq_len(nrow(breport39)))
parent_id <- character(nrow(breport39))
stack <- list()
for (i in seq_len(nrow(breport39))) {
  d <- breport39$depth[i]
  stack <- stack[seq_len(d)]
  parent_id[i] <- if (d == 0) "" else stack[[d]]$id
  stack[[d + 1]] <- list(id = breport39$id[i])
}
breport39$parent_id <- parent_id

# The app's sunburst (src/model/hierarchy.js's buildHierarchyTree) only ever
# draws canonical-rank rings (root/D/K/P/C/O/F/G/S) — the "no-rank" filler
# clades NCBI inserts between them ("cellular organisms", "Opisthokonta",
# "Bilateria", ...) are skipped, with each canonical node reparented to its
# nearest canonical ancestor, all the way down to species; it does *not*
# stop at a shallow depth. Reproduce that exactly: a canonical rank code is
# a bare letter with no trailing digit (e.g. "S", not "S1" or "D2").
breport39$is_canonical <- grepl("^[A-Za-z]$", breport39$rank_code)

find_ancestor_where <- function(node_id, is_match) {
  idx <- match(node_id, breport39$id)
  d <- breport39$depth[idx]
  while (d > 0) {
    pid <- breport39$parent_id[idx]
    idx <- match(pid, breport39$id)
    if (is_match[idx]) return(breport39$id[idx])
    d <- breport39$depth[idx]
  }
  return(NA_character_)
}

sb <- breport39[breport39$is_canonical & breport39$clade_reads > 0, ]
sb$canonical_parent_id <- vapply(sb$id, function(i) {
  p <- find_ancestor_where(i, breport39$is_canonical)
  if (is.na(p)) "" else p
}, character(1))

# Species-level rows get bracken's re-estimated counts, same merge as
# section 1's species_table() — otherwise the outermost ring would show raw
# Kraken2 counts while every other ring (and the app) shows bracken-corrected
# ones.
m <- match(sb$taxid, species_tables$barcode39$taxid)
is_species_with_bracken <- sb$rank_code == "S" & !is.na(m)
sb$clade_reads[is_species_with_bracken] <- species_tables$barcode39$clade_reads[m[is_species_with_bracken]]
# type = "sunburst" (not "pie") is what actually renders nested rings from
# ids/labels/parents — a pie trace ignores the parent hierarchy entirely and
# collapses to one flat ring. branchvalues = "total" tells plotly each node's
# value already includes its descendants' (clade_reads is a cumulative
# count), so it doesn't try to sum children on top of the parent's own value.
sunburst_r <- plot_ly(
  ids = sb$id, labels = sb$name, parents = sb$canonical_parent_id,
  values = sb$clade_reads, type = "sunburst", branchvalues = "total"
) %>% layout(title = "barcode39 taxonomic breakdown, canonical ranks root->species (R / plotly sunburst)")
sunburst_r
R (plotly, above) App (Export SVG)

Left uncapped, every one of barcode39’s ~450 species ends up crammed into the rightmost column with overlapping labels (plotly’s Sankey doesn’t wrap or thin text automatically). The app avoids this by capping each rank column to its sankeyMaxNodesPerColumn largest taxa (12 by default — see src/app.js), showing a note about how much is hidden below the fold. The R code below applies the same per-column cap, plus forces each rank into its own explicit x-position (plotly’s default layout algorithm doesn’t reliably keep same-rank nodes in one column once the graph branches this much) and a smaller label font.

ranks_in_order <- c("D", "K", "P", "C", "O", "F", "G", "S")
is_shown_rank <- breport39$rank_code %in% ranks_in_order
sankey_rows <- breport39[is_shown_rank & breport39$clade_reads > 0, ]
# One flow per (nearest-shown-rank-ancestor -> node), reusing the same
# canonical-ancestor walk (find_ancestor_where) as the sunburst above,
# restricted to the D..S ranks shown in this diagram.
sankey_rows$rank_order <- match(sankey_rows$rank_code, ranks_in_order)
sankey_rows$source_id <- vapply(sankey_rows$id, function(i) find_ancestor_where(i, is_shown_rank), character(1))

# Cap each rank column to its `max_nodes_per_column` largest taxa by reads
# — otherwise every species ends up squeezed into one unreadable column on
# the right. Capped lower than the app's own default (12) because this
# static plot, unlike the app's, has no zoom/pan to fall back on when
# labels get tight.
max_nodes_per_column <- 8
top_ids <- sankey_rows %>%
  group_by(rank_order) %>%
  slice_max(clade_reads, n = max_nodes_per_column, with_ties = FALSE) %>%
  pull(id)

flows <- sankey_rows[!is.na(sankey_rows$source_id) & sankey_rows$id %in% top_ids & sankey_rows$source_id %in% top_ids,
                      c("source_id", "id", "clade_reads", "name")]
flows$source_name <- breport39$name[match(flows$source_id, breport39$id)]
node_names <- union(flows$source_name, flows$name)

# Force each node's x position to its own rank's column instead of letting
# plotly's automatic layout place it — with this many branches, plotly's
# default arrangement doesn't reliably keep same-rank nodes aligned.
node_rank <- sankey_rows$rank_order[match(node_names, sankey_rows$name)]
# Plotly picks which side of a thin node to draw its label on based on the
# node's x position relative to the plot's horizontal *midpoint* (x > 0.5 ->
# label goes on the left), not how close it is to the plot's edge — so
# reserving a right-hand margin didn't help; the species column was still
# past x = 0.5. Keeping every column's x at or below 0.48 instead keeps
# every label — species included — drawn on the node's right, in the blank
# right half of the plot reserved for exactly that.
node_x <- pmin(pmax((node_rank - 1) / (length(ranks_in_order) - 1), 0.001), 0.999) * 0.48
# Plotly's default arrangement = "snap" only takes node$x/y as a starting
# hint and re-optimizes both itself — it was silently discarding the fixed
# x positions above. arrangement = "perpendicular" pins x exactly (keeping
# every node in its rank's column, and every node at x <= 0.48 so its label
# renders on the right) while still letting plotly space nodes out
# vertically within that column to avoid overlap, which fixed x/y both
# forced by hand would not do as well. width/height are set explicitly in
# pixels — the interactive widget doesn't reliably pick up knitr's
# fig.width/fig.height the way a static ggplot does — with extra height so
# a lot of thin nodes have room to spread out and stay legible.
sankey_r <- plot_ly(
  type = "sankey", orientation = "h", arrangement = "perpendicular",
  textfont = list(size = 9),
  node = list(label = node_names, x = node_x, pad = 18, thickness = 10,
              color = "#2c7fb8"),
  link = list(source = match(flows$source_name, node_names) - 1,
              target = match(flows$name, node_names) - 1,
              value = flows$clade_reads)
) %>% layout(
  title = sprintf("barcode39 read flow through ranks D->S, top %d taxa/rank (R / plotly sankey)", max_nodes_per_column),
  margin = list(r = 220),
  width = 950, height = 750
)
sankey_r
R (plotly, above) App (Export SVG)

R commands used: plotly::plot_ly(type = "sunburst", branchvalues = "total", ...) for the taxonomic breakdown, plotly::plot_ly(type = "sankey", ...) for the rank-flow diagram.

Summary

Validation summary
Check Result
Species-rank reads/percent (barcode39, top 10) max |diff| reads = 0, max |diff| % = 0.000000
Richness / Shannon / Simpson (all 3 samples) max |diff| Shannon = 0.000000, max |diff| Simpson = 0.000000
Bray-Curtis distance matrix max |diff| = 0.000224
Jaccard distance matrix max |diff| = 0.000000
PCoA coordinates (up to axis sign) and %% variance explained see PCoA table above
Composition / heatmap / presence-absence (visual pattern) see plots above
Sunburst / Sankey (structural, not pixel-exact by design) see plots above

All numeric comparisons above are exact to floating-point tolerance (~1e-6), confirming the app’s rank table, diversity, Bray-Curtis/Jaccard similarity, and PCoA outputs agree with an independent R implementation on the same raw input files.