bio-genome-assembly-scaffolding

What this skill does

## Version Compatibility

Reference examples tested with: YaHS 1.2+, SALSA2 2.3+, juicer_tools 1.22/2.0+, bwa 0.7.17+, chromap 0.2+, samtools 1.19+, RagTag 2.1+, tidk 0.2.3+, seqkit 2.6+.

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

hifiasm output filenames and YaHS/SALSA2 enzyme and resolution defaults have changed across releases; confirm output names and `-e`/`-r` behaviour against the installed version. The Hi-C-mapping flag set (`bwa mem -5SP` vs the Arima/VGP per-end pipeline) varies by vendor (Arima, Dovetail/Omni-C, Phase Genomics) - check the current mapping guide before pasting flags. If a command errors, introspect the tool and adapt rather than retrying.

# Genome Scaffolding

**"Turn my contigs into chromosomes"** -> Order and orient contigs by long-range linking signal, insert N-gap spacers, then curate the contact map - scaffolding adds no sequence, so the output is a draft chromosome structure, not finished sequence.
- CLI: `yahs contigs.fa hic_to_contigs.bam` (Hi-C default), `run_pipeline.py -a contigs.fa -b sorted.bed -m yes` (SALSA2), `ragtag.py scaffold ref.fa contigs.fa` (reference-guided), `allmaps path` (linkage maps)

## The Single Most Important Modern Insight -- Automated Scaffolding Produces a Draft; the Genome Is Trustworthy Only After Manual Curation of the Contact Map

YaHS finishes in minutes and emits a file literally named `_scaffolds_final.agp`. It is **not final**. Hi-C scaffolding is a contact-frequency *inference* - it orders/orients contigs by the polymer-physics signal that contact frequency decays with 1D distance - and it is error-prone exactly at the joins. Every reference-grade pipeline (VGP, Darwin Tree of Life, Earth BioGenome) treats the automated AGP as a *first pass* that a human curator then corrects in PretextView or Juicebox: breaking misjoins, flipping inversions, removing false duplications, assigning chromosomes by eye. Howe 2021 quantified the hidden labor across 111 VGP/DToL assemblies: on average **221 interventions per Gb (67 breaks, 105 joins, 49 false-duplication removals)**. Three load-bearing consequences:

1. **The contact map is the QC, not a decoration.** A correct chromosome shows one bright diagonal with smooth off-diagonal decay; off-diagonal blocks, anti-diagonal "bowties," and bleeding between chromosomes are errors to inspect and break (see Reading the Contact Map). An assembly paper whose methods say "scaffolded with [tool]" but never mention curation/Juicebox/PretextView is shipping a draft - downgrade any claim that depends on large-scale order (synteny, fusions, structural variants).

2. **Scaffold N50 is not contig N50, and conflating them is the classic reviewer catch.** Scaffolding inserts runs of **N** (estimated, often arbitrary length) and adds zero sequence. A genome can post "scaffold N50 = 60 Mb, chromosome-scale!" while its contig N50 is 200 kb - the contiguity is borne by gaps, not finished sequence. Always report contig N50, scaffold N50, gap count, and total N bases separately.

3. **Reordering cannot rescue a bad contig.** Scaffolders take contigs as given; a chimeric contig drags the wrong sequence into the wrong chromosome. Chimeras must be **broken before scaffolding** (see Chimera-Breaking).

## Linking-Data Modality Taxonomy

| Modality | How it links | Scale | Status |
|----------|--------------|-------|--------|
| Hi-C / Omni-C | proximity ligation; contact frequency ~ 1/(1D distance) | chromosome-scale | Dominant modern method. Omni-C is enzyme-free/sequence-agnostic; Arima uses fixed enzyme motifs |
| Bionano optical maps (DLS) | labeled motif patterns aligned to in-silico contig maps | megabase | Resolves large SVs and bridges complex repeats; separate instrument + library (cost decision) |
| 10x linked reads | barcoded short reads from one long molecule | ~10-100 kb | DISCONTINUED by 10x (2020). Legacy data only; Tigmint/ARCS still run, no new data |
| Genetic / linkage maps | marker recombination order = chromosome order | chromosome-scale, low resolution | ALLMAPS integrates multiple maps; orthogonal validation of Hi-C |
| Reference (homology) | align contigs to a related genome, copy its order | as good as synteny | RagTag. Fast, but imposes the reference's karyotype - see hazard below |
| Long reads as linkers | k-mer pairs spanning junctions | read length | LINKS; mostly superseded by assembling with the long reads directly |
| Mate-pair / jumping | large-insert paired short reads | 2-40 kb | Legacy (SSPACE); chimeric-insert artifacts; never recommend today |

## Decision Tree by Available Linking Data

| Available linking data | Use | Why |
|-------------|-----|-----|
| Hi-C/Omni-C, want fast contiguous default | YaHS | community default; fast, high N90, AGP + Juicebox outputs in one command |
| Hi-C, want graph-aware conservative joins + input-error correction | SALSA2 `-m yes` (`-g graph.gfa`) | iterative; breaks chimeric input; uses assembly graph to avoid orientation errors |
| Hi-C, interactive curation central to workflow | 3D-DNA + Juicer + Juicebox (JBAT) | the Aiden-lab `.assembly` <-> Juicebox round-trip is built around hand-editing |
| Hi-C, vertebrate/reference-grade | YaHS or SALSA2 -> PretextView/JBAT curation | VGP/DToL standard: automate then curate |
| Diploid, Hi-C, want both haplotypes as chromosomes | hifiasm `--h1/--h2` to PHASE first, then YaHS to SCAFFOLD each haplotype | same Hi-C, two jobs - phasing picks the homolog, scaffolding orders along it (see Phasing vs Scaffolding) |
| Closely related reference, karyotype NOT a question | RagTag `correct` then `scaffold` | homology ordering in minutes - but never if karyotype is the biology |
| One or more genetic/linkage maps | ALLMAPS | integrates multiple maps; robust to marker errors; validates/anchors Hi-C |
| Bionano optical maps + sequence assembly | Bionano Solve hybrid scaffold (before Hi-C) | megabase maps bridge complex repeats and large SVs |
| 10x linked reads (legacy data) | Tigmint (break) -> ARCS/ARKS (link) | break misassemblies first; no new 10x data exists |
| Contigs not yet QC'd / haplotigs present | -> assembly-qc, purge_dups | scaffolding haplotigs strings them up as fake chromosomes |
| Hi-C contact map for TADs/loops, not chromosomes | -> hi-c-analysis/matrix-operations | different use of the same assay |

**Modern reference-grade recipe (vertebrate/eukaryote):** HiFi -> hifiasm contigs -> (optional Bionano hybrid scaffold) -> Hi-C scaffold (YaHS/SALSA2) -> **manual curation in PretextView/JBAT** -> gap-fill (TGS-GapCloser) -> QC (tidk telomeres, contact-map diagonal, contig-vs-scaffold N50).

## Hi-C Read Mapping (the step people botch)

Hi-C reads are **not** a normal paired-end library: the two ends come from *different* genomic loci ligated together, so standard PE proper-pair/insert-size logic mis-handles them.

```bash
# Map each end SEPARATELY with no mate rescue/pairing; -5 reports the 5' (junction) portion as primary.
bwa index contigs.fa
bwa mem -5SP -T0 -t 16 contigs.fa hic_R1.fq.gz hic_R2.fq.gz | \
    samtools view -@ 8 -b - > aligned.bam
# Mandatory: mark/remove PCR + optical duplicates (Hi-C is duplicate-rich) before scaffolding.
samtools sort -@ 8 -n aligned.bam | samtools fixmate -m - - | \
    samtools sort -@ 8 - | samtools markdup -@ 8 - hic_to_contigs.bam
samtools index hic_to_contigs.bam
```

- **Each end separately, no mate rescue** (`-5SP`): `-5` = 5' portion of a chimeric junction read as primary; `-S`/`-P` skip mate rescue and pairing. The Arima/VGP per-end pipeline (`bwa mem` on R1 and R2 independently, `filter_five_end.pl`, `two_read_bam_combiner.pl`) is the higher-fidelity equivalent.
- **MAPQ filter** at scaffolding (YaHS `-q`, SALSA2 reads MAPQ) to drop repeat-ambiguous reads.
- **Enzyme awareness:** Arima/Dovetail-DpnII cut at fixed motifs (`GATC`); the scaffolder uses `-e` to model legitimate read starts. **Omni-C / DNa