# Bio Applied Variant Surveillance

> Assign SARS-CoV-2 Pango lineages with Pangolin/Nextclade and deconvolve wastewater into lineage fractions with Freyja (variants/demix). Use when doing variant surveillance, wastewater epidemiology, lineage tracking, or VOC/VOI classification.

- Skill: `pavel-kravchenko/bio-applied-variant-surveillance` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-variant-surveillance`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-variant-surveillance/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```bash
# 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
```

```python
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.

```bash
# 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.

```python
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