Classifying Microbial Species From Metagenomic Sequencing

Metagenomic classification looks like a machine-learning problem and the method that dominates it is not one. It is exact substring matching against a taxonomy, and understanding why that wins tells you exactly where it breaks.

What the classifier is being asked

Shotgun metagenomic sequencing takes DNA from a sample containing many organisms and sequences it without isolating anything. The output is tens of millions of short reads from an unknown mixture. The question is which organism each read came from, and the answer has to come from a reference database of known genomes.

Three constraints shape every method. The read is short, commonly 150 bases, which is not much evidence. The number of reads is enormous, so per-read cost must be tiny. And the correct answer is frequently not a species at all: a read from a conserved gene is genuinely consistent with an entire genus, and a classifier that names a species anyway is producing a confident wrong answer rather than a useful one.

That third constraint is what rules out framing this as flat classification over species labels. The label space is a tree, and the right output is a node in it at whatever depth the evidence supports.

Exact k-mers and lowest common ancestors

The dominant approach builds a map from k-mer to taxonomy node, once, offline. For every k-mer in every reference genome, find every genome it occurs in, and store the lowest common ancestor of those genomes’ taxonomy nodes. A k-mer unique to one strain maps to that strain; one shared across all enterobacteria maps to the family; one shared across all life maps to the root and is uninformative.

Kraken 2 is the reference implementation of this idea, published by Derrick Wood, Jennifer Lu and Ben Langmead in Genome Biology in 2019. Its published defaults are k equal to 35 and a minimiser length ell equal to 31. The minimiser is the mechanism that makes the table fit in memory: rather than storing all 35-mers, it stores, for each 35-mer, the lexicographically smallest canonical 31-mer inside it, applies a spaced seed mask, and hashes that into a compact table. Adjacent 35-mers along a genome usually share a minimiser, so a run of overlapping k-mers collapses to a single table entry, and the database shrinks by roughly an order of magnitude relative to storing every k-mer.

Classification of a read is then: extract its k-mers, look each one up, and collect the taxonomy nodes they hit. The read is assigned by scoring every root-to-leaf path through the taxonomy by the number of the read’s k-mers that map to nodes on that path, and taking the leaf of the highest-scoring path. Because the score accumulates along a path, a read supported by many genus-level k-mers and a few species-level ones lands at the species, and a read whose evidence is spread across two genera lands at their common ancestor.

One read, worked

A 150-base read at k equal to 35 contains 150 minus 35 plus 1 equals 116 k-mers. Suppose the lookups come back like this:

total k-mers                                116

  map to Escherichia coli (species)          62
  map to Escherichia (genus)                 18
  map to Enterobacteriaceae (family)          9
  map to root / uninformative                 2
  no hit in database                         25

Score the paths. The path root to Enterobacteriaceae to Escherichia to E. coli collects 62 plus 18 plus 9 equals 89, because every k-mer at or above a node on the path counts toward it. Any path through a different genus collects only the 9 family-level k-mers plus the 2 root ones. The read is assigned to E. coli.

The 25 no-hit k-mers are worth pausing on. They are consistent with several things: sequencing errors, each of which destroys the 35 k-mers overlapping it, so a single miscalled base accounts for up to 35 of them; strain-specific sequence absent from the reference; or a genuinely different organism. Twenty-five is consistent with one error near the middle of the read.

This is what the confidence threshold controls. It is the fraction of a read’s k-mers that must map to the assigned path, and the default is zero — meaning a read is assigned on any evidence at all, including a single k-mer. Here the fraction is 89 over 116, about 0.77, so the call survives any reasonable threshold. A read with 3 supporting k-mers out of 116, a fraction of 0.026, would also be reported as E. coli at the default and would vanish at a threshold of 0.05. Raising the threshold trades sensitivity for precision, and leaving it at the default is the single most common reason a report contains species that are not in the sample.

Read counts are not abundance

A classification report gives reads per taxon, and converting that to “what fraction of the community is this organism” requires two corrections that are frequently skipped.

• Genome length. Shotgun sequencing samples DNA, not organisms. A 6-megabase genome yields twice the reads of a 3-megabase genome at equal cell counts, so read fraction is a DNA fraction and converting to a cell fraction means dividing by genome length.
• Reads stuck above the species. A large share of reads land at genus or family because they fall in conserved regions. Those reads belong to some species; they just could not be assigned to one. Abundance re-estimation methods redistribute them down to the species level in proportion to the species-level evidence already observed, which is a Bayesian correction and is a separate step from classification. Reporting only the species-level read counts systematically understates every taxon and understates the abundant ones most.

What the database decides for you

Every property of the output is downstream of the reference database, and this is the dominating limitation of the whole approach.

You can only find organisms whose genomes are in it. Reads from an unsequenced organism either miss entirely or, worse, land on its nearest sequenced relative, producing a confident and wrong species call. The database is also enormously biased toward what people have had reason to sequence — clinical isolates and model organisms are represented by thousands of strains, environmental organisms by a handful — so classification rates differ by ecosystem for reasons that have nothing to do with the sample. And host DNA will swamp anything from a tissue or swab sample unless the host genome is either in the database, so reads can be assigned to it and set aside, or removed by alignment beforehand.

Two practical consequences. Include the host genome and common laboratory contaminants in the database deliberately, so those reads are labelled rather than mis-assigned. And when an unexpected species appears in a report, treat it as a hypothesis to be checked — the discriminating evidence is whether reads cover its genome evenly or pile up on a few conserved regions, which is the diagnostic worked through in detecting contamination in sequencing data.