Bio Genome Assembly Contamination Detection

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

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

FCS-adaptor (Run First - Cheap and Unambiguous)

# 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)

# 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 - manual inspection required). A separate INFO action specifically marks sequence known to be integrated into host genomes (e.g. endosymbiont insertions). clean genome acts on EXCLUDE/TRIM/FIX only; it does not touch REVIEW or INFO. Those non-auto-cleaned tiers exist precisely because FCS-GX can flag legitimate HGT/endosymbiont sequence as contaminant - never bulk-delete every flagged contig without reading the report.

The RAM wall (an operational, cloud-forcing fact). The GX database is ~470 GiB on disk and the documented sweet spot is a ~512 GiB-RAM host. Underprovisioned, the run does not fail - it crawls (minutes become days when the DB spills out of RAM). Best practice copies the DB into a tmpfs RAM disk; most labs rent a high-memory cloud VM for the screen. The first question before advising FCS-GX is "is there ~1/2 TB RAM or a cloud budget?".

BlobToolKit (The Blob Plot - Visual Triage)

blobtools create --fasta assembly.fasta ./BlobDir
blobtools add --hits diamond.out --taxrule bestsumorder --taxdump /path/taxdump \
    --cov mapping.bam ./BlobDir
blobtools view --remote ./BlobDir   # interactive viewer; GC(x) vs coverage(y), sized by length, colored by taxonomy

Each contig is plotted by GC fraction (x) vs coverage (y), sized by length, colored by best-hit taxonomy. The target organism forms one tight GC x coverage cloud; a free-living contaminant grows independently and forms its own cloud at a different GC and a different coverage. Inputs: a hit file (BLAST/DIAMOND vs nt/UniProt), a coverage BAM, and an NCBI taxdump. blobtools filter can extract or drop by taxon - but see the failure mode below before doing so.

Merge CheckM2 and GUNC

Goal: Produce the joint CheckM2 x GUNC table that is the field standard for MAG QC, so a chimera invisible to CheckM2 is not reported as clean.

Approach: Load both reports, join on genome name, and apply MIMAG thresholds AND the GUNC pass flag together; never gate on contamination % alone.

import pandas as pd

CONTAM_HQ = 5      # MIMAG high-quality: contamination < 5% (Bowers 2017); above this, gene-content inference is unreliable
COMPLETE_HQ = 90   # MIMAG high-quality: completeness > 90%
RRS_TRUST = 0.5    # GUNC pass is trustworthy only when reference_representation_score > 0.5; below it 'pass' means 'can't tell'

checkm = pd.read_csv('checkm2_out/quality_report.tsv', sep='\t')
gunc = pd.read_csv('gunc_out/GUNC.progenomes_2.1.maxCSS_level.tsv', sep='\t')
merged = checkm.merge(gunc, left_on='Name', right_on='genome', how='left')

merged['gunc_trustworthy'] = merged['reference_representation_score'] > RRS_TRUST
merged['high_quality'] = ((merged['Completeness'] > COMPLETE_HQ) &
                          (merged['Contamination'] < CONTAM_HQ) &
                          (merged['pass.GUNC'] == True) &
                          merged['gunc_trustworthy'])
merged.to_csv('combined_qc.tsv', sep='\t', index=False)

CheckM2 also runs as checkm2 predict -i bins/ -o checkm2_out --threads 16 -x fa (downloads its own DIAMOND DB), and GUNC as gunc run -d bins/ -o gunc_out -t 16 -e .fa against a downloaded reference (gunc download_db). GTDB-Tk taxonomy is gtdbtk classify_wf --genome_dir bins/ --out_dir gtdbtk_out -x fa --cpus 16.

GUNC and the RRS Trap

pass.GUNC = TRUE at **CSS <= 0.45** is not a clean bill of health when reference_representation_score (RRS) is low. Low RRS means the genome is barely represented in GUNC's reference DB - there is not enough signal to *detect* chimerism, so "pass" degenerates into "can't tell, not clean". For genuinely novel lineages (microbial dark matter), a GUNC pass is weak evidence; trust it only when RRS > 0.5, and otherwise lean on manual contig inspection or MDMcleaner. Read the configurations: CheckM2-high-contam + GUNC-pass = duplicative contamination within one lineage; CheckM2-low-contam + GUNC-fail = the dangerous disjoint-marker chimera; both pass with high RRS = the only "clean" verdict.

The Tardigrade Lesson (HGT vs Contamination)

Boothby 2015 (PNAS 112:15976) reported ~17% of tardigrade genes arrived by horizontal gene transfer; Koutsovoulos 2016 (PNAS 113:5053) showed almost all of it was undetected bacterial contamination in the assembly, dropping real HGT to ~1-2%. The original paper was corrected, not retracted - it stands as a cautionary monument. A foreign-looking gene has two explanations the sequence alone cannot distinguish: a contaminant contig, or a real HGT/endosymbiont gene integrated into the host genome. The tell is physical integration, not taxonomy. A real HGT gene sits on a host contig (host-gene flanks, host GC, host coverage, host introns); a contaminant sits on its own contig (foreign GC, its own coverage cloud, prokaryotic gene structure). Deleting everything that looks bacterial manufactures the opposite of Boothby's error - erasing real biology. Require integration evidence (host flanks, host-typical GC/coverage/introns, long reads spanning the host->foreign junction) before keeping or stripping.