src/science/genome.ts

Functions

Imports

It imports src/render/data.ts, src/render/figure.ts, src/render/plot.ts, src/render/series.ts. Imported by src/science/index.ts.

Source

1import * as Plot from "@observablehq/plot";
2import { defaultFormat, markKey, maxOf } from "../render/data.js";
3import { type FigureOptions, figure, type SeriesColor, validateFigure } from "../render/figure.js";
4import { decorate, keyed, renderPlot } from "../render/plot.js";
5import { planSeries } from "../render/series.js";
6
7/** One annotated feature on a genome: a gene, a coding sequence, a primer, a mutation site. */
8export interface GenomeFeature {
9  readonly name: string;
10  /** The first position, 1-based and inclusive, as in GenBank and GFF. */
11  readonly start: number;
12  /** The last position, inclusive. */
13  readonly end: number;
14  /** 1 for the forward strand, -1 for reverse, 0 or omitted for none. */
15  readonly strand?: 1 | -1 | 0;
16  /** The feature type (gene, CDS, UTR), which sets its color and legend entry. */
17  readonly type?: string;
18  /** The track (row) it belongs to, e.g. a segment or an annotation source. */
19  readonly track?: string;
20}
21
22/** Options for {@link genomeTrack}. */
23export interface GenomeTrackOptions extends FigureOptions {
24  readonly features: readonly GenomeFeature[];
25  /** The sequence length in nucleotides; the axis runs from 1 to this. */
26  readonly length: number;
27  /** Tracks in display order. Without it, first-seen order. */
28  readonly trackDomain?: readonly string[];
29  /** Feature types in color and legend order. Without it, first-seen order. */
30  readonly typeDomain?: readonly string[];
31  readonly colors?: Readonly<Record<string, SeriesColor>>;
32  /** The coordinate axis label (default "Position (nt)"). */
33  readonly positionLabel?: string | null;
34  /** Print feature names (default true). Names that do not fit are left to hover and the table. */
35  readonly featureLabels?: boolean;
36}
37
38interface Placed {
39  readonly key: string;
40  readonly feature: GenomeFeature;
41  readonly type: string;
42  readonly track: string;
43  readonly row: string;
44  /** The lane within the track, for keyboard order. */
45  readonly lane: number;
46}
47
48const KIND = "genomeTrack";
49const DEFAULT_TRACK = "Features";
50const DEFAULT_TYPE = "Feature";
51const ROW_SEPARATOR = "\u0000";
52
53/** Greedy interval packing: each feature takes the first lane of its track where it fits. */
54function packLanes(
55  features: readonly GenomeFeature[],
56  trackOf: (f: GenomeFeature) => string,
57): number[] {
58  const laneEnds = new Map<string, number[]>();
59  const order = features
60    .map((f, i) => ({ f, i }))
61    .sort((a, b) => a.f.start - b.f.start || a.i - b.i);
62  const lanes = new Array<number>(features.length).fill(0);
63  for (const { f, i } of order) {
64    const ends = laneEnds.get(trackOf(f)) ?? [];
65    let lane = ends.findIndex((end) => end < f.start);
66    if (lane === -1) lane = ends.length;
67    ends[lane] = f.end;
68    laneEnds.set(trackOf(f), ends);
69    lanes[i] = lane;
70  }
71  return lanes;
72}
73
74/**
75 * A linear genome map: features drawn to scale along a nucleotide axis as arrows pointing in
76 * the direction of their strand, one row per track, with overlapping features stacked into
77 * lanes. Suited to viral genomes, segments, plasmids and amplicon maps.
78 *
79 * @example
80 * genomeTrack({
81 *   length: 1778,
82 *   features: [
83 *     { name: "HA1", start: 52, end: 1035, strand: 1, type: "Mature peptide" },
84 *     { name: "HA2", start: 1036, end: 1698, strand: 1, type: "Mature peptide" },
85 *   ],
86 *   alt: "Influenza A segment 4 ...",
87 * });
88 */
89export function genomeTrack(options: GenomeTrackOptions): string {
90  const { features, length } = options;
91  if (!Number.isInteger(length) || length < 1)
92    throw new Error(`${KIND}: length must be a positive integer, is ${length}`);
93  validateFigure(KIND, options);
94  if (features.length === 0) throw new Error(`${KIND}: features is empty`);
95  features.forEach((f, i) => {
96    if (
97      !Number.isInteger(f.start) ||
98      !Number.isInteger(f.end) ||
99      f.start < 1 ||
100      f.end > length ||
101      f.end < f.start
102    ) {
103      throw new Error(
104        `${KIND}: feature ${i + 1} "${f.name}" runs ${f.start} to ${f.end}; positions must be whole numbers from 1 to ${length}, start before end`,
105      );
106    }
107  });
108  const trackOf = (f: GenomeFeature): string => f.track ?? DEFAULT_TRACK;
109  const typeOf = (f: GenomeFeature): string => f.type ?? DEFAULT_TYPE;
110  const tracks = options.trackDomain
111    ? [...options.trackDomain]
112    : [...new Set(features.map(trackOf))];
113  const types = options.typeDomain ? [...options.typeDomain] : [...new Set(features.map(typeOf))];
114  for (const f of features) {
115    if (!tracks.includes(trackOf(f)))
116      throw new Error(`${KIND}: track "${trackOf(f)}" is missing from trackDomain`);
117    if (!types.includes(typeOf(f)))
118      throw new Error(`${KIND}: type "${typeOf(f)}" is missing from typeDomain`);
119  }
120  const plan = planSeries(KIND, types, options.colors);
121
122  const lanes = packLanes(features, trackOf);
123  const placed: Placed[] = features.map((f, i) => ({
124    key: markKey(trackOf(f), f.name, f.start),
125    lane: lanes[i] ?? 0,
126    feature: f,
127    type: typeOf(f),
128    track: trackOf(f),
129    row: `${trackOf(f)}${ROW_SEPARATOR}${lanes[i]}`,
130  }));
131  // Keyboard order: along the sequence, then down the lanes.
132  const byStart = new Map(
133    [...placed]
134      .sort((a, b) => a.feature.start - b.feature.start || a.lane - b.lane)
135      .map((p, i) => [p, i]),
136  );
137  const rows: string[] = [];
138  for (const track of tracks) {
139    const laneCount = Math.max(
140      0,
141      maxOf(placed.filter((p) => p.track === track).map((p) => p.lane + 1)),
142    );
143    for (let lane = 0; lane < Math.max(1, laneCount); lane += 1)
144      rows.push(`${track}${ROW_SEPARATOR}${lane}`);
145  }
146
147  const strandWord = (s: GenomeFeature["strand"]): string =>
148    s === 1 ? "forward strand" : s === -1 ? "reverse strand" : "no strand";
149  const describe = (p: Placed): string =>
150    `${p.feature.name} (${p.type}): ${defaultFormat(p.feature.start)} to ${defaultFormat(p.feature.end)}, ${defaultFormat(p.feature.end - p.feature.start + 1)} nt, ${strandWord(p.feature.strand)}`;
151
152  const width = options.width ?? 720;
153  const rowHeight = 30;
154  const height = options.height ?? rows.length * rowHeight + 64;
155  const trackLabel = (row: string): string => {
156    const [track, lane] = row.split(ROW_SEPARATOR);
157    return lane === "0" && tracks.length > 1 ? (track ?? "") : "";
158  };
159  const longestTrack = tracks.length > 1 ? maxOf(tracks.map((t) => t.length)) : 0;
160  const plotLeft = tracks.length > 1 ? Math.min(160, longestTrack * 7 + 20) : 16;
161  const pixelsPerNt = (width - plotLeft - 24) / length;
162
163  // Features are drawn as rects by Plot, then reshaped into arrows pointing along their strand.
164  const toArrow = (el: Element, p: Placed): void => {
165    const x = Number(el.getAttribute("x"));
166    const y = Number(el.getAttribute("y"));
167    const w = Number(el.getAttribute("width"));
168    const h = Number(el.getAttribute("height"));
169    const head = p.feature.strand ? Math.min(w * 0.4, h * 0.6) : 0;
170    const r = (n: number): number => Math.round(n * 100) / 100;
171    const points =
172      p.feature.strand === 1
173        ? [
174            [x, y],
175            [x + w - head, y],
176            [x + w, y + h / 2],
177            [x + w - head, y + h],
178            [x, y + h],
179          ]
180        : p.feature.strand === -1
181          ? [
182              [x + head, y],
183              [x + w, y],
184              [x + w, y + h],
185              [x + head, y + h],
186              [x, y + h / 2],
187            ]
188          : [
189              [x, y],
190              [x + w, y],
191              [x + w, y + h],
192              [x, y + h],
193            ];
194    const polygon = el.ownerDocument.createElementNS("http://www.w3.org/2000/svg", "polygon");
195    for (const name of el.getAttributeNames()) {
196      if (!["x", "y", "width", "height"].includes(name))
197        polygon.setAttribute(name, el.getAttribute(name) ?? "");
198    }
199    polygon.setAttribute(
200      "points",
201      points.map(([px, py]) => `${r(px ?? 0)},${r(py ?? 0)}`).join(" "),
202    );
203    polygon.append(...el.childNodes);
204    el.replaceWith(polygon);
205  };
206
207  const showLabels = options.featureLabels ?? true;
208  const fits = (p: Placed): boolean =>
209    (p.feature.end - p.feature.start + 1) * pixelsPerNt >= p.feature.name.length * 7 + 12;
210
211  const marks: Plot.Markish[] = [
212    Plot.ruleY(rows, {
213      y: (r: string) => r,
214      x1: 1,
215      x2: length,
216      stroke: "var(--enarratio-grid)",
217      strokeWidth: 2,
218    }),
219    Plot.rectX(placed, {
220      x1: (p: Placed) => p.feature.start - 0.5,
221      x2: (p: Placed) => p.feature.end + 0.5,
222      y: "row",
223      fill: "type",
224      inset: 0,
225      insetTop: 4,
226      insetBottom: 4,
227      title: describe,
228      render: keyed(
229        placed,
230        (p) => ({
231          key: p.key,
232          filter: { field: "type", value: p.type },
233          column: byStart.get(p) ?? 0,
234          row: p.lane,
235          x: String(p.feature.start),
236          series: p.type,
237        }),
238        toArrow,
239      ),
240    }),
241  ];
242  if (showLabels) {
243    const labelled = placed.filter(fits);
244    marks.push(
245      Plot.text(labelled, {
246        x: (p: Placed) => (p.feature.start + p.feature.end) / 2,
247        y: "row",
248        text: (p: Placed) => p.feature.name,
249        render: decorate(labelled, (el, p) => {
250          el.setAttribute("fill", `var(--enarratio-series-text-${plan.slot(p.type)})`);
251        }),
252      }),
253    );
254  }
255
256  const { svg } = renderPlot(
257    {
258      width,
259      height,
260      marginTop: 16,
261      marginLeft: plotLeft,
262      marginRight: 24,
263      marginBottom: 44,
264      x: {
265        type: "linear",
266        domain: [1, length],
267        label: options.positionLabel === undefined ? "Position (nt)" : options.positionLabel,
268        labelAnchor: "center",
269        labelArrow: "none",
270        grid: true,
271      },
272      y: {
273        type: "band",
274        domain: rows,
275        tickFormat: trackLabel,
276        label: null,
277        tickSize: 0,
278        padding: 0,
279      },
280      color: { domain: types, range: [...plan.range] },
281      marks,
282    },
283    options.alt,
284  );
285
286  return figure(options, {
287    kind: "genome",
288    svg,
289    legend: plan.legend,
290    seriesField: "type",
291    slotColors: plan.slotColors,
292    table: {
293      columns: [
294        "Feature",
295        "Type",
296        ...(tracks.length > 1 ? ["Track"] : []),
297        "Start",
298        "End",
299        "Length (nt)",
300        "Strand",
301      ],
302      rows: placed
303        .slice()
304        .sort(
305          (a, b) =>
306            tracks.indexOf(a.track) - tracks.indexOf(b.track) || a.feature.start - b.feature.start,
307        )
308        .map((p) => [
309          p.feature.name,
310          p.type,
311          ...(tracks.length > 1 ? [p.track] : []),
312          String(p.feature.start),
313          String(p.feature.end),
314          String(p.feature.end - p.feature.start + 1),
315          p.feature.strand === 1 ? "+" : p.feature.strand === -1 ? "-" : "",
316        ]),
317    },
318  });
319}