Bio Applied Variant Surveillance

Variant Surveillance and Wastewater Epidemiology

When to Use

• Assigning Pango lineages (e.g. BA.2, XBB.1.5) to SARS-CoV-2 consensus sequences from clinical or surveillance sequencing
• Deconvolving pooled wastewater sequencing into per-lineage relative abundances with Freyja
• Tracking spike protein mutations of concern (RBD, furin cleavage site) across a lineage time series
• Correlating wastewater viral load with clinical case counts to estimate early-warning lead time
• Classifying a lineage as VOC/VOI/VUM using WHO-style transmissibility/immune-escape/severity criteria

Version Compatibility

• Pangolin ≥4.3 with pangolin-data ≥1.29 (lineage designations update frequently — pin a dataset version for reproducibility)
• Nextclade ≥3.x (CLI) with the sars-cov-2 dataset
• Freyja ≥1.5 (requires usher barcode files, updated via freyja update)
• Python ≥3.10, pandas ≥2.0, numpy ≥1.24, scipy ≥1.11, matplotlib ≥3.8

Prerequisites

• Install: pip install freyja, pip install pangolin (or conda install -c bioconda pangolin nextclade)
• A reference-aligned consensus FASTA (for Pangolin/Nextclade) or a BAM aligned to the Wuhan-Hu-1 reference (for Freyja)
• Familiarity with bio-variant-calling-vcf-basics for interpreting the variants TSV Freyja produces
• Familiarity with bio-genome-intervals-bed-file-basics if masking primer/amplicon dropout regions

Lineage Assignment (Pangolin + Nextclade)

Goal: assign a Pango lineage and clade/QC flags to each consensus sequence in a FASTA. Approach: run Pangolin for the lineage call and Nextclade for clade + amino-acid changes + QC, then join on sequence name.

# Lineage assignment (writes lineage, conflict, ambiguity_score, scorpio_call)
pangolin sequences.fasta --outfile lineages.csv --threads 4

# Clade + quality flags + per-site amino acid substitutions
nextclade run --input-fasta sequences.fasta \
    --input-dataset sars-cov-2 --output-tsv nextclade.tsv

import pandas as pd

def load_lineage_calls(pangolin_csv: str, nextclade_tsv: str) -> pd.DataFrame:
    """Join Pangolin lineage calls with Nextclade QC flags on sequence name.

    Flags any sequence where Nextclade QC != 'good' so low-quality calls
    can be excluded before frequency estimation.
    """
    lin = pd.read_csv(pangolin_csv)
    nc = pd.read_csv(nextclade_tsv, sep="\t")
    merged = lin.merge(nc[["seqName", "clade", "qc.overallStatus"]],
                        left_on="taxon", right_on="seqName", how="left")
    merged["pass_qc"] = merged["qc.overallStatus"] == "good"
    return merged

Wastewater Deconvolution (Freyja)

Goal: estimate the fraction of each circulating lineage from a single pooled wastewater BAM. Approach: Freyja calls variants against the reference, then solves a constrained linear regression (observed mutation frequencies = barcode matrix × lineage fractions) using UShER-derived lineage-defining mutation barcodes.

# 1. Call variants + per-site depth from wastewater BAM
freyja variants wastewater.bam --variants variants.tsv --depths depths.tsv \
    --ref reference.fasta

# 2. Update lineage-defining mutation barcodes (do this regularly!)
freyja update

# 3. Deconvolve lineage fractions (linear mixture model)
freyja demix variants.tsv depths.tsv --output lineages.csv

# 4. Aggregate a directory of per-sample demix outputs into a time series
freyja aggregate --inputdir ./samples/ --output aggregated.tsv

# 5. Stacked-area plot of lineage fractions over time
freyja plot aggregated.tsv --output lineage_plot.pdf

Wastewater Lead-Time Analysis

Goal: quantify how many days wastewater viral load leads clinical case counts. Approach: smooth the noisy wastewater signal, then scan a range of lag offsets and take the lag with maximum Pearson correlation against clinical cases.

import numpy as np
from scipy.stats import pearsonr
from scipy.ndimage import gaussian_filter1d

def estimate_lead_time(viral_load: np.ndarray, cases: np.ndarray,
                        max_lag: int = 14) -> tuple[int, float]:
    """Find the lag (in days) at which wastewater viral load best predicts
    clinical case counts.

    Returns (best_lag_days, pearson_r) where wastewater[t] is compared
    against cases[t + lag] for lag in [0, max_lag].
    """
    ww_smooth = gaussian_filter1d(viral_load, sigma=3)
    best_lag, best_r = 0, -np.inf
    for lag in range(max_lag + 1):
        if lag == 0:
            r, _ = pearsonr(ww_smooth, cases)
        else:
            r, _ = pearsonr(ww_smooth[:-lag], cases[lag:])
        if r > best_r:
            best_lag, best_r = lag, r
    return best_lag, best_r

Reference Facts

Pango nomenclature: root lineages A/B (early Wuhan) → sub-lineages B.1 → B.1.1 → B.1.1.7 (Alpha); after ~3 sub-levels a new alias is issued (BA, BQ, XBB); recombinants get an X prefix. WHO tiers: VOC (Concern) > VOI (Interest) > VUM (Monitoring).

Spike protein regions: NTD (aa 13–305, antibody binding), RBD (aa 319–541, ACE2 contact — E484K/K417N/F486V drive immune escape, N501Y increases ACE2 affinity), furin cleavage site (aa 681–685, P681H/R increases fitness).

Wastewater epidemiology: detects roughly 1 infected person per 100,000 catchment population; typically provides 3–7 days lead time before a clinical case surge; population-level signal is unaffected by clinical testing access gaps; normalize viral RNA copies/L to PMMoV or crAssphage as a fecal-strength indicator before comparing across sites or time.

Pitfalls

• Amplicon dropout: primer mismatches in divergent lineages cause missing mutations and incorrect lineage calls — update primer schemes (e.g. ARTIC) with each major variant wave
• Freyja barcode staleness: lineage barcodes must be updated regularly (freyja update); stale barcodes misclassify new sublineages as "Other"
• Wastewater normalization: viral load varies with flow rate and precipitation — always normalize to PMMoV or crAssphage, never compare raw copies/L across sites
• Consensus vs. quasispecies: Pangolin assigns one lineage per consensus sequence; co-infections or recombinants require manual inspection of per-site variant frequencies (Freyja's variants output) rather than trusting a single consensus call
• Coordinate systems: BED is 0-based half-open; VCF/GFF are 1-based — off-by-one errors are common when merging Freyja/Nextclade output with BED-based masks
• Multiple testing: apply Benjamini-Hochberg FDR when testing associations across many genomic positions or lineages simultaneously

See Also

• bio-epidemiological-genomics-phylodynamics — for transmission tree and Rt estimation from the same sequence data
• bio-epidemiological-genomics-pathogen-typing — for non-SARS-CoV-2 pathogen lineage/strain typing
• bio-variant-calling-vcf-basics — for parsing the variants TSV/VCF that Freyja and Nextclade emit
• bio-population-genetics-selection-statistics — for testing whether a spike mutation shows signatures of positive selection