bio-metagenomics-abundance

What this skill does

## Version Compatibility

Reference examples tested with: Bracken 2.9+, Kraken2 2.1.3+, pandas 2.2+, scikit-bio 0.6+, R zCompositions 1.5+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `bracken -h`, `bracken-build -h` to confirm flags and defaults
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('zCompositions')` then `?cmultRepl` to verify parameters

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

The Bracken `databaseRLENmers.kmer_distrib` is built per database at a fixed read length and MUST match both the Kraken2 database and the actual (post-trim) read length; `-r` is not auto-detected and a mismatch silently biases every species fraction. Bracken needs the default Kraken2 report format, not mpa-style.

# Abundance Estimation

**"How much of each taxon is in my sample?"** -> Re-estimate species reads with Bracken, then treat the table as a composition - because the sequencer fixed the total, so the numbers are relative and a change in one forces apparent changes in the rest.
- CLI: `bracken -d DB -i kraken.kreport -o out.bracken -w out.bracken.kreport -r 150 -l S -t 10`

Scope: Bracken mechanics plus what happens to the numbers afterward - estimand choice, compositional transforms, normalization, reference-frame DA logic, absolute conversion. Read classification -> kraken-classification, metaphlan-profiling. Diversity/ordination/DA tool mechanics -> metagenome-visualization. Generic 16S diversity -> the microbiome category.

## The Single Most Important Modern Insight -- A Relative-Abundance Change Is Not a Change

A shotgun abundance table is a composition: the sequencer fixes the total number of reads, not the sample's microbial load. So every number is relative to every other, and "taxon X went up 2-fold" is undefined without a reference frame or an external load anchor. Without one, a single blooming taxon makes every other taxon look depleted (the blooming-taxon illusion), and a Pearson correlation of two taxa's proportions is biased negative by arithmetic, not biology. Bracken is step one of about six: classify -> re-estimate -> compositional transform -> normalize -> reference-frame test -> (optional) absolute conversion. Two corollaries:

1. **Bracken percent and MetaPhlAn percent are different physical quantities.** Bracken is a fraction of READS (large-genome biased, approximately fraction of DNA); MetaPhlAn is a genome-size-normalized marker abundance (approximately fraction of cells). They legitimately disagree 2-3x for the same sample. Picking one chooses an estimand (DNA vs cells), not the "more accurate tool."
2. **No library-size normalizer fixes the genome-size confound.** TSS/CSS/TMM/GMPR all operate on the count table and inherit the large-genome bias; only a coverage-based estimand (CoverM) or a genome-normalized profiler (MetaPhlAn) removes it.

## What Quantity Is Being Estimated?

| Estimand | Definition | Bias / use |
|----------|------------|-----------|
| Read count | raw reads to a taxon | library-size and genome-size dependent; never compare raw across samples |
| Relative abundance | read count / total | compositional (sums to 1); still genome-size biased |
| Coverage abundance | reads x readlen / genome size | removes large-genome bias; ~ fraction of genomes/cells (CoverM) |
| Cell fraction | fraction of organisms | needs coverage + an assumption of one genome per cell; polyploidy/growth bias it |
| Absolute load | composition x external total | the only basis for "increased/decreased" (flow/qPCR/spike-in) |

## Bracken: Step One, Not the Finish Line

Bracken redistributes reads stranded at the genus/family node (where shared k-mers stopped Kraken) down to species, using a database-derived expectation of where length-L reads classify. It does not classify and it does not add precision.

```bash
bracken-build -d "$KRAKEN_DB" -t 8 -k 35 -l 150   # one-time; -k MUST equal the Kraken2 build k (35)
bracken -d "$KRAKEN_DB" -i kraken.kreport -o out.bracken -w out.bracken.kreport \
    -r 150 \   # MUST equal the bracken-build -l AND the actual post-trim read length (not auto-detected)
    -l S -t 10 # species level; -t drops taxa with fewer than 10 clade-level reads (strict <) before redistribution
```

Output columns: `name, taxonomy_id, taxonomy_lvl, kraken_assigned_reads, added_reads, new_est_reads, fraction_total_reads`. `fraction_total_reads` is a fraction of classified-and-retained reads, not of all input - so two samples with different unclassified (e.g. host) fractions have non-comparable denominators. `combine_bracken_outputs.py --files *.bracken -o matrix.tsv` builds a taxa-by-sample matrix.

### Bracken's defining failure mode
Bracken can only redistribute among species already in the database. A true organism absent from the database has its reads parked by Kraken at the shared genus node, and Bracken hands them to the database-present congeners - confidently fabricating or inflating those species. Distrust any species with high `added_reads` but tiny `kraken_assigned_reads`; gate presence on upstream unique-minimizer evidence (kraken-classification) and a coverage breadth check.

## Treat the Table as a Composition

**Goal:** Make the abundance table valid for multivariate stats and correlation by removing closure with a centered log-ratio, after handling zeros (log of zero is undefined).

**Approach:** Impute count zeros with Bayesian-multiplicative replacement (preserves ratios), then CLR-transform; use Aitchison distance (Euclidean on CLR) downstream, never raw-proportion Pearson or Bray-Curtis for correlation.

```python
import pandas as pd
import numpy as np
from skbio.stats.composition import clr, multi_replace   # multiplicative_replacement was renamed multi_replace in skbio 0.6

counts = pd.read_csv('matrix.tsv', sep='\t', index_col=0)   # taxa x samples (new_est_reads)
mat = counts.T.values.astype(float)                          # samples x taxa for transform
mat_nozero = multi_replace(mat / mat.sum(axis=1, keepdims=True))
clr_mat = clr(mat_nozero)                                    # closure removed; rows are CLR coordinates
clr_df = pd.DataFrame(clr_mat, index=counts.columns, columns=counts.index)
```

For sparse tables prefer R `zCompositions::cmultRepl()` (Bayesian-multiplicative, posterior-imputed) over a fixed +1 pseudocount, which is arbitrary, distorts ratios, and re-opens closure. Structural zeros (taxon genuinely absent) and sampling zeros (present below detection) are usually indistinguishable from the table - disclose the assumption rather than pretend otherwise.

## Library-Size Normalization

| Method | What it does | Shotgun applicability |
|--------|--------------|-----------------------|
| TSS (proportions) | divide by library size | IS the compositional closure; fine for viz, biased for DA/correlation |
| Rarefaction | subsample to common depth | discards data; defensible for diversity, contested for DA (see below) |
| CSS (metagenomeSeq) | scale by a cumulative-sum quantile | robust to dominant taxa; designed for marker surveys, usable on counts |
| TMM (edgeR) | trimmed mean of M-values | "most features unchanged" often violated in microbiome; use with caution |
| RLE (DESeq median-of-ratios) | geometric-mean reference | breaks on zeros (geometric mean -> 0); needs poscounts workaround |
| GMPR | pairwise median ratios then geometric-mean | purpose-built for zero-inflated counts; good default size factor |

None of these fixes the genome-size confound - that needs a coverage estimand or a genome-normalized profiler.

### The rarefaction debate (do not take a side; decide per analysis)
McMurdie & Holmes 2014 (*PLoS Comput Biol* 10:e1003531) showed rarefying for DIFFERENTIAL ABUNDANCE is statistically wasteful versus modeling library size. Schloss 2024 (*mSphere* 9:e00354-23