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 | 2x 2x 1114x 1114x 1114x 2x 2x 2x 369x 369x 2x 4x 4x 2x 5x 5x 5x 5x 5x 5x 2x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 473x 2x 39x 39x 39x 39x 39x 39x 39x 39x 39x 39x 39x 1322x 1322x 1322x 39x 39x 39x 663x 663x 663x 39x 39x 39x 39x 39x 2x 93x 93x 93x 93x 93x 93x 93x 1310x 1310x 1310x 1310x 93x 93x 93x 93x 93x 2x 472x 472x 472x 472x 472x 614x 1114x 1114x 1114x 904x 904x 904x 614x 472x 472x 2x 1322x 1322x 1322x 1322x 2908x 2908x 2908x 1416x 1416x 1416x 1322x 1322x 2x 511x 511x 511x 511x 511x 511x 511x 511x 511x 511x 1556x 1556x 1556x 1556x 511x 511x 511x 511x 511x 2x 604x 604x 604x 604x 604x 604x 604x 604x 604x 604x 604x 604x 604x 2868x 2868x 2868x 2871x 2871x 2871x 2871x 2871x 2871x 2866x 2866x 2868x 2868x 602x 604x 2x 604x 604x 604x 604x 604x 604x 185x 185x 185x 185x 185x 2452x 2452x 2452x 306x 2458x 2458x 306x 306x 306x 306x 2452x 2452x 604x 2x 512x 1x 1x 1x 1x 1x 512x 2x | /**
* A `.hic`'s pixels, as the rest of the reader wants them.
*
* The counterpart to `hic/handle.ts`: that module answers the questions a
* cooler's *metadata* answers, and this one answers the question its pixel
* table answers. Together they are the whole of the container's surface, and
* nothing above either knows which file it is looking at.
*
* ## The one place this is better than a cooler, and why it is not exploited
*
* A cooler's index is one-dimensional: `bin1_offset` locates a row band, and
* whichever columns that band happens to hold come back and are discarded by
* the scatter. A `.hic`'s index is two-dimensional, so an off-diagonal
* rectangle — the pair view, the exact shape this application exists to draw —
* costs the blocks it covers rather than a whole row band across the genome.
*
* That advantage is taken here and is deliberately *not* pushed upward. No
* caller is given a faster path, no view branches on the container, and no
* budget is relaxed on the strength of it. A read that is cheaper is simply
* cheaper; making it visible would make the containers different, which is the
* property this whole subtree is built to avoid.
*
* ## Emitting each stored pixel exactly once
*
* Every accumulator above — `scatterDense`, `pixelAccumulator`,
* `sparseAggregator` — is **additive**, and mirrors what it is given across the
* diagonal under `symmetric-upper`. So a pixel handed over twice is a cell
* counted twice, silently, in a picture that still looks like a contact map.
* Every walk below therefore visits each stored chromosome pair at most once,
* and the pair set is built as a deduplicated set rather than as a nested loop.
*/
import { ViewerError } from '../errors'
import { yieldToEventLoop } from '../async'
import type { Extent, PixelChunk } from '../types'
import type { ReadPlan, SpanScanOptions } from '../metacooler/matrix'
import type { PixelStore, ScanAllOptions } from '../metacooler/store'
import type { HicFile } from './HicFile'
/** A chromosome pair key, so a set can hold one. */
function pairKey(a: number, b: number): number {
return a * 1e6 + b
}
export class HicStore implements PixelStore {
constructor(private readonly file: HicFile) {}
/**
* A `.hic` states no record count, and this does not guess one.
*
* `sumCounts` per matrix is a sum of *values*, not a count of records, and a
* block's record count is only knowable by decompressing it. Every caller
* treats `null` as unknown; the one that used it for a progress denominator
* counts chromosome pairs instead, which is a real number rather than a
* plausible one.
*/
readonly nnz = null
/** Blocks hold int16 counts or float32 values; the float kind covers both. */
readonly countKind = 'f' as const
/** The index is two-dimensional: a rectangle costs its own blocks (`PixelStore.rectangular`). */
readonly rectangular = true
scoreColumns(): string[] {
// One value per record and no room in the format for a second. A `.hic`
// cannot carry the extra per-pixel columns a cooler may (`expected`, a
// p-value), so the score picker offers the one that exists.
return ['count']
}
binOffsets(): Float64Array {
throw new ViewerError(
'unsupported',
'A .hic has no per-bin pixel index, so the stored-pixel density column cannot be ' +
'computed for one. It is a metagenome diagnostic and a .hic has no genomes, so ' +
'nothing that draws asks for it.',
{},
)
}
/**
* The stored pixels a rectangle needs.
*
* `plan` is not consulted, and that is correct rather than lazy: it exists to
* tell a one-dimensional index *which* row band to fetch for a rectangle it
* cannot express, and this index can express the rectangle. What comes back is
* every stored pixel inside `rows × cols` **or** inside its transpose, which
* is the same set under upper-triangular storage and is exactly what
* `scatterDense` expects to be handed for any of the three plan modes.
*/
async forEachWindow(
rows: Extent,
cols: Extent,
_plan: ReadPlan,
onChunk: (chunk: PixelChunk) => void,
options: SpanScanOptions = {},
): Promise<boolean> {
this.refuseOtherColumns(options.value)
const pairs = this.pairsFor(rows, cols)
// A window inside one chromosome is one pair, read through the chunk cache a
// pan comes back to; a window across chromosomes — the whole genome at a
// coarse rung, at its widest — is a walk, and reads exactly (`HicFile`).
const sparse = pairs.length > 1
return this.walk(
pairs,
(a, b) => this.file.pairRecords(a, b, rows, cols, sparse),
onChunk,
options,
() => [rows, cols],
)
}
/**
* The selection view's read: a set of row bands, against every column.
*
* A pair is visited when its **first** chromosome — the lower of the two, and
* therefore the one holding `bin1` — is touched by a band, which mirrors the
* cooler path's semantics exactly: `forEachPixelRangesAsync` walks the pixels
* whose `bin1` lies in the band, and the accumulator mirrors the rest.
*/
async forEachBands(
bands: readonly Extent[],
onChunk: (chunk: PixelChunk) => void,
options: SpanScanOptions = {},
): Promise<boolean> {
this.refuseOtherColumns(options.value)
const touched = new Set<number>()
for (const band of bands) for (const chrom of this.file.chromsSpanning(band)) touched.add(chrom)
const whole: Extent = [0, this.file.nBins]
const pairs = this.file.indexedPairs().filter(([a]) => touched.has(a))
const sparse = pairs.length > 1
// The row window is the bands clipped to `a`; the column window is the file.
// Passing the union of the bands rather than one band at a time is what keeps
// a pair visited once when two bands land on one chromosome.
const windowOf = (a: number): readonly [Extent, Extent] | null => {
const extent = this.clipTo(a, bands)
return extent === null ? null : [extent, whole]
}
return this.walk(
pairs,
(a, b) => {
const window = windowOf(a)
return window === null ? Promise.resolve(null) : this.file.pairRecords(a, b, window[0], window[1], sparse)
},
onChunk,
options,
windowOf,
)
}
/**
* Every stored pixel, pair by pair — the genome × genome overview's read.
*
* **`options.chunkSize` is not honoured, and that is a limitation rather than
* a decision.** `CoolerStore.scanAll` reads a bounded window of the pixel
* table at a time, which is what lets `estimateBlockBytes` charge a constant
* `PIXEL_CHUNK_BYTES` for the read instead of a term in `nnz`. A `.hic`'s
* pixels are reachable only through `HicFile.pairRecords`, whose unit is a
* whole chromosome pair: it holds every decoded block of that pair live and
* then allocates one combined COO buffer over their total, so peak memory is
* roughly 28 B per stored record of the largest pair — hundreds of MB for
* chr1 at a fine rung, against the constant the estimator assumes. The block
* LRU does not bound it: it bounds what is *cached*, and the read holds
* strong references to blocks it has already been handed.
*
* Narrowing it from this side is not available. `HicFile.localSpan` unions
* the two axes for an intra pair, so asking for a sub-rectangle of one asks
* for the bounding square, and two neighbouring sub-rectangles come back
* holding the same blocks — a pixel handed over twice, which is the one thing
* the doctrine at the head of this file forbids. The fix is a block-batched
* read inside `HicFile` that `pairRecords` is then written on top of; until
* that exists the same caveat applies to {@link forEachWindow} and
* {@link forEachBands}, whose `chunkSize` is dropped for the same reason.
*
* **This walk runs behind the others**, and is written to. The census of §6.5
* is one, started the moment a `.hic` opens, on the thread every view reads
* through — see {@link inOrder} for what keeps it from holding them up.
*/
async scanAll(options: ScanAllOptions): Promise<void> {
const whole: Extent = [0, this.file.nBins]
const pairs = this.file.indexedPairs()
const sparse = pairs.length > 1
await this.inOrder(
pairs,
(a, b) => this.file.pairRecords(a, b, whole, whole, sparse),
(records, done) => {
const stopped = records !== null && options.onChunk(records) === false
options.onProgress?.(done, pairs.length)
return !stopped
},
undefined,
options.pause,
() => [whole, whole],
)
}
/**
* The stored pairs a rectangle touches.
*
* Built from the chromosomes on each axis and then filtered against the master
* index, so a file written intra-only costs no request for the inter pairs it
* does not have. Normalised to `(low, high)` and deduplicated, because
* `rows × cols` and `cols × rows` are one stored matrix.
*/
private pairsFor(rows: Extent, cols: Extent): Array<[number, number]> {
const rowChroms = this.file.chromsSpanning(rows)
const colChroms = this.file.chromsSpanning(cols)
const seen = new Set<number>()
const pairs: Array<[number, number]> = []
for (const r of rowChroms) {
for (const c of colChroms) {
const pair: [number, number] = r <= c ? [r, c] : [c, r]
const key = pairKey(pair[0], pair[1])
if (seen.has(key)) continue
seen.add(key)
pairs.push(pair)
}
}
return pairs
}
/** The span of `bands` that falls on one chromosome, as one global extent. */
private clipTo(chrom: number, bands: readonly Extent[]): Extent | null {
const [base, end] = this.file.chromExtent(chrom)
let lo = Number.POSITIVE_INFINITY
let hi = Number.NEGATIVE_INFINITY
for (const band of bands) {
const from = Math.max(band[0], base)
const to = Math.min(band[1], end)
if (to <= from) continue
lo = Math.min(lo, from)
hi = Math.max(hi, to)
}
return Number.isFinite(lo) ? [lo, hi] : null
}
/**
* The common loop of the three reads above: pair by pair, cancellable between
* them, answering `false` when it was cancelled.
*/
private walk(
pairs: ReadonlyArray<[number, number]>,
read: (a: number, b: number) => Promise<PixelChunk | null>,
onChunk: (chunk: PixelChunk) => void,
options: SpanScanOptions,
windowOf: (a: number, b: number) => readonly [Extent, Extent] | null,
): Promise<boolean> {
return this.inOrder(
pairs,
read,
(records, done) => {
if (records !== null && records.length > 0) onChunk(records)
options.onProgress?.(done, pairs.length)
return true
},
() => options.cancelled?.() === true,
undefined,
windowOf,
)
}
/**
* Each pair's records handed to `take` in the order given — `take` answering
* `false` stops the walk, and so does `cancelled` between pairs.
*
* **It yields before each pair, not after.** The reads under a window are
* synchronous on the worker, so whatever a walk reads before it next yields,
* every other request to the reader waits for — and a census is started the
* moment a `.hic` opens, so the request waiting was the loading screen's own.
* Yielding first means nothing is read until what was already queued has been
* answered. The yield is also what lets a cancel arrive at all, the argument
* `forEachPixelSpanAsync` makes for its window loop. A walk given `pause`
* waits on that instead, which is how the census stands aside for as long as
* a view is reading rather than for one turn of the loop (§12).
*
* **And a walk across pairs reads them a plan at a time.** Its reads are
* exact ones that do not hold the worker (`HicFile`), and they are round trips
* and little else, so the pairs are cut into plans (`HicFile.walkPlans`) and
* each plan's records and blocks are read together (`HicFile.prefetchWalk`),
* by the first of its pairs the walk reaches, before any of them is read —
* which then finds them held. Up to `HicFile.walkAhead` pairs are read at
* once, so the next plan is being read while the last of this one's pairs are
* handed over. The pixels are still handed over one pair at a time and in
* order, so a consumer sees exactly what a one-at-a-time walk would give it.
*/
private async inOrder(
pairs: ReadonlyArray<[number, number]>,
read: (a: number, b: number) => Promise<PixelChunk | null>,
take: (records: PixelChunk | null, done: number) => boolean,
cancelled?: () => boolean,
pause?: () => Promise<void>,
windowOf?: (a: number, b: number) => readonly [Extent, Extent] | null,
): Promise<boolean> {
const ahead = pairs.length > 1 ? this.file.walkAhead() : 1
const planned = this.planned(pairs, cancelled, pause, windowOf)
const queue: Array<Promise<PixelChunk | null>> = []
let next = 0
let done = 0
while (done < pairs.length) {
await (pause?.() ?? yieldToEventLoop())
if (cancelled?.() === true) return false
while (queue.length < ahead && next < pairs.length) {
const [a, b] = pairs[next]
// A plan a cancel stopped is read no further; the loop stops at its next turn.
const reading = planned(next).then((ready) => (ready ? read(a, b) : null))
next += 1
// A read the walk stops before reaching is left to finish unawaited,
// and its failure is nobody's to report. The one awaited below still
// throws.
reading.catch(() => undefined)
queue.push(reading)
}
const records = await queue.shift()!
done += 1
if (!take(records, done)) return false
}
return true
}
/**
* The plan each pair of a walk belongs to, read once, by the first of its
* pairs that asks: a function from a pair's place in the walk to whether its
* plan was read — `false` when a cancel stopped it.
*
* A walk of one pair, or one that does not say which rectangle of each pair
* it reads, is not planned, and each pair is read as it comes.
*/
private planned(
pairs: ReadonlyArray<[number, number]>,
cancelled: (() => boolean) | undefined,
pause: (() => Promise<void>) | undefined,
windowOf: ((a: number, b: number) => readonly [Extent, Extent] | null) | undefined,
): (at: number) => Promise<boolean> {
if (pairs.length < 2 || windowOf === undefined) return () => Promise.resolve(true)
const plans = this.file.walkPlans(pairs)
const planOf: number[] = []
plans.forEach((plan, index) => plan.forEach(() => planOf.push(index)))
const reads = new Map<number, Promise<boolean>>()
return (at) => {
const index = planOf[at]
let reading = reads.get(index)
if (reading === undefined) {
const walk = plans[index].flatMap(([a, b]) => {
const window = windowOf(a, b)
return window === null ? [] : [{ a, b, rows: window[0], cols: window[1] }]
})
reading = this.file.prefetchWalk(walk, { pause, cancelled })
reads.set(index, reading)
}
return reading
}
}
/**
* A value column that is not `count` is a question this format cannot answer.
*
* Refused by name rather than silently served as `count`, which is the exact
* failure `SpanScanOptions.value` was added to fix for coolers: one picture
* drawn from two quantities, with no tell.
*/
private refuseOtherColumns(value: string | undefined): void {
if (value === undefined || value === 'count') return
throw new ViewerError(
'unsupported',
`A .hic stores one value per record, so there is no '${value}' column to draw. ` +
'Extra per-pixel columns are a cooler feature.',
{ value },
)
}
}
|