Skip to content

SymbioSeas/primeval

Repository files navigation

primeval: in silico PCR assay validation

primeval evaluates the sensitivity and specificity of PCR assays (including probe-based assays designed for ddPCR/qPCR) against a user-provided set of genome assemblies. For each assay, it identifies valid amplicons using BLAST-based primer alignment and reports detection calls, mismatch counts, and amplicon sizes per assembly.

Two tools in this repository

Tool Location Purpose
primeval repository root (workflow/, config/, …) Validate PCR assays in silico against a set of genome assemblies.
assay-design assay-design/ Identify clade-conserved / clade-specific orthologs from a Panaroo pangenome and extract representative sequences. This is the candidate-gene discovery step used to develop the Vpop assays. See assay-design/README.md.

Both tools share the single conda environment defined in environment.yaml.

Once installed (below), each tool has a command on your PATH: primeval and assay-design. Run either with --help.

Features

  • Evaluates primer and probe binding across a user-provided set of genome assemblies
  • Supports probe-based assays (hydrolysis probes, e.g., ddPCR/qPCR) and probe-free assays (SYBR/dsDNA)
  • IUPAC degenerate base support in all primer and probe sequences
  • Handles primer binding on either strand of an assembly
  • Configurable mismatch tolerances and 3′-exact match requirements
  • Outputs per-assay detection summaries, amplicon details, and detection heatmaps
  • Runs locally or on SLURM HPC clusters via Snakemake profiles

Requirements

conda install -c conda-forge ncbi-datasets-cli

System requirements

Disk and memory scale with the number and size of input assemblies. Estimates below are from the manuscript's Vibrionaceae runs (RefSeq assemblies average ~5 MB each). Raw BLAST output is the dominant transient cost; with the default keep_blast: false it is deleted as the run proceeds.

Component Per assembly Draft set (10,715) Complete set (927)
Input assemblies (.fna) ~5 MB 49 GB 4.3 GB
Cached BLAST DBs (resources/blast_db/) ~1.3 MB 14 GB ~1.2 GB
Raw BLAST output (transient) ~15 MB 158 GB ~14 GB
Reports + amplicons ~0.2 GB small
Peak disk (keep_blast: false) ~65 GB ~6 GB
Peak disk (keep_blast: true) ~220 GB ~20 GB
Peak RAM ~16 GB ~16 GB

Rule of thumb: budget roughly 6–7 MB of transient disk per assembly with keep_blast: false, or ~21 MB per assembly if retaining raw BLAST, plus ~16 GB RAM for the final aggregation step.

Installation

git clone https://github.com/SymbioSeas/primeval.git
cd primeval
conda env create -f environment.yaml
conda activate primeval

The pipeline runs inside this activated environment; the Snakemake profiles set use-conda: false so no per-rule environments are built. Creating the environment also installs three commands onto your PATH: primeval (the pipeline), assay-design (the companion tool), and download-assemblies (the assembly-download helper).

Quick start

1. Download assemblies

Use the download-assemblies command to download RefSeq assemblies for your taxon of interest:

download-assemblies -t "Vibrionaceae" -o assemblies/

This downloads all RefSeq assemblies (complete through contig level) for the specified taxon and writes a metadata.csv to the output directory. See Downloading assemblies for options and HPC usage.

2. Configure the pipeline

Copy the config template into your analysis directory (the folder holding your assemblies/, where results will be written) and edit the copied config to point at your assemblies and set detection thresholds:

cp /path/to/primeval/config/config.yaml ./config.yaml
assembly_dir: "assemblies"          # directory containing .fna files
metadata: "assemblies/metadata.csv" # metadata CSV from download step
results_dir: "results"              # all outputs written here
assay_table: "assay_table.csv"      # your assay definitions

max_primer_mismatches: 2            # mismatches allowed per primer
prime3_exact_nt: 3                  # 3′-terminal bases that must match exactly
max_probe_mismatches: 1             # mismatches allowed in probe
max_amplicon_size: 500              # maximum amplicon size (bp)

# group_by: metadata column(s) used to group assemblies in the reports (a
# detection matrix + heatmap is written per column). Optional: leave unset to
# auto-group by NCBI ANI-derived species (requires metadata from
# download-assemblies). See "Metadata format" below.
# group_by: "species"
# group_by: ["species", "phenotype"]

3. Prepare your assay table

Create a CSV file named assay_table.csv with one row per assay (see Assay table format).

4. Run

primeval is installed on your PATH (see Installation). Run it from any analysis directory containing your assemblies/ directory and an updated config.yaml file:

primeval --run-name Vpop

Results are written to results/Vpop_<date>/ (amplicons/, blast/, reports/, and run.log). Options:

Flag Default Description
--run-name NAME results Names the results directory: results/<NAME>_<date>/.
--directory DIR current dir Analysis directory (assemblies, config, outputs).
--configfile FILE <dir>/config.yaml Pipeline configuration.
--force off Reuse today's <NAME>_<date>/ and resume unfinished work.
--cores N 8 CPU cores.

Re-running the same name on the same day without --force creates <NAME>_<date>_2, _3, … so previous results are never overwritten. Anything after -- is passed straight to Snakemake (e.g. primeval --run-name Vpop -- -n for a dry run). Set keep_blast: true in config.yaml to retain the raw per-assembly BLAST output (see System requirements), or keep_logs: true to retain the per-assembly logs/ directory (both are deleted on success by default; run.log is always kept).

Advanced (direct Snakemake / SLURM): invoke the workflow directly, passing the config explicitly:

snakemake --configfile config/config.yaml --profile workflow/profiles/local \
  --config results_dir=results/Vpop_manual
# SLURM: swap in --profile workflow/profiles/slurm

Assay table format

The assay table is a CSV file with the following columns:

Column Required Description
assay Yes Unique assay name (used in all output files)
fwd Yes Forward primer sequence (5′-3′)
rev Yes Reverse primer sequence (5′-3′, same orientation as fwd - primeval handles RC internally)
probe Column yes, value no Probe sequence (5′-3′). The column must be present, but leave the value empty to declare a probe-free (SYBR/dsDNA-dye) assay. See Probe-free assays.
target_group No (column may be omitted) The metadata group this assay is designed to detect, as column:value (e.g. phenotype:protective, species:Vibrio mediterranei). A bare value (no colon) is matched against the primary group_by column. Drives assay_performance.csv.
target_gene No (column may be omitted) Free-text gene/target label. Not used in detection; carried through to assay_performance.csv if present.

Sequence notation:

  • Standard IUPAC ambiguity codes are supported (R, Y, S, W, K, M, B, D, H, V, N)
  • Modifications can be noted inline using /ModName/ or [ModName] notation; these are stripped before alignment (e.g., /56-FAM/ACGT[BHQ1]ACGT)

Which columns must exist: the assay, fwd, rev, and probe columns must all be present. The run fails with Assay table missing required columns if any is absent. Only probe may carry an empty value (that is how you declare a probe-free assay); assay, fwd, and rev need real values. The target_group and target_gene columns are optional and may be omitted entirely.

Extra columns: you may add any additional columns to assay_table.csv (e.g. reference, notes) to keep your work organized; primeval ignores them. The one exception is target_gene — if you include it, it is carried through to assay_performance.csv alongside fwd, rev, and probe to make your life easier.

Probe-free assays

Leave probe empty for SYBR/dsDNA-dye chemistry. primeval then requires only a valid amplicon to call a detection: no probe oligo is searched, and the Primer Only call is structurally unreachable for that assay. Any valid amplicon is Detected, and Not Detected means no valid amplicon was found. This mirrors the chemistry, where any double-stranded product fluoresces. Note the consequence for cross-assay comparison, described under Reading the numbers.

Example:

assay,probe,fwd,rev,target_group
Assay1,ACGGGACAAAAAGGATGGCGAGTAC,AGCCGAGCGTTACCAGC,CGAACGCAATGATTCTCTGAGC,species:Vibrio mediterranei
Assay2,,GCTACGCCCTCCATCATCC,GCGCGTGATTATCTGATAGC,

Assay1 is probe-based and scored against a target group. Assay2 has an empty probe (probe-free) and an empty target_group, so it is reported but not scored.


Metadata format

metadata.csv describes your assemblies. Only one column is required:

  • accession — must match the assembly filenames (<accession>.fna).

Every other column is optional and carried through to the per-assembly outputs. To group assemblies in the reports, set group_by in config.yaml to one or more metadata columns:

group_by: "species"                 # single grouping
group_by: ["species", "phenotype"]  # a detection matrix + heatmap per column

Grouping-column names must not contain whitespace (use e.g. isolation_source, not isolation source).

If you used download-assemblies, the generated metadata.csv carries NCBI ANI/BioSample fields and primeval auto-groups by ANI-derived species when group_by is unset. If you bring your own metadata.csv without those NCBI columns, you must set group_by. When you set group_by, assemblies missing from metadata.csv or with an empty cell in a grouping column are still analyzed and reported under Ungrouped. In the default ANI-auto mode (no group_by), only species-confident assemblies get their own group; assemblies that are only genus-resolvable, unclassified, or otherwise low-confidence fold into a single Unclassified (low confidence ANI) group. Each assembly's tier is still retained per-assembly as ani_confidence (High/Genus/Low) in detection_by_assembly.csv, and the run manifest reports the high/genus/low counts.

Because grouping columns are independent, an assay can be scored against a different resolution than the one used to lay out a matrix (e.g. group the report by species while scoring a nested clade assay with target_group: phenotype:protective).

Example minimal metadata:

accession,species,phenotype
GCF_000000001.1,Vibrio mediterranei,protective
GCF_000000002.1,Vibrio harveyi,

Detection thresholds

primeval reports a detection call per assay per assembly using thresholds set in config/config.yaml. The defaults reflect PCR biochemistry:

Parameter Default Rationale
max_primer_mismatches 2 A primer with one or two internal mismatches still typically primes efficiently; this tolerates strain-level SNPs while excluding poor binders. Counted IUPAC-aware (a degenerate base matches any of its represented bases).
prime3_exact_nt 3 Mismatches at the 3′ terminus inhibit polymerase extension, so the last several bases must match exactly regardless of max_primer_mismatches.
max_probe_mismatches 1 Used for probe-based assays only. Hydrolysis probes tolerate less mismatch than primers, so the default is stricter.
max_amplicon_size 500 Typical qPCR/ddPCR amplicons are ~70–200 bp; 500 bp captures valid products while rejecting spurious long-range primer pairings.

Detection calls:

  • Detected — valid amplicon found with probe contained within it (probe assays), or valid amplicon found (probe-free assays)
  • Primer Only — valid amplicon found but probe not detected within it
  • Not Detected — no valid amplicon found

Interpreting Primer Only: both primers bind and would amplify, but the probe site is diverged or absent. For a hydrolysis-probe (ddPCR/qPCR) assay this usually means no fluorescent signal despite amplification, so assay_performance.csv counts it as a non-detection; for SYBR/probe-free chemistry any valid amplicon is a detection. pct_detected_or_primer in the summaries lets you see both interpretations.

The BLAST search itself is run with deliberately permissive settings (evalue=1000, perc_identity=70, word_size=7, tuned for short oligo queries) so that no candidate binding site is missed; stringency is enforced downstream by the mismatch, 3′-exact, and amplicon-size filters above.


Outputs

All outputs are written to the run directory (results/<run-name>_<date>/):

results/<run-name>_<date>/
├── amplicons/
│   └── {accession}.csv           # per-assembly, one row per assay: detection call,
│                                  #   mismatch counts, amplicon sizes/contigs/positions,
│                                  #   and (if enabled) amplicon sequences
├── run.log                       # full run log
└── reports/
    ├── {column}_detection_matrix.csv   # one per group_by column: group × assay % detected
    ├── detection_summary_long.csv       # tidy long table: assay × grouping × group (incl. misses)
    ├── assay_summary.xlsx               # same as detection_summary_long, one worksheet per assay
    ├── assay_performance.csv            # per-assay sensitivity / specificity vs target_group
    ├── detection_by_assembly.csv        # one row per assembly: metadata + 0/1 per assay
    ├── run_manifest.txt                 # parameters, tool versions, checksums, input accounting
    ├── per_assay/                       # {assay}_results.csv — every assembly's call for one assay
    └── figures/
        └── {column}_detection_heatmap.pdf/png   # one per group_by column

Notes

  • {column}_detection_matrix.csv and its matching figures/{column}_detection_heatmap.* are written once per group_by column (or once as group_detection_matrix.csv / group_detection_heatmap.* in ANI-auto mode, since the auto-derived column is named group). See Metadata format.
  • detection_summary_long.csv (one tidy table) and assay_summary.xlsx (one worksheet per assay) are the same data in different layouts, not a duplicate file.
  • detection_summary_long.csv includes assay × grouping-column × group combinations an assay misses (pct_detected = 0), so you can see both what an assay detects and what it doesn't, across every group_by column.
  • assay_performance.csv scores each assay against its declared target_group. Because targets are column:value, a nested clade assay (phenotype:protective) and a species-wide assay (species:Vibrio mediterranei) are each scored against the right axis in one run. See Interpreting assay_performance.csv for the full column reference and how to read the numbers.
  • n_multi_amplicon / max_amplicons flag assemblies where an assay produces more than one product. Interpretation is target-dependent: for a single-copy target this can indicate off-target priming or overestimated abundance, whereas for a multi-copy target (e.g. 16S rRNA) multiple amplicons per genome are expected.
  • detection_by_assembly.csv joins the binary detection matrix (1 = Detected) with every column of the input metadata.csv, one row per assembly.
  • Raw BLAST output and the per-assembly logs/ directory are intermediate and are deleted once a run completes successfully. Set keep_blast: true / keep_logs: true in config.yaml to retain them; the top-level run.log is always kept.

Interpreting assay_performance.csv

For every assay carrying a target_group, primeval classifies each assembly in the run two ways and cross-tabulates them:

  • Truth — is this assembly in the assay's intended target group? (i.e. does metadata[target_column] equal the value in target_group?)
  • Call — did the assay produce a detection in this assembly? (detection_call is Detected)
Detected (in silico) Not detected
In target group tp (correct hit) fn (missed target)
Not in target group fp (off-target hit) tn (correct reject)

Columns

Column Meaning
assay Assay name, from assay_table.csv.
fwd, rev, probe Carried over verbatim from assay_table.csv (including any /ZEN/-style modification notation), so the results file is self-contained.
target_gene Carried over from assay_table.csv when you include that optional column.
target_group The target exactly as written in the assay table (e.g. phenotype:protective).
target_column The metadata column the target resolved to (e.g. phenotype). A bare value resolves to the primary group_by column.
n_target Assemblies in the target group (tp + fn).
n_nontarget Assemblies outside the target group (fp + tn).
tp True positives — target assemblies the assay detected.
fn False negatives — target assemblies the assay failed to detect (i.e. sensitivity gaps).
fp False positives — non-target assemblies the assay detected (i.e. cross-reactivity).
tn True negatives — non-target assemblies the assay correctly did not detect.
sensitivity 100 × tp / (tp + fn) — the % of intended targets the assay detects (assay inclusivity).
specificity 100 × tn / (tn + fp) — the % of non-targets the assay correctly rejects (assay exclusivity).
n_detected_total Every assembly this assay detected in the run (tp + fp).

Control/reference assays with a blank target_group still get a row, with blank metrics and only n_detected_total filled, so the file remains a complete roster of every assay. Any rate whose denominator is zero is left blank rather than reported as 0 — e.g. a target_group that matches no assemblies yields n_target = 0, a blank sensitivity, and a warning in run.log.

Reading the numbers

  • These are in silico predictions, not wet-lab performance. They report whether the primers and probe have acceptable binding sites in each genome under the configured detection thresholds. Primeval does NOT test amplification efficiency, Tm, secondary structure, or partitioning behaviour. Read them as the sequence-level expectation that empirical validation is tested against.
  • Every denominator is your input assembly set. specificity = 100 means "no off-target detections among the genomes you supplied," not "no off-targets exist." See Assay specificity validation.
  • Primer Only counts as a non-detection (see Detection thresholds), so a probe-based assay that amplifies off-target without probe binding is correctly not counted as a false positive.
  • Probe-based and probe-free assays are not directly comparable. A probe-free (SYBR/dsDNA-dye) assay scores a detection on two binding sites; a probe-based assay needs three. All else equal that gives probe-free assays systematically higher apparent sensitivity and lower apparent specificity. The probe column carried into this file tells you which is which (blank = probe-free), so compare like with like.

Downloading assemblies

download-assemblies is installed with the environment and runs from any directory.

Basic usage

download-assemblies -t "Taxon name" -o assemblies/

Options:

Flag Description Default
-t TAXON Taxon name or NCBI tax ID (required; repeatable — see below)
-o OUTDIR Output directory assemblies
-l LEVELS Assembly levels (comma-separated) complete,chromosome,scaffold,contig
-s SOURCE Assembly source: refseq, genbank, or all refseq
-e EMAIL NCBI e-mail (or set NCBI_EMAIL env var)
-k API_KEY NCBI API key for higher rate limits (or set NCBI_API_KEY env var)

Multiple taxa

Pass -t more than once to download the de-duplicated union of several taxa into a single output directory:

download-assemblies \
    -t "Vibrio jasicida" \
    -t "Vibrio owensii" \
    -o assemblies/

Each taxon is queried in turn; assemblies shared between taxa (e.g. when one query nests inside another) are downloaded and written to metadata.csv only once. Mixing names and tax IDs is fine, e.g. -t "Vibrio owensii" -t 661487.

Setting your NCBI API key once

An NCBI API key raises your download rate limit from 3 to 10 requests/sec. This is worth setting if you're downloading large sets of assemblies (i.e., hundreds or thousands of assemblies). Rather than passing -k every time, save it once:

cp config/ncbi_credentials.example.sh config/ncbi_credentials.sh
# edit config/ncbi_credentials.sh and paste your key into NCBI_API_KEY

download-assemblies sources this file automatically on every run. To keep the key elsewhere, point PRIMEVAL_CREDENTIALS at your own file.

The key is resolved as: -k flag → NCBI_API_KEY environment variable → credentials file (first one set wins).

HPC / SLURM

Local available storage requirements for primeval are directly scaled by the assembly dataset provided (i.e., you need space to store the downloaded assemblies you provide primeval!). If needed, primeval runs can easily be submitted in a SLURM environment using the wrapper below.

Wrap the command in an sbatch job for large downloads:

sbatch --time=24:00:00 --mem=8G \
  --wrap="download-assemblies -t Vibrionaceae -o assemblies/"

download-assemblies is resume-aware: if interrupted, re-running it will skip assemblies already successfully downloaded.


Assay specificity validation

primeval tests assay sensitivity and specificity against your input assembly dataset. The scope of specificity testing is therefore determined by which assemblies you provide.

Recommended workflow for specificity screening:

  1. Single-primer BLAST screen (NCBI web interface): Individually BLAST each primer and probe sequence against the NCBI nt database, excluding your target taxon. This identifies any off-target binding sites outside your group of interest. If no hits are returned for any oligo, off-target amplification outside the taxon is extremely unlikely (a primer must bind for any amplicon to form).

  2. Expand the input dataset if needed: If step 1 returns hits in a non-target taxon, download assemblies from that taxon and add them to your assembly_dir. primeval will then determine whether those single-primer hits form complete, detectable amplicons.

This two-stage approach is computationally efficient, such that you only download and evaluate assemblies in taxa where off-target primer binding is possible.

Step 1 is a manual pre-screen performed through the NCBI web interface; it is not part of the reproducible primeval pipeline. The reproducible specificity evidence for a manuscript comes from primeval's assay_performance.csv over your assembled dataset (Step 2 onward).

In the future: A future release will support BLASTing directly against NCBI pre-built reference databases (e.g., ref_prok_rep_genomes) as a single-step broader specificity check, without requiring manual assembly downloads.


Test dataset

A small set of 5 assemblies is provided for validating your installation:

bash test_data/download_test_data.sh

Then update config/config.yaml:

assembly_dir: "test_data/assemblies"
metadata: "test_data/assemblies/metadata.csv"

And run from the repo root (using the repo's config directly):

primeval --run-name test --configfile config/config.yaml

The test dataset covers all detection scenarios: Detected (including via minus-strand primer binding), Primer Only, and Not Detected.

Example output

With group_by: "species" set and a target_group of species:Vibrio mediterranei declared for assay VmedA, the report format looks like this:

# species_detection_matrix.csv
group,n_assemblies,VhPath,VmedA
Vibrio harveyi,2,100.0,0.0
Vibrio mediterranei,3,0.0,66.67

# assay_performance.csv
assay,fwd,rev,probe,target_group,target_column,n_target,n_nontarget,tp,fn,fp,tn,sensitivity,specificity,n_detected_total
VmedA,GCTACGCCC…,GCGCGTGAT…,ACGACCTTC…,species:Vibrio mediterranei,species,3,2,2,1,0,2,66.67,100.0,2

assay-design — clade-specific target discovery

assay-design is the companion tool used to design the assays that primeval validates. It parses a Panaroo pangenome to find orthologs that are conserved within a clade and specific to it (absent elsewhere), then extracts representative protein/nucleotide sequences — the candidate targets from which the Vpop dPCR assays were built.

It installs its own command, assay-design, on your PATH:

assay-design \
    --matrix          gene_presence_absence.csv \
    --isolates-dir    isolate_groups \
    --gene-data       gene_data.csv \
    --representatives representatives.tsv \
    --output-dir      results

A tiny, worked example lives in assay-design/example/. For inputs, thresholds, and full usage, see the assay-design README.


Citation

If you use primeval, please cite the archived release (TBD!):

Smith S, et al. (2026) [manuscript title]. [journal]. doi:[doi] primeval [version] (2026). Zenodo. doi:[zenodo-doi]


License

MIT — see LICENSE.

About

Batch in silico PCR assay validation against genome assemblies. BLAST-based sensitivity/specificity screening for ddPCR/qPCR/SYBR primers and probes, with a companion tool for clade-specific target discovery.

Topics

Resources

License

Stars

0 stars

Watchers

0 watching

Forks

Packages

 
 
 

Contributors