Profiling and feature tables¶

Estimate who is in a community and how much of each is present, and turn those estimates into a feature table (a sample × feature abundance matrix). Two complementary methods: woltka_ogu builds an OGU feature table from reads you have already aligned to a reference, while sylph_profile profiles raw shotgun reads directly against a pre-built reference database. Both emit long-form output that feeds straight into diversity analysis and can be written to BIOM with COPY ... (FORMAT BIOM).

Table of Contents¶

woltka_ogu - OGU (Operational Genomic Unit) feature table from alignments, global or per-sample.
sylph_profile - FracMinHash relative-abundance profiling of shotgun reads against a .syldb database.
sylph_index_create - Build a sylph .syldb or two-stage .syl2db reference database from a table.

`woltka_ogu`¶

Compute Woltka OGU (Operational Genomic Unit) counts over SAM-like alignment data. Implements Woltka's classification algorithm, assigning reads to taxonomic units while fractionally distributing multi-mapped reads. When the optional sample_id named parameter is supplied, the aggregation runs in parallel across distinct sample values — one DuckDB query per sample on a dedicated per-thread connection — which bounds memory to a single sample's footprint.

Function signature:

woltka_ogu(relation, sequence_id_field [, sample_id])

Parameters: - relation (VARCHAR): Name of a table or view containing SAM-like alignment data (must be a resolvable catalog name; pass as a string literal) - sequence_id_field (VARCHAR): Name of the column holding sequence identifiers — typically read_id, or a numeric index column for better hash performance - sample_id (VARCHAR, named, optional): Name of the column holding sample identifiers. When supplied, adds that column to the output and runs the aggregation per distinct sample value in parallel.

Required columns in relation: - Column named by sequence_id_field: read/sequence identifier - reference (VARCHAR, BIGINT, or UUID): reference sequence name (becomes feature_id, preserving the storage type). A non-identifier type is rejected at bind time. - flags (USMALLINT): SAM alignment flags - When sample_id is supplied: the named column (any comparable type) — NULLs are rejected at bind time

Output schema: - When sample_id is omitted: (feature_id <reference_type>, value DOUBLE) — feature_id preserves the reference column's storage type (VARCHAR/BIGINT/UUID), mirroring align_minimap2/alignment_slice id-type preservation. - When sample_id is supplied: (<sample_id_column> <its_type>, feature_id <reference_type>, value DOUBLE) — the first column's name matches the value you passed to sample_id, and both id columns preserve their input storage types.

This long-form (sample_id, feature_id, value) shape is exactly what the diversity functions consume and what COPY ... (FORMAT BIOM) writes.

Algorithm: 1. Orients reads using alignment flags (forward/reverse via alignment_is_read1). 2. For each read orientation, divides 1 by the number of unique features aligned to. 3. Aggregates fractional counts per feature (and per sample, when sample_id is used).

Correctness assumption for per-sample mode: read IDs are unique across samples. The per-sample subset then yields the same distribution as the global aggregation.

Examples:

-- Global aggregation across the whole relation
SELECT * FROM woltka_ogu('my_alignments', 'read_id');

-- Per-sample aggregation — one aggregation per distinct sample value, in parallel
SELECT * FROM woltka_ogu(
    'my_alignments',
    'read_id',
    sample_id := 'sample_id'
);

-- Filter high-quality alignments via a view, then classify
CREATE OR REPLACE VIEW primary_alignments AS
    SELECT * FROM read_alignments('alignments.bam')
    WHERE alignment_is_primary(flags) AND mapq >= 20;

SELECT * FROM woltka_ogu('primary_alignments', 'read_id');

-- Multi-sample: union sources under a sample_id column, then classify per sample
CREATE OR REPLACE VIEW all_samples AS
    SELECT *, 'sample1' AS sample_id FROM read_alignments('sample1.bam')
    UNION ALL
    SELECT *, 'sample2' AS sample_id FROM read_alignments('sample2.bam')
    UNION ALL
    SELECT *, 'sample3' AS sample_id FROM read_alignments('sample3.bam');

SELECT * FROM woltka_ogu('all_samples', 'read_id', sample_id := 'sample_id')
ORDER BY sample_id, feature_id;

-- Export per-sample results to BIOM for downstream analysis
COPY (
    SELECT * FROM woltka_ogu('my_alignments', 'read_id', sample_id := 'sample_id')
) TO 'ogu_table.biom' (FORMAT BIOM, COMPRESSION 'gzip');

-- Export global results to BIOM (add a sample_id column)
COPY (
    SELECT feature_id, 'MySample' AS sample_id, value
    FROM woltka_ogu('my_alignments', 'read_id')
) TO 'ogu_single.biom' (FORMAT BIOM);

Behavior: - Arguments relation and sequence_id_field (and sample_id if used) must be quoted string literals; they are catalog/column names resolved at bind time. - Multi-mapped reads are fractionally assigned: each mapping receives weight 1/N where N is the number of unique references the read maps to, within the same orientation. - Read orientation (forward/reverse) is considered separately via the alignment_is_read1() flag — paired-end reads are handled automatically. - For better performance on large datasets, add a numeric index column and pass it as sequence_id_field instead of read_id. - Output row order is non-deterministic when sample_id is used (parallel per-sample execution). Use an explicit ORDER BY if stable ordering is required. - sample_id mode cannot read a TEMP table or view. Every other relation parameter in miint accepts CREATE TEMP TABLE/CREATE TEMP VIEW output, but the per-sample path builds a fixed-name temporary view on a private per-thread connection, which must stay isolated so parallel workers do not collide on that name. Materialize the source as an ordinary table (or a view in a schema) first. Tracked as #207.

`sylph_index_create`¶

Builds a sylph reference database from a DuckDB table or view of reference sequences, then writes it to disk. The in-SQL counterpart to the upstream sylph sketch CLI: reference genomes come from a table (e.g. the output of read_fastx over reference FASTAs) rather than file paths, so the reference set can be assembled with ordinary SQL (joins, filters, unions). The resulting file is consumable by sylph_profile and by the upstream sylph CLI.

The output format is chosen by the two_stage parameter, or, when it is not given, by the output_path suffix. Two-stage output (two_stage := true, or a .syl2db suffix) writes sylph 1.0's two-stage database (what sylph convert-db-two-screen produces): a small stage-1 screen index that is loaded in full, plus per-genome compressed dense blocks that are read on demand only for genomes passing the screen. Profiling against it scales with the sample rather than the reference, so it is the right choice for large (atlas-scale) references of typical-size microbial genomes. Upstream advises against it for small genomes such as viruses and plasmids, which the sparse stage-1 screen leaves with too few k-mers. Otherwise a plain .syldb is written.

Sketching is done from in-memory sequence bytes via the same sylph 1.0.0-miint Rust static library that backs sylph_profile; a genome built here is byte-identical to one from sylph sketch given the same contigs in the same order. Linux and macOS only; not registered on WASM or Windows builds.

Function signature:

sylph_index_create(source_table, output_path, genome_id='col', [options])

Parameters:

source_table (VARCHAR, positional): Name of a table or view with columns read_id (VARCHAR or BIGINT) and sequence1 (VARCHAR). Each row is one contig; sequence1 is the contig sequence and read_id is its name. If a comment column is also present (as in read_fastx output), the contig name is reconstructed as the full FASTA header read_id || ' ' || comment — matching how sylph sketch (via needletail) stores first_contig_name. Without a comment column, the contig name is read_id as-is.
output_path (VARCHAR, positional): Destination path. Overwritten if it exists. When two_stage is not given, a .syl2db suffix selects the two-stage format.
two_stage (BOOLEAN, default = inferred from the output_path suffix): Write a two-stage .syl2db database instead of a plain .syldb. An explicit value wins over the suffix, as FORMAT does for COPY; readers detect the format from the file, not the name.
genome_id (VARCHAR, required): Name of the grouping column. One GenomeSketch is produced per distinct value; all contigs sharing a value are merged into that genome (matching sylph sketch's default whole-file behavior). A natural choice is the filepath column from read_fastx(..., include_filepath := true) — one genome per reference FASTA.
order_by (VARCHAR, default sequence_index): Integer column that orders the pieces of a contig stored across several rows before they are joined, and picks each genome's first_contig_name: its contig with the lowest value. The default is the column read_fastx adds, so for a genome from one FASTA file the first contig is the first in the file, as in sylph sketch; a source without it names another column. Among contigs that share the lowest value (a genome built from several files, where sequence_index restarts in each), the one whose name sorts last is picked. Values are cast to BIGINT as each genome is read, and one that does not cast (a non-numeric read_id, for example) fails the build with an Invalid Input Error naming the genome and giving DuckDB's Conversion Error.
k (INTEGER, default 31): k-mer size. Only 21 and 31 are supported. Must match the k used at profile time.
c (INTEGER, default 200): FracMinHash subsampling rate. Must be ≥ the c of any sample sketch profiled against the database.
min_spacing (INTEGER, default 30): Minimum k-mer spacing; thins densely-seeded regions.
pseudotax (BOOLEAN, default true): Track min-spacing-dropped k-mers so the database supports pseudotax/profiling mode. Set false for query-only databases. Not allowed with a two-stage output, which always needs profiling k-mers.
threads (INTEGER, default 0): Number of genomes sketched at once, each on its own thread. 0 uses DuckDB's threads setting; another value sets the count for this call only, above or below that setting, without changing it. The count is capped at the number of genomes.
screen_c (INTEGER, default 3000, two-stage only): FracMinHash rate of the stage-1 screen index. Must be ≥ c; larger is a smaller, faster, more permissive screen. Matches --screen-c on sylph convert-db-two-screen.
min_sparse_kmers (INTEGER, default 50, two-stage only): Minimum screen k-mers per genome; genomes that would fall short at screen_c get a denser, genome-specific screen rate. Matches --min-sparse-kmers.

Output schema: a single status row.

Column	Type	Description
`output_path`	VARCHAR	The path the database was written to
`k`	INTEGER	k-mer size used
`c`	INTEGER	FracMinHash subsampling rate used
`num_genomes`	UBIGINT	Number of genomes written to the database
`status`	VARCHAR	`'ok'` on success

Behavior:

Grouping: genomes are the distinct non-NULL genome_id values. With more than one thread, the order of the genomes in the written file can differ from run to run; their sketches do not. Within each genome, rows sharing a contig name are joined into one contig in order_by order.
Memory: each thread reads and sketches one genome's contigs at a time; only the finalized (compact) genome sketches accumulate before the single write at the end.
Synchronous: the build runs as a side effect of the query; the status row is emitted once it completes.

Examples:

-- Build a .syldb from a directory of reference FASTAs, one genome per file.
CREATE TABLE refs AS
  SELECT * FROM read_fastx('refs/*.fasta.gz', include_filepath := true);

SELECT num_genomes
FROM sylph_index_create('refs', 'my_refs-c200.syldb',
                        genome_id := 'filepath',
                        order_by  := 'sequence_index');

-- Then profile reads against the freshly built database.
SELECT genome_name, taxonomic_abundance
FROM sylph_profile('reads', 'my_refs-c200.syldb')
ORDER BY taxonomic_abundance DESC;

Error conditions:

Error if source_table does not exist or is missing required columns (read_id, sequence1).
Error if the genome_id parameter is not supplied, or if the genome_id / order_by column does not exist.
Error if an order_by value does not cast to BIGINT, naming the genome being read (when several genomes fail, the first in genome_id order); error if threads is negative.
Error if k is not 21 or 31.
Error if the output is two-stage and pseudotax := false, or screen_c < c; error if screen_c / min_sparse_kmers are given for a plain .syldb output.
Error if no non-NULL genome_id values are found.
IO error surfaces the underlying sylph diagnostic string (e.g., failure to write the output path).

`sylph_profile`¶

FracMinHash-based relative-abundance profiling of paired-end shotgun metagenomic reads against a pre-built reference database, using sylph (Shaw & Yu 2024, Nature Biotechnology). Sequences come from a DuckDB table or view — there is no FASTQ path argument — and the database is loaded once per call. The result is streamed back via the Arrow C Data Interface (zero-copy).

Both sylph database formats are accepted and detected from the file's contents, not its name: a plain .syldb (all genome sketches resident) or sylph 1.0's two-stage .syl2db (screen index resident, dense blocks decoded on demand for genomes that pass the stage-1 screen). For the genomes it reports, a two-stage database gives the same numbers as the .syldb it was converted from.

Embedded as a Rust static library (sylph 1.0.0-miint fork; MIT). Linux and macOS only; the function is not registered on WASM or Windows builds.

Function signature:

sylph_profile(source_table, syldb_path, [sample_id='col'], [options])

Parameters:

source_table (VARCHAR, positional): Name of a table or view with columns read_id (VARCHAR), sequence1 (VARCHAR), and optionally sequence2 (VARCHAR). When sequence2 is present and non-empty per row, the read pair is processed paired-end; otherwise the call is single-end.
syldb_path (VARCHAR, positional): Path to a sylph .syldb or .syl2db reference database, built either via sylph_index_create or the upstream sylph CLI (sketch, convert-db-two-screen). The file is opened read-only and shared across the call.
sample_id (VARCHAR, optional): Name of a column on source_table to partition by. When set, the function fans out per-sample (parallelized via the per-sample helper used by deblur / align_mafft / detect_chimera_uchime_denovo) and prepends a sample_id column to the output. Without this option, the entire table is processed as a single sample.
min_ani (DOUBLE, default = sylph default): Minimum adjusted ANI (percent, 0..100) for a genome to be reported. Negative or unset = use sylph's built-in default.
min_number_kmers (UINTEGER, default = 10): Minimum number of sampled k-mers a reference genome must have to be considered (a size floor: ≈ genome length / c). sylph 1.0 lowered the default from 50 so genomes under ~10 kb are no longer excluded outright; pass 50 for the 0.9 behavior.
min_contain (UINTEGER, default = 7): Minimum number of contained k-mers (observed in the sample) for a genome to count as a hit — the evidence floor sylph 1.0 added alongside the lower min_number_kmers. 0 = sylph default.
min_count_correct (DOUBLE, default = 3.0): Minimum corrected k-mer count for genome detection. Lower values increase sensitivity at the cost of false positives.
screen_ani (DOUBLE, default = 0.85, .syl2db only): Minimum adjusted ANI, as a fraction in [0, 1], for a genome to pass the two-stage stage-1 screen. Deliberately permissive; the dense stage restores specificity. Ignored for .syldb databases.
estimate_unknown (BOOLEAN, default false): Renormalize taxonomic_abundance to fraction-of-reads-explained, accounting for unknown / unmatched fraction.
dedup_paired_reads (BOOLEAN, default true): Enable approximate paired-read deduplication during sketching. Matches sylph CLI behavior.
dedup_fpr (DOUBLE, default 0.0001): Bloom-filter false-positive rate used during sketch-time deduplication. Matches --fpr on the upstream CLI.
threads (UINTEGER, default 0 = auto): Rayon thread count for the inner sylph compute. In per-sample mode, the inner count auto-balances against DuckDB's worker count so that outer × inner ≈ db_threads — same total core utilization as sylph profile -t db_threads.

Output schema:

Column	Type	Description
`sample_id`	(matches source column type)	Only present when `sample_id` is set; the partition value for the row
`genome_index`	UINTEGER	0-based row offset of the genome within the `.syldb`
`genome_name`	VARCHAR	Genome identifier from the sketch (typically the reference FASTA path or accession)
`contig_name`	VARCHAR	First contig name from the genome's sketch — disambiguator when multiple genomes share `genome_name`
`sequence_abundance`	DOUBLE	Fraction of the sample's reads explained by this genome (0..1)
`taxonomic_abundance`	DOUBLE	Coverage-corrected relative abundance (sums to ~1.0 across all reported genomes, or to "explained fraction" when `estimate_unknown := true`)
`adjusted_ani`	DOUBLE	Lambda-adjusted average nucleotide identity (percent) between sample reads and the genome sketch
`eff_cov`	DOUBLE	Effective sequencing coverage of the genome by the sample
`naive_ani`	DOUBLE	Unadjusted ANI (percent) — useful for comparison against `adjusted_ani` to spot low-coverage / high-correction cases
`kmers_reassigned`	UBIGINT	Number of k-mers reassigned away from this genome by the winner-take-all step

Adding columns at the end is non-breaking; reordering is.

Behavior:

Single-vs-paired: auto-detected per-row. Rows whose sequence2 is NULL or empty are sketched single-end; rows with both sequences set are sketched paired. Mixing within one sample is allowed.
Database lifecycle: the database is loaded once per sylph_profile call into the GlobalState and reused across all execute invocations on that call. Back-to-back calls against the same DB pay the load cost each time (a future session-scoped cache is on the roadmap). For a .syl2db the fixed cost is the screen index only; dense blocks are decoded per call for the genomes that pass the screen.
Ordering: order_preservation_type = NO_ORDER — rows may interleave across genomes when running per-sample. Sort downstream if needed.
Numerical parity: bit-identical to upstream sylph profile --reads on the same data, modulo a small documented set of divergences in the embedded library.

Examples:

-- Profile a metagenome against the GTDB sylph DB (single sample).
CREATE TABLE reads AS
  SELECT read_id, sequence1, sequence2
  FROM read_fastx('sample_R1.fq.gz', sequence2 := 'sample_R2.fq.gz');

SELECT genome_name, taxonomic_abundance, adjusted_ani, eff_cov
FROM sylph_profile('reads', '/refs/gtdb-r220-c200-dbv1.syldb')
ORDER BY taxonomic_abundance DESC
LIMIT 20;

-- Per-sample profiling: one row stream per sample_accession, parallelized.
CREATE TABLE reads AS
  SELECT sample_accession, read_id, sequence1, sequence2
  FROM 'PRJNA000000.parquet';

SELECT *
FROM sylph_profile('reads', '/refs/gtdb-r220-c200-dbv1.syldb',
                   sample_id := 'sample_accession',
                   estimate_unknown := true);

Error conditions:

Error if source_table does not exist or is missing required columns (read_id, sequence1).
Error if syldb_path cannot be opened, is corrupt, or is neither a sylph .syldb nor .syl2db file.
Error if sample_id is set but the named column does not exist on source_table or collides with an output column name.
IO error surfaces the underlying sylph diagnostic string (e.g., truncated database, unsupported version).
sample_id mode cannot read a TEMP table or view. Every other relation parameter in miint accepts CREATE TEMP TABLE/CREATE TEMP VIEW output, but the per-sample path builds a fixed-name temporary view on a private per-thread connection, which must stay isolated so parallel workers do not collide on that name. Materialize the source as an ordinary table (or a view in a schema) first. Tracked as #207.

Profiling and feature tables¶

Table of Contents¶

woltka_ogu¶

sylph_index_create¶

sylph_profile¶

See also¶

`woltka_ogu`¶

`sylph_index_create`¶

`sylph_profile`¶