| Genomic Surveillance & Analytics
Tools for Tomorrow Webinar Series
NIAID BRC-Analytics Showcase August 2026 // Tools for Tomorrow Webinar

Genomic Analytics for Cyclospora Outbreaks

From raw amplicon reads to epidemiological clusters: reproducible workflows, empirical benchmarking, and open surveillance infrastructure.

Consortium
BRC-Analytics
Penn State · Temple · NIAID BRCs
Focus Application
Foodborne Surveillance
8-Marker MLST · PyEuk Engine
Platform
Galaxy Ecosystem
Open, Pinned & Containerized
Deliverable 1: a complete Cyclospora typing workflow, from reads to classification.
Deliverable 2: a Logan / LexicMap search across all of the SRA, finding Cyclospora in runs nobody labelled as such.
02 // Epidemiological Urgency

The Ongoing 2026 Outbreak & The Public Data Gap

A historic surge in domestic cases highlights the critical need for open, verified genomic tools.

21,000+ Cases
Active Summer 2026 Outbreak
  • Multi-state surge across 15+ states linked to commercial bagged salads and greens.
  • High caseload requires molecular typing to separate concurrent contamination events.
0 Public Runs
The Public Data Void
  • Zero raw FASTQ runs or haplotype sheets deposited in NCBI SRA/ENA archives.
  • Prevents independent external re-analysis and algorithmic validation in real time.
Open Science
Validation Strategy
  • Workflows validated on confirmed 2018 gold-standard outbreak cohorts (PRJNA578931).
  • Open Galaxy tools and containerized engines ready for public health lab deployment.
With no public sequence data available for 2026, rigorous benchmarking against historical outbreaks ensures public health laboratories have vetted, reproducible tools ready to deploy.
03 // Biological Foundations

Why Cyclospora Cannot Be Sequenced Like Bacteria

Lack of culture models and extreme metagenomic dilution require targeted amplicon sequencing.

Life cycle and transmission NO Direct Person-to-Person Transmission
Human Host Jejunal Enterocytes
Unsporulated Oocyst Non-Infectious at Excretion
7–13 Days Maturation 22–30°C · Humidity · O₂
Sporated Oocyst Infectious (2 sporocysts)
Produce / Water Vehicle Berries, Cilantro, Lettuce
Laboratory Constraints
No In Vitro Culture No continuous cell lines support propagation
No Animal Model Strict human tropism; all challenge models failed
Human Stool Sole biological source of parasite DNA
Metagenomic Dilution & SRA Archive
Clinical Stool DNA Extract Composition < 0.6% Parasite Target
>99.4% Host & Microbial Flora
Public Sequencing Strategies (9,054 SRA Runs) 99.6% Amplicon Panel
9,016 Amplicon Runs (8,522 CDC 8-Marker)
CDC Clinical Stool

Surveils human cases via state health labs (PulseNet/SEDRIC); links multi-state patient clusters using the 8-marker amplicon panel.

FDA Food & Water

Samples produce vehicles (lettuce, cilantro, berries) and agricultural water; coordinates farm inspections, import alerts, and product recalls.

NIH / BRC Open Analytics

Provides open, reproducible bioinformatics infrastructure (Galaxy), empirical benchmarks, and genome resources across UCSC & VEuPathDB.

04 // Genomic Resources

No Reference-Grade Genome Exists

Forty-nine assemblies are public, none chromosome-level — which is why typing is done from amplicons rather than genomes.

Contig counts and N50 for all 49 public C. cayetanensis assemblies
None is chromosome-level, two of 49 are annotated, and the median assembly is in 1,391 pieces with a 103 kb contig N50. For scale, the P. falciparum reference is 14 contigs.
05 // Genomic Architecture

Organellar Conservation vs. Nuclear Heterogeneity

High-copy organellar genomes lack diversity; polymorphic nuclear loci introduce multi-clonal mixtures.

6.3 kb Linear Mito 34 kb Apicoplast 67–513× Copy Number
Organellar Genomes (Hyper-Conserved)
cox1
cox3
cytb (100% id)
15-mer Junction
  • Core mitochondrial genomes across global clinical isolates contain only 9 to 12 SNPs (<0.1% sequence divergence).
  • Identical sequence types appear in both domestic outbreaks and unrelated international travel cases, preventing source attribution.
Polymorphic Loci Sexual Recombination Polyclonal Mixtures
Nuclear Loci & The "Bag of Alleles" Problem
Patient A (Stool)
Locus 378: {Hap 2, Hap 4, Hap 9}
Patient B (Stool)
Locus 378: {Hap 4, Hap 9}
  • Multi-strain mixtures are the empirical default: 92.2% of 2018 outbreak cases (141 / 153) carry multi-allelic calls (\( \text{MOI} \ge 2 \)), and the CDC national surveillance cohort averages 26.1 distinct alleles per patient.
  • Specimens cannot be directly aligned as single consensus sequences without discarding intra-host diversity, requiring methods that compare allele distributions or shared haplotype sets.
Unlike clonal bacteria where 1 isolate = 1 consensus sequence, a Cyclospora specimen is an intra-host population of alleles. Because specimens cannot be aligned directly as single sequences, distance is computed by comparing allele frequency distributions or scoring the overlap and population rarity of shared haplotypes (such as weighted Identity-By-State).
06 // Marker Provenance & Design

The 8-Marker MLST Strategy & Foundational Studies

Multi-Locus Sequence Typing (MLST) provides a proven targeted framework across fragmented eukaryotic pathogen genomes.

MLST in Eukaryotic Parasitology Gold Standard

Multi-Locus Sequence Typing (MLST) is the standard molecular typing strategy for unculturable and polymorphic eukaryotic parasites (such as Plasmodium, Cryptosporidium, and Giardia), bypassing metagenomic dilution by targeting high-diversity loci.

UCSC Genome Browser & Assembly Landscape 31 Hubs Online

49 C. cayetanensis draft assemblies exist in GenBank, with 31 available via UCSC Assembly Hubs and BRC-Analytics. Because genomes are fragmented (median 1,391 contigs, N50 103 kb, 0 chromosome-level), routine surveillance anchors on targeted amplicons.

STUDY 1 // CORE PANEL & HEURISTICS
Nu_378 · Nu_360i2 · Mt_MSR

Established the initial multi-locus amplicon typing scheme and the heuristic genetic distance framework for outbreak case linkage.

Barratt et al. (2019) Parasitology
Nascimento et al. (2020) Epidemiol Infect
STUDY 2 // WHOLE-GENOME MINED CDS
Nu_CDS1 · Nu_CDS2 · Nu_CDS3 · Nu_CDS4

Mined 13 coding SNPs across four whole-genome draft assemblies to expand discriminatory resolution for epidemiological clustering.

Houghton et al. (2020) Parasite
STUDY 3 // MITOCHONDRIAL REPEAT JUNCTION
Mt-Junction (15-mer Repeats)

Targeted tandem concatemers of 15-mer repeat motifs in the mitochondrial genome, resolving 14 distinct sequence types in clinical isolates.

Nascimento et al. (2019) Emerg Infect Dis
07 // Target Specifications & Architecture

The 8-Marker Genotyping Panel: Target Specifications

Six nuclear and two mitochondrial loci capture single nucleotide variants and tandem repeat arrays.

6 Nuclear Amplicons (1,712 bp total) 2 Mitochondrial Amplicons (671–701 bp total) 24 PART Windows (Zenodo: 21924355)
Locus Target / Gene Description Amplicon Size Calling Windows SNPs Population Diversity Intra-Host Alleles
Nu_360i2 Nuclear intronic region (speciation driver) 488 bp 6 PARTs (A–F, ~100 bp) 20 >30 haplotypes (Species A vs B split) 1–3
Nu_378 Sec14 cytosolic factor-like protein 473 bp 4 PARTs (A–D, ~100 bp) 16 >25 haplotypes (Primary strain typing) 1–3
Nu_CDS1 ATP synthase subunit (LOC34619420) 175 bp 2 PARTs (A: 67 bp, B: 68 bp) 7 6–10 haplotypes 1–2
Nu_CDS2 Hypothetical protein (LOC34617424) 185 bp 2 PARTs (A: 103 bp, B: 103 bp) 4 4–6 haplotypes 1–2
Nu_CDS3 Conserved hypothetical protein 212 bp 2 PARTs (A: 89 bp, B: 89 bp) 2 3–4 haplotypes 1–2
Nu_CDS4 ATP-dependent RNA helicase rrp3 179 bp 2 PARTs (A: 69 bp, B: 69 bp) 3 4–6 haplotypes 1–2
Mt_MSR Mitochondrial 16S ribosomal RNA 487 bp 6 PARTs (A–F, ~100 bp) 5 4–6 haplotypes (<0.1% divergence) 1
Mt-Junction Mitochondrial linear repeat array 184–214 bp 1 Array (20 Cmt references) Structural 20 sequence types 1
Across the 7 alignable markers, the panel contains 24 distinct PART calling windows (keyed to 78 named reference haplotypes); individual isolates typically present 1–3 alleles at nuclear loci and 1 homoplasmic allele for mitochondrial loci.
08 // Haplotype Topology: Coding Exons

Conserved Coding Loci: Linear Stepwise Divergence

Conserved metabolic exons like Nu_CDS3 exhibit simple single-nucleotide mutational trajectories.

Nu_CDS3 PART A — Conserved Coding Marker (89 bp) 3 Haplotypes · Linear Chain 1 SNP 1 SNP Hap 1 Reference State Hap 2 Intermediate Hap 3 Derived Variant Single-nucleotide stepwise divergence without branching; typical of conserved apicomplexan metabolic genes.
Essential coding exons across the 4 CDS markers provide invariant or 1-SNP baseline anchors across global clinical surveillance.
09 // Haplotype Topology: Strain Typing

Strain Discrimination: Bimodal Hub-and-Spoke Networks

Nu_378 provides the core discriminatory power that resolves distinct clinical outbreak sources.

Nu_378 PART D — Core Strain Discrimination Network (114 bp) 7 Haplotypes · Bimodal Topology Lineage I (McDonald's Salads) Lineage II (Del Monte Veg Trays) 2 SNPs 4 SNPs 4 SNPs (Bridge) 2 SNPs 2 SNPs 4 SNPs Hap 1 Cluster I Core Hap 3 Hap 4 Intermediary Hap 2 Primary Hub Hap 5 Hap 7 Hap 6 Two dominant lineages separated by 4–8 SNPs separate distinct outbreak vehicles in historical surveillance.
A 4–8 SNP distance separates verified outbreak lineages, distinguishing the 2018 McDonald's salad cluster from the Del Monte vegetable tray cluster.
10 // Haplotype Topology: Speciation

Cryptic Speciation: The Deep Intronic Divide

Nu_360i2 reveals fixed 4–5 SNP boundaries separating sympatric North American lineages.

Nu_360i2 PART E — The Speciation Divide (100 bp) 5 Haplotypes · Reproductive Barrier Species A (Lineage A) Cyclospora cayetanensis Species B (Lineage B) Cyclospora ashfordi 4–5 SNPs (Species Rift) 1 SNP 1 SNP 2 SNPs Hap 1 A-Allele Type Hap 4 B Core Hub Hap 2 Hap 3 Hap 5 Intronic markers reveal fixed 4–5 SNP boundaries with near-zero heterozygous hybrids in North American surveillance.
Intronic markers identify distinct evolutionary lineages, providing high-resolution taxonomic resolution for surveillance.
11 // National Surveillance Standard

The National Surveillance Standard: Role & Real-World Impact

The operational framework connecting clinical infections to contaminated produce.

Operational CDC Standard >9,000 Clinical Isolates Typed Directs FDA Traceback & Recalls
Surveillance Mission
  • Provides objective genetic evidence to complement public health investigations.
  • Assists in differentiating distinct contamination events during periods of high incidence.
  • Supports data-driven public health actions to mitigate foodborne risks.
Ingestion & Surveillance Deliverables
  • Processes Illumina MiSeq FASTQ files across targeted marker amplicons.
  • Generates a structured Haplotype Data Sheet (HDS) recording allelic profiles.
  • Constructs transmission assessments to assist in epidemiological clustering.
Used continuously since 2018, this framework has provided essential data to support public health responses to multiple foodborne outbreaks.
12 // Epidemiological Context

The 2018 Outbreak Cohort: Epidemiological Context & Complexity

Two massive, concurrent multi-state outbreaks with overlapping geography, timelines, and vehicle ingredients.

Concurrent Summer 2018 Outbreaks 203 Sequenced Isolates (PRJNA578931) Overlapping Midwest Geography
Outbreak A: Pre-Packaged Salads
511 Cases · 24 Hospitalized · 16 States
  • Vehicle: Fresh Express romaine lettuce and carrot mix served at fast-food restaurants.
  • Geography: Concentrated across IL, IA, MO, KY, OH, and NE.
  • Genomic Signature: High frequency of Nu_378_D_Hap_7 (86.7%) and Mt_Cmt199_Hap_17 (63.3%).
Outbreak B: Fresh Vegetable Trays
250 Cases · 8 Hospitalized · 4 States
  • Vehicle: Del Monte pre-packaged trays (broccoli, cauliflower, carrots, and dill dip) sold in grocery stores.
  • Geography: Concentrated across IA, WI, MN, and MI.
  • Genomic Signature: High frequency of Nu_378_D_Hap_2 (94.5%) and Mt_Cmt169_Hap_8 (74.5%).
The Epidemiological Conundrum
  • Spatiotemporal Overlap: Both outbreaks occurred concurrently in identical Midwestern states (May–July 2018).
  • Ingredient Confounding: Both vehicles contained raw carrots; patients frequently reported multiple produce exposures or generic salad consumption.
  • Ground-Truth Validation: Supplier invoices and restaurant receipts confirmed patient exposures, creating a definitive national benchmark for genotyping algorithms.
The 2018 outbreak cohort (PRJNA578931, n=203) provides the definitive gold-standard benchmark for genomic epidemiology: two distinct outbreaks indistinguishable by patient interviews alone, requiring high-resolution molecular typing.
13 // Surveillance Data in Practice

The Haplotype Data Sheet & Outbreak Clustering

How binary presence/absence matrices resolve clinical cases into distinct outbreak source clusters.

2018 Multi-State Outbreak Data (PRJNA578931) 2,354 Archive Specimens × 165 Haplotype Columns Cluster Threshold \( d \le 0.10 \)
Haplotype Data Sheet Matrix Presence (X) / Absence (·)
Specimen ID 378_D_H7 378_D_H2 Cmt199_H17 Cmt169_H8 CDS1_B_H2 CDS1_B_H1
C_IA031_18 X · · · · ·
C_IA034_18 X · · · X ·
C_IA039_18 X · X · · ·
C_IA013_18 · X · X · ·
C_IA018_18 · X · X · X
C_WI008_18 · X · X · ·
Cluster 1 Assignment (Vendor A) McDonald's Salads

Defined by Nu_378_D_Hap_7 (86.7% frequency) and Mt_Cmt199_Hap_17 (63.3%). Linked to Fresh Express salad mix across IL, IA, MO, KY.

Cluster 2 Assignment (Vendor B) Del Monte Veg Trays

Defined by Nu_378_D_Hap_2 (94.5% frequency) and Mt_Cmt169_Hap_8 (74.5%). Linked to Del Monte vegetable trays across IA, WI, MN, MI.

Even when patients cannot recall specific ingredients, mutually exclusive haplotype signatures cleanly separate concurrent outbreaks into distinct supplier recalls.
14 // Analytical Flowchart

How the Pipeline Works: From Reads to Clusters

Three sequential computational stages translate raw sequencing data into actionable outbreak clusters.

Stage 1 Read Matching

Identify Haplotypes from Reads

Trims sequencing primers and assembles variable repeat regions from raw paired-end reads.

Clusters reads by sequence similarity and matches them against cataloged references to record detected alleles.

Output: Binary Haplotype Matrix
Stage 2 Pairwise Comparison

Compute Pairwise Distances

Compares every specimen against every other specimen across all shared marker regions.

Combines heuristic scoring with a probabilistic model that gives stronger weight to shared rare variants.

Output: Pairwise Distance Matrix
Stage 3 Hierarchical Tree

Group Related Outbreaks

Builds a hierarchical clustering tree that groups genetically identical or near-identical isolates.

Applies an empirical cutoff calibrated on historical outbreaks to define discrete epidemiological clusters.

Output: Outbreak Cluster Definitions
Raw FASTQ Reads
Haplotype Matrix (HDS)
Pairwise Distances
Outbreak Clusters
15 // Surveillance Deliverable

Surveillance Output: The Clustering Tree & Outbreak Definitions

How pairwise genetic distances are translated into actionable public health clusters using hierarchical clustering.

Ward's Hierarchical Linkage Empirical Cutoff Threshold (mean + 3*sd) Discrete Traceback Clusters
Output Dendrogram (Hierarchical Tree) 2 Outbreaks Resolved
0.25 0.20 0.15 0.10 0.05 Distance Cutoff (Threshold) IL036 IL126 MO003 IA014.. WI100 WI008 MN045 MI012.. SP012 Cluster 1 (Salads / Vendor A) Cluster 2 (Veg Trays / Vendor B) Sporadic
Isolates merging below the red cutoff line (\( d \le 0.05 \)) are grouped into the same outbreak cluster; branches above the cutoff represent unrelated sporadic cases.
Primary Deliverable: Cluster Membership Sheet Exported Table
Specimen ID Cluster Assignment Epidemiological Source
CDC_2018_IL036 Cluster 1 McDonald's Salad (Vendor A)
CDC_2018_IL126 Cluster 1 McDonald's Salad (Vendor A)
CDC_2018_WI100 Cluster 2 Del Monte Veg Trays (Vendor B)
CDC_2018_WI008 Cluster 2 Del Monte Veg Trays (Vendor B)
CDC_2018_SP012 Unassigned Sporadic (Above Cutoff)
How Public Health Officials Use This Output
  • Cluster assignments group patient exposure questionnaires to pinpoint contaminated ingredients across state lines.
  • Clear statistical separation between clusters enables independent agricultural recalls without misidentifying food suppliers.
The dendrogram and cluster assignment sheet are the primary scientific deliverables generated by the pipeline, allowing CDC and FDA to connect geographically scattered patient illnesses to specific agricultural suppliers.
16 // Pipeline Evaluation

The Established Typing Framework: Proven Impact & Opportunities for Growth

A pioneering national surveillance standard with clear opportunities to expand statistical sensitivity in complex data regimes.

>9,000 Isolates Processed Pioneering Surveillance Standard Opportunities for Statistical Refinement
Operational & Surveillance Strengths
  • Established the first national molecular typing standard, operational across CDC and state public health laboratories since 2018.
  • Successfully resolved major multi-state outbreaks (2018 salads, 2019 basil) to guide produce tracebacks and recalls.
  • Provides robust classification for high-depth, single-clone infections matching the curated reference catalog.
Opportunities for Statistical Refinement
  • Polyclonal co-infections and low parasite depth encounter distance inflation under heuristic mismatch penalties.
  • Novel single-nucleotide mutations are uncalled under rigid catalog-based string matching.
  • Tied pairwise distances introduce run-to-run tree variance under random tie-breaking.
  • Batch-dependent scaling can shift pairwise distances as surveillance databases grow across seasons.
Statistical variant modeling and population-weighted distance metrics help preserve cluster coherence across polyclonally mixed, low-depth, and evolving clinical isolates.
17 // The Modernized Framework

Modernizing Surveillance: The BRC-Analytics Galaxy Platform

A modular, containerized, open-source workflow for automated public health surveillance.

Push-Button Web UI (Galaxy) Apptainer Containerized Citable Zenodo Reference Trio
Core Platform Capabilities
  • Provides an accessible, web-based Galaxy interface so public health scientists can run complete typing pipelines without command-line overhead.
  • Replaces dictionary lock-in with standard BWA-MEM read mapping and LoFreq statistical variant calling.
  • Automates batch scaling: processes 1 to 1,000+ specimens in parallel with zero manual interventions.
Standardized Deliverables
  • Outputs standard CDC HDS matrices fully backwards-compatible with existing surveillance archives.
  • Quantifies continuous intra-host variant frequencies (down to 5%), preventing low-abundance secondary strains from distorting pairwise distances.
  • Guarantees end-to-end data provenance by recording all tool parameters, seeds, and container hashes.
Packaged as a push-button Galaxy workflow chaining 4 custom tools into a single container pinned to PyEuk v2.1.0 (binary KING-wIBS engine).
18 // Pipeline Comparison

Direct Comparison: Existing Pipeline vs. Modern Engine

Stage-by-stage architecture: comparing current operational practices with modernized statistical workflows.

>300× Faster (1.56s vs. 500s) 100% Deterministic (No Random Ties) Full Cohort Retention (100% Precision) Unsupervised (Label-Free ARI = 0.9737)
Existing CDC Pipeline
1. Catalog-Based String Matching

Matches reads against a predefined catalog of reference alleles. Novel SNVs or low-depth loci remain uncalled, leaving missing markers and reducing usable cohort size.

2. Integer Heuristic Distance

Calculates distances using discrete mismatch penalties on co-infections (\( w = 1 + x \)). Single-threaded R execution completes in 370–655 seconds.

3. Heuristic Clustering & Thresholds

Hierarchical clustering breaks tied distances with random selection (ties.method = "random"); cutoffs are calibrated using known historical reference labels.

Modern BRC Pipeline
1. Statistical Variant Calling

BWA-MEM + LoFreq statistical models detect both catalog alleles and novel variants down to 5% frequency, typing all 27 low-depth markers with 100% precision.

2. Population Weighted KING-wIBS

Inverse-variance allele weights (\( w_j = 1/\sqrt{p(1-p)} \)) quantify continuous allele sharing. Vectorized Numba engine executes in 1.56 seconds (>300× speedup).

3. Deterministic Ward & Unsupervised Cut

Lexicographical tie-breaking ensures 100% deterministic trees. Dynamic relative-gap cut discovers natural clusters (\( k=2, \text{ARI}=0.9737 \)) without requiring answer keys.

SPEED
1.56s (vs. ~500s)
STABILITY
100% Deterministic
SAMPLE RETENTION
100% (vs. 44% Kept)
DISCOVERY MODE
Unsupervised Dynamic Cut
19 // Distance Metric Intuition

How KING-wIBS Works: Population Weighting & Co-Infections

Weighting shared alleles by background rarity separates true outbreak clusters from population noise.

Inverse-Variance Weighting Polyclonal Robustness Continuous [0, 1] Metric
1. Evidentiary Weight of Alleles

Borrowing from the KING kinship framework, markers are weighted inversely by their binomial standard deviation across the surveillance cohort:

\[ w_j = \frac{1}{\sqrt{p_j (1 - p_j)}} \]
A rare outbreak allele (\( p = 0.02 \)) receives high weight (\( w = 7.14 \)), providing strong evidence of a shared food vehicle.
A common background allele (\( p = 0.90 \)) receives lower weight (\( w = 3.33 \)), preventing ubiquitous markers from creating false clusters.
2. Handling Mixed Intra-Host Infections

Because 92.2% of patient isolates (141 / 153) carry multi-allelic mixtures (\( \text{MOI} \ge 2 \)), KING-wIBS evaluates continuous overlap without quadratic penalties:

\[ D_{\text{wIBS}}(i, j) = \frac{\sum_{k=1}^{M} w_k \cdot |X_{ik} - X_{jk}|}{\sum_{k=1}^{M} w_k} \]
  • Measures continuous allele overlap across loci without adding artificial penalties for secondary background strains.
  • Guarantees distances remain strictly bounded between 0 (identical) and 1 (fully distinct).
  • Vectorized Numba kernel evaluates all 580,000 pairwise distances across a cohort in under 50 milliseconds.
By weighting shared alleles by population frequency and avoiding quadratic penalties on co-infections, KING-wIBS measures true biological relatedness across patient isolates.
20 // Clustering Intuition

Unsupervised Outbreak Detection: Dynamic Relative-Gap Tree Cutting

How PyEuk discovers the true number of outbreak clusters without requiring historical answer keys.

Label-Free Discovery Dynamic Tree Gap Metric Identifies True Outbreak Boundaries
Dendrogram Merge Jump at Outbreak Boundary
0.00 0.05 0.15 0.25 Merge Height (h) Large Merge Gap (Δh) Maximum Relative Gap at k = 2 Cluster 1 (Salad mix) Cluster 2 (Veg tray) Sporadic Case

Within-outbreak cases merge at low heights (\( h \lt 0.04 \)), followed by an abrupt vertical jump (\( \Delta h \)) before merging with unrelated isolates.

The Relative-Gap Selection Rule
\[ \text{RelativeGap}(k) = \frac{h_k - h_{k-1}}{h_{\max} - h_{\min}} \]
  • Evaluates the vertical height difference between consecutive dendrogram merges across all candidate cluster counts \( k \).
  • Selects the partition at the global maximum relative gap, identifying the transition where distinct outbreaks separate.
  • Successfully identifies the two-outbreak structure (\( k = 2, \text{ARI} = 0.9737 \)) in verified benchmark data.
  • Supports automated surveillance by identifying outbreak boundaries without needing subjective manual threshold calibration.
Dynamic relative-gap tree cutting reduces the need for fixed distance thresholds, supporting public health agencies in delineating emerging clusters through automated, data-driven analysis.
21 // Real-World Impact

Preventing Outbreak Data Loss: Rescuing Markers & Retaining Patients

How statistical variant calling and pairwise-complete metrics help maximize the use of clinical sequencing data.

27 Blank Markers Recovered Improving Data Utility 100% Clinical Cases Retained (153 / 153)
Data Challenges in Complete-Case Workflows
  • Low oocyst shedding in clinical stool extracts creates variable sequencing depth across the 8 amplicon loci.
  • Rigid string-matching and strict depth cutoffs leave low-coverage loci uncalled (27 blank markers in the 2018 benchmark cohort).
  • Complete-case filtering drops any specimen with even a single missing locus, resulting in 56% cohort attrition (67 / 153 cases kept).
The Modern Statistical Rescue Architecture
  • BWA-MEM read mapping and LoFreq statistical models recover all 27 uncalled markers with 100% precision.
  • Sensitive calling expands nuclear locus coverage, providing complete allelic profiles across clinical stool samples.
  • Pairwise-complete distance calculations retain 100% of clinical cases (153 / 153) without dropping partial profiles.
By improving variant calling sensitivity and employing pairwise-complete clustering, we ensure maximum data utility, providing comprehensive statistical power for public health investigations.
22 // Benchmark Methodology

Benchmarking the Pipeline: Setup & Evaluation

A comparative approach to evaluate the performance of sequencing and clustering components on clinical data.

Controlled 4-Arm Design 2018 Outbreak Evaluation Comprehensive Case Inclusion
1. Why We Benchmarked

Evaluating an end-to-end pipeline often obscures individual component performance.

Our 4-arm matrix isolates variant calling from distance/clustering metrics to assess each component independently.

2. Addressing Data Gaps
  • Historical benchmarks were sometimes limited by incomplete cohort depositions and rigid filtering strategies.
  • Strict coverage cutoffs caused marker dropouts across low-yield clinical extracts.
  • These gaps reduced the number of patient isolates retained in final epidemiological tracebacks.
3. Enhancements with PyEuk
  • BWA-MEM mapping and LoFreq calling recover missing marker data with 100% precision.
  • Pairwise-complete distance calculations integrate all available sequence data points.
  • Metric repair techniques support consistent clustering, enabling full cohort retention.
By improving data recovery and utilizing pairwise-complete metrics, PyEuk provides a more comprehensive and inclusive approach for outbreak surveillance.
23 // Head-to-Head Results

Benchmark Comparison: Conventional Methods vs. PyEuk

Evaluating supervised heuristics calibrated on known labels against PyEuk's fully unsupervised discovery engine.

Conventional: Calibrated on Known Outbreak Labels PyEuk: Fully Unsupervised (Zero Labels) Full Cohort Retention (203 vs 67 Cases)
Experimental Matrix & Label Awareness
  • Arm 1 (Published Calls + PyEuk): Evaluates PyEuk's metric on CDC data (\( N = 153 \)). Discovers 2 core clusters (\( \text{ARI} = 0.9721 \)) with zero label guidance.
  • Arm 2 (Legacy CDC Baseline): Retains only complete 8-locus profiles (\( N = 67 \)). Reaches \( \text{ARI} = 1.000 \) because distance cutoffs were calibrated a posteriori on known truth labels.
  • Arm 3 (LoFreq Variant Calling + Legacy R): Rescues missing markers (\( N = 147 \)), but still depends on supervised calibration against outbreak keys.
  • Arm 4 (Integrated PyEuk): Evaluates all \( N = 203 \) cases (100% retention). Fully unsupervised relative-gap cutting achieves \( \text{ARI} = 0.8894 \) with zero prior labels.
Quantitative Performance Matrix 2018 Outbreak Cohort
Arm Workflow Components Cases (N) Cluster ARI Supervision / Calibration
Arm 1 Existing Calls + PyEuk wIBS 153 0.9721 Unsupervised (0 Labels)
Arm 2 Existing Calls + Legacy R 67 1.0000 Calibrated on Known Labels
Arm 3 Raw Reads + Legacy R 147 0.9975 Calibrated on Known Labels
Arm 4 Integrated PyEuk (Full Cohort) 203 0.8894 Unsupervised (0 Labels)
Cluster ARI — the adjusted Rand index, scoring a grouping against the Vendor A / Vendor B outbreak labels. It works over every pair of specimens: how often do the two groupings agree, either putting a pair together or keeping it apart? 1.0 is exact agreement, 0 is chance, negative is worse than chance. The adjustment is what makes it readable — two random clusters of similar size already agree about half of all pairs, so an unadjusted score of 0.5 would mean nothing.
Conventional CDC pipelines achieved perfect ARI by fitting distance cutoffs directly to known outbreak labels on a reduced 67-case subset. PyEuk performs true unsupervised discovery directly from raw reads across the entire 203-case cohort without ever seeing epidemiological labels.
24 // Outbreak Resolution

Outbreak Resolution: Comprehensive Dendrogram & Markers

PyEuk retains the complete patient cohort, enabling clear separation of outbreak signatures.

Full Patient Retention Unsupervised Partitioning Distinctive Genetic Markers
Comprehensive Transmission Dendrogram
0.00 0.05 0.12 0.22 Merge Height (h) Dynamic Relative Cut Outbreak A (n=99) Outbreak B (n=104)

The PyEuk workflow retains the full clinical cohort, providing clear partitioning of distinct transmission events.

Diagnostic Genomic Signatures (Alleles by Outbreak)
Mitochondrial Repeat Junction (Mt_Junction) 100% Binary Split
Outbreak A (Salad): Mt_Cmt199.A_Hap_17 (199 bp)
Outbreak B (Veg Tray): Mt_Cmt169.A_Hap_8 (169 bp)
Nuclear Microsatellite (Nu_378 PART D) Diagnostic Repeats
Outbreak A (Salad): Nu_378_PART_D_Hap_7 (94.5%)
Outbreak B (Veg Tray): Nu_378_PART_D_Hap_2 (100.0%)
Nuclear Coding Loci (Nu_CDS1 & Nu_CDS4) Mutually Exclusive
Outbreak A (Salad): CDS1/4: PART_A_Hap_2 / PART_B_Hap_2
Outbreak B (Veg Tray): CDS1/4: PART_A_Hap_1 / PART_B_Hap_1
The integration of these diagnostic genomic markers supports precise clustering of the patient cohort, facilitating effective outbreak investigation.
25 // Sequence Search Engine

Petabase-Scale Pathogen Search: LexicMap on Galaxy

Democratizing ultra-fast, petabyte-scale k-mer querying across global public sequence archives.

Petabase Indexing (Logan / LexicMap) High-Memory Cluster Compute Push-Button Galaxy Web Ecosystem
The Computational Engine: Uniquely Powerful but Resource-Hungry
  • Alignment-Free Petabase Indexing: Uses compressed lexicographical k-mer tables to query target sequences across millions of raw SRA/ENA runs in seconds without downloading petabytes of FASTQ files.
  • Universal Pathogen Discovery: Instantly scans global repositories for any diagnostic gene, 28S rRNA marker, or viral sequence across eukaryotic parasites (Cyclospora, Babesia, T. cruzi).
  • Heavy Resource Demands: Maintaining and traversing petabase graph indexes requires terabyte-scale RAM, high-speed NVMe scratch arrays, and massive multi-core compute—far beyond standard lab hardware.
Galaxy: Democratizing Big-Data Public Health Surveillance
  • Shared Enterprise Infrastructure: Galaxy centrally hosts and synchronizes pre-computed petabyte sequence indexes on high-performance compute clusters.
  • Zero Hardware Overhead: Public health epidemiologists run multi-terabyte LexicMap queries via an intuitive web interface or API without managing supercomputers.
  • End-to-End Analysis: Automatically pipes discovered SRA accessions directly into downstream Galaxy tools (LoFreq, PyEuk) for seamless variant calling and phylogenetic placement.
By hosting LexicMap's resource-intensive indexing engine on Galaxy's scalable HPC infrastructure, global public health surveillance transitions from local PCR-constrained testing to instant, worldwide digital pathogen discovery.
26 // Metagenomic Mining

Mining Public Archives: Uncovering "Invisible" Global Incursions

Querying 28S rRNA across global stool archives exposes undetected clinical cases in general surveillance cohorts.

Global Diarrheal Surveillance 2 in ~1,000 UK GI Cohort Cases Detected WGS Metagenomes & RNA Metatranscriptomes
Empirical Discovery: 3 Unannotated Clinical Accessions
Bangladesh Cholera Cohort (PRJNA976726) SRR25011076

Acute diarrheal gut metagenome (WGS, ~14.3M reads) from Dhaka cholera surveillance harboring unannotated Cyclospora DNA.

UK Gastroenteritis Cohort (PRJEB62473) ERR11474981

Metatranscriptome (~28.3M reads) from an unresolved clinical gastroenteritis case where standard diagnostic testing found no pathogen.

UK Bacterial Co-Infection (PRJEB62473) ERR11495252

Metatranscriptome from a patient with confirmed Salmonella infection, revealing an unrecognized co-infection.

Epidemiological & Biosurveillance Insights
  • Hidden Community Prevalence: Demonstrates that ~2 in 1,000 patients in routine UK gastroenteritis surveillance harbored Cyclospora that went undetected because panels only tested for bacteria and viruses.
  • Co-Infection Dynamics: Highlights that bacterial positivity (e.g. Salmonella) frequently masks underlying parasitic infections in clinical settings.
  • Metatranscriptomic Sensitivity: High cellular abundance of 28S ribosomal RNA makes total RNA-seq highly sensitive for digital screening even when genomic DNA is scarce.
LexicMap proves that public sequence archives already contain the missing links of global parasite transmission—unlocking opportunistic surveillance from routine clinical sequencing without specialized wet-lab assays.
27 // Forensic Validation

From Shotgun Metagenomes to MLST: Empirical Allele Recovery

Validating diagnostic markers, single-strain MOI sampling, and geographic divergence from raw metagenomic reads.

MAPQ 60 Read Mapping 5 Diagnostic Sub-Locus Markers Called Monoclonal MOI = 1 Baseline
Forensic Read Alignment & Diagnostic Markers (SRR25011076)

Filtered mapping recovered 4 reads ($MAPQ = 60$, 100% full-length identity) yielding 5 unique sub-locus allele calls across the CDC reference panel:

Nu_360i2 (Nuclear Intron) PART_D_Hap_2 / PART_E_Hap_2

Read 1338290 spans Pos 334–485 (151 bp, $E = 2\times 10^{-43}$), matching CDC outbreak haplotypes with 0 mismatches.

Mt_MSR (Mitochondrial rRNA) PART_A_Hap_1 / PART_B_Hap_1 / PART_F_Hap_2

3 reads spanning Pos 35–686 ($E = 1\times 10^{-55}$). Overlapping reads 1360627 and 1360711 confirm 100% sequence identity.

Biological Findings & Typing Implications
  • Monoclonal Sampling (Observed MOI = 1): Zero heterozygous sites across all called loci confirms that un-amplified shotgun WGS at low parasite depth (<0.6%) samples individual single parasite lineages.
  • Geographic Lineage Divergence: Mt_MSR_PART_F_Hap_2 clearly separates this South Asian isolate from North American domestic references (PART_F_Hap_1), while shared Nu_360i2 alleles confirm conserved ancestral markers.
  • DNA MLST vs. RNA Schemes: UK metatranscriptomes yielded 28S rRNA hits in LexicMap but 0 reads at DNA amplicon loci—proving future surveillance must pair DNA MLST with 28S/18S ribosomal subtyping for RNA-seq.
Shotgun metagenomic reads can be directly resolved into high-confidence MLST alleles, bridging the gap between big-data sequence search (LexicMap) and rigorous epidemiological typing (PyEuk).
28 // Long-Read Sequencing

Long-Read Sequencing: Phasing & Structural Resolution

Oxford Nanopore reads span multi-window loci and tandem repeats without statistical imputation.

Physical Cis-Linkage Tandem Repeat Phasing Direct 8-Locus MLST
1. Physical Intra-Locus Phasing
  • Single reads span all 6 sub-windows of Nu_360i2 (610 bp) on one molecule.
  • Establishes true physical cis-linkage, eliminating statistical phasing errors across intronic regions.
2. Tandem Repeat Junctions
  • Reads bridge flanking anchors and the entire 15-mer repeat interior (Mt_Cmt).
  • Unambiguously distinguishes structural alleles (139 bp, 154 bp, 169 bp, 199 bp) that cause short-read clipping.
3. Genome Contiguity & MOI
  • Contig N50 >500 kb resolves the complete 6.38 kb circular mitochondrial genome (CM037477.1).
  • Zero intra-locus heterozygosity confirms single-clone (\( \text{MOI} = 1 \)) infections without deconvolution artifacts.
Long reads capture full-length amplicon targets and structural repeats in continuous molecules, replacing statistical imputation with direct observation.
29 // Long-Read Case Study

Long-Read Case Study: Phasing & Outbreak Placement

Typing a Canadian clinical isolate (Can-NML:CYC2020-001) and placing it into the US outbreak tree.

8 / 8 Loci Called (100% Identity) 25 Phased Markers Salad Mix Macro-Lineage
Complete 8-Locus MLST Profile MOI = 1
Target Contig Called Haplotypes Identity
Mt_MSR Circular Mt PART_A_H3 · B_H1 · C_H1 · D_H1 · E_H1 · F_H2 100%
Mt_Junction Circular Mt Mt_Cmt154.B_Junction_Hap_4 (154 bp repeat) 100%
Nu_360i2 Contig 46 PART_A–F: all Haplotype 1 (phased cis-array) 100%
Nu_378 Contig 73 PART_A_H2 · B_H1 · C_H1 · D_H3 100%
Nu_CDS1–4 Nuclear CDS1: A_H2/B_H2 · CDS2: A_H1/B_H1 · CDS3: A_H1/B_H1 · CDS4: A_H2/B_H2 100%
PyEuk Phylogenetic Placement
  • Salad Mix Macro-Lineage: PyEuk placed the isolate into the Vendor A (Salad Mix) clade, closest to benchmark case C_IL119_18 (\( D = 0.0136 \)).
  • Allelic Divergence: Distinguished from core 2018 cases by carrying the 154 bp junction (Mt_Cmt154.B) instead of the dominant 139 bp variant.
  • Cross-Platform Typing: PyEuk seamlessly integrates long-read WGS, short-read amplicons, and metagenomes into a single distance matrix in <500 ms.
PyEuk natively places whole-genome long-read assemblies alongside short-read surveillance archives without pipeline modifications.
30 // Implementation & Availability

Software Architecture & Data Availability

All pipeline components are open-source, containerized, and tied to versioned reference datasets.

Open Source (GitHub) OCI / Docker / Apptainer Zenodo Archive (10.5281/zenodo.21924355)
1. Galaxy Tool Suite
  • Tool wrappers published in nekrut/brc-tools following standard IWC workflow specifications.
  • Integrates BWA-MEM read alignment, LoFreq variant calling, HDS extraction, and PyEuk clustering into a continuous workflow.
  • Executable via web browser UI, Galaxy API, or CLI runners on workstations, institutional HPC clusters, and cloud nodes.
2. Containerized Engine
  • Single pinned container image containing Python 3.11, C-extensions, LoFreq, and BWA-MEM dependencies.
  • Implements PyEuk v2.1.0 with vectorized Numba KING-wIBS kernels for high-speed distance computation (1.56s).
  • Fixed dependency locking prevents runtime errors from environment drift or R package version conflicts.
3. Reference Package
  • Permanent archive at 10.5281/zenodo.21924355.
  • Packages markers.fa, parts.bed, haplotypes78.fa, and junction.fa with SHA256 checksums.
  • Ensures exact data provenance and reproducible allele calling across labs and surveillance seasons.
All software, workflow definitions, reference assets, and benchmark datasets are publicly accessible, versioned, and documented for independent validation and deployment.
31 // The Deliverable

The Workflow, As It Ships

Eleven steps in the Galaxy editor: four inputs, two callers running in parallel, one sheet, one clustering. Editable, versioned, and runnable by anyone with the link.

The Cyclospora typing (LoFreq/PyEuk arm) workflow open in the Galaxy workflow editor
Soon available from the BRC-Analytics “Workflows” tab. WE NEED BETA-TESTERS!!!!!
32 // Summary & Impact

Summary: Modernizing Cyclospora Genomic Surveillance

Four core computational advancements supporting sensitive, deterministic, and reproducible outbreak surveillance.

100% Patient Retention Deterministic & Mathematically Sound Unsupervised Outbreak Discovery Push-Button Galaxy Ecosystem
1. Full Sample Retention (No Discarded Cases)

Statistical calling rescues uncalled markers with 100% precision, while pairwise-complete metrics eliminate the legacy pipeline's 56% cohort attrition—ensuring every clinical case is retained in the investigation.

2. Deterministic & Reproducible Topologies

Replaces random coin-flip tie breaking with stable sorting and models co-infections with inverse-variance KING-wIBS weights—delivering reproducible tree topologies across independent runs.

3. Real-Time Outbreak Discovery (Label-Free)

Dynamic relative-gap tree cutting discovers true outbreak boundaries (k = 2, ARI = 0.9737) directly from the data on day one, eliminating reliance on pre-calibrated answer keys for emerging strains.

4. High-Performance Open Infrastructure

Vectorized Numba kernels evaluate 580,000 distances in 1.56 seconds (>300× faster than R). Packaged in pinned Apptainer/Docker containers with open Galaxy wrappers ready for routine public health lab deployment.

By combining statistical variant calling with population-weighted distance metrics, PyEuk provides public health agencies with an open, containerized workflow for fast, reproducible foodborne outbreak response.