Large and repeat-rich genomes (CLI)#

Community experience from a ~21.7 Gbp conifer assembly (GitHub issue #34), plus CLI behaviour in v4.5.1 for query -k.

This page is for the command-line MFEprimer (mfeprimer index / spec / qc / search) when the background database is very large and/or highly repetitive (for example plant / conifer assemblies of tens of Gbp). The web Specificity Check uses curated databases and is a different workflow.

Why default k=9 can hurt#

MFEprimer seeds binding-site search from the primer 3′ k-mer. Chance hits scale roughly as genome size / 4k:

Rough genome size Expected chance hits at k=9
~3 Gbp (human) ~104
~21.7 Gbp ~8×104

On a short genome-templated primer (e.g. 16–17 bp KASP-style cores), reported Plus/Minus binding counts (PBN/MBN) can reach the same order. Then:

  • Amplicon tables become unstable near the Tm cutoff (-t): salt / ion changes can swing reported amplicon counts (e.g. 0 vs 30+) even though the seed-level burden (PBN/MBN) is unchanged.
  • A plain -b (full binding workup) is often not viable; use --bind-amp-only (see below).

Raising index -k (e.g. 13–15; maximum 15) is the main lever. Query -k alone cannot invent a longer seed than the index stores.

Auto-detect query -k from the index#

You no longer need to keep query -k in sync by hand:

Situation Behaviour
Omit -k Read k from each -d index (.primerqc.bin header; legacy: first .primerqc.fai length)
Set -k explicitly Must equal the index k for every -d, or the run errors (avoids silent 0 sites)
Multiple -d All part indexes must share the same k 
# Build with a genome-scale k
mfeprimer index -i genome.fa -k 13 -c 32 -f

# Query: omit -k (auto from index)
mfeprimer spec -i primers.fa -d genome.fa -o out.spec \
  --bind-amp-only -b -j

Choose k when indexing; at query time, omit -k unless you want an explicit match check.

Split the genome and pass multiple -d#

A single index of a 20+ Gbp FASTA may be impractical. Split the FASTA, index each part with the same -k, then pass every part with repeated -d. Amplicons and binding sites are aggregated across parts.

Sequence IDs in the graphical report:

  • One -d: bare names (e.g. chr11a) — easy to join to a whole-genome .fai.
  • More than one -d: names are prefixed with the database path (/path/parts/part.5.fa:chr11a) so hits from different parts stay distinguishable. Strip the DB: prefix if your downstream joins expect chromosome-only IDs.

--bind-amp-only#

On large assemblies, prefer:

mfeprimer spec ... --bind-amp-only -b

--bind-amp-only implies -b, but only computes detailed binding patterns inside predicted amplicons. Genome-wide hit counts (PBN/MBN) are still reported. A plain -b that expands every seed hit is usually too expensive. (Related: issue #33.)

How to read specificity (not just amplicon count)#

Prefer Why
PBN / MBN Seed-level binding burden; large values warn that the primer is shared across many loci
Sequence-first Does the predicted amplicon match the intended design? Marker name / coordinates are corroborating evidence
3′ mismatch from the alignment Read the binding alignment graphic; do not rely only on summarised mismatch counts
-t 0 as a diagnostic Separates “many candidates near the Tm threshold” from a true change in seed burden

Ion / oligo flags are --diva and --oligo (not --divalent / --dna). Also --mono, --dntp.

Suggested workflow#

1. Pick index -k from genome size (keep expected chance hits bounded)
2. Optionally split FASTA; index every part with the same -k
3. Run spec/qc with --bind-amp-only; omit query -k
4. Judge with PBN/MBN + intended product sequence + 3′ alignment
5. If amplicon counts swing with salt/-t, try -t 0
6. If using multiple -d, strip DB: prefixes before joining to a global .fai

See also#

Backward Human pan-genome check MP-Ref Forward