bio-metagenomics-strain-tracking

What this skill does

## Version Compatibility

Reference examples tested with: inStrain 1.8+, StrainPhlAn/MetaPhlAn 4.1+, dRep 3.4+, skani 0.2+, Bowtie2 2.5+, samtools 1.19+, pandas 2.2+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `inStrain profile -h`, `strainphlan -h`, `skani dist -h` 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 mapping REFERENCE defines the answer. Map to genomes present in the sample set (dRep-dereplicated MAGs from this dataset), not generic database genomes - a distant reference inflates apparent SNVs and corrupts popANI. Record the reference set, the popANI/nGD threshold, the minimum coverage and breadth, and the co-detection rate; the "strain" is defined by these, not by nature.

# Strain Tracking

**"Is the same strain in samples A and B?"** -> Compare per-position SNV populations (not consensus genomes) at adequate depth - because a strain is defined by the chosen threshold, and ANI cannot resolve the difference that matters.
- CLI: `inStrain profile sample.bam reps.fasta -o sample.IS -s reps.stb -p 8` then `inStrain compare`

Scope: in-situ strain resolution, sharing, and deconvolution from a community. Pure-culture isolate outbreak SNP/cgMLST trees -> epidemiological-genomics. MAG assembly/binning -> genome-assembly/metagenome-assembly. Species presence/abundance -> kraken-classification, metaphlan-profiling.

## The Single Most Important Modern Insight -- A Strain Is a Threshold, Not a Thing

There is no universal definition of a metagenomic strain. A strain is an operational construct fixed by the reference mapped to, the genome fraction comparable at adequate depth, and the cutoff drawn. inStrain's popANI >= 99.999% over >= 50% of the genome IS the strain definition; Valles-Colomer's per-species nGD threshold IS the strain definition. Changing the cutoff changes how many strains exist. Two corollaries:

1. **ANI answers "same genome?"; popANI/nGD answer "same population, in situ?"** Genome-to-genome ANI (MASH/skani/fastANI) saturates - two genuinely distinct, separately transmissible strains routinely share >99.9% ANI, and ANI's resolution floor sits above the difference that matters. Strain sharing cannot be done with an ANI number; it needs microdiversity-aware, position-level concordance.
2. **Detecting a shared strain is a narrow statement:** across the genome fraction both samples covered at >= 5x, their SNV populations were >= 99.999% concordant. That is not proof an organism was transmitted between two people - the threshold certifies genomic identity within a coverage window, nothing more.

## The Three Tasks (Do Not Conflate)

| Task | Question | Tools |
|------|----------|-------|
| Identification | which known reference strain is present? | StrainGST, sourmash gather, MIDAS |
| Tracking / sharing | is the SAME strain in samples A and B? | inStrain compare, StrainPhlAn, MIDAS2, metaSNV, SameStr |
| Deconvolution | how many strains coexist in ONE sample, what are their haplotypes? | DESMAN, Strainberry, strainFlye, Strainy |

Genome-to-genome ANI (MASH/skani/fastANI) is a fourth, orthogonal task - "are these two assembled genomes the same?" - isolate/MAG comparison and dereplication, NOT in-situ strain resolution.

## Tool Taxonomy

| Tool | Citation | Mechanism / role | When |
|------|----------|------------------|------|
| inStrain | Olm 2021 *Nat Biotechnol* 39:727 | popANI/conANI microdiversity from reads mapped to MAGs | the reference standard for shared-strain detection |
| StrainPhlAn | Truong 2017 *Genome Res* 27:626 | marker-SNV consensus -> phylogeny -> nGD | large cross-sample marker surveys, no assembly needed |
| MIDAS2 | Zhao 2023 *Bioinformatics* 39:btac713 | UHGG pan-genome SNV + gene CNV | accessory-genome strain signal at scale |
| StrainGE | van Dijk 2022 *Genome Biol* 23:74 | k-mer search + low-coverage variant calling | low-abundance strains down to 0.5x coverage |
| metaSNV v2 | Van Rossum 2022 *Bioinformatics* 38:1162 | SNV distances + subspecies clustering | subspecies structure across samples |
| skani | Shaw 2023 *Nat Methods* 20:1661 | sparse-chaining ANI | genome-vs-genome ANI; robust on fragmented MAGs (prefer over fastANI) |
| Strainberry / strainFlye | Vicedomini 2021 *Nat Commun* 12:4485; Fedarko 2022 *Genome Res* 32:2119 | long-read haplotype separation | deconvolute co-occurring strains (-> genome-assembly) |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Is a strain shared between two metagenomes? | inStrain compare (popANI) | microdiversity-aware; the field standard |
| Cross-sample transmission survey, many samples | StrainPhlAn (per-species nGD) | marker-based, scalable, no assembly |
| Low-abundance pathogen (< 1% / < 5x) | StrainGE | detects/compares down to 0.5x |
| Accessory-genome / pan-genome strain signal | MIDAS2 | adds gene-content axis SNV tools miss |
| Separate co-occurring strains into haplotypes | DESMAN (many samples) or long-read Strainberry/strainFlye | SNV tools do not partition a mixture |
| Compare two assembled genomes / dereplicate | skani (or fastANI) | genome-vs-genome ANI, not in-situ strains |
| Pure-culture isolate outbreak tree | -> epidemiological-genomics | cgMLST/SNP-distance on one genome per sample |

## inStrain: popANI vs conANI

**Goal:** Decide whether two metagenomes share a strain without being fooled by which allele happens to be the majority.

**Approach:** dRep the dataset's MAGs into representative genomes, map reads to the concatenated references, profile each sample, then compare on popANI (microdiversity-aware) over the co-covered genome fraction.

```bash
# 1. dRep -> representative genomes (97-99% ANI); concatenate; build scaffold-to-bin (.stb).
# 2. Map reads to the concatenated reps - your OWN MAGs, not database genomes.
bowtie2 -x reps -1 r1.fq.gz -2 r2.fq.gz | samtools sort -o sampleA.bam
inStrain profile sampleA.bam reps.fasta -o sampleA.IS -s reps.stb -g genes.fna -p 8
inStrain profile sampleB.bam reps.fasta -o sampleB.IS -s reps.stb -g genes.fna -p 8
inStrain compare -i sampleA.IS sampleB.IS -o compare.out -s reps.stb -p 8
```

conANI calls a difference whenever the consensus base differs - confounded by within-sample microdiversity (a minor-allele flip fakes a difference). popANI calls a difference only if the two samples share NO alleles at all, including minor ones, so popANI >= conANI always and is what detects shared strains consensus tools miss. Read genome-level calls from `genomeWide_compare.tsv` (breadth column `percent_compared`); the per-scaffold `comparisonsTable.tsv` uses `percent_genome_compared`.

## StrainPhlAn: Marker SNVs and nGD

```bash
metaphlan sample.fq.gz --input_type fastq -s sample.sam.bz2 --bowtie2out sample.bz2 -o profile.tsv  # need the SAM (-s)
sample2markers.py -i sams/*.sam.bz2 -o consensus_markers -n 8
extract_markers.py -c t__SGB1877 -o clade_markers/
strainphlan -s consensus_markers/*.json -m clade_markers/t__SGB1877.fna \
    -r reference_genomes/*.fna.bz2 -o output -c t__SGB1877 \
    --marker_in_n_samples_perc 80 --sample_with_n_markers 20 --nproc 8  # 4.0 named this --marker_in_n_samples
```

The output tree gives a pairwise nGD (normalized genetic distance). There is no universal nGD strain cutoff - derive a per-species threshold from the data (same-individual-different-timepoint pairs fall below it, unrelated pairs above), as in Valles-Colomer 2023. Low coverage means too few markers pass the filters and the sample is dropped from the species tree silently - so a missing shared-strain call is not evidence of no shared strain.

## Genome-vs-Genome ANI (Isolate/MAG Comparison, NOT In-Situ Strains)

```bash
skani dist genomeA.fasta genomeB.fasta   # prefer skani over fastANI: r