Primery: turning BLAST mismatches into species-specific primers

Manuscript in preparation

Primer-BLAST checks a primer you have already designed. Primer3 places a primer without knowing how the target differs from its relatives. Primery fills the step between them: locate the discriminatory positions, score them, and constrain design to the regions that score well.

Pipeline and parameters

Primery runs as a PySide6 desktop application in Python 3.12, with the computational core separated from the interface so the same stages run from a script. Each stage writes a structured artefact that the next stage reads.

Stage 1: retrieval

The downloader queries NCBI Nucleotide through Entrez for an organism name, restricts hits to 200–10,000 bp, and excludes whole-genome records. A heuristic parser groups the results by gene name from the sequence title, which is where the locus list for a run comes from. Output is one FASTA per locus plus a manifest.tsv holding accession, organism label, locus assignment and QC status.

Check the manifest before continuing. Mislabelled organism annotations in public records propagate through the whole run, and the cheapest place to catch them is here.

Stage 2: mismatch analysis

Each reference goes to remote BLASTN against core_nt with megablast enabled, HITLIST_SIZE = 200, EXPECT = 10.0 and a minimum query coverage of 0.5. Results come back as compressed JSON2 and parse locally.

Hits belonging to the target species drop out by taxonomy ID or scientific name. For each surviving non-target hit, the high-scoring pair projects onto query coordinates to build a query-anchored alignment. Positions outside the HSP carry N for uncovered; deletions relative to the query carry a gap character. That distinction matters: an uncovered position is missing evidence, and treating it as a match would inflate the apparent conservation of the locus.

At each query position i the analyser records three quantities: the count of non-target accessions carrying a substitution or deletion, the set of distinct non-target species those accessions belong to, and the fraction of accessions per species showing the mismatch.

Stage 3: clustering and scoring

A greedy pass groups positions into clusters. Two positions join the same cluster when they lie within L base pairs of each other, where L is the primer length for the chosen chemistry: 25 for PCR, 35 for RPA. Each cluster receives a score equal to the number of distinct non-target species with mismatches inside it, multiplied by the mean accession fraction across those species. Breadth times consistency.

The second quantity is local density, defined as the count of mismatch positions divided by cluster span in base pairs. Density predicts the free-energy penalty a primer laid over the cluster will impose on a non-target template. A companion note works through how to use it and what threshold to hold.

Stage 4: design

The designer enumerates cluster pairs whose outer separation falls inside the product-size range for the chemistry, then calls Primer3 constrained to those regions with the SantaLucia nearest-neighbour model and 1998 salt corrections. Candidate pairs rank by a composite of four terms: 10 points per non-target species discriminated by both primers together, 5 points per mismatch position covered, a 3′ proximity bonus, and a thermal penalty at half weight.

The 3′ bonus awards 15 points for a mismatch at the ultimate base, 5 at the penultimate, 5 at the antepenultimate and nothing further in. The thermal penalty sums the maximum hairpin melting temperature, the maximum self-dimer melting temperature and the pair-complementarity temperature. In the benchmark, a hairpin at 44.9 °C on one reverse primer reduced a pair from a potential 425.0 to 402.5, moving it from first to second rank. Read the penalty column before you read the rank.

Stage 5: validation

Each designed pair returns to BLAST through Biopython's NCBIWWW.qblast with word size 7 and EXPECT = 1000, tuned for short queries. The analyser reconstructs putative amplicons wherever forward and reverse primers align to the same subject on opposite strands within 50–3,000 bp, then labels the pair species-specific, genus-specific, non-specific or no-hits.

Scoring, as code

The two functions that decide where a primer goes. Cluster scoring runs on the mismatch table from stage 2; pair scoring runs on Primer3 output plus the clusters each primer covers:

def cluster_score(cluster, species_at, accession_fraction):
    """Distinct non-target species in the cluster, times mean consistency."""
    species = {s for pos in cluster for s in species_at[pos]}
    if not species:
        return 0.0
    mean_fraction = sum(accession_fraction[s] for s in species) / len(species)
    return len(species) * mean_fraction


THREE_PRIME_BONUS = {0: 15, 1: 5, 2: 5}   # distance from the 3' terminus


def pair_score(n_species, mismatch_positions, primer, thermo):
    """10 per species discriminated, 5 per mismatch covered, positional
    bonus, half-weight thermal penalty."""
    bonus = sum(
        THREE_PRIME_BONUS.get(len(primer) - 1 - i, 0)
        for i in mismatch_positions
    )
    penalty = (thermo.hairpin_tm
               + thermo.self_dimer_tm
               + thermo.pair_complementarity_tm)
    return 10 * n_species + 5 * len(mismatch_positions) + bonus - 0.5 * penalty

Stage 2 submission parameters, kept explicit because changing any of them changes which species end up in the exclusion set:

BLAST_PARAMS = dict(
    database="core_nt",
    program="blastn",
    megablast=True,
    hitlist_size=200,
    expect=10.0,
    min_query_coverage=0.5,
    format_type="JSON2",
)

# Stage 5 uses different parameters: short queries need a small word size
VALIDATION_PARAMS = dict(program="blastn", word_size=7, expect=1000)
AMPLICON_WINDOW = (50, 3000)   # bp between opposing primer alignments

Stack

Language

Python 3.12, with the computational core free of GUI imports

Interface

PySide6, dockable panels, bilingual labels, background threads via QThread

Sequence handling

Biopython for Entrez, sequence I/O and NCBIWWW.qblast

Thermodynamics

primer3-py, SantaLucia nearest-neighbour model with 1998 salt corrections

Remote services

NCBI Entrez, BLAST Common URL API, Primer-BLAST; requests for transport

Statistics

NumPy and SciPy for descriptive statistics and one-way ANOVA over cluster scores

Outputs

JSON, FASTA, XLSX per stage; PDF and DOCX reports with annotated sequence maps

Chemistry presets

Two consequences follow for design work. RPA primers at 30–36 nt span more mismatch positions than PCR primers at 18–25 nt, so a wide cluster that yields nothing useful under PCR constraints can still support an RPA assay. LAMP outer primers are short and depend on the 3′ bonus, so for LAMP, sort the cluster list by 3′-proximate mismatches. Cluster score is the wrong key there.

Benchmark

Four organisms across twenty-one loci covered two diagnostic domains with different phylogenetic structure. Fusarium holds over 300 species, and sister species inside the F. graminearum complex share more than 95% identity at common loci. Equine cyathostomin nematodes comprise more than fifty morphologically similar species that resist visual identification.

Reading the cluster profile before designing anything

The decay shape across ranked clusters tells you more about planning than the top score does. C. nassatus COX1 produced 25 clusters scoring from 17.9 down to 8.8, a shallow decline in which most clusters remain usable. F. culmorum RPB2 produced 34 clusters scoring from 8.44 down to 0.80, so the usable set stops after the first few ranks.

Use the profile as a go/no-go step. A shallow decay across many clusters means the locus supports several independent assay designs, which gives room to redesign after a wet failure. A steep decay means one shot, and a second locus belongs in the plan from the start.

When a locus returns zero clusters

Ten of the fourteen F. bulbicola loci produced no clusters. This is a database limitation. Secondary-metabolism genes of the FUM cluster and single-accession references have no sequenced relatives to discriminate against. Work through four checks in this order before discarding a locus.

• Count non-target hits in the BLAST result. Fewer than about ten distinct species means the locus cannot support a discrimination claim at any threshold.
• Confirm the reference is not a whole-genome record or a chimeric submission, which pushes query coverage below the 0.5 floor and drops most hits.
• Widen the length window and re-run retrieval. Locus names in sequence titles vary between submitters, and the heuristic parser misses variants.
• Substitute a locus with broader taxonomic sampling. HIS3, TEF1 and LAE1 yielded clusters; the FUM series did not.

Primer outcomes

Composite scores across all designed pairs ranged from 97.5 to 410.0. Melting temperatures for the top-ranked pairs fell between 59.8 and 61.2 °C against a 60 °C target, a mean absolute deviation of 0.7 °C. Product sizes ranged from 106 to 316 bp inside the 100–500 bp PCR window, which confirms that cluster-pair enumeration preserves amplicon geometry.

Validation outcomes

The top C. nassatus COX1 pair returned 29 target hits, each a C. nassatus accession, at a 160 bp product with zero mismatches in either primer and zero off-target amplicons.

The top F. culmorum RPB2 pair, holding the highest composite score in the benchmark, returned genus-specific. It amplified F. culmorum together with F. graminearum, F. asiaticum and F. meridionale at a 122 bp product with zero mismatches, while excluding all non-Fusarium taxa. RPB2 sequences inside the F. graminearum complex share more than 97% identity, so no primer pair anchored in that locus can separate those species. The classification step reports it as genus-specific, which is the correct answer.

A high composite score therefore promises the best discrimination the locus supports, which is sometimes a genus. Plan the assay around that outcome: a genus-specific amplicon followed by restriction digestion or sequencing still delivers species-level identification, at the cost of a second step.

Taking a pair to the bench

The report includes the thermodynamic parameters needed to set up the reaction.

• Set the initial annealing temperature 3–5 °C below the lower of the two reported T_m values, then run a gradient across ±5 °C. The cluster-based designs hold a target/off-target differential in the 6–8 °C range, so the gradient has a window to find.
• Include the three nearest non-target species from the cluster report as template controls. Those are the species the primers were scored against, and they are where cross-amplification will appear.
• Reject any pair whose thermal penalty exceeds roughly 40 before ordering oligos, even at a high composite rank. A 44.9 °C hairpin means the primer folds on itself at annealing temperature.
• For a genus-specific pair, size the amplicon so that a diagnostic restriction site inside it separates the complex members, or plan for Sanger sequencing of the product.

Cost

Remote NCBI operations dominate runtime. BLAST submission and retrieval averaged 8.2 ± 3.4 min. Mismatch analysis and clustering took 2.1 ± 0.8 s. Primer design took 0.4 ± 0.1 s. Primer-BLAST validation took 6.5 ± 2.8 min per primer pair. End to end that gives 15 to 45 min per locus depending on server load. Local computation accounts for under three seconds of it, so batch size depends on how many BLAST submissions the session can hold, not on CPU.

Limits

• Discrimination inherits the sampling of core_nt. A sparsely sequenced genus yields empty clusters, not weak ones. The pipeline cannot tell an absent species from a conserved one.
• Results are in silico throughout. Specificity classification predicts amplification and does not measure it.
• LAMP support covers outer primers; inner primers need a separate design step.
• The clustering pass is greedy, and the composite score is a heuristic that tracked in-silico specificity across four organisms. Four organisms demonstrate the method and do not validate it.

Reports export to PDF and DOCX with cluster composition tables, thermodynamic parameters, specificity results and an optional sequence map carrying a ruler, mismatch heatmap, cluster regions and primer arrows.

References

1. Ye J. et al. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. doi:10.1186/1471-2105-13-134
2. Untergasser A. et al. Primer3: new capabilities and interfaces. doi:10.1093/nar/gks596
3. Camacho C. et al. BLAST+: architecture and applications. doi:10.1186/1471-2105-10-421
4. Cock P. J. A. et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. doi:10.1093/bioinformatics/btp163
5. Sayers E. W. et al. Database resources of the National Center for Biotechnology Information. doi:10.1093/nar/gkab1112
6. O'Donnell K. et al. Phylogenetic analyses of RPB1 and RPB2 support a middle Cretaceous origin for a clade comprising all agriculturally and medically important fusaria. doi:10.1016/j.fgb.2012.09.002
7. Hoffmann M. et al. PriSeT: efficient de novo primer discovery. doi:10.1186/s12859-023-05175-6
8. SantaLucia J., Hicks D. The thermodynamics of DNA structural motifs. doi:10.1146/annurev.biophys.32.110601.141800