src/science/titer.ts

Functions

Imports

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

Source

1import * as Plot from "@observablehq/plot";
2import { domainOf, type KeysOfType, maxOf, minOf, readLabel, readNumber } from "../render/data.js";
3import { type FigureOptions, figure, type SeriesColor, validateFigure } from "../render/figure.js";
4import { keyed, renderPlot } from "../render/plot.js";
5import { planSeries } from "../render/series.js";
6import { DEFAULT_WIDTH } from "../render/size.js";
7import { type GeometricSummary, geometricSummary } from "./stats.js";
8
9/** A serial dilution series: the first dilution and the factor between steps. */
10export interface DilutionSeries {
11  /** The reciprocal of the first dilution tested, e.g. 10 for 1:10. */
12  readonly start: number;
13  /** The factor between steps, e.g. 2 for two-fold or 10 for ten-fold dilutions. */
14  readonly factor: number;
15}
16
17/** Options for {@link titerPlot}. */
18export interface TiterPlotOptions<T extends object> extends FigureOptions {
19  readonly data: readonly T[];
20  /** The field naming the group each sample belongs to: a cohort, a timepoint, a vaccine arm. */
21  readonly group: KeysOfType<T, string | number>;
22  /**
23   * The field holding each sample's titer as a reciprocal dilution (40 for 1:40). Values below
24   * the limit of detection are shown at half the limit, hollow, as is conventional.
25   */
26  readonly titer: KeysOfType<T, number>;
27  /** The dilution series the assay used, which sets the axis ticks (default 1:10, two-fold). */
28  readonly dilution?: DilutionSeries;
29  /** The limit of detection as a reciprocal dilution (default: the series start). */
30  readonly limitOfDetection?: number;
31  /** Groups in display order. */
32  readonly groupDomain?: readonly (string | number)[];
33  /** Color each group as its own series (default), or all in the first series color. */
34  readonly colorByGroup?: boolean;
35  readonly colors?: Readonly<Record<string, SeriesColor>>;
36  /** Tick labels as ratios ("1:40", default) or reciprocal titers ("40"). */
37  readonly tickLabels?: "ratio" | "reciprocal";
38  readonly groupLabel?: string | null;
39  readonly titerLabel?: string | null;
40}
41
42interface Sample {
43  readonly key: string;
44  readonly group: string;
45  readonly groupIndex: number;
46  readonly titer: number;
47  readonly plotted: number;
48  readonly belowLimit: boolean;
49  offset: number;
50}
51
52const KIND = "titerPlot";
53
54const formatTiter = (value: number): string =>
55  value >= 100 ? String(Math.round(value)) : String(Math.round(value * 10) / 10);
56
57/**
58 * Antibody or viral titers by group on a log dilution axis: every sample as a point, the
59 * geometric mean titer (GMT) with its 95% confidence interval, and the assay's limit of
60 * detection. Identical titers are spread side by side rather than drawn on top of each other,
61 * so the number of samples at each dilution is visible. The GMT counts a value below the limit
62 * of detection as half the limit, the usual convention, and the data table says how many there were.
63 *
64 * @example
65 * titerPlot({
66 *   data: sera,
67 *   group: "timepoint",
68 *   titer: "hai",
69 *   dilution: { start: 10, factor: 2 },
70 *   alt: "HAI titers rose from a GMT of 14 before vaccination to 160 at day 28 ...",
71 * });
72 */
73export function titerPlot<T extends object>(options: TiterPlotOptions<T>): string {
74  const { data, group, titer } = options;
75  validateFigure(KIND, options);
76  if (data.length === 0) throw new Error(`${KIND}: data is empty`);
77  const dilution = options.dilution ?? { start: 10, factor: 2 };
78  if (!(dilution.start > 0) || !(dilution.factor > 1)) {
79    throw new Error(`${KIND}: dilution needs start above 0 and factor above 1`);
80  }
81  const lod = options.limitOfDetection ?? dilution.start;
82  const groupLabel = options.groupLabel === undefined ? group : options.groupLabel;
83  const titerLabel = options.titerLabel === undefined ? "Titer" : options.titerLabel;
84  const asRatio = (options.tickLabels ?? "ratio") === "ratio";
85  const tick = (v: number): string => (asRatio ? `1:${formatTiter(v)}` : formatTiter(v));
86
87  const rawGroups = data.map((row, i) => readLabel(KIND, row, group, i));
88  const groups = domainOf(KIND, group, rawGroups, options.groupDomain?.map(String));
89  const samples = data.map((row, i): Sample => {
90    const value = readNumber(KIND, row, titer, i);
91    if (value <= 0)
92      throw new Error(`${KIND}: row ${i + 1} titer is ${value}; titers are positive reciprocals`);
93    const g = rawGroups[i] as string;
94    const belowLimit = value < lod;
95    return {
96      key: String(i),
97      group: g,
98      groupIndex: groups.indexOf(g),
99      titer: value,
100      plotted: belowLimit ? lod / 2 : value,
101      belowLimit,
102      offset: 0,
103    };
104  });
105
106  // Spread samples that share a group and titer across the group's width, centered.
107  const stacks = new Map<string, Sample[]>();
108  for (const s of samples) {
109    const id = `${s.group}|${s.plotted}`;
110    stacks.set(id, [...(stacks.get(id) ?? []), s]);
111  }
112  const widest = maxOf([...stacks.values()].map((s) => s.length));
113  const spacing = Math.min(0.08, 0.7 / Math.max(1, widest - 1));
114  for (const stack of stacks.values()) {
115    stack.forEach((s, k) => {
116      s.offset = (k - (stack.length - 1) / 2) * spacing;
117    });
118  }
119
120  // Keyboard order within a group: lowest titer first, then left to right.
121  const withinGroup = new Map<Sample, number>();
122  for (const g of groups) {
123    samples
124      .filter((s) => s.group === g)
125      .sort((a, b) => a.plotted - b.plotted || a.offset - b.offset)
126      .forEach((s, i) => {
127        withinGroup.set(s, i);
128      });
129  }
130  const summaries = new Map<string, GeometricSummary>(
131    groups.map((g) => [
132      g,
133      geometricSummary(samples.filter((s) => s.group === g).map((s) => s.plotted)),
134    ]),
135  );
136  const colorByGroup = options.colorByGroup ?? true;
137  const plan = planSeries(KIND, colorByGroup ? groups : ["Samples"], options.colors);
138  const colorOf = (s: Sample): string => (colorByGroup ? s.group : "Samples");
139
140  const all = samples.map((s) => s.plotted);
141  const top = Math.max(maxOf(all), lod);
142  const ticks: number[] = [];
143  for (let v = dilution.start; v <= top * dilution.factor; v *= dilution.factor) {
144    ticks.push(v);
145    // A step this fine is not a dilution series, and would draw thousands of ticks (F15).
146    if (ticks.length > 40) {
147      throw new Error(
148        `${KIND}: dilution factor ${dilution.factor} needs over 40 steps from ${dilution.start} to ${top}`,
149      );
150    }
151  }
152  if (lod / 2 < dilution.start) ticks.unshift(lod / 2);
153
154  const summaryRows = groups.map((g, i) => ({ g, i, s: summaries.get(g) as GeometricSummary }));
155  const intervalRows = summaryRows.filter(
156    (r) => r.s.lower !== undefined && r.s.upper !== undefined,
157  );
158  const describe = (s: Sample): string =>
159    `${groupLabel ?? group} ${s.group}: ${s.belowLimit ? `below the limit of detection (<${tick(lod)})` : tick(s.titer)}`;
160
161  const marks: Plot.Markish[] = [
162    Plot.ruleY([lod], { strokeDasharray: "4,3" }),
163    Plot.text([lod], {
164      y: (d: number) => d,
165      text: () => "LOD",
166      frameAnchor: "right",
167      textAnchor: "end",
168      dy: -6,
169    }),
170    Plot.link(intervalRows, {
171      x1: (r: (typeof intervalRows)[number]) => r.i,
172      x2: (r: (typeof intervalRows)[number]) => r.i,
173      y1: (r: (typeof intervalRows)[number]) => r.s.lower ?? r.s.mean,
174      y2: (r: (typeof intervalRows)[number]) => r.s.upper ?? r.s.mean,
175      stroke: "var(--enarratio-text)",
176      strokeWidth: 1.5,
177    }),
178    Plot.link(summaryRows, {
179      x1: (r: (typeof summaryRows)[number]) => r.i - 0.28,
180      x2: (r: (typeof summaryRows)[number]) => r.i + 0.28,
181      y1: (r: (typeof summaryRows)[number]) => r.s.mean,
182      y2: (r: (typeof summaryRows)[number]) => r.s.mean,
183      stroke: "var(--enarratio-text)",
184      strokeWidth: 2.5,
185    }),
186    Plot.dot(samples, {
187      x: (s: Sample) => s.groupIndex + s.offset,
188      y: "plotted",
189      r: 4,
190      stroke: colorOf,
191      fill: colorOf,
192      strokeWidth: 1.5,
193      fillOpacity: 0.85,
194      title: describe,
195      render: keyed(
196        samples,
197        (s) => ({
198          key: s.key,
199          filter: { field: group, value: s.group },
200          column: s.groupIndex,
201          row: withinGroup.get(s) ?? 0,
202          x: s.group,
203          series: s.group,
204        }),
205        // Below the limit of detection is hollow: a shape difference, not only a color one.
206        (el, s) => {
207          if (s.belowLimit) el.setAttribute("fill", "var(--enarratio-background)");
208        },
209      ),
210    }),
211  ];
212
213  const { svg } = renderPlot(
214    {
215      width: options.width ?? DEFAULT_WIDTH,
216      height: options.height ?? 360,
217      marginTop: 20,
218      marginLeft: titerLabel === null ? 56 : 72,
219      marginBottom: 44,
220      x: {
221        type: "linear",
222        domain: [-0.5, groups.length - 0.5],
223        ticks: groups.map((_, i) => i),
224        tickFormat: (i: number) => groups[i] ?? "",
225        label: groupLabel,
226        labelAnchor: "center",
227        labelArrow: "none",
228      },
229      y: {
230        type: "log",
231        base: dilution.factor,
232        domain: [minOf(ticks), maxOf(ticks)],
233        ticks,
234        tickFormat: tick,
235        label: titerLabel,
236        labelAnchor: "center",
237        labelArrow: "none",
238        grid: true,
239      },
240      color: { domain: plan.names, range: [...plan.range] },
241      marks,
242    },
243    options.alt,
244  );
245
246  return figure(options, {
247    kind: "titer",
248    svg,
249    legend: [],
250    slotColors: plan.slotColors,
251    table: {
252      columns: [groupLabel ?? group, "n", "GMT", "95% CI", "Below LOD", "Titers"],
253      rows: summaryRows.map(({ g, s }) => {
254        const members = samples.filter((x) => x.group === g);
255        return [
256          g,
257          String(s.n),
258          formatTiter(s.mean),
259          s.lower === undefined || s.upper === undefined
260            ? ""
261            : `${formatTiter(s.lower)} to ${formatTiter(s.upper)}`,
262          String(members.filter((x) => x.belowLimit).length),
263          members
264            .map((x) => (x.belowLimit ? `<${formatTiter(lod)}` : formatTiter(x.titer)))
265            .join(", "),
266        ];
267      }),
268    },
269  });
270}