app/lib/primers.mjs

Functions

Imports

It imports no other module of this repo. Imported by app/kb/registry/align.mjs, app/kb/registry/catalog.mjs, app/kb/registry/reference.mjs.

Source

1// PCR primer arithmetic: reverse complements, primer sites on a template, and product sizes.
2//
3// Dustin's protocol rule: primer sequences and reverse complements are handled by tested code and
4// diffed against the cited paper, never worked out by hand on a page. Melting temperature is not
5// computed here: the registry computes one under stated conditions (app/kb/registry/tm.mjs), because
6// the papers do not all give salt, primer and Mg2+ concentrations, so no paper's figure is comparable.
7//
8// COORDINATES are 1-based and inclusive, the way GenBank numbers a record, so a site or product
9// reported here can be checked against the accession's own feature table without an off-by-one.
10//
11// IUPAC codes are honoured on both sides. A primer base and a template base MATCH when the sets of
12// nucleotides they stand for share at least one member, so an `N` in either matches anything.
13
14/** @type {Record<string, string>} each IUPAC code and the nucleotides it stands for */
15const IUPAC = {
16  A: "A",
17  C: "C",
18  G: "G",
19  T: "T",
20  R: "AG",
21  Y: "CT",
22  S: "CG",
23  W: "AT",
24  K: "GT",
25  M: "AC",
26  B: "CGT",
27  D: "AGT",
28  H: "ACT",
29  V: "ACG",
30  N: "ACGT",
31};
32
33/** @type {Record<string, string>} */
34const COMPLEMENT = {
35  A: "T",
36  T: "A",
37  C: "G",
38  G: "C",
39  R: "Y",
40  Y: "R",
41  S: "S",
42  W: "W",
43  K: "M",
44  M: "K",
45  B: "V",
46  V: "B",
47  D: "H",
48  H: "D",
49  N: "N",
50};
51
52/**
53 * Uppercase, whitespace removed, RNA `U` read as `T`. Anything that is not an IUPAC nucleotide code
54 * THROWS, naming the character and its position: a stray digit or dash from a pasted table would
55 * otherwise shift every coordinate after it without a word.
56 *
57 * @param {string} seq
58 * @returns {string}
59 */
60export function normalize(seq) {
61  const out = seq.replace(/\s+/g, "").toUpperCase().replace(/U/g, "T");
62  for (let i = 0; i < out.length; i += 1) {
63    if (!(out.charAt(i) in IUPAC)) {
64      throw new Error(`not an IUPAC nucleotide code: ${JSON.stringify(out.charAt(i))} at position ${i + 1}`);
65    }
66  }
67  return out;
68}
69
70/**
71 * @param {string} seq
72 * @returns {string} the reverse complement, IUPAC-aware (R <-> Y, K <-> M, B <-> V, D <-> H; S, W, N map to themselves)
73 */
74export function reverseComplement(seq) {
75  const s = normalize(seq);
76  let out = "";
77  for (let i = s.length - 1; i >= 0; i -= 1) out += COMPLEMENT[s.charAt(i)];
78  return out;
79}
80
81/**
82 * The first record of a FASTA file.
83 *
84 * @param {string} text
85 * @returns {{ header: string, sequence: string }}
86 */
87export function parseFasta(text) {
88  const lines = text.split(/\r?\n/);
89  const start = lines.findIndex((line) => line.startsWith(">"));
90  if (start === -1) throw new Error("no FASTA header line (a line starting with >)");
91  const body = [];
92  for (const line of lines.slice(start + 1)) {
93    if (line.startsWith(">")) break;
94    body.push(line);
95  }
96  return { header: (lines[start] ?? "").slice(1).trim(), sequence: normalize(body.join("")) };
97}
98
99/**
100 * @param {string} a one IUPAC code
101 * @param {string} b one IUPAC code
102 * @returns {boolean}
103 */
104function compatible(a, b) {
105  if (a === b) return true;
106  const set = IUPAC[b] ?? "";
107  for (const n of IUPAC[a] ?? "") if (set.includes(n)) return true;
108  return false;
109}
110
111/**
112 * @typedef {{ type: "extra-primer-base" | "missing-primer-base", primerPosition: number }} Indel
113 *   `extra-primer-base`: the primer carries a base the template lacks, at `primerPosition` (1-based
114 *   from the primer's 5′ end). `missing-primer-base`: the template carries a base the primer lacks,
115 *   immediately 3′ of primer position `primerPosition` (0 means before the first base).
116 */
117
118/**
119 * @typedef {object} Site
120 * @property {"+" | "-"} strand `+`: the primer's own sequence reads on the template as given, so it
121 *   primes synthesis rightward. `-`: its reverse complement reads there, so it primes leftward.
122 * @property {number} start 1-based first template base the primer covers
123 * @property {number} end 1-based last template base the primer covers, inclusive
124 * @property {number} mismatches substitutions, not counting an indel
125 * @property {number[]} mismatchPositions 1-based positions IN THE PRIMER, counted from its 5′ end
126 * @property {Indel | null} indel
127 * @property {string} templateSite the template bases under the primer, written 5′ to 3′ in the
128 *   primer's own orientation, so it lines up with the primer character for character when there is no indel
129 */
130
131/**
132 * Hamming comparison of a primer (in template orientation) against the template at `at`.
133 *
134 * @param {string} template
135 * @param {string} probe
136 * @param {number} at 0-based
137 * @param {number} limit give up past this many mismatches
138 * @returns {number[] | null} 0-based probe indices that mismatch, or null past the limit
139 */
140function hamming(template, probe, at, limit) {
141  /** @type {number[]} */
142  const miss = [];
143  for (let i = 0; i < probe.length; i += 1) {
144    if (!compatible(probe.charAt(i), template.charAt(at + i))) {
145      miss.push(i);
146      if (miss.length > limit) return null;
147    }
148  }
149  return miss;
150}
151
152/**
153 * Every place a primer binds a template within the allowed number of edits, on both strands.
154 *
155 * An INDEL is opt-in (`maxIndels: 1`) and at most one is considered: a primer with one extra base
156 * relative to the reference, or one base short of it. It is reported, never silently absorbed.
157 *
158 * Where two candidate sites on the same strand overlap, only the better is kept (fewer edits, and
159 * an ungapped fit over a gapped one on a tie), because one stretch of template is one binding site.
160 *
161 * @param {string} template
162 * @param {string} primer 5′ to 3′
163 * @param {{ maxMismatches?: number, maxIndels?: number }} [options]
164 * @returns {Site[]} sorted by start
165 */
166export function findSites(template, primer, options = {}) {
167  const maxMismatches = options.maxMismatches ?? 0;
168  const maxIndels = options.maxIndels ?? 0;
169  if (maxIndels !== 0 && maxIndels !== 1) throw new Error("maxIndels must be 0 or 1");
170  const t = normalize(template);
171  const p = normalize(primer);
172  if (p.length === 0) throw new Error("empty primer");
173
174  /** @type {Site[]} */
175  const found = [];
176  for (const strand of /** @type {const} */ (["+", "-"])) {
177    // The primer written as it appears on the + strand of the template.
178    const onTemplate = strand === "+" ? p : reverseComplement(p);
179    const L = onTemplate.length;
180    // Map an index in `onTemplate` back to a 1-based position in the primer from its 5′ end.
181    /** @param {number} i */
182    const primerPos = (i) => (strand === "+" ? i + 1 : L - i);
183
184    /**
185     * @param {number} at 0-based template start
186     * @param {number} span template bases covered
187     * @param {number[]} miss 0-based indices into `onTemplate` that mismatch
188     * @param {Indel | null} indel
189     */
190    const record = (at, span, miss, indel) => {
191      const covered = t.slice(at, at + span);
192      found.push({
193        strand,
194        start: at + 1,
195        end: at + span,
196        mismatches: miss.length,
197        mismatchPositions: miss.map(primerPos).sort((a, b) => a - b),
198        indel,
199        templateSite: strand === "+" ? covered : reverseComplement(covered),
200      });
201    };
202
203    for (let at = 0; at + L <= t.length; at += 1) {
204      const miss = hamming(t, onTemplate, at, maxMismatches);
205      if (miss) record(at, L, miss, null);
206    }
207    if (maxIndels === 0) continue;
208
209    // One extra base in the primer: drop each primer base in turn and fit what is left.
210    for (let drop = 0; drop < L; drop += 1) {
211      const probe = onTemplate.slice(0, drop) + onTemplate.slice(drop + 1);
212      for (let at = 0; at + probe.length <= t.length; at += 1) {
213        const miss = hamming(t, probe, at, maxMismatches);
214        if (!miss) continue;
215        // Indices past the dropped base belong one further along in the full primer.
216        const full = miss.map((i) => (i >= drop ? i + 1 : i));
217        record(at, probe.length, full, { type: "extra-primer-base", primerPosition: primerPos(drop) });
218      }
219    }
220    // One base missing from the primer: the template carries one more base than the primer does.
221    for (let gap = 1; gap < L; gap += 1) {
222      for (let at = 0; at + L + 1 <= t.length; at += 1) {
223        const window = t.slice(at, at + gap) + t.slice(at + gap + 1, at + L + 1);
224        const miss = hamming(window, onTemplate, 0, maxMismatches);
225        if (!miss) continue;
226        // `gap` template bases precede the skipped one. On the - strand the primer runs the
227        // other way, so the skipped base sits 3′ of primer position L - gap.
228        const after = strand === "+" ? gap : L - gap;
229        record(at, L + 1, miss, { type: "missing-primer-base", primerPosition: after });
230      }
231    }
232  }
233  return keepBest(found);
234}
235
236/** @param {Site} s */
237function edits(s) {
238  return s.mismatches + (s.indel ? 1 : 0);
239}
240
241/**
242 * @param {Site} a
243 * @param {Site} b
244 * @returns {number} negative when `a` is the better fit
245 */
246function better(a, b) {
247  return edits(a) - edits(b) || (a.indel ? 1 : 0) - (b.indel ? 1 : 0) || a.start - b.start;
248}
249
250/**
251 * @param {Site[]} sites
252 * @returns {Site[]}
253 */
254function keepBest(sites) {
255  /** @type {Site[]} */
256  const kept = [];
257  for (const s of [...sites].sort(better)) {
258    const overlaps = kept.some((k) => k.strand === s.strand && k.start <= s.end && s.start <= k.end);
259    if (!overlaps) kept.push(s);
260  }
261  return kept.sort((a, b) => a.start - b.start || (a.strand < b.strand ? -1 : 1));
262}
263
264/**
265 * @typedef {object} Product
266 * @property {number} start 1-based template position of the 5′ end of the + strand primer
267 * @property {number} end 1-based template position of the 5′ end of the - strand primer
268 * @property {number} templateLength end - start + 1: the span on the reference
269 * @property {number} length the product as synthesised. Each primer is copied into the product
270 *   whole, so a primer with an extra base adds one and a primer missing a base removes one. With no
271 *   indel it equals `templateLength`.
272 * @property {number} forwardMismatches
273 * @property {number} reverseMismatches
274 * @property {Site} forward
275 * @property {Site} reverse
276 */
277
278/**
279 * @param {Site} site
280 * @returns {number} product length change the primer causes relative to the template
281 */
282function indelShift(site) {
283  if (!site.indel) return 0;
284  return site.indel.type === "extra-primer-base" ? 1 : -1;
285}
286
287/**
288 * Every product a forward and a reverse primer make on a template: a + strand site of one primer
289 * upstream of a - strand site of the other, the two pointing at each other. Both pairings are
290 * tried, so a reference filed in the opposite orientation to the primers still gives its product.
291 *
292 * Measured from the forward primer's 5′ end to the reverse primer's 5′ end, inclusive.
293 *
294 * @param {string} template
295 * @param {string} forward 5′ to 3′
296 * @param {string} reverse 5′ to 3′
297 * @param {{ maxMismatches?: number, maxIndels?: number, maxLength?: number }} [options]
298 *   `maxLength` (default 10000) drops products longer than a PCR run would make.
299 * @returns {Product[]} best first: fewest edits, then shortest
300 */
301export function amplicons(template, forward, reverse, options = {}) {
302  const maxLength = options.maxLength ?? 10000;
303  const f = findSites(template, forward, options);
304  const r = findSites(template, reverse, options);
305  /** @type {Product[]} */
306  const out = [];
307  /**
308   * @param {Site[]} plusFrom
309   * @param {Site[]} minusFrom
310   * @param {boolean} forwardIsPlus
311   */
312  const pair = (plusFrom, minusFrom, forwardIsPlus) => {
313    for (const plus of plusFrom) {
314      if (plus.strand !== "+") continue;
315      for (const minus of minusFrom) {
316        if (minus.strand !== "-" || minus.end < plus.start) continue;
317        const templateLength = minus.end - plus.start + 1;
318        if (templateLength > maxLength) continue;
319        const fwd = forwardIsPlus ? plus : minus;
320        const rev = forwardIsPlus ? minus : plus;
321        out.push({
322          start: plus.start,
323          end: minus.end,
324          templateLength,
325          length: templateLength + indelShift(plus) + indelShift(minus),
326          forwardMismatches: fwd.mismatches,
327          reverseMismatches: rev.mismatches,
328          forward: fwd,
329          reverse: rev,
330        });
331      }
332    }
333  };
334  pair(f, r, true);
335  pair(r, f, false);
336  return out.sort(
337    (a, b) =>
338      edits(a.forward) + edits(a.reverse) - (edits(b.forward) + edits(b.reverse)) ||
339      a.length - b.length ||
340      a.start - b.start,
341  );
342}
343
344/**
345 * The best product (fewest edits, then shortest). THROWS when the pair makes none within the
346 * allowed edits: a missing product is an answer the caller must see, not a length of zero.
347 *
348 * @param {string} template
349 * @param {string} forward
350 * @param {string} reverse
351 * @param {{ maxMismatches?: number, maxIndels?: number, maxLength?: number }} [options]
352 * @returns {Product}
353 */
354export function amplicon(template, forward, reverse, options = {}) {
355  const best = amplicons(template, forward, reverse, options)[0];
356  if (!best) {
357    throw new Error(
358      `no product within ${options.maxMismatches ?? 0} mismatch(es) and ${options.maxIndels ?? 0} indel(s) per primer`,
359    );
360  }
361  return best;
362}