---
title: "Apache Parquet-based Mass Spectrometry Data Backend"
output:
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{Apache Parquet-based Mass Spectrometry Data Backend}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\VignettePackage{MsBackendParquet}
---

```{r style, echo = FALSE, results = "asis", message = FALSE}
BiocStyle::markdown()
```

**Package**: `r Biocpkg("MsBackendParquet")`<br />
**Authors**: `r packageDescription("MsBackendParquet")[["Author"]] `<br />
**Compiled**: `r date()`

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

# Introduction

The `r Biocpkg("MsBackendParquet")` package provides a backend for
`r Biocpkg("Spectra")` objects that stores mass spectrometry (MS) data
in an on-disk dataset built on
[Apache Parquet](https://parquet.apache.org/).
Datasets are *written* with the
[Apache Arrow](https://arrow.apache.org/docs/r) R bindings and *queried*
with [DuckDB](https://duckdb.org/), keeping the spectra metadata
together with the `mz`/`intensity` peak values in a columnar format that
supports compression, dictionary encoding and Hive-style partitioning.

A dataset is a directory identified by an `mzStack.json`
manifest, holding one or more *runs*. A run is one source of spectra
and owns a contiguous block of the dataset's spectrum ids.
A run is either
`native` (converted from a raw MS data file, with peaks stored alongside
the metadata) or `mzpeak`, an index over external
[HUPO-PSI mzPeak](https://github.com/HUPO-PSI/mzPeak) archives that are
read where they are and never modified (see
*Querying collections of mzPeak archives* below). The manifest records
which kind each run is, so the two never have to be told apart by
inspecting the directory.

To keep access fast, the (small) non-peak spectra variables are held in
a bounded in-memory index while the peak data is read from disk on
demand. As a result, metadata filters cost about as much as with an
in-memory backend, while the bulk peak storage stays out of core.

The backend extends `MsBackendCached` from `r Biocpkg("Spectra")` so it
benefits from caching of user-defined spectra variables while keeping
the on-disk dataset read-only.

## Relation to other Spectra backends

`MsBackendMzR` from `r Biocpkg("Spectra")` keeps the spectra metadata in
memory and reads peaks from the raw mzML files on demand: nothing has to be
converted, but every peak access parses the files again.
`r Biocpkg("MsBackendSql")` stores spectra in an SQL database, which suits
data shared from a database server. `r Biocpkg("SpectraStash")` saves and
restores whole `Spectra` objects in interoperable formats such as HDF5. It
is for exchanging an object, not for querying it on disk.

`MsBackendParquet` converts once, like `MsBackendSql`, and then queries a
database engine, DuckDB. But its storage is plain Parquet files that need
no server and that other tools, such as Python, DuckDB or Spark, read
directly. It can also index mzPeak archives where they are, without
converting them.

# Installation

The package is installed from Bioconductor:

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("MsBackendParquet")
```

Reading raw mzML files needs the `r Biocpkg("mzR")` package, which
`BiocManager::install("MsBackendParquet", dependencies = TRUE)` installs
too. Datasets made from mzPeak archives or from in-memory data do not
need it.

# Creating a Parquet dataset

The examples below use a DDA data file of a pesticide mix, provided by
the `r Biocpkg("MsDataHub")` package. The file is downloaded once and
cached locally.

```{r data-file, message = FALSE}
library(Spectra)
library(MsBackendParquet)
library(MsDataHub)

fl <- MsDataHub::PestMix1_DDA.mzML()
```

`MsDataHub` returns the file under its cache name, and the id of a run is
taken from the name of the file it was imported from. Copying the file to
a named path first therefore gives the dataset readable run ids:

```{r name-file}
fl_named <- file.path(tempdir(), "PestMix1_DDA.mzML")
invisible(file.copy(fl, fl_named, overwrite = TRUE))
```

`mzMLToParquet()` converts raw MS data files into a Parquet dataset. It
validates the input files, writes the dataset and returns an initialised
backend, ready to be wrapped in a `Spectra` object:

```{r mzml-to-parquet}
path <- tempfile()
be <- mzMLToParquet(fl_named, path = path, verbose = FALSE)

sps <- Spectra(be)
sps
```

Each input file becomes its own run, written to
`<path>/spectra/run_id=<id>/`, so a filter that selects whole files
prunes at the directory level:

```{r run-dirs}
list.files(file.path(path, "spectra"))
```

`mzMLToParquet()` imports MS data via `MsBackendMzR()` and writes the
spectra metadata together with the `mz`/`intensity` values as Parquet
`list<double>` columns.

A `Spectra` object that is already loaded, with any backend, can instead be
moved onto Parquet with `setBackend()`. The destination is given as `path`:

```{r set-backend}
sps_mem <- Spectra(fl_named, source = MsBackendMzR())
sps_pq <- setBackend(sps_mem, MsBackendParquet(), path = tempfile())
sps_pq
```

# Accessing data

`spectraData()`, `peaksData()`, `mz()`, `intensity()` and related
methods work the same way as for any other `MsBackend`. The spectra
variables of the dataset are:

```{r spectra-variables}
spectraVariables(sps)
```

Values are returned as usual, with the peak data read from disk on
demand:

```{r access}
spectraData(sps, c("msLevel", "rtime", "precursorMz"))

head(peaksData(sps)[[1L]])

head(tic(sps))
```

# Filtering

Filter methods are translated into DuckDB SQL predicates that are pushed
down to the Parquet scan. DuckDB prunes whole row groups using their
min/max statistics, and skips whole partitions for Hive-partitioned
datasets, so a range or set filter reads only the data it needs. When
the filter column is held in the in-memory metadata index (see *The
in-memory metadata cache* below), the predicate is evaluated in memory
instead.

Selecting spectra by MS level and by retention time, which is given in
seconds:

```{r filter-basic}
filterMsLevel(sps, 2L)

filterRt(sps, c(200, 300))
```

Filters return a `Spectra` object and can be chained. Each call narrows
the predicate rather than reading data, so the combination is answered
by a single scan:

```{r filter-chained}
sps |>
    filterRt(c(200, 300)) |>
    filterMsLevel(2L) |>
    length()
```

MS2 spectra can be selected by precursor m/z, either over a range or
against a set of values with a mass tolerance in ppm:

```{r filter-precursor}
filterPrecursorMzRange(sps, c(278, 279))

ms2 <- filterPrecursorMzValues(sps, 278.093, ppm = 30)
ms2

range(rtime(ms2))
```

`filterPrecursorMzValues()` accepts several values at once, which become
one OR-ed predicate and therefore one scan rather than one scan per
value:

```{r filter-precursor-values}
filterPrecursorMzValues(sps, c(278.093, 124.086), tolerance = 0.01)
```

None of the filters above read peak data. They are answered from the
spectra variables, and the resulting predicate is carried into any later
peak access.

# Working with peaks

Peak data is read when it is asked for, and only for the spectra that
survived the filters. The nine MS2 spectra selected above share one
precursor, so their fragments can be compared directly:

```{r peaks}
lengths(ms2)

pk <- peaksData(ms2)[[which.max(lengths(ms2))]]
head(pk[order(-pk[, "intensity"]), ])
```

`containsMz()` from `r Biocpkg("Spectra")` reports which of those
spectra contain a peak at a given m/z:

```{r contains-mz}
which(containsMz(ms2, mz = 164.054, tolerance = 0.01))
```

Because the peaks are stored columnwise, DuckDB reads only the `mz` and
`intensity` columns, and only the rows belonging to the selected
spectra.

# Controlling the dataset layout

`createMsBackendParquetDataset()` writes a dataset without opening it
and exposes the layout options: the reader backend, chunk size,
compression codec, Parquet row-group size, and Hive partitioning.

```{r create-ms-backend}
path_part <- tempfile()
createMsBackendParquetDataset(path = path_part, x = fl,
                              partitioning = "msLevel",
                              compression = "zstd")

be_part <- backendInitialize(MsBackendParquet(), path = path_part)
length(be_part)
```

`partitioning` adds a further Hive level *below* the run.
The partition directories carry mzPeak column names (see *How the
dataset interfaces with mzPeak's model* below):

```{r partition-dirs}
list.files(list.files(file.path(path_part, "spectra"), full.names = TRUE))
```

The same function writes an existing `DataFrame` of spectra data when it
is passed as `data` instead of `x`.

# Querying collections of mzPeak archives

[mzPeak](https://github.com/HUPO-PSI/mzPeak) is the HUPO-PSI successor to
mzML. An mzPeak *archive* holds **one MS run** as a folder (or an
uncompressed ZIP) of Parquet files plus a JSON manifest naming them. That
structure is well suited to distributing and archiving individual runs,
but inefficient for queries spanning thousands of runs.
Answering such a query would require opening and parsing each archive and
its manifest.

`MsBackendParquet` bridges the two. `createMzPeakDataset()` registers a
set of archives as one dataset. The package ships two small archives for
this purpose:

```{r create-mzpeak-dataset}
archives <- system.file("extdata", c("QC01.mzpeak", "QC02.mzpeak"),
                        package = "MsBackendParquet")

qc_path <- tempfile()
createMzPeakDataset(archives[1L], path = qc_path, verbose = FALSE)
```

Further archives join the same dataset later. Only the index grows:

```{r add-mzpeak-archives}
addMzPeakArchives(qc_path, archives[2L], verbose = FALSE)

be_qc <- backendInitialize(MsBackendParquet(), path = qc_path)
sps_qc <- Spectra(be_qc)
sps_qc

table(basename(dataOrigin(sps_qc)))
```

With two runs in the dataset, `filterDataOrigin()` selects one of them:

```{r filter-data-origin}
filterDataOrigin(sps_qc, grep("QC02", dataOrigin(sps_qc), value = TRUE)[1L])
```

The archives are neither modified nor copied. Ingest reads each
archive's manifest, its spectrum metadata table and its scan, precursor
and selected-ion side tables, and writes a small derived index next to
them, around 3% of the size of the archives. It never reads a byte of
signal data, so the cost scales with the number of spectra rather than
the size of the data. Peaks are then read from the archives themselves,
on demand. Pass `link = "copy"` to place a copy of each archive inside
the dataset directory instead of referencing it in place.

Because the archives stay exactly as they were, they remain readable by
every other mzPeak tool, and exporting a run back out of the dataset is a
directory copy rather than a conversion.

## How the dataset interfaces with mzPeak's model

Several aspects of the data are translated or transformed on input, so it
is worth being precise about what changes.

- **Units and names.** mzPeak stores retention time as `time` in *minutes*;
  `Spectra` expects `rtime` in *seconds*. The conversion, and every other
  column rename and re-encoding (`scan_polarity` ±1, `spectrum_representation`
  as a CURIE, the isolation window as a target plus offsets, and so on),
  happens in a single SQL view. The same view serves a natively converted
  dataset: its `spectra/` files are written with mzPeak's column names and
  encodings too, so both kinds are directly comparable and either can be
  exported back to an archive mechanically. `spectrum_index` is always the
  dataset's unique 0-based key; a spectrum's position within its source file
  (`scanIndex`, which restarts per file) is kept in `scan_index`. The package
  adds only a handful of reserved columns (`spectrum_id_`, `acquisition_num_`,
  the `n_*` counts).
- **Flattening.** mzPeak keeps scans, precursors and selected ions in
  separate tables linked back to the spectrum by a foreign key. The index
  takes the first of each, and records how many there were in `n_scans`,
  `n_selected_ions` and `n_precursors`. The archive still holds the full
  detail, but a spectrum with several scans is summarised in the index.
- **Profile and centroid.** mzPeak can hold both representations of the same
  spectrum, in separate files. Choose with `representation`, which defaults
  to `"auto"`:

```{r mzpeak-repr}
be_centroid <- backendInitialize(MsBackendParquet(), path = qc_path,
                                 representation = "centroid")
be_centroid
```

- **Differing schemas.** Two perfectly valid archives often have different
  columns, because writers are encouraged to promote whichever parameters
  they have into columns. They still combine into one queryable dataset;
  columns missing from an archive read as `NA`.

Archives that use mzPeak's *chunked* signal layout are rejected with a clear
message; only the *point* layout is supported so far.

## Searching for peaks across runs

Filters on spectra variables work as they do for any other dataset.
Finding spectra that *contain a peak* at a given m/z is different: it has
to look at the signal, and mzPeak stores signal ordered by spectrum,
which means a filter on m/z cannot skip anything, because every block of
rows spans the whole mass range.

`filterContainsMz()` runs that search inside the storage engine rather
than pulling peaks into R:

```{r filter-contains-mz}
hits <- filterContainsMz(be_qc, c(201, 10101), ppm = 20)
length(hits)

basename(dataOrigin(hits))
```

`filterContainsMz()` works on the backend, and the result is wrapped in a
`Spectra` object like any other backend:

```{r filter-contains-mz-spectra}
Spectra(hits)
```

For datasets where this is a common query, build an m/z-ordered copy of
the signal. Sorting is what makes the block statistics useful, so the
same search then reads a handful of blocks instead of everything:

```{r build-projection}
buildProjection(qc_path, verbose = FALSE)

length(filterContainsMz(be_qc, 201, ppm = 20))
```

On a 5-million-point test dataset this made `filterContainsMz()` around
eight times faster. A projection is a cache: it is a derived copy of the
signal and so costs roughly as much disk again, it can be dropped with
`dropProjection()`, and results are identical with or without it. Only
the performance changes. Build it per run, for the runs you query often.

```{r drop-projection}
dropProjection(qc_path, verbose = FALSE)
```

# Sample metadata

An experiment records factors about the *sample* a run came from, such as
subject, timepoint, treatment group or injection order, that are not in
the raw data. `runData()` stores them in the dataset, one row per run,
keyed by `run_id`.

```{r run-data}
runData(be_qc) <- data.frame(run_id = c("QC01", "QC02"),
                             subject = c("S1", "S2"),
                             timepoint = c(0L, 6L))
runVariables(be_qc)
```

The columns then behave like any other spectra variable:

```{r run-data-access}
spectraData(be_qc, c("rtime", "subject", "timepoint"))
```

and `filterSampleData()` selects on them:

```{r filter-sample-data}
Spectra(filterSampleData(be_qc, timepoint == 6))
```

The sample metadata is held once per run, not copied onto every spectrum.
Each run owns a
*contiguous block* of `spectrum_id_`, so `timepoint == 6` is exactly a
range of ids, and `filterSampleData()` resolves it against a table of a
few thousand rows rather than reading a column of a few million. The
resulting predicate is a range on the column the files are ordered by,
which is the case Parquet's row-group statistics prune best, and it is
carried into any later `peaksData()`.

Keeping it out of the per-spectrum data also makes it correctable. Sample
metadata is the one thing that reliably changes after ingest. Rewriting it
touches a few kilobytes and leaves the signal, and any projections built
over it, untouched.

Runs with no row get `NA`, so metadata covering only some runs is fine. To
go further, and link samples to spectra in an experiment-wide object, pass
the same table to `r Biocpkg("MsExperiment")` and link the two on
`run_id`:

```{r link-sample-data, message = FALSE}
library(MsExperiment)

samples <- S4Vectors::DataFrame(run_id = c("QC01", "QC02"),
                                subject = c("S1", "S2"),
                                timepoint = c(0L, 6L))
mse <- MsExperiment(sampleData = samples, spectra = Spectra(be_qc))
mse <- linkSampleData(mse, with = "sampleData.run_id = spectra.run_id")
mse

## The spectra of the second sample.
spectra(mse[2])
```

# The in-memory metadata cache

The non-peak spectra variables are typically small compared to the peak
data, so `MsBackendParquet` keeps them in a bounded, read-through
in-memory index. The cache is held at the package level, keyed by
dataset path and shared by every backend over that dataset, so it is
independent of the individual, serialisable backend objects. It is
warmed on `backendInitialize()` and dropped automatically whenever the
dataset on disk is (re-)written. It is entirely separate from the
user-modified variables tracked by `MsBackendCached`, which always take
precedence.

The cache is controlled through three options, shown here with their
defaults:

```{r options}
## Enable ("auto", the default) or disable ("off") the metadata cache.
options(MsBackendParquet.cacheMetadata = "auto")

## Columns held by default (the ones the optimised filters need).
options(MsBackendParquet.cacheColumns =
            c("msLevel", "rtime", "precursorMz", "dataOrigin"))

## Budget in cells (rows x columns); 5e7 is roughly 400 MB of doubles.
## Requests over budget quietly fall back to reading from disk.
options(MsBackendParquet.cacheMaxCells = 5e7)
```

Two further options tune the underlying DuckDB session.
`MsBackendParquet.threads` sets how many threads DuckDB uses, and Arrow
too. When it is unset, both use at most two threads, and one when the
registered `r Biocpkg("BiocParallel")` backend is serial, so the package
does not take every core of a shared machine. DuckDB reads the option when
the package first opens its connection, and Arrow when the package is
loaded, so set it before either.
`MsBackendParquet.memoryLimit` maps to DuckDB's `memory_limit`, and is
left at DuckDB's default when unset.

# Streaming large mzML files

For very large mzML files, the default conversion path goes through a
`Spectra` object built with `MsBackendMzR()`, which materialises a
whole chunk of peaks in memory before writing it to Parquet. The
`engine = "mzr"` option bypasses `Spectra` entirely and streams
batches of spectra directly from `mzR::openMSfile()` into Apache
Arrow's incremental `ParquetFileWriter`:

```{r mzml-to-parquet-streaming}
be_stream <- mzMLToParquet(fl, path = tempfile(),
                           engine = "mzr", batch_size = 1000L,
                           verbose = FALSE)
length(be_stream)
```

Memory use stays bounded by `batch_size * mean(peaks per spectrum)`,
regardless of total file size. The streaming engine processes files
serially and requires the `mzR` package.

# Performance

Storing the peaks as columnar `list<double>` values lets DuckDB skip
whole row groups when a query touches only part of the data, while
projection and predicate push-down keep both the columns and rows read
to a minimum. Combined with the in-memory metadata index, this makes
metadata filters roughly as fast as an in-memory backend while the peak
storage stays on disk.

Across a suite of seven queries over 380,100 spectra (50 mzML files),
the Parquet backend is the fastest of the on-disk `Spectra` backends on
5 of 7 queries (medians in ms, lower is better):

| query | MsBackendParquet | MsBackendMzR | MsBackendHdf5Peaks | MsBackendOfflineSql |
|---|---|---|---|---|
| `filterRt()`, wide range | **12.24** | 15.44 | 24.05 | 380.19 |
| `filterRt()`, narrow range | 11.23 | **5.33** | 14.76 | 348.17 |
| `filterPrecursorMzValues()`, ppm | **14.09** | 21.66 | 28.66 | 336.33 |
| `filterRt()` and `filterPrecursorMzValues()` | **12.56** | 17.25 | 27.75 | 389.22 |
| `peaksData()` of MS1 spectra | **314.82** | 1510 | 1440 | 919.04 |
| `tic()` of MS1 spectra | **506.79** | 1390 | 1490 | 890.82 |
| `spectraData()` of MS1 spectra | 26.09 | **17.27** | 22.30 | 718.91 |

The benchmark scripts live in
[`inst/benchmarks/`](https://github.com/computational-metabolomics/MsBackendParquet/tree/devel/inst/benchmarks);
see `inst/benchmarks/performance.Rmd` for the full numbers, the
configuration sweeps over the tunable options, and the design
considerations behind them.

# Saving and reloading the backend

Because `MsBackendParquet` objects only contain a file system path,
they are fully serialisable and can be sent to parallel workers
without re-opening connections:

```{r save-load}
rds <- tempfile(fileext = ".rds")
saveRDS(be, rds)
be2 <- readRDS(rds)
identical(spectraNames(be), spectraNames(be2))
```

# Session information

```{r si}
sessionInfo()
```