All files / src/metacooler range.ts

96.41% Statements 188/195
91.95% Branches 80/87
92.3% Functions 12/13
96.41% Lines 188/195

Press n or j to go to the next uncovered block, b, p or k for the previous block.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441                2x 2x 2x 2x                                                             2x 2x 2x 2x 2x 2x 2x 2x     25x 25x 25x               2x 745x     745x 744x 745x 736x 736x 736x 8x 341x 8x 8x 341x 7x 7x           377x 377x 377x 377x 1116x 1116x 1116x 1116x 743x 5x 1116x 377x 377x                                   2x 34x 34x 34x 342x 342x 15x 15x 15x 15x 327x 327x 34x 34x 34x   43x 43x 43x                                   2x 503x 503x 503x 503x 503x 503x 503x       377x 377x 377x   377x 1115x 1115x 1111x 1115x   377x 377x 377x 503x 365x 503x         5x 5x   5x 5x 5x 360x 360x                         2x 331x 331x 331x 331x 331x 331x     331x 331x 331x 331x 12x 12x 315x 315x                                                             2x 369x 369x 369x 369x 369x   369x     369x 369x   369x 53x 18x 18x 18x       18x 18x 18x 53x 11x 11x 11x 11x 11x 11x 11x 24x 24x   358x 272x 272x                                   369x 369x               5x 5x 5x 5x         5x 5x 5x 5x                     233x 5x 5x 5x 5x 5x 5x 5x 34x 34x 34x 34x 34x 34x 34x 34x                                                 277x 277x 277x 277x 277x 277x 277x 8x 277x       4x 4x 4x 4x 4x                                                                         2x 315x 315x 315x 315x 315x 315x 315x 315x 315x 315x 315x 315x  
/**
 * Range resolution — a mirror of `metacooler.ops.range` and, underneath it,
 * cooler's `parse_region` / `Cooler.extent` (spec §5.4).
 *
 * The same string must mean the same bins in the browser and in the CLI, so
 * the parsing is ported rather than reimplemented, humanised units and all.
 */
 
import { UNBINNED_NAME } from '../constants'
import { NameIndex, type NameContext } from '../genome/naming'
import { truncateNames, ViewerError } from '../errors'
import { pyRepr } from '../format'
import type { Extent, ParsedRange, Table } from '../types'
 
/**
 * The slice of the reader {@link parseMetacoolerRange} needs.
 *
 * Structural rather than the concrete `MetaCooler` for two reasons: it breaks
 * what would otherwise be an import cycle between the class and this module,
 * and it lets the worker's `ReaderApi` — the same object seen through a
 * narrower interface — be passed straight in.
 */
export interface RangeTarget {
  /** `null` for a plain cooler: there are no genome names to resolve. */
  readonly genomes: Table | null
  readonly contigs: Table
  /**
   * This file's sequences under every convention that might name one.
   *
   * Optional because a synthesised target need not carry one, and because it is
   * only ever reached on the path where the strict answer was already an error
   * — see {@link parseMetacoolerRange}. Absent, one is built from the contig
   * table, which is the same index a moment later and costs the same.
   */
  readonly names?: NameIndex
  /** Row indices of {@link genomes} that may be shown (§7.4). */
  displayGenomes(): number[]
  genomeExtent(genome: string): Extent
  extent(region: string): Extent
}
 
/** Multipliers cooler's `parse_humanized` understands, after upper-casing. */
const UNITS: Record<string, number> = {
  K: 1e3,
  KB: 1e3,
  M: 1e6,
  MB: 1e6,
  G: 1e9,
  GB: 1e9,
}
 
/** Throw a malformed-range refusal carrying cooler's own wording. */
function badRange(message: string, detail: Record<string, unknown> = {}): never {
  throw new ViewerError('bad-range', message, detail)
}
 
/**
 * Port of cooler's `parse_humanized`: `"1,000"`, `"1.5kb"`, `"2M"` → an integer.
 *
 * The comma strip happens first, so thousands separators are legal everywhere;
 * the result is truncated toward zero exactly as Python's `int(float)` is.
 */
export function parseHumanized(text: string): number {
  const parts = text.replace(/,/g, '').split(/([0-9.]+)/)
  // Python unpacks the split into exactly three names, so anything with a
  // second numeric run (`"1k2"`) is an error there and must be one here.
  if (parts.length !== 3 || parts[0] !== '') badRange(`Unexpected token "${text}"`, { token: text })
  const [, digits, rawUnit] = parts
  if (rawUnit.length === 0) {
    if (!/^[0-9]+$/.test(digits)) badRange(`Unexpected token "${text}"`, { token: text })
    return Number(digits)
  }
  const value = Number.parseFloat(digits)
  if (!Number.isFinite(value)) badRange(`Unexpected token "${text}"`, { token: text })
  const unit = rawUnit.toUpperCase().trim()
  const multiplier: number | undefined = UNITS[unit]
  if (multiplier === undefined) badRange(`Unknown unit '${unit}'`, { unit })
  return Math.trunc(value * multiplier)
}
 
/** One token of a coordinate suffix, as cooler's `_tokenize` yields them. */
type Token = { kind: 'HYPHEN' | 'COORD' | 'OTHER'; text: string }
 
/** Port of cooler's `_tokenize`: hyphens, coordinates, and everything else. */
function tokenize(text: string): Token[] {
  const re = /\s*(-)|\s*([0-9,]+(?:\.[0-9]*)?(?:[a-z]+)?)|\s*(.+)/gi
  const tokens: Token[] = []
  for (const match of text.matchAll(re)) {
    const hyphen: string | undefined = match[1]
    const coord: string | undefined = match[2]
    const other: string | undefined = match[3]
    if (hyphen !== undefined) tokens.push({ kind: 'HYPHEN', text: hyphen })
    else if (coord !== undefined) tokens.push({ kind: 'COORD', text: coord })
    else tokens.push({ kind: 'OTHER', text: other ?? '' })
  }
  return tokens
}
 
/**
 * Split a list of regions on its separators, leaving thousands separators alone.
 *
 * A box that takes several sequences is comma-separated, and every coordinate in
 * this application is written with thousands separators (§26.1) — so a plain
 * `split(',')` reads the app's own spelling of a locus, `IV:1,000,001-1,400,000`,
 * as four regions called `IV:1`, `000`, `001-1` and `400`. The reconstruction
 * task came back *Expected HYPHEN token missing* for a locus copied out of the
 * heading directly above its box.
 *
 * The rule is the one a reader already has in their head: a comma is part of a
 * number when a digit is on its left and exactly three digits are on its right.
 * Everything else separates. So `I, II` and `1,2` are two regions each — a pair
 * of Ensembl-style contigs is not a mis-grouped number — while
 * `IV:1,000-2,000, V` is the locus and the sequence a reader meant.
 */
export function splitRegionList(text: string): string[] {
  const out: string[] = []
  let part = ''
  for (let at = 0; at < text.length; at++) {
    const char = text[at]
    if (char === ',' && !thousandsComma(text, at)) {
      out.push(part)
      part = ''
      continue
    }
    part += char
  }
  out.push(part)
  return out.map((one) => one.trim()).filter((one) => one !== '')
}
 
function thousandsComma(text: string, at: number): boolean {
  return /[0-9]$/.test(text.slice(0, at)) && /^[0-9]{3}(?![0-9])/.test(text.slice(at + 1))
}
 
/**
 * Port of cooler's `parse_region_string`.
 *
 * Splits on the first `:` only, so a contig whose name contains a colon is
 * unaddressable — a cooler limitation this viewer inherits deliberately rather
 * than diverging from.
 *
 * `reversed` lifts the one refusal this port makes that a **circle** does not
 * deserve: on a sequence that joins its own ends, `chr:2,900,000-100,000` is not
 * a typo, it is the 200 kb window across the origin — the most local contacts
 * there are, and the only window a reader cannot otherwise ask for. It is opt-in
 * and reached from one place ({@link parseMetacoolerRange}, for a sequence the
 * reader has declared), so every linear range keeps cooler's own error.
 *
 * @returns `[chrom, start, end]` with `null` for coordinates the string omits.
 */
export function parseRegionString(
  text: string,
  options: { reversed?: boolean } = {},
): [string, number | null, number | null] {
  const parts = text.split(':')
  const chrom = parts[0].trim()
  if (chrom.length === 0) badRange('Chromosome name cannot be empty', { region: text })
  if (parts.length < 2) return [chrom, null, null]
 
  // Only `parts[1]` is tokenised, as in cooler: a second colon and anything
  // after it is discarded rather than being an error.
  const tokens = tokenize(parts[1])
  let i = 0
  const next = (): Token | undefined => tokens[i++]
 
  const expect = (token: Token | undefined, kinds: Token['kind'][]): Token => {
    if (token === undefined) badRange(`Expected ${kinds.join(' or ')} token missing`, { region: text })
    if (!kinds.includes(token.kind)) badRange(`Unexpected token "${token.text}"`, { region: text })
    return token
  }
 
  const start = parseHumanized(expect(next(), ['COORD']).text)
  expect(next(), ['HYPHEN'])
  const tail = next()
  if (tail === undefined) return [chrom, start, null]
  const end = parseHumanized(expect(tail, ['COORD']).text)
  if (end < start && options.reversed !== true) {
    // cooler's own wording, with the one thing cooler has no way to say after
    // it: on a sequence the reader has declared circular this is a legal window
    // (§8.15), and the reader who typed it is usually a base pair away from
    // meaning it.
    badRange(
      'End coordinate less than start — unless the sequence joins its own ends, ' +
        'which is a tick in the circular column of the Chromosomes table',
      { region: text, start, end },
    )
  }
  return [chrom, start, end]
}
 
/**
 * Port of cooler's `parse_region`: fill in the open ends and bounds-check.
 *
 * Note what this is *not*: despite the UCSC-looking syntax, cooler applies no
 * 1-based-to-0-based conversion anywhere — `parse_region_string` returns the
 * digits as written and `parse_region` uses them as a 0-based half-open
 * interval. The docstrings in both cooler and metacooler say "1-based, fully
 * closed"; the code does not, and the code is what the CLI actually does, so
 * it is what this mirrors. Diverging would silently shift every region view by
 * one bin relative to `metacooler show --range`.
 */
export function parseRegion(
  region: string,
  chromsizes: Map<string, number>,
): [string, number, number] {
  const [chrom, rawStart, rawEnd] = parseRegionString(region)
  const clen = chromsizes.get(chrom)
  if (clen === undefined) {
    throw new ViewerError('unknown-label', `Unknown sequence label: ${chrom}`, { label: chrom })
  }
  const start = rawStart === null ? 0 : rawStart
  const end = rawEnd === null ? clen : rawEnd
  if (end < start) badRange('End cannot be less than start', { region, start, end })
  if (start < 0 || end > clen) {
    badRange(`Genomic region out of bounds: [${start}, ${end})`, { region, start, end, length: clen })
  }
  return [chrom, start, end]
}
 
// ---------------------------------------------------------------------------
// metacooler-level range resolution (ops/range.py)
// ---------------------------------------------------------------------------
 
/**
 * Resolve a range string to a half-open bin interval, mirroring
 * `parse_metacooler_range` including every one of its refusals.
 *
 * Resolution order, and it matters: the label (everything before the first
 * `:`) is tested against **genome** names first. Genome and contig namespaces
 * can collide, and Python lets the genome win; matching contigs first would
 * make the same string mean different things in the two implementations.
 *
 * Three refusals, each distinct on purpose:
 * - `__unbinned__` → `unbinned-not-displayable`. Unassigned contigs are not a
 *   genome (§7.4); folding this into "unknown label" would leave the user
 *   doubting a name they read off the file.
 * - a coordinate suffix on a genome name → `genome-coordinate-suffix`. Python
 *   raises rather than silently ignoring the suffix, so a user who asks for
 *   part of a genome is told that genome ranges are whole-genome.
 * - anything else → `unknown-label`, listing both namespaces truncated to 20.
 *
 * ...except that the last of those three is now tried twice. Before refusing, a
 * label is offered to {@link NameIndex}, which knows that `chr1`, `1`, `01` and
 * `I` are one chromosome under four conventions (§5.9). `context` is the naming
 * context of the file the label came *out of*, where the caller has one — a
 * companion panel following a lead of another nomenclature — because a Roman
 * numeral is Roman by virtue of the list it was written in.
 */
export function parseMetacoolerRange(
  mc: RangeTarget,
  rangeStr: string,
  context?: NameContext,
  circular?: readonly string[],
): ParsedRange {
  // No trim: Python splits the raw string, so `' MAG_1'` is genuinely unknown.
  const label = rangeStr.split(':')[0]
  // Membership by scanning the name columns, exactly as `ops/range.py` tests
  // `label in mc.genomes["name"].tolist()`.
  const genomeNames = mc.genomes === null ? [] : mc.genomes.strings('name')
  const contigNames = mc.contigs.strings('name')
 
  if (genomeNames.includes(label)) {
    if (label === UNBINNED_NAME) {
      throw new ViewerError(
        'unbinned-not-displayable',
        `Range ${pyRepr(rangeStr)}: ${pyRepr(UNBINNED_NAME)} holds every contig that binning ` +
          `failed to assign. It is a bookkeeping row, not a genome — there is no coherent ` +
          `contact map to draw for it — so it cannot be displayed. Its share of the file is ` +
          `reported as a QC statistic in the file info panel.`,
        { label },
      )
    }
    if (rangeStr.includes(':')) {
      throw new ViewerError(
        'genome-coordinate-suffix',
        `Range ${pyRepr(rangeStr)}: ${pyRepr(label)} is a genome name — coordinate suffixes ` +
          `are not allowed for genome ranges. Use just ${pyRepr(label)}.`,
        { label, region: rangeStr },
      )
    }
    return { kind: 'genome', label, extent: mc.genomeExtent(label) }
  }
 
  if (contigNames.includes(label)) {
    return { kind: 'contig', label, extent: extentOf(mc, label, rangeStr, circular) }
  }
 
  // Nothing under that spelling. **This, and only this, is where the naming
  // conventions get a hearing** (§5.9) — after both namespaces have been asked
  // literally and the answer would otherwise be a refusal.
  //
  // Putting it here rather than at the top is the whole of how this stays a
  // faithful port. Every string the file *does* carry resolves exactly as
  // `ops/range.py` resolves it, in the same order, with the same genome-wins
  // rule; a session of one file cannot tell this clause exists. What changes is
  // only the error case: a deck holding a UCSC `.hic` beside an Ensembl `.cool`
  // used to answer one of them `unknown-label` and show a blank panel, and now
  // answers both.
  //
  // Contigs only, deliberately. A genome name in a metacool is a MAG somebody
  // binned, not a chromosome anybody has a convention for, and §8.5's argument
  // — two independently assembled organisms agree about a name by coincidence —
  // applies to those names hardest of all.
  const matched = (mc.names ?? indexOf(mc)).find(label, { from: context })
  if (matched !== null) {
    // Answered here rather than by re-entering this function with the rewritten
    // string. A second pass would start again at the *genome* test, and in a
    // metacool whose MAGs are numerically named a contig matched across
    // conventions would come back as a genome — or, with a coordinate suffix on
    // it, as `genome-coordinate-suffix`. What was matched is a contig; this
    // says so, and it is the only way `parseMetacoolerRange` can be reached
    // twice for one string.
    const suffix = rangeStr.slice(label.length)
    return {
      kind: 'contig',
      label: matched.name,
      // Under **this file's** spelling, which is what the wrap is looked up by:
      // the declaration is a set of this file's own names (§8.15), and a reader
      // who typed `chr1:2.9M-100k` against a file that writes `1` is asking
      // about the same molecule.
      extent: extentOf(mc, matched.name, `${matched.name}${suffix}`, circular),
      matched,
    }
  }
 
  // The genome list is the *displayable* one: `__unbinned__` is not selectable
  // in any view (§7.4), so offering it as an available name would contradict
  // the refusal a user gets the moment they take the suggestion.
  //
  // A file with no genome layer has no genome to name: `ops/range.py` is
  // metacooler's and only ever sees a metacool, so its sentence — and the empty
  // `Genomes: []` it would print — is kept exactly for the files that have one,
  // and a cooler is told about its sequences alone. That is also every file a
  // page with the multi-genome format switched off can open (§29).
  if (mc.genomes === null) {
    throw new ViewerError(
      'unknown-label',
      `Range ${pyRepr(rangeStr)}: ${pyRepr(label)} is not a sequence of this file.\n` +
        `  Sequences: ${truncateNames(contigNames)}`,
      { label, region: rangeStr },
    )
  }
  throw new ViewerError(
    'unknown-label',
    `Range ${pyRepr(rangeStr)}: ${pyRepr(label)} is not a known genome or contig.\n` +
      `  Genomes: ${truncateNames(mc.displayGenomes().map((i) => genomeNames[i]))}\n` +
      `  Contigs: ${truncateNames(contigNames)}`,
    { label, region: rangeStr },
  )
}
 
/**
 * A contig range's bin extent — **unwrapped** where the sequence is a circle and
 * the coordinates run backwards over its origin (§8.15).
 *
 * `chr:2,900,000-100,000` on a 3 Mb circle is two stored runs, the last 100 kb
 * and the first 100 kb, and it is the window a reader most wants on a bacterial
 * chromosome: the junction is where the origin is, and on a line it is the one
 * place that cannot be looked at. What comes back is written in the same
 * **unwrapped** bin space `worker/wrap.ts` consumes — the first bin is where the
 * window opens and the last runs off the end of the sequence by however far the
 * window carries on past it — because that is the space the view pans in and the
 * only one in which panning across the junction is monotonic.
 *
 * Built out of two ordinary extents rather than fresh arithmetic, so the bin
 * edges are cooler's own on both sides of the junction: `chr:start-` for the
 * head and `chr:0-end` for the tail, with the sequence's own bin count between
 * them. A `.mcool` rung whose last bin absorbs a remainder (see
 * {@link regionToExtent}) is therefore handled here by being handled there.
 *
 * Every other range — a linear sequence, a window that reads forwards, a bare
 * contig name — goes straight to {@link RangeTarget.extent} and keeps cooler's
 * refusal.
 */
function extentOf(
  mc: RangeTarget,
  name: string,
  region: string,
  circular: readonly string[] | undefined,
): Extent {
  if (circular === undefined || !circular.includes(name)) return mc.extent(region)
  const [, start, end] = parseRegionString(region, { reversed: true })
  if (start === null || end === null || end >= start) return mc.extent(region)
  // Bounds-checked by the two reads below, which are ordinary ranges: a start
  // past the end of the sequence, or an end below zero, refuses there in
  // cooler's own words rather than being wrapped into something plausible.
  const whole = mc.extent(name)
  const head = mc.extent(`${name}:${start}-`)
  const tail = mc.extent(`${name}:0-${end}`)
  return [head[0], whole[1] + (tail[1] - whole[0])]
}
 
/** A name index over a target that did not bring one. */
function indexOf(mc: RangeTarget): NameIndex {
  const names = mc.contigs.strings('name')
  const lengths = mc.contigs.numbers('length')
  return NameIndex.of(names.map((name, i) => ({ name, length: lengths[i] })))
}
 
/**
 * The bin extent of a `contig` / `contig:start-end` region, as
 * `cooler.Cooler.extent` computes it for a fixed-bin file —
 * `chrom_offset[cid] + floor(start / binsize)` to
 * `chrom_offset[cid] + ceil(end / binsize)` — **clamped to the contig's own
 * bins**.
 *
 * The clamp is not cosmetic, and it is the one place this port deliberately
 * answers something cooler's arithmetic does not.
 *
 * A coarsened `.mcool` level does not carry `binnify(chromsizes, binsize)`'s
 * bin table. `cooler zoomify` merges the fine bins in groups of the factor, and
 * where a chromosome's fine bin count is not a multiple of it, the remainder is
 * absorbed into the last coarse bin rather than given one of its own — yeast
 * chromosome I at 2 kb ends in a single `228,000–230,218` bin, so the level has
 * 115 bins for it where `ceil(230218 / 2000)` says 116. The extra one is the
 * *first bin of chromosome II*, and the arithmetic hands it back for every
 * chromosome whose fine bin count was odd: I, VI, XII, XIV and XVI in this
 * assembly.
 *
 * What that costs is not one stray column. A block that spans two chromosomes
 * has two axis groups, so the view stops being a single-coordinate one: the
 * base-pair ticks along the edge disappear, the label switches to two names,
 * and the picture quietly gains a row and a column of another chromosome's
 * contacts. `chrom_offset` is the authority on which bins belong to a
 * chromosome — it is what the pixels are indexed by — so it, not the division,
 * decides where the block ends.
 */
export function regionToExtent(
  chromOffset: ArrayLike<number>,
  contigId: number,
  start: number,
  end: number,
  binsize: number,
): Extent {
  const base = chromOffset[contigId]
  const limit = chromOffset[contigId + 1] ?? base
  const lo = Math.min(base + Math.floor(start / binsize), limit)
  const hi = Math.min(base + Math.ceil(end / binsize), limit)
  return [lo, Math.max(lo, hi)]
}