bio-genome-assembly-contamination-detection

What this skill does

## Version Compatibility

Reference examples tested with: FCS-GX 0.5+, FCS-adaptor 0.5+, BlobToolKit 4.3+, CheckM2 1.0+, GUNC 1.0+, GTDB-Tk 2.4+, pandas 2.2+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures

Database release drives results more than the binary version. Record: the **FCS-GX GX database** release (it grows each release; ~470 GiB and rising), the **CheckM2 DIAMOND DB** version, the **GUNC reference** (proGenomes 2.1 default vs GTDB), and the **GTDB-Tk reference package release** (e.g. R220) - GTDB-Tk fails loudly if the DB release does not match the binary. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Contamination Detection

**"Is my assembly contaminated?"** -> First decide which of two unrelated questions is being asked - "which contigs are not my organism?" (foreign-sequence screen) or "how complete/clean/chimeric is this bin?" (MAG quality) - then run the matching toolset; crossing them returns confidently wrong numbers.
- CLI (foreign): `run_fcsadaptor.sh` then `fcs.py screen genome --fasta asm.fa --gx-db <db> --tax-id <taxid>`, `blobtools create/add/view`
- CLI (MAG): `checkm2 predict -i bins/ -o out/` AND `gunc run -d bins/ -o out/`, `gtdbtk classify_wf`

## The Single Most Important Modern Insight -- "Contamination" Is Two Disjoint Problems With Disjoint Toolsets

There is no single "contamination" measurement. There are two questions that share a word and almost nothing else, and applying one toolset to the other problem runs to completion and prints a plausible, meaningless number.

1. **Foreign sequence in a single-organism assembly** (a eukaryotic nuclear genome or a cultured isolate with bacterial/human/vector contigs in it). The question is *which contigs are not my target organism?* Answer = a **set of contigs to remove, trim, or split out**. Tools: **NCBI FCS-GX** + **FCS-adaptor** (now GenBank-submission-mandatory), **BlobToolKit** (the GC x coverage x taxonomy blob plot). There is no "contamination %".

2. **MAG/bin quality** (a bin recovered from a metagenome). The question is *how complete is this bin and how much is from other organisms?* Answer = a **percentage pair** (completeness %, contamination %) plus a chimerism verdict, judged against **MIMAG**. Tools: **CheckM2** + **GUNC** together, **GTDB-Tk** for taxonomy.

**The cardinal category error:** running CheckM2/CheckM on a eukaryotic nuclear assembly (its bacterial/archaeal marker sets are meaningless for a eukaryote - eukaryote completeness is BUSCO; see assembly-qc), or running FCS-GX on a single MAG and reading the flagged-length fraction as "the MIMAG contamination %". Both wrong applications complete silently and emit a number. The only defense is knowing which question each tool answers.

**The second load-bearing fact - CheckM2 and GUNC are blind to different things, so report the pair, never the % alone.** CheckM2 estimates contamination from **marker-gene redundancy** (a single-copy marker seen twice = contamination signal). A **chimera of two organisms with disjoint marker complements** - organism A contributed markers 1-60, organism B markers 61-120 - shows *no* duplicated markers, so CheckM2 reports low contamination while the bin is biological nonsense. GUNC catches exactly this, because it scores whether taxonomic signal is *consistent across contigs* (clade separation score), not marker counts. Neither is a superset of the other.

## Tool Taxonomy

| Tool | Citation | Role | Problem |
|------|----------|------|---------|
| FCS-adaptor | Astashyn 2024 *Genome Biol* | adaptor/vector screen (VecScreen successor); tiny DB, trivial RAM | foreign (run FIRST) |
| FCS-GX | Astashyn 2024 *Genome Biol* | cross-taxon foreign-sequence screen vs a declared tax-id; GenBank-mandatory | foreign |
| BlobToolKit | Challis 2020 *G3* | GC x coverage x taxonomy blob plot for visual triage; the maintained successor to BlobTools | foreign |
| tiara | Karlicki 2022 *Bioinformatics* | deep-learning euk/prok/organelle/plastid/mito contig classifier (alignment-free) | foreign / organelle partitioning |
| CheckM2 | Chklovski 2023 *Nat Methods* | ML completeness + contamination (lineage-agnostic; handles novel/reduced lineages) | MAG (default) |
| CheckM (legacy) | Parks 2015 *Genome Res* | lineage-specific collocated markers on a placement tree; reports Strain heterogeneity | MAG (legacy; slow, ~40 GB RAM) |
| GUNC | Orakov 2021 *Genome Biol* | taxonomic-chimerism detection via clade separation score (CSS) | MAG (run WITH CheckM2) |
| GTDB-Tk | Chaumeil 2020 *Bioinformatics* | GTDB taxonomic classification of bacterial/archaeal genomes | MAG taxonomy |

BlobTools (Laetsch & Blaxter 2017 *F1000Res*) and CheckM are the legacy predecessors; BlobToolKit and CheckM2 are the maintained defaults. CAT/BAT (contig taxonomy) and MAGpurify (reference-free MAG cleaning) exist but are secondary; RefineM is deprecated by its author (use GUNC + manual curation, or MDMcleaner for SAGs/dark-matter lineages).

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Eukaryotic nuclear assembly headed for GenBank | FCS-adaptor -> FCS-GX (mandatory) -> BlobToolKit for triage | submission bounces without a clean FCS screen |
| Cultured bacterial/archaeal isolate | FCS-GX + FCS-adaptor; CheckM2 as a completeness sanity check | an isolate should be ~100% complete, low contam |
| MAG(s) binned from a metagenome | CheckM2 AND GUNC together; score vs MIMAG; GTDB-Tk for taxonomy | the % pair plus the orthogonal chimerism verdict |
| Unbinned metagenome | bin first (-> metagenome-assembly), THEN CheckM2 + GUNC per bin | CheckM2/GUNC consume bins, do not make them |
| Eukaryotic contigs mixed into a metagenome | tiara / Whokaryote / EukRep to partition euk from prok | CheckM2 markers are prokaryote-only |
| Organellar (mito/chloroplast) contigs in a nuclear assembly | tiara organelle classes + coverage spike -> separate, do NOT delete | organelles are real biology; submit as own record |
| Foreign-looking gene embedded in a host scaffold | investigate integration (-> comparative-genomics) before removing | could be real HGT/endosymbiont, not contamination |
| Read-level host removal before assembly | -> read-qc/contamination-screening, metagenomics/kraken-classification | this skill screens assembled sequence, not reads |
| ANI/species placement of the cleaned genome | -> comparative-genomics/genome-distance-and-species-delineation | taxonomy after decontamination |

## FCS-adaptor (Run First - Cheap and Unambiguous)

```bash
# Adaptor/vector screen; tiny DB, trivial RAM. --euk or --prok for the lineage.
run_fcsadaptor.sh --fasta-input assembly.fa.gz --output-dir ./adaptor_out --euk
```

Adaptor and vector hits are unambiguous - always trim or exclude. Running adaptor first removes the obvious junk and shrinks the input before the expensive GX pass.

## FCS-GX (The GenBank-Mandatory Foreign Screen)

```bash
# Screen against the GX database; --tax-id is the NCBI taxid of the SOURCE organism.
python3 ./fcs.py screen genome --fasta assembly.fa.gz --out-dir ./gx_out/ \
    --gx-db "$GXDB_LOC/gxdb" --tax-id 9606

# Apply ONLY the auto-clean actions (EXCLUDE/TRIM/FIX) to produce a cleaned FASTA.
zcat assembly.fa.gz | python3 ./fcs.py clean genome \
    --action-report ./gx_out/assembly.fa.9606.fcs_gx_report.txt \
    --output clean.fasta --contam-fasta-out contam.fasta
```

The action report (`<name>.<taxid>.fcs_gx_report.txt`, where `<name>` keeps the input name minus only `.gz` - e.g. `assembly.fa.gz` -> `assembly.fa.9606.fcs_gx_report.txt`) assigns each flagged region an action: **EXCLUDE** (drop the whole sequence), **TRIM** (remove a contaminated end), **FIX** (hard-mask an internal span), or **REVIEW** (flagged but NOT auto-cleaned - manu