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.

# 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

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

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.

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

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