What each step does, and the exact rule it applies. The step numbers match the subfolders of the output.

kmc -k99 -t<threads> -m<memory> -fm -ci<N> -cx<N × d>  @inclusion_list.txt  inclusion  work/
kmc -k99 -t<threads> -m<memory> -fm -ci1 -cx1000000000 @exclusion_list.txt  exclusion  work/
kmc_tools -t<threads> simple inclusion exclusion kmers_subtract inclusion_specific
kmc_tools -t<threads> transform inclusion_specific dump kmers.txt

KMC counts the kmers of all the files of a group together, so with N inclusion genomes a kmer that occurs at least N times (-ci N) is, in practice, a kmer every inclusion genome carries. -p/--min-inclusion lowers that threshold to a fraction of N, rounded up, for the regions that most but not all of the inclusion genomes share. When N × d passes 255, -cs N × d is added as well: KMC stores a count in one byte by default and applies -ci/-cx to that saturated value, so without it -ci 300 matches nothing at all — measured with kmc 3.2.4 on one kmer repeated 300 times, where -ci 250 reports a count of 255 and -ci 256 reports nothing. -cx N × d (-d, default 1) drops the kmers that occur more often than that, which are the repeated regions: they would give an assay that amplifies several places at once. The trade-off is that a kmer which happens to occur N times inside fewer genomes also passes this step — which is why step 4 checks the presence of every candidate in every inclusion genome with blast.

For the exclusion group, every kmer counts, however rare: kmers_subtract then removes from the inclusion set every kmer seen anywhere in the exclusion group. What is left is the inclusion-specific kmers, written as a fasta file (kmer_0, kmer_1, ...).

A shorter -k finds fewer regions, not more, which is the opposite of what one might expect and worth being clear about. A kmer survives the subtraction only if no exclusion genome holds it, so for a region that does exist in the exclusion genomes the kmer has to span a difference. A long kmer spans one more often, so more kmers survive and the specific stretches they assemble into are longer. A short kmer survives only sitting almost on top of a difference.

Measured on the Listeria panel, nine clonal complexes against the other 78:

All four were measured under the span rule of the time — a 22-base pair rather than the 25 of step 3 below — so the comparison between kmer sizes holds, but the absolute counts are two releases old: they run a little under what 1.7.0 gives, and 1.8.0's move from megablast to -task blastn moved them again. At -k 99 it reports 332 regions over those nine targets rather than 306.

Which assays are recovered does not change, so the kmer size buys candidates rather than answers. 99 is a good starting point for bacteria: it is near the top of the yield, and below it whole runs start coming back empty. 127 adds 4% more candidates and nothing else. A longer kmer is also more sensitive to assembly errors, since one wrong base spoils every kmer covering it.

skesa --cores <threads> --mem <memory> --fasta kmers.fasta --contigs_out assembly.fasta \
      --min_contig <kmer size>
# or
spades.py --s 1 kmers.fasta --isolate --only-assembler --threads <threads> --memory <memory> -o spades/
# or, with no external program at all
primer-finder ... -a unitig

The kmers overlap one another by k-1 bases, so assembling them turns thousands of kmers into a handful of contigs: the specific regions themselves. This is only a way to have fewer and longer sequences to examine — the sequences are not new information. SKESA is the default and is faster here; SPAdes (-a spades) sometimes joins regions that SKESA leaves apart. --only-assembler is used because kmers carry no quality values to correct.

-a unitig assembles the kmers without an assembler, and keeps more of them. SKESA and SPAdes are built for sequencing reads: they estimate coverage, treat what they see once as error, and reconcile disagreement between reads. None of that applies to a kmer set, where every kmer is an exact substring of real genomes and is present or absent with nothing in between — and a coverage model therefore throws the signal away. What is left once the noise handling goes is the de Bruijn graph on k-1-mers and a rule for walking it: extend while the step is forced, which means exactly one way on and exactly one way in. A stretch so built is one the kmers determine, never a guess between alternatives, and the result does not depend on the order the kmers arrive in.

Measured on the Brucella panel, where the two differ most:

On the bundled example it reports the same region as SKESA and 52 bases more of it. It is not the default: the two assemblers have years of use behind their output and this does not, so it is there to be chosen. Choose it when the specific stretches are short or the kmer count is small — and note that a kmer set says nothing about which strand it came from, so a stretch may be reported as the reverse complement of what another assembler would give.

A run does not have to ask for it to be rescued by it. An assembler built for reads can get nothing at all out of a kmer set, and the error that used to end such a run advised a smaller -k, which runs into the same assembler from the other end. So when the chosen assembler exits with an error, or assembles less than 0.1 bases per kmer, the kmers are assembled directly instead, with a warning, and run_info.json records assembler_used beside assembler.

Two kinds of failure are not answered this way, because they say the machine could not run the assembler rather than that the assembler could make nothing of the kmers, and the walk honours no memory budget at all: an assembler killed by a signal, which is what a kernel out-of-memory kill looks like, and SPAdes exiting 68 (its own --memory check) or 255 (GeneralError, which is what a SPAdes whose child was killed reports). Those end the run with the error. The threshold comes from 21 paired runs on three validation panels: SKESA assembled 0.000, 0.000 and 0.019 bases per kmer on the three where it returned nothing usable, and 0.702 to 1.388 on the other eighteen — a 37-fold gap with nothing in it (the record).

The fallback deliberately does not try to catch an assembler that is merely leaving a lot behind. On those same panels SKESA recovered between 10.7% and 69.8% of the bases the walk found while staying inside the healthy band, so a run cannot tell from its own numbers that it is losing most of the signal — and a rule that guessed would either fire on healthy runs or miss the worst case.

--min_contig is set to the kmer size, which is not the assembler's own default. Both assemblers stop at 200 bases if left alone, and a kmer is 99: an inclusion-specific stretch shorter than 200 bases would then never leave this step, and no later step could recover it — it would simply be missing, with nothing in the output to say so. Those short stretches are exactly what is left when the two groups are closely related. In the Listeria validation, where the groups are clonal complexes of one species, one run went from 67 assembled contigs to 398 and from 25 reported regions to 66, of which 41 are shorter than 200 bases (that run used -p 0.93, which CC121 needs). --min-region sets the floor to something else.

Only SKESA and the unitig walk can be asked for less. --min_contig tells SKESA what to report, and -a unitig builds every stretch the kmers determine and then drops the short ones, so its floor is exactly what is asked for and the shortest it can report is a single kmer. SPAdes has no such option and its own floor is in graph simplification, not in what it writes: assembled from every 99-mer of a synthetic stretch, stretches of 150, 180 and 199 bases are already missing from its simplified graph, while 210 and longer come through. So -a spades does not report the short regions however low the floor is set, and says so in a warning rather than quietly dropping them; the floor is applied to its output only so that it does not report contigs the run asked to leave out.

Lifting the floor only puts the short stretches in front of the later steps; it does not excuse them from anything. The 145-base stretch carrying a published Listeria probe now reaches the candidate stage and is then turned down by the exclusion check in step 5, which is a decision about the sequence rather than about its length.

minimap2 -t <threads> -a --eqx <one exclusion genome> assembly.fasta > mapping.sam

--eqx makes minimap2 write = and X in the cigar string instead of M, which is what lets primer-finder tell a match from a mismatch. The contigs are mapped to one exclusion genome: the first of the exclusion folder in alphabetical order, or the one given with -r. Only the primary alignment of each contig is read.

A contig is kept when, judged on that cigar string:

• it does not map at all — the whole contig is missing from that exclusion genome, the best case; or
• a mismatch run, an insertion or a deletion is longer than one base; or
• two differences (mismatch, insertion or deletion) are close enough to sit under one primer — the pair, counting both differences and the matches between them, spans at most 25 bases, which is the widest primer Primer3 may return. Two that touch count, as do two separated only by a deletion, since a deletion takes up no base of the contig; or
• an end of the contig is clipped over at least 21 bases, which means the aligner could not place a whole primer's worth of sequence anywhere in that genome.

Drawn, with a stretch of contig against the exclusion genome it mapped to — · where they agree, X a mismatch, I a base the contig has and the exclusion genome does not, - the other way round, and S an end the aligner could not place at all. Each of these is what is_candidate actually answers:

  kept    ··········XX··········          a run of two: one oligo covers both
  kept    ·····X·······················X···  23 between them, a span of 25: the widest primer covers both
  no      ·····X························X···  24 between them, a span of 26: no primer reaches both
  kept    ··········X-X·········          either side of a deleted base, which takes up none of the contig
  kept    ··········III·········          three bases the exclusion genome does not have
  kept    SSSSSSSSSSSSSSSSSSSSSSS·······  23 bases unplaceable: a whole primer's worth is simply absent
  no      ··········X···········          one mismatch, and nothing else near it
  no      ··········I···········          one inserted base, likewise
  no      SSSSSSSSSSSSSSSSSSSS··         20 clipped: not quite a primer's worth

The last three are the point of the rule. A contig that differs by a single base is specific — but no oligo placed on it discriminates well enough to be worth designing, so it is not a candidate.

Unless nothing else is left. A run where no contig meets that rule used to stop with an error. It now tries once more accepting any single difference, and reports what it finds marked last_resort in the header of every region, in run_info.json (last_resort: true, min_differences: 1), and in a warning from both find and design. The two answers are different claims and the mark is what keeps them apart: a region that only the second pass accepted rests on one mismatch, which is the kind of assay that is hardest to make work, so it is a lead to look at rather than a result to act on. --no-last-resort asks for the strict answer and the error. Nothing changes for a run that finds something under the strict rule: the second pass only runs when the first reported nothing at all, which also means no region it reports would have passed strictly.

The relaxation is only to the number of differences. A contig identical to the exclusion genome is still not a candidate, and the 21-base clipped end is unchanged — that is how much unplaceable sequence is enough, not how many differences there are.

25 bases is the longest primer Primer3 is asked for here, so two differences that close can sit in the same primer, which is what makes an assay selective. A single mismatch somewhere in the contig is not enough, because a primer carrying one mismatch still amplifies the exclusion template often enough to be useless.

The span and the clipped length are different quantities and are not tied together: the first is how far apart two differences may be, the second is how much unplaceable sequence is worth keeping on its own. The span follows the design step's widest primer, which is what it means for both to sit under one. Making it 25 rather than a span of 22 gained 8.5% more candidate regions on the Listeria panel and changed no verdict — the record has the numbers. The clipped length keeps the 21 it has always had, since raising it would only throw regions away.

primer-finder design is built to spend them that way: it aims its Primer3 requests at the stretch of at most one oligo's width holding the most differences, and by default sets aside an assay that covers fewer than two of them under a single primer. The probe does not count towards that, since it is longer and hotter than a primer and two mismatches under it are no evidence it will not bind. See Designing assays.

The kept contigs are written to best_kmers.fasta, the ones with the most differing bases first — which puts a contig the exclusion genome does not hold at all at the top — with those bases in lower case. Mismatched, inserted and clipped bases are marked; a deletion is in the cigar string only, since those bases are not in the contig. When a contig maps to the minus strand, the marks are put back on the forward sequence, which is the one written out.

Because only one exclusion genome is used here, the differences found may be particular to that genome. Step 5 is what checks them against the whole exclusion group.

makeblastdb -in <each inclusion genome> -dbtype nucl -out <output folder>/...
blastn -task blastn -db <genome> -query best_kmers.fasta -evalue 1e-10 -max_target_seqs 1 \
       -num_threads 1 -outfmt "6 qseqid evalue"

Each inclusion genome is searched for each candidate; a hit with an e-value of 1e-10 or better counts as present. The presence of every candidate in every genome is written to inclusion_blast_hits.tsv, and a candidate missing from even one inclusion genome is dropped: an assay must amplify all of them. This is where the kmer counting of step 1 is checked properly, since "occurs at least N times" is not quite "is in N genomes". With -p/--min-inclusion, a candidate has to reach the same number of genomes as the kmers did; the ones that survive say how many hold them (inclusion=11/12) and the widest are listed first. The genomes are searched in parallel (-t), and the blast databases are built inside the output folder, never next to the input genomes.

Gzipped genomes are decompressed into the output folder for this step: blast cannot read gzip.

blastn -task blastn -db <each exclusion genome> -query all_inclusion_contigs.fasta -evalue 1e-10 \
       -max_target_seqs 500 -num_threads 1 -outfmt "6 qseqid qstart qend evalue qseq sseq"

Every matching subject is asked for rather than the ten best: a draft assembly can put a short stretch of a contig on a sequence of its own, which at ten could be crowded out by the sequences holding the rest — and that would read as the genome lacking a stretch it has.

Every hit of a contig in an exclusion genome is read, not just the best one, and each gives the positions of the contig that differ from that copy: mismatches, bases the contig has and the genome does not, and the position where bases the genome has are missing from the contig. The aligned sequences blast returns are walked base by base, so the positions are positions in the contig, counted from 0.

A position counts as different in a genome only when every copy of that region in the genome differs there: one matching copy is enough for the assay to amplify that genome. A position is then kept when it differs in at least 90% of the exclusion genomes that align it — a few exclusion genomes carrying the inclusion allele do not spoil an assay, and the threshold leaves room for them. A genome whose hits do not reach that position is left out of the count rather than counted as identical: it does not hold that region, which only makes the assay more selective.

A stretch that no exclusion genome aligns is a different matter: the whole of it is missing from the group, which is the strongest evidence of specificity there is, so it counts as differing. It has to run to at least 40 bases — not because that is what a primer needs, but because of what this search can see. Measured against a 3 Mb genome at this step's own -evalue 1e-10, an exact island of 25 or 30 bases is not reported at all while one of 33 is, so a shorter run could be present and merely invisible. The shape this is for is a novel block sitting in shared backbone — a repeat-targeted marker, most often — where every exclusion genome matches the flanks either side and nothing in between (the record).

A contig is kept when:

• no exclusion genome hits it at all (it is specific on its own), in which case it is written as it came out of step 3; or
• more than one kept position remains, and two of them span at most 25 bases, as in step 3.

The survivors are final_kmers.fasta, with those positions listed in the header and in lower case in the sequence.

The last-resort pass of step 3 reaches this step too. At that stringency a single kept position anywhere is enough here, where two close together are normally wanted, and every region the pass reports is marked last_resort in its header — including one no exclusion genome hits at all, which got through step 3 on one difference even though step 5 has nothing to hold against it.

primer-finder design does this part: Primer3 on the regions, a check of what would make each assay selective, and the primer files for in silico PCR. It aims its requests at the differences a single oligo could cover — which is what step 3 selected the region for — and by default carries no further an assay that covers fewer than two of them under one primer (the probe does not count). See Designing assays, which also explains how the assays are ordered and why the total number of differences under the primers matters more than a single one at the 3' end.

• Perfect matches only. A kmer must be in every inclusion genome without a single mismatch, and in no exclusion genome. One bad assembly in either group can remove every candidate. Draft assemblies with missing regions are the usual reason a run finds nothing. -p/--min-inclusion relaxes the inclusion side of this, at the price of an assay that does not amplify every inclusion genome; nothing relaxes the exclusion side, where a single matching genome makes an assay useless.
• Assembled genomes only. Reads are not an input; assemble them first.
• One exclusion genome in step 3. The candidate list depends on which genome that is (-r). The differences are then re-checked against the whole group, so this affects what gets looked at, not what survives.
• find designs nothing. Its output is the region to design on, not a primer pair; design is the step that proposes assays, and neither command tests one. Only a laboratory does that.
• Plasmids and repeats. -d 1 discards anything repeated inside a genome; a kmer carried by a plasmid that moves between the two groups will not be specific.
• Where the output folder may be. It must be outside both input folders — they are searched recursively, so results written inside one would come back as input genomes — and its path must not contain a space, which BLAST cannot handle in a database path. The input folders themselves may contain spaces: every genome is linked into the output folder before blast sees it.