Virology Bioinformatics

virology-bioinformatics

When to Use

• Assembling a SARS-CoV-2/influenza/RSV consensus genome from ARTIC-tiled amplicon Illumina reads
• Calling intra-host low-frequency variants (minor allele deconvolution, coinfection/reinfection detection)
• Assigning Pango lineages or Nextclade clades to a consensus FASTA
• Building a time-resolved (molecular clock) phylogeny for outbreak investigation
• Deconvolving pooled wastewater sequencing into per-lineage proportions for genomic surveillance

Version Compatibility

• iVar ≥1.4.2, minimap2 ≥2.26, samtools ≥1.19
• LoFreq ≥2.1.5
• Nextclade CLI ≥3.8, pangolin ≥4.3 (pangolin-data ≥1.25)
• Nextstrain Augur ≥24.4 (TreeTime ≥0.11)
• Freyja ≥1.5
• Python ≥3.10, pandas ≥2.2

Prerequisites

• Bioconda env with ivar, minimap2, samtools, lofreq, nextclade, pangolin, augur, freyja installed (conda install -c bioconda -c conda-forge <tool>)
• Reference genome + ARTIC primer BED (e.g. NC_045512.2 + nCoV-2019/V4.1 scheme); wastewater workflow additionally needs a UShER barcode set (bundled with Freyja)
• Familiarity with BAM/VCF fundamentals (bio-applied-ngs-fundamentals) and general tree building (bio-core-phylogenetics)

Quick Reference

SARS-CoV-2 ARTIC Assembly Pipeline

# 1. Align reads to Wuhan-Hu-1 reference (NC_045512.2)
minimap2 -ax sr NC_045512.2.fa sample_R1.fastq.gz sample_R2.fastq.gz \
    | samtools sort -o sample.bam && samtools index sample.bam

# 2. Trim ARTIC v4.1 primers
ivar trim -i sample.bam -b nCoV-2019.primer.bed -p sample_trimmed -e
samtools sort -o sample_trimmed.sorted.bam sample_trimmed.bam
samtools index sample_trimmed.sorted.bam

# 3. Generate consensus (N-mask positions below 20x coverage)
samtools mpileup -A -d 0 -Q 0 sample_trimmed.sorted.bam \
    | ivar consensus -p consensus -n N -m 20 -t 0.5

# 4. Quality check
samtools flagstat sample_trimmed.sorted.bam
samtools depth -a sample_trimmed.sorted.bam > sample.depth.tsv

Goal: turn the samtools depth table into a pass/fail genome-completeness QC metric before trusting the consensus. Approach: read the per-base depth file and compute the fraction of the genome above the consensus (20x) and variant-calling (200x) thresholds.

import pandas as pd


def genome_completeness(depth_tsv: str, genome_len: int = 29903,
                         thresholds: tuple[int, ...] = (20, 200)) -> dict[int, float]:
    """Compute % of genome positions at or above each depth threshold.

    depth_tsv: output of `samtools depth -a ref.bam` (chrom, pos, depth; no header).
    genome_len: expected reference length, used to catch missing/truncated positions.
    Returns {threshold: pct_genome_covered}.
    """
    depth = pd.read_csv(depth_tsv, sep="\t", header=None,
                         names=["chrom", "pos", "depth"])
    if len(depth) < genome_len:
        # missing rows means samtools never emitted those positions (no coverage at all)
        missing = genome_len - len(depth)
        depth = pd.concat([depth, pd.DataFrame({"depth": [0] * missing})], ignore_index=True)
    return {t: round(100 * (depth["depth"] >= t).mean(), 1) for t in thresholds}


if __name__ == "__main__":
    print(genome_completeness("sample.depth.tsv"))

Intra-Host Variant Calling with LoFreq

# Recalibrate base qualities (optional but recommended)
samtools calmd -b sample_trimmed.sorted.bam NC_045512.2.fa > sample_calmd.bam

# Call variants at >=1% frequency
lofreq call-parallel --pp-threads 8 \
    -f NC_045512.2.fa \
    -o sample_lofreq.vcf \
    --sig 0.01 --bonf dynamic \
    sample_calmd.bam

# Filter: min AF 1%, min depth 100x
lofreq filter -i sample_lofreq.vcf -o sample_filtered.vcf \
    --af-min 0.01 --cov-min 100

Goal: load the filtered LoFreq VCF into a DataFrame to inspect minority variants (e.g. flag possible coinfection/reinfection). Approach: parse the VCF INFO field (LoFreq stores AF, DP, SB as key=value pairs) without a heavyweight VCF library.

import pandas as pd


def parse_lofreq_vcf(vcf_path: str, af_min: float = 0.01) -> pd.DataFrame:
    """Parse a LoFreq VCF into a tidy variant table.

    Extracts AF (allele frequency), DP (depth), SB (strand bias) from INFO.
    Returns rows with AF >= af_min, sorted by position.
    """
    rows = []
    with open(vcf_path) as fh:
        for line in fh:
            if line.startswith("#"):
                continue
            fields = line.rstrip("\n").split("\t")
            chrom, pos, _id, ref, alt, qual, flt, info = fields[:8]
            info_dict = dict(kv.split("=", 1) for kv in info.split(";") if "=" in kv)
            af = float(info_dict.get("AF", 0.0))
            if af < af_min:
                continue
            rows.append({
                "chrom": chrom, "pos": int(pos), "ref": ref, "alt": alt,
                "af": af, "depth": int(info_dict.get("DP", 0)),
                "strand_bias": float(info_dict.get("SB", 0)),
                "filter": flt,
            })
    return pd.DataFrame(rows).sort_values("pos").reset_index(drop=True)


if __name__ == "__main__":
    variants = parse_lofreq_vcf("sample_filtered.vcf")
    print(f"{len(variants)} variants >=1% AF; "
          f"{(variants['af'] < 0.05).sum()} are minority (<5%) variants")

Nextclade + pangolin Lineage Assignment

nextclade run --dataset-name sars-cov-2 --output-tsv nextclade.tsv consensus.fa
pangolin consensus.fa --outfile lineage_report.csv

import subprocess
import pandas as pd


def assign_lineage(consensus_fa: str) -> pd.DataFrame:
    """Run Nextclade against a consensus FASTA and return clade/QC columns.

    Requires the `sars-cov-2` Nextclade dataset (`nextclade dataset get --name sars-cov-2 --output-dir data/`).
    """
    subprocess.run(
        ["nextclade", "run", "--dataset-name", "sars-cov-2",
         "--output-tsv", "nextclade.tsv", consensus_fa],
        capture_output=True, text=True, check=True,
    )
    results = pd.read_csv("nextclade.tsv", sep="\t")
    return results[["seqName", "clade", "Nextclade_pango", "qc.overallStatus"]]


if __name__ == "__main__":
    print(assign_lineage("consensus.fa"))

Nextstrain Augur Timed Phylogeny

# Filter and subsample
augur filter \
    --sequences sequences.fasta \
    --metadata metadata.tsv \
    --output filtered.fasta \
    --group-by country year month \
    --sequences-per-group 10 \
    --min-date 2020-01-01

# Align to reference
augur align \
    --sequences filtered.fasta \
    --reference-sequence NC_045512.2.gbk \
    --output aligned.fasta --fill-gaps

# Build tree (IQ-TREE under the hood)
augur tree \
    --alignment aligned.fasta \
    --output tree_raw.nwk \
    --nthreads 8

# Timed tree with TreeTime
augur refine \
    --tree tree_raw.nwk \
    --alignment aligned.fasta \
    --metadata metadata.tsv \
    --timetree --coalescent skyline \
    --output-tree tree.nwk \
    --output-node-data branch_lengths.json

Freyja Wastewater Deconvolution

# Call variants + depths directly from a mixed-population wastewater BAM
freyja variants wastewater_sample.bam \
    --variants variants.tsv \
    --depths depths.tsv \
    --ref NC_045512.2.fasta

# Deconvolve variant proportions against the UShER barcode matrix (NNLS)
freyja demix variants.tsv depths.tsv --output demixed.tsv

# Aggregate multiple samples into a time series
freyja aggregate output_dir/ --output aggregated.tsv
freyja plot aggregated.tsv --output wastewater_variants.pdf --lineages

Real-world reference: Karthikeyan et al. 2022 (Nature 609:101-108, PRJNA819090) built this exact iVar-trim → minimap2 → freyja variants/demix/aggregate pipeline on UC San Diego wastewater and clinical amplicon data, detecting VOCs up to 14 days before they appeared in clinical genomic surveillance.

Pitfalls

• Amplicon dropout: some ARTIC amplicons fail; expect coverage gaps at pool boundaries
• Consensus masking: mask with N (not a random base) at positions below the depth threshold
• Primer contamination: for amplicon data, primers must be trimmed BEFORE variant calling or every read end looks like a fixed "variant"
• Phylogenetic signal: check root-to-tip regression (TempEst) before trusting augur refine's clock rate
• Recombinants: SARS-CoV-2 XBB-lineage genomes show phylogenetic incongruence and can break a single-tree model
• Freyja barcode lag: the UShER barcode library lags newly designated lineages by days-to-weeks, so brand-new variants can be under- or mis-called

Key Databases

• GISAID EpiCoV: SARS-CoV-2 genomes (>15M sequences)
• NCBI Virus / SRA: all viral genomes and raw reads with metadata
• Nextstrain: real-time phylogenies for flu, COVID, mpox
• Pango designations: lineage nomenclature repository

See Also

• bio-applied-variant-calling-and-snp-analysis — general VCF/SNP calling and filtering concepts
• bio-applied-phylodynamics — deeper molecular-clock and coalescent modeling
• bio-applied-variant-surveillance — lineage/mutation tracking beyond a single genome
• bio-core-phylogenetics — tree building/manipulation fundamentals underlying Augur