bio-epidemiological-genomics-variant-surveillance

Version Compatibility

Reference examples tested with: Nextclade 3.3+, ggplot2 3.5+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Variant Surveillance

"Classify my viral sequences into lineages" → Assign pathogen lineages and track variants of concern using Nextclade or pangolin for real-time genomic surveillance.

• CLI: nextclade run -d sars-cov-2 -i sequences.fasta
• CLI: pangolin sequences.fasta for SARS-CoV-2 Pango lineage assignment

Nextclade CLI

# Install Nextclade
npm install -g @nextstrain/nextclade

# Or download binary
curl -fsSL "https://github.com/nextstrain/nextclade/releases/latest/download/nextclade-x86_64-unknown-linux-gnu" -o nextclade
chmod +x nextclade

# List available datasets
nextclade dataset list

# Download dataset (e.g., SARS-CoV-2)
nextclade dataset get --name sars-cov-2 --output-dir data/sars-cov-2

# Run analysis
nextclade run \
    --input-dataset data/sars-cov-2 \
    --output-tsv results.tsv \
    --output-json results.json \
    sequences.fasta

Pangolin for SARS-CoV-2

# Install pangolin
pip install pangolin

# Update lineage definitions
pangolin --update

# Run lineage assignment
pangolin sequences.fasta -o pangolin_results.csv

# With specific version
pangolin sequences.fasta --analysis-mode accurate -o results.csv

Parse Nextclade Results

import pandas as pd

def parse_nextclade(results_file):
    '''Parse Nextclade TSV output

    Key columns:
    - seqName: Sequence identifier
    - clade: Nextstrain clade (e.g., 21L for Omicron BA.2)
    - Nextclade_pango: Pangolin lineage
    - qc.overallStatus: Quality control status
    - substitutions: List of mutations
    - aaSubstitutions: Amino acid changes
    '''
    df = pd.read_csv(results_file, sep='\t')

    # Filter by QC status
    df['pass_qc'] = df['qc.overallStatus'].isin(['good', 'mediocre'])

    return df


def summarize_lineages(results_df, lineage_col='Nextclade_pango'):
    '''Summarize lineage distribution'''
    # Filter passed QC
    passed = results_df[results_df['pass_qc']]

    summary = {
        'total_sequences': len(results_df),
        'passed_qc': len(passed),
        'unique_lineages': passed[lineage_col].nunique(),
        'lineage_counts': passed[lineage_col].value_counts().to_dict()
    }

    return summary

Track Variants of Concern

Goal: Classify viral sequences into WHO-defined variants of concern and track their prevalence over time.

Approach: Map Pango lineages to VOC labels using pattern matching, then group by time period and compute proportional representation of each VOC.

# WHO Variants of Concern/Interest definitions
VOC_DEFINITIONS = {
    'Alpha': {'lineages': ['B.1.1.7', 'Q.*'], 'key_mutations': ['N501Y', 'P681H']},
    'Beta': {'lineages': ['B.1.351'], 'key_mutations': ['K417N', 'E484K', 'N501Y']},
    'Gamma': {'lineages': ['P.1'], 'key_mutations': ['K417T', 'E484K', 'N501Y']},
    'Delta': {'lineages': ['B.1.617.2', 'AY.*'], 'key_mutations': ['L452R', 'P681R']},
    'Omicron': {'lineages': ['B.1.1.529', 'BA.*', 'XBB.*'], 'key_mutations': ['G339D', 'N501Y']}
}

def classify_voc(lineage):
    '''Classify lineage as VOC'''
    for voc, definition in VOC_DEFINITIONS.items():
        for pattern in definition['lineages']:
            if pattern.endswith('*'):
                if lineage.startswith(pattern[:-1]):
                    return voc
            elif lineage == pattern:
                return voc
    return 'Other'


def track_voc_prevalence(results_df, date_col='collection_date'):
    '''Track variant of concern prevalence over time'''
    results_df = results_df.copy()
    results_df['VOC'] = results_df['Nextclade_pango'].apply(classify_voc)

    # Group by week
    results_df['week'] = pd.to_datetime(results_df[date_col]).dt.to_period('W')

    prevalence = results_df.groupby(['week', 'VOC']).size().unstack(fill_value=0)
    prevalence_pct = prevalence.div(prevalence.sum(axis=1), axis=0) * 100

    return prevalence_pct

Mutation Analysis

def parse_mutations(mutation_string):
    '''Parse Nextclade mutation string

    Format: 'A123T,C456G' (nucleotide) or 'S:N501Y,S:D614G' (amino acid)
    '''
    if pd.isna(mutation_string) or mutation_string == '':
        return []
    return mutation_string.split(',')


def find_mutation_prevalence(results_df, mutation_col='aaSubstitutions'):
    '''Calculate prevalence of each mutation'''
    all_mutations = []
    for muts in results_df[mutation_col].dropna():
        all_mutations.extend(parse_mutations(muts))

    mutation_counts = pd.Series(all_mutations).value_counts()
    mutation_prevalence = mutation_counts / len(results_df) * 100

    return mutation_prevalence


def detect_emerging_mutations(results_df, date_col='collection_date', threshold=5):
    '''Detect mutations increasing in frequency'''
    # Goal: Flag mutations showing rapid prevalence increase between
    # early and recent time periods (>2-fold increase above threshold).
    # Approach: Split data at median date, compute per-mutation prevalence
    # in each half, and report mutations with significant frequency gains.
    results_df = results_df.copy()
    results_df['date'] = pd.to_datetime(results_df[date_col])

    # Split into early and recent
    midpoint = results_df['date'].median()
    early = results_df[results_df['date'] < midpoint]
    recent = results_df[results_df['date'] >= midpoint]

    early_prev = find_mutation_prevalence(early)
    recent_prev = find_mutation_prevalence(recent)

    # Find emerging (low->high)
    emerging = []
    for mut in recent_prev.index:
        early_freq = early_prev.get(mut, 0)
        recent_freq = recent_prev[mut]

        if recent_freq > threshold and recent_freq > early_freq * 2:
            emerging.append({
                'mutation': mut,
                'early_prevalence': early_freq,
                'recent_prevalence': recent_freq,
                'fold_change': recent_freq / max(early_freq, 0.1)
            })

    return sorted(emerging, key=lambda x: -x['fold_change'])

Surveillance Report

def generate_surveillance_report(results_df, period='week'):
    '''Generate variant surveillance report'''
    passed = results_df[results_df['pass_qc']]

    report = {
        'period': period,
        'total_sequences': len(results_df),
        'passed_qc': len(passed),
        'qc_pass_rate': f"{len(passed)/len(results_df)*100:.1f}%"
    }

    # Lineage distribution
    lineage_counts = passed['Nextclade_pango'].value_counts()
    report['dominant_lineage'] = lineage_counts.index[0]
    report['dominant_lineage_pct'] = f"{lineage_counts.iloc[0]/len(passed)*100:.1f}%"
    report['top_5_lineages'] = lineage_counts.head(5).to_dict()

    # VOC tracking
    passed['VOC'] = passed['Nextclade_pango'].apply(classify_voc)
    voc_counts = passed['VOC'].value_counts()
    report['voc_distribution'] = voc_counts.to_dict()

    return report

Related Skills

• epidemiological-genomics/phylodynamics - Time-scaled analysis of variants
• variant-calling/variant-annotation - Mutation annotation
• data-visualization/ggplot2-fundamentals - Visualize variant dynamics