Genome Scaffolding
Hi-C Data Preprocessing
# Align Hi-C reads to draft assembly
bwa index draft_assembly.fa
bwa mem -5SP -t 16 draft_assembly.fa hic_R1.fq.gz hic_R2.fq.gz | \
samtools view -@ 8 -bhS - > aligned.bam
# Filter for Hi-C contacts (pairtools)
pairtools parse --min-mapq 40 --walks-policy 5unique --max-inter-align-gap 30 \
--nproc-in 8 --nproc-out 8 --chroms-path draft_assembly.fa.fai aligned.bam | \
pairtools sort --nproc 8 | \
pairtools dedup --nproc 8 --mark-dups | \
pairtools split --output-pairs contacts.pairs.gz
YaHS Scaffolding (Recommended)
# Index assembly
samtools faidx draft_assembly.fa
# Convert BAM to BED
bedtools bamtobed -i aligned.bam | sort -k4 > aligned.bed
# Run YaHS
yahs draft_assembly.fa aligned.bed -o scaffolds
# Output files:
# scaffolds_scaffolds_final.fa - final scaffolds
# scaffolds_scaffolds_final.agp - AGP file
# scaffolds.bin - contact matrix
YaHS with Error Correction
# Run with error correction
yahs draft_assembly.fa aligned.bed -o scaffolds --no-contig-ec
# Generate contact map for juicebox
juicer pre scaffolds.bin scaffolds_scaffolds_final.agp draft_assembly.fa.fai | \
sort -k2,2d -k6,6d -T ./ --parallel=8 -S 50G | \
awk 'NF' > scaffolds.pre.txt
# Create .hic file
java -Xmx48G -jar juicer_tools.jar pre scaffolds.pre.txt scaffolds.hic \
<(cut -f1,2 scaffolds_scaffolds_final.fa.fai)
3D-DNA Pipeline
# Prepare input (requires Juicer aligned data)
# Run Juicer first to get merged_nodups.txt
# Run 3D-DNA
run-asm-pipeline.sh -r 2 draft_assembly.fa merged_nodups.txt
# Output: draft_assembly.final.fasta
# Generate review assembly for Juicebox
run-asm-pipeline-post-review.sh -r draft_assembly.final.review.assembly \
draft_assembly.final.fasta merged_nodups.txt
SALSA2 Scaffolding
# Run SALSA2
python run_pipeline.py -a draft_assembly.fa -l draft_assembly.fa.fai \
-b aligned.bed -e GATC -o salsa_output -m yes
# With multiple restriction enzymes
python run_pipeline.py -a draft_assembly.fa -l draft_assembly.fa.fai \
-b aligned.bed -e GATC,GANTC -o salsa_output -m yes -p yes
Generate Contact Map
# Using cooler
cooler cload pairs -c1 2 -p1 3 -c2 4 -p2 5 \
draft_assembly.fa.fai:10000 contacts.pairs.gz scaffolds.cool
# Balance matrix
cooler balance scaffolds.cool
# Multi-resolution (mcool)
cooler zoomify scaffolds.cool -o scaffolds.mcool
Visualize with HiGlass
# Convert to higlass format
clodius aggregate bedfile --chromsizes-filename chrom.sizes \
--output-file scaffolds.beddb scaffold_boundaries.bed
# Load into higlass server
docker run --detach --publish 8888:80 \
--volume ~/hg-data:/data \
higlass/higlass-docker:latest
Manual Curation (Juicebox)
# Load .hic file in Juicebox Assembly Tools (JBAT)
# Perform manual corrections:
# - Break misjoins
# - Order/orient scaffolds
# - Merge scaffolds
# Export corrected assembly
# File -> Export Assembly -> FASTA
Post-Scaffolding Gap Filling
# TGS-GapCloser for long-read gap filling
tgsgapcloser --scaff scaffolds.fa --reads ont_reads.fq.gz \
--output filled --thread 16 --ne
# LR_Gapcloser alternative
LR_Gapcloser.sh -i scaffolds.fa -l ont_reads.fq.gz -t 16 -o gapclosed.fa
Validate Scaffolding
# Check chromosome-scale contiguity
seqkit stats scaffolds.fa
# BUSCO on scaffolds
busco -i scaffolds.fa -l eukaryota_odb10 -o busco_scaffolds -m genome -c 16
# N50/L50 statistics
assembly-stats scaffolds.fa
# Compare pre/post scaffolding
quast.py draft_assembly.fa scaffolds.fa -o quast_comparison
Check Telomeres
# Find telomeric repeats (vertebrate TTAGGG)
seqkit locate -i -p 'TTAGGG{10,}' scaffolds.fa > telomeres_forward.bed
seqkit locate -i -p 'CCCTAA{10,}' scaffolds.fa > telomeres_reverse.bed
# Count chromosomes with telomeres on both ends
awk '$2 < 1000' telomeres_forward.bed | cut -f1 | sort -u > left_telomeres.txt
awk -v OFS='\t' 'NR==FNR{len[$1]=$2;next} $3 > len[$1]-1000' \
scaffolds.fa.fai telomeres_reverse.bed | cut -f1 | sort -u > right_telomeres.txt
comm -12 left_telomeres.txt right_telomeres.txt > complete_chromosomes.txt
Rename to Chromosomes
# After manual curation, rename scaffolds to chromosomes
awk '/^>/{print ">chr" ++i; next}{print}' scaffolds.fa > chromosomes.fa
# Or with mapping file
seqkit replace -p '(.+)' -r '{kv}' -k name_mapping.tsv scaffolds.fa > chromosomes.fa
Related Skills
- genome-assembly/long-read-assembly - Generate initial contigs
- genome-assembly/assembly-polishing - Polish before scaffolding
- genome-assembly/assembly-qc - Validate final assembly
- hi-c-analysis/hic-data-io - Hi-C data processing
bg-szy