A Comprehensive Guide to Hybrid Assembly Pipeline for Genomic Sequencing

Long reads tell you how a microbial genome is spanned. Short reads tell you which bases to trust. A hybrid assembly uses each read set for the job it is good at: span first, then correct.

The pipeline below is written for Oxford Nanopore long reads and Illumina paired short reads. The commands are the ones this note was built around. Several flags are tied to the chemistry and to the tool version, and those are called out where a copy-paste would be the wrong experiment.

What to install

Anaconda or another conda distribution is the practical way to pin these tools. Medaka and BUSCO pull different stacks, so this note activates them in their own environments rather than forcing one environment to hold every package.

conda create --name myenv
conda activate myenv

Quality control:

• NanoPlot plots length and quality for long reads.
• FastQC reports the Illumina reads.
• Filtlong drops long reads that are short or low quality, and it can use the Illumina reads as an external reference while it does so.

Assembly and polishing:

• Flye builds the long-read assembly.
• Medaka polishes that assembly with a neural network trained for Nanopore reads.
• BWA maps the Illumina reads back to the polished contigs.
• Polypolish uses those alignments for a second, short-read polish.

Assessment and annotation come after the sequence exists:

• BUSCO asks whether expected single-copy genes are present.
• QUAST reports length, contiguity, GC, and, when you give it a reference, misassemblies. This note runs it without a reference.

The file names below (LR_input.fastq, SR_input_1.fastq, SR_input_2.fastq) are placeholders. Thread count is 8 throughout. Change both.

Quality control

Long reads:

mkdir -p nanoplot_lr_raw
NanoPlot --fastq LR_input.fastq --N50 --verbose --outdir nanoplot_lr_raw/ -t 8

Short reads:

mkdir -p fastqc_reports
fastqc SR_input_1.fastq SR_input_2.fastq -o fastqc_reports/ -t 8

Read the length distribution and the per-base quality before filtering. A filter cannot repair a library that failed. It can only discard reads.

Filter the long reads

Filtlong keeps long reads of at least 1000 bases and retains 90 percent of them by its ranking. The Illumina pair is the external reference for that ranking.

filtlong -1 SR_input_1.fastq -2 SR_input_2.fastq --min_length 1000 --keep_percent 90 LR_input.fastq > LR_filtered.fastq
mkdir -p nanoplot_lr_filtered
NanoPlot --fastq LR_filtered.fastq --N50 --verbose --outdir nanoplot_lr_filtered/ -t 8

Compare the second NanoPlot report with the first. You want the N50 to hold or improve because short fragments left, not because the filter kept a handful of reads.

Assemble with Flye

flye --nano-raw LR_filtered.fastq --out-dir Flye/ --threads 8 --scaffold -g 6m

Two flags are specific to this example. --nano-raw is Flye’s preset for that class of Nanopore reads. If your reads were basecalled with a high-quality model, use the matching Flye preset instead of copying this one. -g 6m tells Flye to expect a genome of about 6 megabases, which is a plausible bacterial size and a poor guess for anything else. Replace 6m with the size you actually expect.

Polish with the long reads

conda activate medaka
medaka_consensus -i LR_filtered.fastq -d Flye/assembly.fasta -o Polish1/ -m r941_min_fast_g303 -t 8

The model string r941_min_fast_g303 belongs to an older R9.4.1 and Guppy combination. Medaka will apply whatever model you name. Name the one that matches the flow cell and the basecaller that produced LR_filtered.fastq. The input assembly is Flye’s assembly.fasta. The polished sequence is written under Polish1/.

Polish with the short reads

Polypolish wants the Illumina alignments, including secondary and supplementary ones, which is why bwa mem is run with -a. The insert-filter helper is part of the Polypolish interface these commands were written for. A newer release may have renamed that step. Use the flags your installed version documents.

mkdir -p Polish2
bwa index Polish1/consensus.fasta
bwa mem -t 8 -a Polish1/consensus.fasta SR_input_1.fastq > Polish2/alignments_1.sam
bwa mem -t 8 -a Polish1/consensus.fasta SR_input_2.fastq > Polish2/alignments_2.sam

polypolish_insert_filter.py --in1 Polish2/alignments_1.sam --in2 Polish2/alignments_2.sam --out1 Polish2/filtered_1.sam --out2 Polish2/filtered_2.sam

polypolish Polish1/consensus.fasta Polish2/filtered_1.sam Polish2/filtered_2.sam > final_assembly.fasta

The file you assess from here on is final_assembly.fasta.

Did the assembly hold together?

BUSCO here uses the bacterial lineage bacteria_odb10. That lineage is the right check for the genome this example assumes. A different microbe needs its own lineage.

conda activate busco
busco -i final_assembly.fasta -l bacteria_odb10 -m genome -o busco_final_assembly

quast -o quast_final_assembly -t 8 final_assembly.fasta

Read the two reports together. BUSCO can look complete on an assembly whose contigs are still in the wrong order, and QUAST’s N50 can look excellent on a sequence that is missing genes. Without a reference, QUAST cannot call misassemblies.

Annotation

Prokka predicts genes and attaches the annotation it can infer. RAST is a second, web-based functional annotation. They will disagree on names. That is a reason to keep both, not a reason to average them.

conda install -c conda-forge -c bioconda prokka
prokka --outdir prokka_annotation --prefix final_assembly final_assembly.fasta

RAST takes the same assembly through its upload page: rast.nmpdr.org.

Bottom line

The long reads build the contigs. Medaka polishes them with the same long reads. Polypolish polishes again with the Illumina pair that already helped choose which long reads to trust. BUSCO and QUAST tell you whether that sequence is complete and contiguous. Prokka and RAST name the genes. The numbers in the commands that are not about your sample — genome size, Medaka model, Flye preset, BUSCO lineage — have to be changed to match it.