# Bio Applied Assembly Sv

> Assemble ONT/HiFi reads with Flye/Hifiasm, polish with Medaka, QC with QUAST/BUSCO, call SVs (DEL/INS/INV/DUP/BND) with Sniffles2. Use for long-read assembly, N50/BUSCO QC, or nanopore/HiFi SV calling to VCF.

- Skill: `pavel-kravchenko/bio-applied-assembly-sv` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-assembly-sv`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-assembly-sv/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-assembly-sv

---


# Long-Read Genome Assembly and Structural Variant Calling

## When to Use

- Assembling a bacterial, eukaryotic, or metagenomic genome de novo from ONT or PacBio HiFi long reads
- Choosing between Flye (versatile: ONT/HiFi/CLR/metagenome) and Hifiasm (haplotype-resolved HiFi, optional trio binning)
- Polishing an ONT draft assembly with Medaka and grading it with QUAST (contiguity/reference accuracy) and BUSCO (gene completeness)
- Calling structural variants — deletions, insertions, inversions, duplications, breakends — from an aligned long-read BAM with Sniffles2
- Jointly genotyping SVs across a cohort from per-sample `.snf` files without re-running alignment

## Version Compatibility

Flye ≥2.9, Hifiasm ≥0.19, Medaka ≥1.11, Sniffles2 ≥2.2, QUAST ≥5.2, BUSCO ≥5.7 (odb10 lineages), bcftools ≥1.19, Python ≥3.10 with pysam ≥0.22, pandas ≥2.0, numpy ≥1.26.

## Prerequisites

- Conda/mamba env with `flye`, `hifiasm`, `medaka`, `sniffles` (Sniffles2), `quast`, `busco`, `bcftools`, `samtools` on PATH
- Reads already basecalled and QC'd (`bio-long-read-sequencing-basecalling`, `bio-long-read-sequencing-long-read-qc`)
- For SV calling: reads aligned to a reference with minimap2 (`bio-long-read-sequencing-long-read-alignment`), sorted/indexed BAM
- Python: `pysam`, `pandas` for parsing `assembly_info.txt` and VCF output

## De Novo Assembly (Flye / Hifiasm)

**Goal:** produce contigs from raw long reads.
**Approach:** pick the assembler by read type — Flye for ONT (any mode) or when metagenomics is needed; Hifiasm for PacBio HiFi when haplotype phasing matters. `--genome-size` in Flye only needs to be within ±50% of truth.

```bash
# Flye: ONT R10.4.1 (Q20+) bacterial genome
flye --nano-hq reads.fastq.gz --genome-size 5m --out-dir flye_out/ --threads 8
# Flye modes: --nano-raw (R9.4.1, Q<15) | --nano-corr (pre-corrected ONT)
#             --pacbio-hifi (CCS/HiFi) | --pacbio-raw (CLR)
grep -c ">" flye_out/assembly.fasta          # contig count
cat flye_out/assembly_info.txt               # per-contig length/cov/circularity

# Hifiasm: PacBio HiFi, haplotype-resolved by default
hifiasm -o sample.asm -t 8 hifi_reads.fastq.gz
awk '/^S/{print ">"$2"\n"$3}' sample.asm.bp.p_ctg.gfa > assembly_primary.fasta

# Trio binning (parental Illumina reads -> fully phased hap1/hap2)
yak count -b37 -t 8 -o pat.yak paternal.fastq.gz
yak count -b37 -t 8 -o mat.yak maternal.fastq.gz
hifiasm -o sample.trio -t 8 -1 pat.yak -2 mat.yak hifi_reads.fastq.gz
```

```python
import numpy as np
import pandas as pd


def assembly_n_stats(lengths, genome_size=None, thresholds=(50, 90)):
    """Compute N-stats (and NG-stats if genome_size given) from contig lengths.

    lengths: iterable of contig lengths (bp)
    genome_size: expected genome size in bp, enables NGxx metrics
    """
    lengths = np.sort(np.asarray(lengths, dtype=int))[::-1]
    total = lengths.sum()
    reference = genome_size if genome_size else total
    cumsum = np.cumsum(lengths)
    out = {"n_contigs": len(lengths), "total_bp": int(total), "largest_bp": int(lengths[0])}
    for t in thresholds:
        idx = min(np.searchsorted(cumsum, reference * t / 100), len(lengths) - 1)
        out[f"N{t}"] = int(lengths[idx])
        if genome_size:
            out[f"NG{t}"] = int(lengths[idx])
    return out


# Parse Flye's assembly_info.txt (tab-separated, header row)
info = pd.read_csv("flye_out/assembly_info.txt", sep="\t")
stats = assembly_n_stats(info["length"], genome_size=5_000_000)
print(stats)
```

## Polishing and Quality Assessment

**Goal:** correct systematic ONT homopolymer-indel errors, then verify contiguity and gene completeness.
**Approach:** Medaka model MUST match flow cell chemistry (encoded in the model name); skip Medaka for HiFi assemblies (already Q20+). Follow with QUAST (needs a reference for full metrics) and BUSCO (reference-free, gene-space).

```bash
medaka tools list_models                              # find your chemistry's model
medaka_consensus -i reads.fastq.gz -d flye_out/assembly.fasta \
    -o medaka_out/ -t 8 -m r1041_e82_400bps_hac_v4.2.0  # must match flow cell
# Output: medaka_out/consensus.fasta

quast.py flye_out/assembly.fasta medaka_out/consensus.fasta \
    --reference reference.fasta --output-dir quast_report/ --threads 8

busco -i medaka_out/consensus.fasta -l bacteria_odb10 \
    -o busco_out/ -m genome --cpu 8    # eukaryota_odb10 / mammalia_odb10 for eukaryotes
```

A good bacterial assembly: BUSCO Complete >95%, Duplicated <2%; QUAST # misassemblies ≈0. A good eukaryotic assembly: Complete >90%, Duplicated <5%.

## Structural Variant Calling (Sniffles2)

**Goal:** call DEL/INS/INV/DUP/BND ≥50 bp directly from long reads, which span most SVs (unlike short-read discordant-pair/split-read inference).
**Approach:** call per-sample with `--snf` saved for joint genotyping; filter on `FILTER=PASS` and `SUPPORT`; parse the VCF with pysam rather than shelling out to bcftools for downstream analysis.

```bash
sniffles --input aligned_sorted.bam --vcf svs.vcf --reference ref.fa \
    --threads 8 --snf sample.snf        # .snf enables joint genotyping later

# Multi-sample joint genotyping (no re-alignment needed)
sniffles --input sample1.snf sample2.snf sample3.snf \
    --vcf joint_svs.vcf --reference ref.fa

bcftools view -i 'FILTER="PASS" && INFO/SUPPORT>=5' svs.vcf > svs_filtered.vcf
```

```python
import pysam
import pandas as pd


def parse_sv_vcf(vcf_path, min_support=5):
    """Load a Sniffles2 VCF into a tidy DataFrame of PASS calls.

    Returns columns: chrom, pos, svtype, svlen, af, support.
    """
    rows = []
    with pysam.VariantFile(vcf_path) as vf:
        for rec in vf:
            info = rec.info
            support = info.get("SUPPORT", 0)
            if rec.filter.keys() and "PASS" not in rec.filter.keys():
                continue
            if support < min_support:
                continue
            rows.append({
                "chrom": rec.chrom,
                "pos": rec.pos,
                "svtype": info.get("SVTYPE"),
                "svlen": abs(info.get("SVLEN", 0)),
                "af": info.get("AF"),
                "support": support,
            })
    return pd.DataFrame(rows)


hits = parse_sv_vcf("svs_filtered.vcf")
print(hits["svtype"].value_counts())
print(hits.loc[hits.svtype == "DEL", "svlen"].median())          # median DEL size
alu_ins = hits[(hits.svtype == "INS") & hits.svlen.between(250, 320)]  # Alu-sized INS
```

## Pitfalls

- **Wrong Flye input mode** (`--nano-hq` vs `--nano-raw` vs `--pacbio-hifi`) degrades assembly quality — match it to the actual basecaller/chemistry
- **Medaka model/chemistry mismatch** (e.g. using an R9.4.1 model on R10.4.1 data) silently produces worse polishing, not an error
- **Polishing HiFi assemblies** is usually unnecessary and can waste compute — HiFi consensus is already Q20+
- **GFA vs FASTA**: Hifiasm outputs GFA, not FASTA — convert primary contigs with `awk '/^S/{print ">"$2"\n"$3}'`
- **Coordinate systems**: BED is 0-based half-open; VCF/GFF is 1-based inclusive — mixing them causes off-by-one errors in SV coordinates
- **SNF vs VCF merging**: joint-genotyping `.snf` files in Sniffles2 is far more efficient than VCF merging (e.g. SURVIVOR) for cohort-scale SV analysis
- **Multiple testing**: apply FDR correction when comparing SV burden across many samples/genes

## See Also

- `bio-long-read-sequencing-structural-variants`
- `bio-long-read-sequencing-long-read-alignment`
- `bio-genome-assembly-hifi-assembly`
- `bio-genome-assembly-assembly-polishing`