---
title: "Using BOLDNODE functionality"
author: "Sameer Padhye, Spencer K. Monckton, Dirk Steinke, Jireh Agda, Teresita M. Porter"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: cosmo
    highlight: pygments
    toc: true
    toc_depth: 5
    toc_float: true
    number_sections: false
    fig_width: 8
    fig_height: 5
    self_contained: false
vignette: >
  %\VignetteIndexEntry{Using BOLDNODE}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE,
  warning = FALSE,
  message = FALSE
)
```

```{=html}
<style>
body {
  font-size: 18px;
  line-height: 1.7;
}

h1 {
  font-weight: 600;
  color: #2c3e50;
  border-bottom: 2px solid #3498db;
  padding-bottom: 0.4em;
}

h2 {
  font-weight: 600;
  color: #34495e;
  margin-top: 2em;
}

h3 {
  color: #2980b9;
}

pre {
  background: #f6f8fa;
  border: 1px solid #e1e4e8;
  border-radius: 6px;
  padding: 1em;
}

code {
  background: #f1f3f5;
  color: #c7254e;
  padding: 0.15em 0.4em;
  border-radius: 3px;
  font-size: 0.88em;
}

pre code {
  background: transparent;
  color: inherit;
  padding: 0;
}

.note {
  background: #e8f4fd;
  border-left: 4px solid #3498db;
  padding: 1em 1.2em;
  border-radius: 0 6px 6px 0;
  margin: 1.5em 0;
}

.warning {
  background: #fef9e7;
  border-left: 4px solid #f39c12;
  padding: 1em 1.2em;
  border-radius: 0 6px 6px 0;
  margin: 1.5em 0;
}

table {
  border-collapse: separate;
  border-spacing: 0;
  border-radius: 6px;
  overflow: hidden;
  box-shadow: 0 1px 4px rgba(0,0,0,0.06);
}

thead th {
  background: #3498db;
  color: white;
  border: none;
  padding: 0.8em 1em;
  font-size: 16px;
  line-height: 1.2;
}

tbody tr:nth-child(even) {
  background: #f8f9fa;
}

tbody td {
  padding: 0.7em 1em;
  border-top: 1px solid #ecf0f1;
}

#TOC {
  border-right: 1px solid #ecf0f1;
  padding-right: 1em;
}

.dataTables_wrapper {
  max-width: 1000px;
  margin: 0 auto;
}

table.dataTable {
  width: 100% !important;
}

div.dataTables_scrollBody {
  overflow-x: auto !important;
}
</style>
```

# Introduction

`BOLDNODE` (**N**o-API **O**ffline **D**ata **E**xplorer) is a sister package to `BOLDconnectR` that provides local access to BOLD data packages through DuckDB, enabling exploration and analysis of public BOLD data without API limitations. `BOLDNODE` therefore, offers an offline solution for working with large scale BOLD datasets locally. Please refer to the package documentation <https://github.com/sameerpadhye/BOLDNODE> for more details about the package.

This vignette demonstrates a complete workflow integrating the *search*, *collect*, *get\_* and *bcdm_to* functions in `BOLDNODE` for efficient and scalable exploration of BOLD data packages.

As a case study, we use Canadian **Cerambycidae** (longhorn beetles) to illustrate how these functions can be combined to search, retrieve, and transform the BCDM data.

# Objectives

1.  Search all public Canadian Cerambycidae records from the BOLD data package.
2.  Collect (download) the search.
3.  Summarize the retrieved data using a concise overview.
4.  Determine the lowest supported taxonomic identification for each Cerambycidae BIN through BIN-level taxonomic concordance.
5.  Construct a `BIN × province` occurrence matrix, calculate the species richness estimates,pairwise beta diversity using the `vegan` package and visualize the species X province association via Non Metric Multidimensional Scaling (nMDS) using `vegan` and `ggplot2`.
6.  Obtain representative records for each BIN based on sequence length and identification method.
7.  Extract *Clytus* sequences as a `DNAStringSet` object, calculate nucleotide base proportions using the `Biostrings` package, align the sequences using `Muscle` and visualize the alignment using `ape` and `phangorn`.
8.  Convert *Monochamus* occurrence records into an `sf` object and visualize their geographic distribution using `ggplot2` package.

# Practical Workflow

## Installing and importing `BOLDNODE` and other packages

The latest (in development) version of `BOLDNODE` can be installed via `pak`. Packages like `Biostrings` require `BiocManager` since they are hosted by **Bioconductor**.

```{r install packages,message=F,warning=F,echo=FALSE,include=F}
# CRAN packages
library(BOLDNODE)
library(dplyr)
library(ggplot2)
library(tibble)
library(tidyr)
```

```{r package_install_for_vignette,message=F,warning=F,echo=FALSE}
# library(BOLDNODE)
```
## Data packages

BOLD data packages can be downloaded in the `parquet` format from the BOLD website (<https://boldsystems.org/data/data-packages/>). Historical Data of quarterly BOLD snapshots (each assigned a stable DOI to support reproducible research) can also be accessed by scrolling down the same webpage.

## Data package import

The path where the file is downloaded is saved as a variable that can be used as the input for the `search` function.

```{r setup-parquet_toshow}
# Path to the downloaded BOLD Parquet release file
# parquet_file <- "G:/usr/path/file.parquet"
```

```{r setup-parquet,include=F}
# Path to the downloaded BOLD Parquet release file

parquet_file <- system.file(
  "extdata",
  "test_data.parquet",
  package = "BOLDNODE"
)
```

### 1. Search from the data package

`bold_parquet_search` queries the parquet data package using user-defined search terms such as **COI-5P** for marker type, **Cerambycidae** for taxonomic group, and **Canada** for geographic region.

```{r search,message=F,warning=F}
# Search the BOLD dataset
cerambycidae_search <- bold_parquet_search(
  input.parquet = parquet_file,
  taxonomy = "Cerambycidae",
  geography = "Canada",
  marker = "COI-5P"
)
```

```{r search-records,include=F}
tot_records <- cerambycidae_search %>%
  dplyr::summarise(Total_records = dplyr::n()) %>%
  dplyr::collect() %>%
  dplyr::pull(Total_records)
```

The total number of records is shown in the console.

```{r print_records,echo=FALSE}
cat("The search has", tot_records, "records in the dataset")
```

### 2. Collect the dataset

The search result is then collected into memory (i.e., downloaded in the R session) as a data frame using `bold_search_collect` (*Please note that only first 100 results and a few columns are shown here*).

```{r collect,warning=F,message=F}
# Collect the search results into memory
cerambycidae_data <- bold_search_collect(
  cerambycidae_search,
  chunk.size = 50000,
  export = FALSE
)
```

```{r display-tabular_output,message=F,warning=F}
# Inspect the collected data
DT::datatable(
  head(cerambycidae_data, 100) %>% select(processid, sampleid, bin_uri, family, genus, species),
  options = list(pageLength = 10, scrollX = TRUE)
)
```

### 3. Generate a concise summary

`get_concise_summary` provides a summary of the search results, including total records, unique species, BINs, countries, institutes, markers, and marker length range.

```{r summary,message=F,warning=F}
# Generate concise summary
cerambycidae_summary <- get_concise_summary(cerambycidae_search)
DT::datatable(cerambycidae_summary)
```

### 4. Compute BIN consensus taxonomy

`get_bin_consensus` computes the consensus taxonomic identifications for each BIN in the Cerambycidae dataset. The function works upwards through taxonomic ranks (i.e., from subspecies to kingdom) to find the lowest concordant identification that satisfies the provided threshold (*Please note that only first 50 results are shown here*).

```{r step4-bin-consensus,message=F,warning=F}
# Get BIN consensus with strict consensus threshold (1.0)
bin_consensus <- get_bin_consensus(
  cerambycidae_search,
  threshold = 1.0,
  min.ids = 1
)

# View consensus at different taxonomic ranks
table(bin_consensus$concordant_rank)

# Examine the results
DT::datatable(head(bin_consensus, 50),
  options = list(
    pageLength = 10,
    scrollX = TRUE
  )
)
```

### 5. Obtain a BIN-by-province occurrence matrix for all cerambycids

`bcdm_to_occmatrix` creates a `site-by-taxon` occurrence matrix from the full Cerambycidae dataset, using **province.state** as the geographic grouping variable and **species** as the taxonomic rank (*Please note only 20 species are shown here*).

```{r step5-occurrence-matrix,message=F,warning=F}
# Generate occurrence matrix at species level, grouped by province/state
occ_matrix <- bcdm_to_occmatrix(
  cerambycidae_search,
  kingdom = "Animalia",
  taxon.rank = "species",
  site.cat = "province.state"
)

# Visualize the output
DT::datatable(
  head(occ_matrix, 20),
  options = list(
    pageLength = 10,
    scrollX = TRUE
  )
)
```

*The resulting occurrence matrix can be directly used as input data for many functions in widely used packages like* `vegan`.

#### 5a. Species richness estimate (using external package `vegan`)

```{r vegan-example-richness,message=F,warning=F}
library(vegan)
species_richness_est <- poolaccum(occ_matrix)
species_richness_est
```

*Species based richness estimates*

#### 5b. Beta diversity between the provinces (using external package `vegan`)

```{r column-to-rowname, include=F,message=F,warning=F}
occ_matrix <- occ_matrix %>%
  column_to_rownames("province.state")
```

```{r vegan-example-beta-diversity,message=F,warning=F}
# Calculate pairwise beta diversity
beta_diversity <- vegdist(occ_matrix, method = "bray") |> round(2)

# Visualize
DT::datatable(
  data.frame(as.matrix(beta_diversity)),
  options = list(
    pageLength = 10,
    scrollX = TRUE
  )
)
```

*Pairwise Bray Curtis dissimilarity between the provinces*

#### 5c. Using Non Metric Multidimensional Scaling for visualizing the pairwise distances between provinces (using external packages `vegan` & `ggplot2`)

```{r nMDS,,message=F,warning=F,fig.width=12, fig.height=12, dpi=300}
library(ggrepel)
# Using the occurrence matrix generated above for nMDS
nmds <- metaMDS(
  occ_matrix,
  distance = "bray",
  k = 2,
  trymax = 100,
  trace = FALSE
)

# Site scores
sites <- as.data.frame(scores(nmds, display = "sites"))
sites$Site <- rownames(sites)

# Species scores
species <- as.data.frame(scores(nmds, display = "species"))
species$Species <- rownames(species)

# Keep only the 40 species furthest from the origin
species$dist <- sqrt(species$NMDS1^2 + species$NMDS2^2)
species20 <- species[order(species$dist, decreasing = TRUE), ][1:40, ]

# Plot
ggplot() +
  geom_point(
    data = sites,
    aes(NMDS1, NMDS2),
    size = 4
  ) +
  geom_text(
    data = sites,
    aes(NMDS1, NMDS2, label = Site),
    size = 4,
    alpha = 0.6,
    vjust = -0.4
  ) +
  geom_text_repel(
    data = species20,
    aes(NMDS1, NMDS2, label = Species),
    colour = "red",
    fontface = "bold",
    size = 3,
    max.overlaps = Inf
  ) +
  coord_equal() +
  theme_bw() +
  labs(
    title = "nMDS (Bray-Curtis)",
    x = "NMDS1",
    y = "NMDS2"
  )
```

*A nMDS plot showing the associations between cerambycid species and provinces (Only 40 species shown here for representation)*

### 6. Select representative records for each BIN

`get_bin_reps` sub-samples the search results by BIN and selects representative records according to the provided criteria, yielding a representative sample for downstream analyses. As we saw above, some BINs include multiple species, therefore we will select representatives for all BIN-taxon combinations. Here, we select up to three records per BIN, prioritizing 658-bp sequences from specimens identified via morphological examination. With `enforce.scientific = TRUE`, interim or provisional names are ignored as distinct taxa (e.g. "Sternidius sp. A" is treated the same as "Sternidius").

```{r step5-bin-reps,message=F,warning=F}
# Sample 3 records per BIN-taxon combination
bin_reps <- get_bin_reps(
  bold.search.res = cerambycidae_search,
  Nreps = 3,
  by.tax = TRUE,
  enforce.scientific = TRUE,
  criteria = list(
    seq_length = 658,
    id_method = "Morphology",
    vouchered = TRUE
  )
)

# Compare sample to full dataset
n_bins <- length(unique(bin_reps$bin_uri))
cat(
  "Sampled", nrow(bin_reps), "representatives from", n_bins, "BINs",
  "(out of", tot_records, "total records)."
)
```

The actual number of sampled records varies according to BIN size and taxonomy, as we can see in the collected results (*Please note only 50 records are shown here*). In the case of BINs with multiple species, further analysis of barcode sequences often reveals small but consistent differences between them.

```{r step5-bin-rep-DT,message=F,warning=F}
DT::datatable(
  head(bin_reps, 50),
  options = list(
    pageLength = 10,
    scrollX = TRUE
  )
)
```

### 7. Generate a DNAStringset object for the genus *Clytus*

`bcdm_to_dnastringset` converts the search results into a `DNAStringSet` object. First, `bold_parquet_search` is used to filter the COI-5P *Clytus* records from Canada. Additional search parameters are put in to return only full length barcodes having very low (or none) ambiguous bases (**Only top three results displayed**) (*Please note: The* `bcdm_to_dnastringset` *function requires* `Biostrings` *to be installed before use*).

```{r DNAStringSet-object,message=F,warning=F}
library(Biostrings)
# Filter for genus Clytus
clytus_search <- bold_parquet_search(
  input.parquet = parquet_file,
  taxonomy = "Clytus",
  geography = "Canada",
  marker = "COI-5P",
  basecount = 658,
  ambi.base.cutoff = "<1%"
)

# Collect the filtered data
clytus_data <- bcdm_to_dnastringset(clytus_search,
  cols_for_seq_names = c("processid", "bin_uri", "species")
)
head(clytus_data, 3)
```

*The* `DNAStringSet` *object can then be used in various downstream analyses using packages like* `Biostrings`, `muscle`, `ape` *etc.*

#### 7a. Obtain base proportions from the sequences (using external package `Biostrings`)

```{r basepair-freq,message=F,warning=F}
library(Biostrings)
# Count bases across all sequences
base_counts <- colSums(alphabetFrequency(clytus_data, baseOnly = TRUE))
# Proportions
base_props <- base_counts / sum(base_counts)
base_props
```

#### 7b. Multiple sequence alignment using 'Muscle' (using external packages `muscle` and `Biostrings`)

```{r sequence_alignment,message=F,warning=F}
library(Biostrings)
library(muscle)
alignment_muscle <- muscle(clytus_data)
```

##### 7c. NJ tree visualization of the alignment (using external package `ape` and `phangorn`)

```{r NJ tree,message=F,warning=F,fig.width=12, fig.height=12, dpi=300}
library(ape)
library(phangorn)
# Convert the alignment to DNABin
dna_bin <- as.DNAbin(alignment_muscle)

# Distance matrix using K80
dist_matrix <- dist.dna(
  dna_bin,
  model = "K80"
)

# Neighbor Joining tree
nj_tree <- nj(dist_matrix)

# midpoint rooting
nj_tree <- midpoint(nj_tree)

# Plot tree
plot(
  nj_tree,
  cex = 0.6,
  main = "Neighbor-Joining Tree"
)
```

*A simple NJ tree visualization of the aligned data*

### 8. Generate a `sf` object for spatial mapping of occurrences of the genus *Monochamus*

```{r sf_object}
# Filter for genus Monochamus
monochamus_search <- bold_parquet_search(
  input.parquet = parquet_file,
  taxonomy = "Monochamus",
  geography = "Canada",
  marker = "COI-5P"
)

# Collect the filtered data
monochamus_data <- bcdm_to_sf(monochamus_search)
```

The `sf` object can be used directly for mapping occurrences with external packages like `ggplot2`. (Background maps can be created in many ways. The `sf` library has been used in this case).

```{r occurrence-map,message=F,warning=F}
library(sf)
# Creating a background map using the maps package; some map_data country names (ID column) are changed to suit the BCDM country.ocean names
map_data <- st_as_sf(maps::map("world",
  plot = FALSE,
  fill = TRUE
)) %>%
  filter(ID == "Canada")
# Convert the data to WGS84
map_data <- st_transform(
  map_data,
  4326
)
# Plot
map_plot <- ggplot() +
  geom_sf(
    data = map_data,
    alpha = 0.3,
    linewidth = 0.4
  ) +
  geom_point(
    data = monochamus_data,
    mapping = aes(
      x = lon,
      y = lat
    ),
    colour = "#011B26",
    fill = "#F78E1E",
    size = 3,
    pch = 21
  ) +
  theme_bw(base_size = 15) +
  theme(panel.grid.major = element_line(
    colour = "grey50",
    size = 0.3,
    linetype = 3
  )) +
  xlab("Longitude") +
  ylab("Latitude") +
  coord_sf(expand = FALSE) +
  ggtitle("Distribution map")

map_plot
```

*Occurrence map of Monochamus*

# Performance and Integration Notes

## Benchmarking

Benchmarking used three representative queries:

1.  \~1 million records (taxonomy = "Hemiptera")
2.  \~4.5 million records (taxonomy = "Diptera", geography = "Costa Rica")
3.  \~10 million records (taxonomy = "Diptera")

Search (`bold_parquet_search`) and data retrieval (`bold_search_collect`) were benchmarked separately using the `rbenchmark` package with three replications on the *30th June 2026 data releaseß*(<https://doi.org/10.5883/DP-BOLD_Public.30-Jun-2026>). Tests were run on a *MacBook Air M2 (8-core CPU, 16 GB RAM)*. Across the three datasets, `bold_parquet_search` completed in 1.1 - 1.3s, while `bold_search_collect` required approximately 4 - 12 min, with runtime increasing as dataset size increased.

`BOLDNODE` can efficiently search and retrieve very large datasets; however, the data retrieval (`collect`) step is constrained by the user's available system memory, and collecting very large query results may exceed the machine's memory capacity.

## Integration with `BOLDconnectR`

Functions from BOLDNODE can be combined with the analysis functions from the `BOLDconnectR` package because BOLD data packages are based on the **Barcode Core Data Model (BCDM)** (<https://github.com/boldsystems-central/BCDM>). This shared data structure enables direct interoperability between packages without requiring additional data conversion.