# Virology Bioinformatics

> Assemble viral genomes with iVar/minimap2, call intra-host SNVs with LoFreq, assign Nextclade/pangolin lineages. Use when trimming ARTIC primers, building a consensus FASTA, calling minority variants, or assigning a Pango lineage.

- Skill: `pavel-kravchenko/virology-bioinformatics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/virology-bioinformatics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/virology-bioinformatics/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/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

| Task | Tool | Notes |
|------|------|-------|
| Trim ARTIC primers | iVar trim | `-b primer.bed -e` |
| Align to reference | Minimap2 / BWA-MEM | `-ax sr` for Illumina |
| Consensus genome | iVar consensus | `-t 0.5` (majority base) |
| Low-freq variants | LoFreq | SNV at ≥1% frequency |
| Clade assignment | Nextclade | SARS-CoV-2, flu, RSV |
| Lineage classification | pangolin | Pango lineage |
| Timed phylogeny | Nextstrain Augur | TreeTime integration |
| Wastewater deconvolution | Freyja | Variant proportions |

## SARS-CoV-2 ARTIC Assembly Pipeline

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

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

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

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

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

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

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

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