bio-genome-assembly-genome-profiling

What this skill does

## Version Compatibility

Reference examples tested with: GenomeScope2 2.0+, KMC 3.2+, Jellyfish 2.3+, meryl 1.4+ (Merqury 1.3+), Smudgeplot 0.2.5+, KAT 2.4+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

Smudgeplot changed its backend: classic releases (0.2.x) run `smudgeplot.py hetkmers`/`smudgeplot.py plot` on a KMC dump; newer releases (0.3+) run `smudgeplot hetmers`/`smudgeplot all` on a FastK database. Confirm which interface is installed (`smudgeplot.py --version` or `smudgeplot --version`) before scripting. GenomeScope2 ships as `genomescope2` and as `genomescope.R`; both take the same flags. meryl `best_k.sh` lives in the Merqury install. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Genome Profiling

**"What am I about to assemble, and what should I expect?"** -> Build a k-mer spectrum from raw accurate reads and model it to estimate genome size, heterozygosity, repeat content, and ploidy, which set every downstream assembly expectation and parameter.
- CLI: `kmc -k21 ... reads kmc_db tmp/ && kmc_tools transform kmc_db histogram reads.histo` (count) then `genomescope2 -i reads.histo -o gs_out -k 21 -p 2` (model); `smudgeplot.py hetkmers` / `smudgeplot.py plot` (ploidy)

## The Single Most Important Modern Insight -- Profile the Genome Before Assembling, or the Assembler Guesses For Itself

A k-mer spectrum built from the raw reads -- reference-free, before a single contig exists -- estimates genome size, heterozygosity, repeat content, and ploidy, and those four numbers set the rest of the project: the NG50 denominator (assembly-qc), Flye's `-g`/`--genome-size`, hifiasm's `--hom-cov`/purge level, and whether short reads can produce the assembly being asked for at all. Skipping it leaves the assembler to infer the homozygous-coverage peak itself; when it mis-estimates (heterozygosity, odd ploidy, contamination, a bimodal coverage spectrum) it over- or under-purges, and that is exactly why people publish genomes inflated 1.5-2x by uncollapsed haplotigs. Three load-bearing moves:

1. **The estimate is the denominator AND the sanity check.** NG50 is N50 against the *expected* genome size, not the assembly size, so without the estimate NG50 cannot be reported honestly. After assembly, the same number is the haplotig test: an assembly 1.5-2x the GenomeScope size with high BUSCO-Duplicated is uncollapsed haplotypes, not a big genome -- purge before believing the size (see hifi-assembly, assembly-qc).
2. **The spectrum reads as a diagnostic, not just a size estimate.** A single homozygous peak ~= haploid/inbred; a distinct half-coverage (AB) peak *left* of the homozygous (AA) peak is heterozygous diploid, and the het peak's area gives the heterozygosity rate. A heavy high-multiplicity tail is repeat content; a spike at very high multiplicity is organelle or high-copy repeat; a left-shoulder near multiplicity 1 is sequencing error or a low-coverage contaminant. Read it before trusting any number from it.
3. **Count with accurate reads only.** GenomeScope's negative-binomial mixture model assumes errors are rare and Poisson-like. Noisy ONT (raw/HAC) injects so many unique error k-mers that the error shoulder swamps the real peaks and the fit fails. Count from Illumina or PacBio HiFi; use the ONT reads to *assemble*, never to *profile*.

## Tool Taxonomy

| Tool | Citation | Role | When |
|------|----------|------|------|
| KMC | Kokot 2017 *Bioinformatics* | disk-based k-mer counter -> histogram | default counter; frugal on RAM, fast on large genomes |
| Jellyfish | Marcais & Kingsford 2011 *Bioinformatics* | in-memory k-mer counter -> histogram | alternative counter; classic GenomeScope input |
| meryl | Rhie 2020 *Genome Biol* | k-mer counter + `best_k.sh` | derives k from genome size; feeds Merqury QV downstream |
| GenomeScope2 | Ranallo-Benavidez 2020 *Nat Commun* | model: size, het, repeat, ploidy from the histogram | the profiling model for diploids and polyploids |
| Smudgeplot | Ranallo-Benavidez 2020 *Nat Commun* | ploidy from het k-mer-pair coverage ratios | unknown ploidy; cross-check GenomeScope's `-p` |
| KAT | Mapleson 2017 *Bioinformatics* | spectra plots, reads-vs-assembly k-mer comparison | contamination triage; post-assembly completeness/spectra-cn |

## Decision Tree by Scenario

| Scenario | What the profile indicates | Path |
|----------|---------------------------|------|
| Single sharp peak, size ~ expected, low het | haploid/inbred or clonal; clean | proceed -> short-read-assembly or hifi-assembly with default purge |
| Two peaks (AB at ~half AA) | heterozygous diploid; het rate from AB area | HiFi+phasing best; if short reads only, expect haplotigs -> short-read-assembly (Platanus) |
| High het + Illumina only | short-read DBG will fragment and inflate 1.5-2x | get long reads, or plan purge_dups; do not report inflated size |
| GenomeScope `-p 2` fits poorly; Smudgeplot shows AAB/AABB | triploid/tetraploid; ploidy not 2 | re-run GenomeScope with correct `-p`; -> hifi-assembly haplotype expectations |
| Coverage peak < ~15-20x | too shallow for a stable model fit | sequence more, or treat size/het as lower-confidence; -> read-qc/quality-reports |
| Extra peak at odd multiplicity, or bimodal spectrum | contamination / organelle / mixed sample | KAT spectra triage; screen reads -> read-qc/quality-reports before assembling |
| Only noisy ONT available | cannot profile reliably from error-dominated spectrum | assemble first, then estimate size from the assembly + Merqury -> assembly-qc |
| Need the genome-size denominator for QC | GenomeScope haploid length | feed as NG50 expected size and Flye `-g` -> assembly-qc, long-read-assembly |

## Choosing k

**Goal:** Pick a k that is large enough that most k-mers are genomically unique but small enough to keep per-k-mer depth high.

**Approach:** Derive k from the expected genome size with Merqury's `best_k.sh` (formula `k = log4(G(1-p)/p)`, default tolerable collision rate `p=0.001`); a vertebrate-scale ~3 Gb genome returns k=21 (the de facto GenomeScope default), ~1 Gb returns k=20, and a ~12 Mb yeast returns k=17.

```bash
sh $MERQURY/best_k.sh 3100000000          # ~3.1 Gb vertebrate -> k=21 (the common GenomeScope default)
sh $MERQURY/best_k.sh 1000000000          # ~1 Gb -> k=20
sh $MERQURY/best_k.sh 12000000            # ~12 Mb yeast -> k=17
```

Too small a k saturates: nearly every k-mer recurs across the genome by chance, the unique/repeat peaks merge, and size is overestimated. Too large a k loses depth (per-k-mer coverage is `c*(L-k+1)/L`, so it drops as k rises) and the peaks blur into the error shoulder. k=21 is the long-standing default at vertebrate (~3 Gb) scale because it sits in this window; smaller genomes want smaller k (best_k.sh returns ~20 at 1 Gb, ~17 at 12 Mb). Use the SAME k for counting and for the `-k` passed to GenomeScope2.

## Counting K-mers and Running GenomeScope2

```bash
# KMC: -ci1 keeps singletons (the error shoulder GenomeScope models), -cs10000 caps the histogram tail
kmc -k21 -t16 -m64 -ci1 -cs10000 @fastq_list.txt kmc_db tmp/
kmc_tools transform kmc_db histogram reads.histo -cx10000

# GenomeScope2: -p 2 = diploid; raise for known/Smudgeplot-suggested polyploidy
genomescope2 -i reads.histo -o gs_out -k 21 -p 2

# Jellyfish alternative (-C canonical k-mers, mandatory for unstranded WGS)
jellyfish count -C -m 21 -s 4G -t 16 reads_*.fastq -o reads.jf
jellyfish histo -t 16 reads.jf > reads_jf.histo
```

GenomeScope2 writes a model fit (`model.txt`), the linear/log spectrum plots, and a summary with haploid genome length, heterozygosity, and a "% unique" (inverse of repeat content). `-cs10000`/`-cx10000` cap the histogram so a single organelle/repeat spike at multiplicity 100000+ does not dominate the file; raise the cap only for very high-coverage data.

## K-mer Coverage vs Sequencing Coverage (the co