Bio Applied Viral Genome Assembly

Viral Genome Assembly and Variant Analysis

When to Use

• Assembling a consensus genome from ARTIC-tiled amplicon sequencing (SARS-CoV-2, Mpox, Influenza, RSV)
• Calling intra-host minority/low-frequency variants to characterize quasispecies structure
• Deciding between reference-guided vs. de novo assembly for a novel or divergent virus
• QC-gating sequencing runs before GISAID/NCBI submission (completeness, N-content, depth)
• Distinguishing amplicon-boundary artifacts from real low-frequency variants

Version Compatibility

• iVar ≥1.4.2, LoFreq ≥2.1.5, minimap2 ≥2.26, samtools ≥1.19, fastp ≥0.23
• Nextclade CLI ≥3.0 (dataset schema v2)
• Python ≥3.10, pandas ≥2.0, scipy ≥1.11, numpy ≥1.26

Prerequisites

conda install -c bioconda ivar=1.4.2 lofreq=2.1.5 minimap2 samtools fastp nextclade
pip install pandas scipy numpy

Prior concepts: FASTQ/read QC (bio-applied-ngs-fundamentals) and VCF fundamentals (bio-applied-variant-calling-and-snp-analysis).

Sequencing Strategy Selection

Key distinction from bacterial/human assembly: viruses exist as quasispecies — a cloud of related variants, not a single clone. Consensus ≠ population.

Coverage Thresholds

Reference-Guided Assembly Pipeline

Goal: turn raw ARTIC amplicon FASTQ into a consensus genome plus intra-host variant calls. Approach: trim adapters and primers before alignment-based consensus/variant calling — primer sequences left in place masquerade as high-frequency variants at amplicon boundaries.

# 1. Quality + adapter trim
fastp -i R1.fastq -I R2.fastq -o R1.trim.fastq -O R2.trim.fastq --thread 4

# 2. Align (Illumina short reads; use -ax ont2d for Nanopore)
minimap2 -ax sr ref.fa R1.trim.fastq R2.trim.fastq | samtools sort -o sorted.bam
samtools index sorted.bam

# 3. Trim ARTIC primers (soft-clips primer regions) -- required before any variant call
ivar trim -i sorted.bam -b primer_scheme.bed -p trimmed -m 20 -q 20
samtools sort -o trimmed.sorted.bam trimmed.bam && samtools index trimmed.sorted.bam

# 4. Consensus (min alt-allele freq 0.75, min depth 10 else N-mask)
samtools mpileup -aa -A -d 0 -Q 0 trimmed.sorted.bam | \
    ivar consensus -p consensus -t 0.75 -m 10

# 5. Intra-host variant calling
lofreq indelqual --dindel -f ref.fa -o trimmed.iq.bam trimmed.sorted.bam
lofreq call --call-indels -f ref.fa -o variants.vcf trimmed.iq.bam
# OR ivar variants (ARTIC standard)
samtools mpileup -aa -A -d 0 --reference ref.fa -Q 0 trimmed.sorted.bam | \
    ivar variants -p variants -q 20 -t 0.03 -m 10

# 6. Clade assignment + QC
nextclade run --input-fasta consensus.fa --input-dataset sars-cov-2 --output-tsv qc.tsv

Intra-Host Variant Filtering (Quasispecies)

Goal: separate real minority variants from sequencing/PCR noise in an iVar/LoFreq variant table. Approach: apply all four filters together — depth, allele frequency, strand bias, base quality — none is sufficient alone.

import pandas as pd
from scipy import stats


def apply_variant_filters(df_var, min_depth=100, min_af=0.01, sb_p_cutoff=0.001):
    """Flag PASS/FAIL for each row of a variant table (iVar/LoFreq output).

    df_var must have columns: Depth, Alt_freq, Strand_bias_p.
    Returns df_var with boolean pass_* columns and a combined PASS column.
    """
    df_var = df_var.copy()
    df_var['pass_depth'] = df_var['Depth'] >= min_depth
    df_var['pass_af'] = df_var['Alt_freq'] >= min_af
    df_var['pass_strand'] = df_var['Strand_bias_p'] > sb_p_cutoff
    df_var['PASS'] = df_var['pass_depth'] & df_var['pass_af'] & df_var['pass_strand']
    return df_var


def strand_bias_pvalue(depth, alt_freq, fwd_alt_frac):
    """Fisher exact test p-value for strand bias at a variant site.

    depth: total read depth at the site.
    alt_freq: alternate allele frequency (0-1).
    fwd_alt_frac: fraction of forward-strand reads among all reads (0-1),
    used as the null expectation for how alt/ref reads should split by strand.
    """
    fwd = int(depth * alt_freq * fwd_alt_frac)
    rev = int(depth * alt_freq * (1 - fwd_alt_frac))
    fwd_ref = int(depth * (1 - alt_freq) * fwd_alt_frac)
    rev_ref = int(depth * (1 - alt_freq) * (1 - fwd_alt_frac))
    _, p = stats.fisher_exact([[fwd, rev], [fwd_ref, rev_ref]])
    return p

Per-Sample QC Gating

Goal: decide which samples in a sequencing run are submission-ready. Approach: gate on mapped-read percentage and genome completeness at the consensus depth threshold (20×); anything failing either should be flagged for repeat sequencing, not submitted.

import pandas as pd

QC_THRESHOLDS = {
    'mapped_pct': 80.0,           # % reads mapping to reference
    'genome_completeness': 90.0,  # % genome covered at >=20x (non-N bases)
    'median_depth': 200,          # target depth for reliable variant calling
}


def sample_qc_pass(df_samples, thresholds=QC_THRESHOLDS):
    """Mark each sample PASS/FAIL against submission QC thresholds.

    df_samples must have columns: Mapped_pct, Completeness_pct.
    Returns df_samples with an added boolean QC_pass column.
    """
    df_samples = df_samples.copy()
    df_samples['QC_pass'] = (
        (df_samples['Mapped_pct'] >= thresholds['mapped_pct'])
        & (df_samples['Completeness_pct'] >= thresholds['genome_completeness'])
    )
    return df_samples

Minority Variant Tool Comparison

Pitfalls

• Primer dimers inflate low-quality variants: always trim ARTIC primers with ivar trim before variant calling — un-trimmed primer sequences appear as high-frequency variants at amplicon boundaries
• Amplicon dropout from primer mismatches: divergent lineages can fail to amplify → dropout regions show false N-masking in consensus; update the primer scheme each variant wave
• Coverage uniformity: ARTIC alternating-pool design means some regions get systematically lower depth — check per-amplicon coverage, not just genome average
• Quasispecies ≠ co-infection: intermediate frequencies (30–70%) can be drift, technical noise, or genuine mixed infection — require independent validation before claiming co-infection
• LoFreq requires realigned indels: run lofreq indelqual before lofreq call --call-indels for accurate indel variant calls
• iVar consensus Ns: positions below min_depth are written as N — a 90%-complete genome still has ~3kb masked; downstream tools (Pangolin, Nextclade) tolerate ≤30% Ns
• Coordinate systems: BED primer files are 0-based, iVar is 0-based internally, VCF output is 1-based — watch for off-by-one errors when cross-referencing positions

See Also

• bio-applied-ngs-fundamentals — FASTQ quality scores and read trimming fundamentals
• bio-applied-variant-calling-and-snp-analysis — VCF fields and genotype decoding shared with LoFreq/iVar output
• bio-applied-genome-assembly — de novo assembly algorithms for novel/divergent viruses
• bio-applied-variant-surveillance — Pangolin/Nextclade lineage assignment downstream of consensus calling