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 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 | 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 500x 500x 476x 476x 290x 290x 290x 290x 27x 27x 252x 252x 290x 27x 290x 194x 194x 194x 194x 194x 194x 194x 194x 194x 194x 194x 194x 194x 362x 362x 362x 362x 362x 362x 362x 362x 362x 362x 6839x 6839x 6839x 6839x 6839x 6839x 6839x 6839x 6839x 6839x 6839x 6839x 279430x 279430x 279430x 279430x 279430x 279430x 6824x 6839x 347x 347x 347x 347x 347x 347x 347x 347x 347x 362x 7831x 7831x 7831x 345x 345x 177x 177x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 58x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 97x 7567x 7567x 7567x 7567x 7567x 7567x 91x 91x 91x 91x 48x 48x 7829x 7829x 7829x 1x 6x 6x 6x 6x 6x 6x 6x 6x 6x 6x 39x 39x 39x 39x 39x 39x 6x 6x 405x 405x 405x 405x 405x 405x 2442x 2442x 2442x 37x 2442x 2405x 3798x 3798x 3798x 2405x 2442x 2442x 2442x 405x 405x 405x 405x 319x 319x 317x 317x 361x 86x 86x 86x 106x 106x 106x 21x 21x 1x 1x 2548x 2548x 2548x 2548x 2548x 405x 405x 405x 405x 2442x 2442x 2442x 2442x 15812x 15812x 1717x 11498x 11498x 473x 2442x 2442x 405x 405x 1x 1x 406x 406x 406x 406x 406x 406x 406x 406x 405x 405x 405x 406x 4127x 121762x 121762x 4127x 405x 402x 402x 405x 120602x 120602x 120602x 120602x 120602x 120602x 120132x 120132x 120602x 470x 470x 120602x 402x 402x 1x 38x 38x 38x 38x | /**
* BAI and CSI: which bytes hold the reads over a region.
*
* ## Why this is not `tabix.ts`
*
* It very nearly is. A `.tbi` and a `.bai` are the same binning index with a
* different header, and `parseTabix` already reads one. What it does with it is
* the difference, and the difference is the whole of this module.
*
* `tabix.ts` collapses every bin of a reference down to **one span per
* sequence** and throws the low sixteen bits of every virtual offset away. That
* is the right trade there and its reasoning is written out at
* `annotations/tabix.ts`: a feature lane holds a whole `FeatureSet` on the main
* thread and draws whichever part of it is on screen (§8.4), so a finer index
* would be carried and never asked. A contig's worth of GTF is a few hundred
* kilobytes.
*
* A contig's worth of BAM is not. There is no zoom level to fall back on and no
* summary layer in the file, so the unit of work has to be the **window**, and
* a window is what a bin index resolves and a per-sequence span does not. So
* this keeps what that one discards: bin-level chunks, both halves of every
* virtual offset, the linear index, and a byte forecast made from all three
* before a single byte of alignment data is fetched.
*
* ## The two formats
*
* BAI is the fixed case — 14-bit minimum interval, six levels — and is stored
* uncompressed. CSI generalises both numbers and is BGZF-compressed like a
* `.tbi`, which is why {@link parseBamIndex} sniffs rather than trusts the
* suffix. CSI matters here for one specific reason: **BAI cannot index a
* reference longer than 2^29 bp**, so any genome with a chromosome past ~512 Mb
* is CSI or nothing. It is not the metagenomic case — a co-assembly is tens of
* thousands of *short* contigs and indexes fine under BAI — but it is the human
* chromosome 1 case, which is a more common thing to point this at.
*
* ## The forecast, and why it is the interesting output
*
* `chunksFor` answers in compressed byte spans, and their total is knowable
* before anything is read. That number is what a coverage lane consults to
* decide whether it can afford a window at all — so a window too deep to draw
* is refused **from the index**, at zero data cost, rather than by starting a
* fetch and giving up part-way through it (§12).
*
* ## The geometry is the file's, and is checked
*
* Every count and both of CSI's shape parameters come out of the same bytes as
* the offsets, and each sizes something: a slot count an allocation, a depth
* the walk every query makes. So each goes through `sources/untrusted.ts` first
* (§10). The depth matters most, because a query's cost is quadratic in it:
* depths of two, four and eight thousand took 12, 31 and 89 ms a query, so a
* twenty-four-byte `.csi` declaring one near a million would hold the
* annotations worker for some twenty minutes a query. And the bins a query
* names are never enumerated past the number the index holds, so a `min_shift`
* of nought cannot turn a window into a billion candidates either.
*/
import { ViewerError } from '../../errors'
import { declaredCount } from '../../sources/untrusted'
import { BGZF_MAX_BLOCK, inflateBgzf, isBgzf } from '../tabix'
/** `BAI\1`. */
const BAI_MAGIC = [0x42, 0x41, 0x49, 0x01]
/** `CSI\1`. */
const CSI_MAGIC = [0x43, 0x53, 0x49, 0x01]
/** BAI's fixed geometry: the smallest bin covers 2^14 bp, over six levels. */
const BAI_MIN_SHIFT = 14
const BAI_DEPTH = 5
/**
* The bin a BAI reserves for its own statistics rather than for reads.
*
* `37450` is one past the last real bin of a six-level index. It holds two
* "chunks" that are not byte ranges at all — the unmapped-read span, and a pair
* of counts — so reading it as data would fetch a region of the file chosen by
* arithmetic on a read count.
*/
const BAI_PSEUDO_BIN = 37450
/** The format as a refusal names it. */
const INDEX = 'BAM index'
/**
* The most bins a BAI reference holds: its six levels number 37,450 bins and
* the pseudo-bin is one more, so 65,536 is past any writer, a duplicated bin
* or two included.
*/
const BAI_MAX_BINS = 1 << 16
/**
* The most linear-index slots a BAI reference holds.
*
* One slot a 16 kb window, and BAI's bins cannot address past 2^29 bp — which
* is why htslib refuses to write one for a longer reference and asks for CSI —
* so 2^29 / 2^14 = 32,768 is the format's own maximum (Picard writes one more).
* Twice it is the ceiling.
*/
const BAI_MAX_INTERVALS = 1 << 16
/**
* The deepest a CSI's binning can go: 10 levels below the root.
*
* A bin number is a `uint32`, and depth `d` numbers (8^(d+1) − 1) / 7 bins —
* 1.2 × 10^9 at ten, 9.8 × 10^9 at eleven — so an eleventh level could not be
* written down. htslib computes its first-bin offsets in an `int` and could not
* read one either. With htslib's default `min_shift` of 14, ten levels cover a
* reference of 2^44 bp; the longest chromosomes sequenced, a few gigabases,
* need seven.
*/
const CSI_MAX_DEPTH = 10
/**
* The widest a CSI's smallest bin can be, as `min_shift + 3 × depth`: htslib
* shifts 64-bit positions by it, so it stops at 63.
*/
const CSI_MAX_SHIFT = 63
/** An `int32` count's own bound, where the format sets none tighter than its bytes. */
const INT32_MAX = 2 ** 31 - 1
/** A span of the file, as the index stores it: two virtual offsets. */
export interface Chunk {
begin: bigint
end: bigint
}
/** One reference's index. */
export interface RefIndex {
/** Chunks by bin number. */
bins: Map<number, Chunk[]>
/**
* The linear index: the lowest virtual offset of any record overlapping each
* 2^`minShift` bp window. BAI only — CSI carries the same information as a
* `loffset` per bin instead.
*/
linear: BigUint64Array | null
/** CSI's per-bin lowest offset, serving the linear index's purpose. */
loffset: Map<number, bigint> | null
}
export interface BamIndex {
format: 'bai' | 'csi'
minShift: number
depth: number
refs: RefIndex[]
}
/** `true` when these bytes are a BAI or a (BGZF-wrapped) CSI. */
export function isBamIndex(bytes: Uint8Array): boolean {
return matches(bytes, BAI_MAGIC) || matches(bytes, CSI_MAGIC) || isBgzf(bytes)
}
function matches(bytes: Uint8Array, magic: number[]): boolean {
if (bytes.length < magic.length) return false
return magic.every((byte, index) => bytes[index] === byte)
}
/**
* Parse a `.bai` or a `.csi`, whichever these bytes are.
*
* The index is read whole, which is safe for a CSI and is a real cost for a
* BAI: one is a bin per populated region, the other is additionally a linear
* index entry per 16 kb of every reference, so a human WGS `.bai` is single-
* digit megabytes. That is a one-off on opening the file and it is why §12
* separates the open cost from the per-window cost rather than holding both to
* one number.
*/
export async function parseBamIndex(raw: Uint8Array): Promise<BamIndex> {
// A CSI is BGZF like a `.tbi`; a BAI is stored flat. Deciding by the bytes
// rather than the suffix costs one comparison and means a `.csi` named `.bai`
// still opens — which happens, because the two are handed around together.
const bytes = isBgzf(raw) ? await inflateBgzf(raw) : raw
const parse = matches(bytes, CSI_MAGIC) ? parseCsi : matches(bytes, BAI_MAGIC) ? parseBai : null
if (parse === null) {
throw new Error('This is not a BAM index: the magic bytes are neither BAI\\1 nor CSI\\1.')
}
try {
return parse(bytes)
} catch (err) {
// Every count is checked before it sizes anything, so the one `RangeError`
// left is a `DataView` read past the last byte: the index stops part-way
// through what its own counts promised.
if (!(err instanceof RangeError)) throw err
throw new ViewerError(
'format',
`This ${INDEX} ends part-way through its own tables, ${bytes.byteLength.toLocaleString('en-US')} ` +
'bytes in. It is truncated; rebuild it with `samtools index`.',
{ bytes: bytes.byteLength },
)
}
}
function parseBai(bytes: Uint8Array): BamIndex {
const view = new DataView(bytes.buffer, bytes.byteOffset, bytes.byteLength)
/** Bytes after `at`. */
const after = (at: number): number => bytes.byteLength - at
// Each reference is at least its two counts.
const refCount = declaredCount(view.getInt32(4, true), {
format: INDEX,
field: 'reference count',
ceiling: INT32_MAX,
available: after(8),
bytesPerItem: 8,
})
let at = 8
const refs: RefIndex[] = []
for (let ref = 0; ref < refCount; ref++) {
const bins = new Map<number, Chunk[]>()
// A bin number and a chunk count apiece, with the linear index's count after them.
const binCount = declaredCount(view.getInt32(at, true), {
format: INDEX,
field: `bin count of reference ${ref}`,
ceiling: BAI_MAX_BINS,
available: after(at + 4) - 4,
bytesPerItem: 8,
})
at += 4
for (let b = 0; b < binCount; b++) {
const bin = view.getUint32(at, true)
at += 4
const chunkCount = declaredCount(view.getInt32(at, true), {
format: INDEX,
field: `chunk count of bin ${bin}`,
ceiling: INT32_MAX,
available: after(at + 4),
bytesPerItem: 16,
})
at += 4
const chunks: Chunk[] = []
for (let c = 0; c < chunkCount; c++) {
chunks.push({
begin: view.getBigUint64(at, true),
end: view.getBigUint64(at + 8, true),
})
at += 16
}
// Read past it, then drop it: the pseudo-bin's "chunks" are counts.
if (bin !== BAI_PSEUDO_BIN) bins.set(bin, chunks)
}
// Sizes an allocation, and was the one count here read on trust: a slot
// count of 2^31 − 1 asked for sixteen gigabytes of `BigUint64Array`.
const intervalCount = declaredCount(view.getInt32(at, true), {
format: INDEX,
field: `linear index size of reference ${ref}`,
ceiling: BAI_MAX_INTERVALS,
available: after(at + 4),
bytesPerItem: 8,
})
at += 4
const linear = new BigUint64Array(intervalCount)
for (let i = 0; i < intervalCount; i++) {
linear[i] = view.getBigUint64(at, true)
at += 8
}
refs.push({ bins, linear, loffset: null })
}
return { format: 'bai', minShift: BAI_MIN_SHIFT, depth: BAI_DEPTH, refs }
}
function parseCsi(bytes: Uint8Array): BamIndex {
const view = new DataView(bytes.buffer, bytes.byteOffset, bytes.byteLength)
const after = (at: number): number => bytes.byteLength - at
// The two numbers every query is shaped by, and the one a twenty-four-byte
// file once made a twenty-minute query of: `reg2bins`, `lowestOffset` and
// `highestOffset` all walk every level, and the level offsets are summed per
// level.
const depth = declaredCount(view.getInt32(8, true), {
format: INDEX,
field: 'CSI depth',
ceiling: CSI_MAX_DEPTH,
})
const minShift = declaredCount(view.getInt32(4, true), {
format: INDEX,
field: 'CSI min_shift',
ceiling: CSI_MAX_SHIFT - 3 * depth,
})
// The auxiliary data is skipped, but the reference count follows it.
const auxBytes = declaredCount(view.getInt32(12, true), {
format: INDEX,
field: 'auxiliary data length',
ceiling: INT32_MAX,
available: after(16) - 4,
})
let at = 16 + auxBytes
const refCount = declaredCount(view.getInt32(at, true), {
format: INDEX,
field: 'reference count',
ceiling: INT32_MAX,
available: after(at + 4),
bytesPerItem: 4,
})
at += 4
const refs: RefIndex[] = []
for (let ref = 0; ref < refCount; ref++) {
const bins = new Map<number, Chunk[]>()
const loffset = new Map<number, bigint>()
// A bin number, its lowest offset and a chunk count apiece.
const binCount = declaredCount(view.getInt32(at, true), {
format: INDEX,
field: `bin count of reference ${ref}`,
ceiling: INT32_MAX,
available: after(at + 4),
bytesPerItem: 16,
})
at += 4
for (let b = 0; b < binCount; b++) {
const bin = view.getUint32(at, true)
at += 4
loffset.set(bin, view.getBigUint64(at, true))
at += 8
const chunkCount = declaredCount(view.getInt32(at, true), {
format: INDEX,
field: `chunk count of bin ${bin}`,
ceiling: INT32_MAX,
available: after(at + 4),
bytesPerItem: 16,
})
at += 4
const chunks: Chunk[] = []
for (let c = 0; c < chunkCount; c++) {
chunks.push({
begin: view.getBigUint64(at, true),
end: view.getBigUint64(at + 8, true),
})
at += 16
}
bins.set(bin, chunks)
}
refs.push({ bins, linear: null, loffset })
}
return { format: 'csi', minShift, depth, refs }
}
/**
* `value >> bits`, for a `bits` that may reach 32.
*
* JavaScript takes a shift count **modulo 32**, so `x >> 32` is `x >> 0` —
* which is `x`, not zero. That is not a corner case here: htslib chooses a
* CSI's depth by `for (n_lvls = 0, s = 1 << min_shift; max_len > s; ++n_lvls,
* s <<= 3)`, so a reference longer than 2^29 bp gets `min_shift 14, depth 6`
* and the top level's shift is exactly 32 — on precisely the genomes CSI exists
* for, and the ones this module's own header names.
*
* Division rather than a special case, because it is right for every width and
* a coordinate here is far below 2^53.
*/
function shiftRight(value: number, bits: number): number {
return Math.floor(value / 2 ** bits)
}
/**
* Every bin that can hold a record overlapping `[start, end)`.
*
* The generalised form, which BAI is the `minShift: 14, depth: 5` case of. Each
* level of the tree is a contiguous run of bin numbers starting at `offset`,
* and a region touches a run at every level — so the answer is `depth + 1`
* ranges rather than a search.
*/
export function reg2bins(
start: number,
end: number,
minShift: number,
depth: number,
): number[] {
const out: number[] = []
const last = Math.max(start, end - 1)
let shift = minShift + depth * 3
let offset = 0
for (let level = 0; level <= depth; level++) {
const from = offset + shiftRight(start, shift)
const to = offset + shiftRight(last, shift)
for (let bin = from; bin <= to; bin++) out.push(bin)
offset += 1 << (level * 3)
shift -= 3
}
return out
}
/**
* The chunk lists of every populated bin {@link reg2bins} names for `[start, end)`.
*
* Level by level, as `reg2bins` walks — but a level whose run of bin numbers is
* longer than the bins the reference holds is answered by looking through those
* instead. The run's length is set by the file: a window of a megabase at a
* `min_shift` of nought is a million bins at the deepest level, and of a whole
* chromosome a quarter of a billion, every one of them a lookup that finds
* nothing. Looking through the index costs its size, at most `depth + 1` times
* over. The chunks are sorted by `chunksFor` afterwards, so the order they
* arrive in is not part of the answer.
*/
function binsOverlapping(ref: RefIndex, start: number, end: number, index: BamIndex): Chunk[][] {
const out: Chunk[][] = []
const last = Math.max(start, end - 1)
let shift = index.minShift + index.depth * 3
let offset = 0
for (let level = 0; level <= index.depth; level++) {
const from = offset + shiftRight(start, shift)
const to = offset + shiftRight(last, shift)
if (to - from + 1 > ref.bins.size) {
for (const [bin, chunks] of ref.bins) if (bin >= from && bin <= to) out.push(chunks)
} else {
for (let bin = from; bin <= to; bin++) {
const chunks = ref.bins.get(bin)
if (chunks !== undefined) out.push(chunks)
}
}
offset += 1 << (level * 3)
shift -= 3
}
return out
}
/**
* The lowest virtual offset any record overlapping `start` can be at.
*
* This is the whole value of the linear index, and it is what makes a query
* over a deep file bounded: a bin near the root of the tree holds chunks
* spanning a whole chromosome, and without this floor a 2 kb window would fetch
* every one of them. Chunks ending below the floor cannot contain an
* overlapping record and are dropped.
*
* CSI has no linear index, and its `loffset` per bin is the same fact stored
* per node instead: the lowest offset in the smallest bin containing `start`.
*/
function lowestOffset(ref: RefIndex, start: number, index: BamIndex): bigint {
if (ref.linear !== null) {
const slot = shiftRight(start, index.minShift)
if (ref.linear.length === 0) return 0n
return ref.linear[Math.min(slot, ref.linear.length - 1)]
}
if (ref.loffset === null) return 0n
// Walk the containing bins from the deepest level up, and take the first that
// is populated: the deepest is the tightest floor available.
let shift = index.minShift
let offset = 0
for (let level = index.depth; level >= 0; level--) {
offset = binOffsetOfLevel(level)
const held = ref.loffset.get(offset + shiftRight(start, shift))
if (held !== undefined) return held
shift += 3
}
return 0n
}
function binOffsetOfLevel(level: number): number {
let offset = 0
for (let l = 0; l < level; l++) offset += 1 << (l * 3)
return offset
}
/**
* The offset of the first record that starts at or after `end` — the ceiling.
*
* The counterpart to {@link lowestOffset}, and leaving it out is a mistake that
* is invisible in every content assertion and enormous in the forecast. The
* linear index prunes chunks *before* a window. Nothing prunes chunks *after*
* it, and a coarse bin holds records from anywhere in the reference — so a bin
* one level from the root, holding a handful of reads scattered the length of a
* chromosome, drags the span for a 10 kb window out to the end of the file.
* Measured on a small fixture: 1.4 kB of actual records, and a forecast of
* 371 kB. Every read is correct; the *gate* is not, and a gate that refuses
* windows it could draw is the whole feature failing quietly.
*
* The bound is available because a record goes in the smallest bin that wholly
* contains it, so every record in a bin starts within that bin's range. A bin
* whose range begins at or after `end` therefore holds only records that start
* at or after `end`, none of which can overlap — and since the file is
* coordinate-sorted, the earliest offset among them is where our records must
* already have finished.
*
* The walk is up to eight bins per level, coarsening as it goes, and it is
* **best-effort by design**. A record lives in exactly one bin, so a short read
* far past the query sits in a deep bin no coarse scan will reach, and the walk
* may find nothing — in which case there is no ceiling and every chunk is a
* candidate, which is exactly the behaviour before this existed. That is the
* safe direction: a ceiling that is too *high* costs bytes, and one that is too
* *low* silently drops reads. Bounded at ~8 × (depth + 1) probes rather than a
* pass over every populated bin, which on a human chromosome is a hundred
* thousand of them, per frame.
*/
function highestOffset(ref: RefIndex, end: number, index: BamIndex): bigint {
let best: bigint | null = null
let shift = index.minShift
for (let level = index.depth; level >= 0; level--) {
const offset = binOffsetOfLevel(level)
// Bins are numbered consecutively *within* a level and the levels run on
// from one another, so a scan that walks off the end of one level walks
// into the next — where bin numbers mean something else entirely. Level 1
// holds eight bins, 1 to 8, and bin 9 is the first of level 2: reading it
// as "the ninth bin of level 1" produces a ceiling from a bin that is not
// after the query at all, and every chunk beyond it is then dropped. That
// loses reads, silently, and only on a query near the end of a reference.
const levelBins = 1 << (level * 3)
const after = shiftRight(Math.max(0, end - 1), shift) + 1
for (let i = 0; i < 8 && after + i < levelBins; i++) {
const chunks = ref.bins.get(offset + after + i)
if (chunks === undefined) continue
for (const chunk of chunks) {
if (best === null || chunk.begin < best) best = chunk.begin
}
}
shift += 3
}
// Nothing after the query: every chunk found is a candidate, which is what
// the absence of a ceiling has always meant.
return best ?? BigInt(Number.MAX_SAFE_INTEGER) << 16n
}
/** A byte range to fetch, in the shape `inflateBgzf` and the sources want. */
export interface ByteSpan {
from: number
to: number
/**
* Where the last block holding a wanted record begins.
*
* The same field, for the same reason, as `tabix.ts`'s `ByteSpan`: a fetch
* has to over-ask by up to a block because a block's length is not knowable
* without reading its header, and the decode should not pay for the slack.
*/
lastBlock: number
/**
* Where the first wanted record begins *inside* the first inflated block.
*
* The other half of a virtual offset, and the half `tabix.ts` throws away —
* which is safe there and is not safe here. A BGZF block boundary is not a
* record boundary: `bgzf_write` copies bytes into a 64 KiB buffer and flushes
* it when full, so `samtools` splits records across blocks routinely, and any
* read longer than a block must be split. Starting the walk at inflated
* offset 0 then begins in the middle of whatever record the previous block
* ran out on, reads its bytes as a record length, and yields either nothing
* or garbage — silently, since a nonsense record fails the position filter
* like any other.
*
* Merged spans take the *first* chunk's offset, which is all that is needed:
* once the walk is aligned to a record it stays aligned, and the extra
* records between merged chunks are filtered by position like any others.
*/
dataFrom: number
}
/** How far apart two spans may be before reading across the gap costs more than skipping it. */
const COALESCE_GAP = 1 << 20
/**
* The byte spans holding every read overlapping `[start, end)` on `refId`.
*
* Chunks are pruned against the linear index, sorted, and merged where they
* touch or come within a megabyte — the same policy, and the same constant, as
* `tabix.ts`'s `spansFor`, because the sources below read in 1 MiB chunks and a
* smaller gap costs more to skip than to read across.
*
* Returns an empty array for a reference the index does not cover, which is the
* honest answer for a contig with no reads on it — not an error.
*/
export function chunksFor(
index: BamIndex,
refId: number,
start: number,
end: number,
fileSize: number,
): ByteSpan[] {
const ref = index.refs[refId]
if (ref === undefined) return []
const floor = lowestOffset(ref, Math.max(0, start), index)
const ceiling = highestOffset(ref, Math.max(1, end), index)
const wanted: Chunk[] = []
for (const chunks of binsOverlapping(ref, Math.max(0, start), Math.max(1, end), index)) {
// Both ends. A chunk ending before the floor holds nothing that reaches the
// window; one beginning at or after the ceiling holds nothing that starts
// before the window ends.
for (const chunk of chunks) {
if (chunk.end > floor && chunk.begin < ceiling) wanted.push(chunk)
}
}
if (wanted.length === 0) return []
wanted.sort((a, b) => (a.begin < b.begin ? -1 : a.begin > b.begin ? 1 : 0))
const spans: ByteSpan[] = []
for (const chunk of wanted) {
// A virtual offset is `(block offset << 16) | offset within the block`.
// The high 48 bits say which bytes to fetch; the low 16 say where the
// record actually starts once they are inflated, and both are needed —
// see {@link ByteSpan.dataFrom}.
const from = Number(chunk.begin >> 16n)
const dataFrom = Number(chunk.begin & 0xffffn)
const lastBlock = Number(chunk.end >> 16n)
const to = Math.min(fileSize, lastBlock + BGZF_MAX_BLOCK)
const previous = spans[spans.length - 1]
if (previous !== undefined && from <= previous.to + COALESCE_GAP) {
previous.to = Math.max(previous.to, to)
previous.lastBlock = Math.max(previous.lastBlock, lastBlock)
// `dataFrom` is *not* widened: it belongs to the span's first block, and
// the walk from there runs continuously through everything merged in.
} else {
spans.push({ from, to, lastBlock, dataFrom })
}
}
return spans
}
/** Compressed bytes these spans would fetch — the forecast §12's gate is spent against. */
export function forecastBytes(spans: readonly ByteSpan[]): number {
let total = 0
for (const span of spans) total += span.to - span.from
return total
}
|