Bio Metagenomics Kraken

Version Compatibility

Reference examples tested with: Kraken2 2.1.3+, Bracken 2.9+, KrakenTools 1.2+, pandas 2.2+.

Before using code patterns, verify installed versions match. If versions differ:

• CLI: kraken2 --version, bracken -h, kraken2 --help to confirm flags and defaults
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

The DATABASE is the version that matters most. A taxon absent from the database is invisible no matter how the binary is configured. Record the Kraken2 database build (Standard / PlusPF / PlusPFP / nt / GTDB, and whether capped to 8/16 GB) and the Bracken databaseRLENmers.kmer_distrib read length, which must match both the Kraken2 database and the actual read length. --minimum-hit-groups enforcement varied across older 2.x point releases; confirm defaults with kraken2 --help on the installed build.

Kraken Classification

"What's in my metagenome?" -> Match each read's k-mers to a reference database by lowest-common-ancestor, then re-estimate abundance with Bracken - because the database, not the algorithm, decides what can be found.

• CLI: kraken2 --db DB --paired R1.fq.gz R2.fq.gz --report out.kreport --confidence 0.1 --output out.kraken

Scope: read-based, assembly-free taxonomic classification of shotgun reads, plus the Bracken handoff. Marker-gene profiling -> metaphlan-profiling. Bracken command mechanics -> abundance-estimation. Genome/MAG recovery -> genome-assembly/metagenome-assembly. Host removal and read QC -> read-qc/contamination-screening, contamination-controls. Amplicon/16S -> the microbiome category.

The Single Most Important Modern Insight -- A Kraken Report Is a Database-Conditioned Similarity Ledger, Not a Sample Inventory

Kraken reports what each read most resembles in THIS database - never what is truly present, and never how much. Every number is hostage to three choices made before the run: the database, the confidence threshold, and the assumption that read count means abundance. Three corollaries each common misuse violates:

1. Classification is not presence. At --confidence 0 with the LCA rule, one shared k-mer labels a read with its nearest database relative even when the true organism is absent. Absence-from-database becomes a confident wrong species.
2. Read count is not abundance. Counts scale with genome length and copy number, so the report percentage is a fragment fraction, not a cell fraction. Bracken fixes the wrong-rank problem; it does not fix this.
3. A taxon at the bottom of the report is a hypothesis, not a finding. Single-region hits, hash collisions, and contaminated references populate the long tail. Unique-minimizer coverage separates a real low-abundance organism from a phantom.

Organize the analysis around defending against these three, not around listing flags. Kraken2 at defaults over-classifies; Kraken2 tuned (right database + confidence + hit-groups + a unique-k-mer floor + host removal) is competitive with any classifier.

Why "Exact K-mer" Is Wrong for Kraken2

Kraken1 (Wood & Salzberg 2014 Genome Biol 15:R46) stored every exact k-mer. Kraken2 (Wood 2019 Genome Biol 20:257) replaced that with three ideas that make it fast and lean but PROBABILISTIC: (a) minimizers collapse each k-mer (default k=35) to the smallest hashed l=31-mer in its window; (b) a spaced seed (s=7 masked positions) tolerates errors at "don't care" positions; (c) a compact hash table stores only high bits of each key. The compact hash can return a wrong or spurious LCA on collision - which is exactly why the precision levers below exist. Calling Kraken2 "exact k-mer" hides where false positives come from.

Tool Taxonomy

Decision Tree by Scenario

Basic Classification

# Paired-end, with the precision levers that defaults omit
kraken2 --db "$KRAKEN_DB" \
    --paired --gzip-compressed --threads 8 \
    --confidence 0.1 \            # raise from default 0 to suppress single-k-mer false positives
    --minimum-hit-groups 2 \      # require >=2 distinct hit regions (default 2; raise to 3 for clinical)
    --report out.kreport \
    --output out.kraken \
    R1.fq.gz R2.fq.gz

--paired joins mates with a k-mer-breaking N and classifies the pair as one fragment, raising specificity. The per-read --output (large) can be dropped to /dev/null once the .kreport is what feeds Bracken. --memory-mapping runs without loading the database into RAM (slower; for low-memory hosts).

False-Positive Control: Unique Minimizers

Goal: Separate a real low-abundance organism from a single-region phantom before believing any tail taxon.

Approach: Enable --report-minimizer-data so the report carries distinct-minimizer counts; a taxon with many reads but few distinct minimizers is hitting one conserved region and is a red flag. KrakenUniq's HyperLogLog unique-k-mer count is the heavier-weight version of the same signal.

kraken2 --db "$KRAKEN_DB" --paired --confidence 0.1 \
    --report-minimizer-data \    # inserts 2 columns: total + DISTINCT minimizers (shifts later columns)
    --report out.kreport --output /dev/null \
    R1.fq.gz R2.fq.gz
# A species with high reads but low distinct-minimizers = false positive (one region lit up repeatedly).

Calibrate a unique-k-mer floor against negative controls rather than hard-coding one; the clinical ">=1024 unique k-mers" cutoff is dataset/database-specific folklore, not a constant.

Build a Custom Database

Goal: Build a database whose contents define exactly the detectable universe (and include human for host capture).

Approach: Download taxonomy, add the libraries the question needs (including human), build the minimizer index, then build the matching Bracken distributions at the actual read length.

kraken2-build --download-taxonomy --db custom_db
for lib in bacteria archaea viral human UniVec_Core; do
    kraken2-build --download-library "$lib" --db custom_db
done
kraken2-build --build --db custom_db --threads 16    # writes hash.k2d, opts.k2d, taxo.k2d
kraken2-build --clean --db custom_db                 # drop library/ + taxonomy/ to shrink
bracken-build -d custom_db -t 16 -k 35 -l 150        # -k MUST equal the Kraken2 k (35); -l = read length

kraken2-build --special gtdb builds a GTDB-taxonomy database (curated; the greengenes/silva/rdp special downloads have rotted). --max-db-size randomly downsamples k-mers to fit a cap - this is how the prebuilt 8gb/16gb databases are made, and the reason confidence collapses classification on them.

Hand Off to Bracken

Kraken strands reads at the shared genus when species share k-mers; Bracken redistributes them down using genome-derived priors. It fixes the wrong-rank problem only - never genome-size bias, and never false positives (it can amplify or even invent a species by reassigning an absent organism's reads to its nearest congener). Run FP control first. Command mechanics live in abundance-estimation:

bracken -d "$KRAKEN_DB" -i out.kreport -o out.bracken -w out.bracken.kreport \
    -r 150 \   # MUST match a built databaseRLENmers.kmer_distrib AND the actual read length
    -l S -t 10 # species level; -t is a redistribution floor (drops taxa with fewer than 10 clade-level reads, strict <), not a confidence

Per-Method Failure Modes

Over-classification at default confidence

Trigger: running --confidence 0 and reporting the species list. Mechanism: one shared k-mer can classify a read; the LCA labels it with its nearest database relative. Symptom: hundreds of low-abundance species, many biologically implausible. Fix: --confidence 0.1-0.4 (database-dependent) plus --minimum-hit-groups >=2; verify the tail with unique minimizers.