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 | 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 46x 46x 46x 46x 46x 46x 46x 46x 126x 126x 46x 46x 2x 4x 4x 2x 2x 2x 2x 2x 2x 2x 10x 10x 10x 16x 16x 16x 16x 10x 10x 10x 2x 30x 30x 30x 30x 30x 30x 30x 30x 86x 86x 86x 86x 30x 30x 2x 2x 2x 2x 2x 41x 15x 15x 41x 2x 2x 2x 2x 2x 2x 2x 2x 13x 13x 13x | /**
* The single-cell container: many coolers in one file, read as one map.
*
* An `.scool` (Wolff et al., *Bioinformatics* 37:2053) is a set of coolers over
* **one shared bin space** — one per cell — and that shared bin space is the
* whole of why this file is short. Every cell indexes the same bins, so a bin
* id means the same thing in all of them, and summing a set of cells is not an
* operation at all: it is the same read performed several times into the same
* accumulator.
*
* ## Why the sum needs no arithmetic
*
* `scatterDense` is **additive and per-pixel** (`metacooler/matrix.ts`), which
* it has to be anyway — a read arrives in chunks and a bounded rectangle may be
* covered by several. So a store that hands the same callback every cell's
* chunks in turn produces their sum, exactly, with no summing code, no
* intermediate matrix and a peak footprint of one chunk. The aggregation *is*
* the fan-out.
*
* The same is true of the two other things a store answers. `nnz` is a total
* over the set. `binOffsets` is a cumulative count, and a sum of cumulative
* counts is the cumulative count of the sum — so the density column, the
* band planner and the selection read all work over a cell set without
* knowing one is there.
*
* ## What this deliberately does not do
*
* **Balance per cell.** Weights are applied once, to the summed block, from the
* shared `/bins` — which is cooler's own definition of a balanced matrix
* applied to the pseudo-bulk the reader asked for. Balancing each cell before
* summing would make the result something other than a contact count, and the
* files that most want summing (raw single-cell coolers, this format's whole
* reason for existing) carry no weights to do it with.
*
* **Hold a cooler per cell.** A cell is a group path and an offsets array, and
* the sibling-rung machinery that opens a whole `MetaCooler` per entry does not
* scale to the several thousand cells a real file has. What is held is one
* borrowed handle per *selected* cell, over the one mount.
*/
import { yieldToEventLoop } from '../async'
import { ViewerError } from '../errors'
import type { Hdf5Handle } from '../h5/api'
import { pyList } from '../format'
import type { Extent, Layout, PixelChunk } from '../types'
import { cellPath } from './detect'
import { DEFAULT_COUNT_COLUMN, type ReadPlan, type SpanScanOptions } from './matrix'
import { CoolerStore, type PixelStore, type ScanAllOptions } from './store'
/**
* How many cells a single aggregate may sum over before the caller must say so.
*
* Not a limit on the picture — the block a sum produces is the same size
* whether it came from one cell or four thousand — but on the *reads*: every
* cell in the set is a walk of its own pixel bands, so the cost is linear in
* cells where every other read in this application is bounded by the rectangle
* asked for. Above this the caller gates, the way the whole-file view gates on
* its pixel count.
*/
export const CELL_SUM_GATE = 64
/** One cell of an scool, resolved to something that can be read. */
interface Cell {
name: string
store: CoolerStore
}
export class ScoolStore implements PixelStore {
readonly countKind: 'i' | 'u' | 'f'
private readonly cells: Cell[]
private offsets: Float64Array | null = null
/**
* Told how far the fan-out has got, per cell.
*
* Every other read in this application is bounded by the rectangle asked for
* and needs no such thing. This one is bounded by the rectangle *times the
* cell count*, so the same window that draws instantly over one cell is a
* minute of small reads over four thousand — and the only honest denominator
* is cells, because that is what the work is counted in.
*
* A field rather than an argument because the callers are the view layer's
* dozen matrix requests, which have no business knowing that this file has
* cells at all: the store is per selection, so whoever mounted it can ask to
* be told and everything above stays as it was.
*/
onCells: ((done: number, total: number) => void) | null = null
/**
* @param handle The file's root handle; each cell is a borrowed view of it.
* @param layout The detected layout, which says where the cells live.
* @param names The selected cells, in the order the reader ticked them.
* @param nBins Bins in the shared space — the same for every cell.
*/
constructor(
handle: Hdf5Handle,
layout: Layout,
names: readonly string[],
private readonly nBins: number,
countKind: 'i' | 'u' | 'f',
) {
this.countKind = countKind
this.cells = names.map((name) => ({
name,
store: new CoolerStore(handle.at(cellPath(layout, name)), nBins, countKind),
}))
}
/** The cells this store is summing, in the order they were given. */
get names(): string[] {
return this.cells.map((cell) => cell.name)
}
/**
* Stored pixels across the set — a total, not a distinct count.
*
* Two cells that both record a contact between the same pair of bins
* contribute two stored pixels and one cell of the summed block, so this
* overstates the block's occupancy and understates nothing. It is used to
* forecast *read cost*, which is what it measures exactly.
*/
get nnz(): number {
let total = 0
for (const cell of this.cells) total += cell.store.nnz
return total
}
/**
* The value columns every selected cell offers.
*
* An intersection rather than the first cell's list: a column one cell has
* and another does not cannot be summed across the set, and offering it would
* put a picker entry on screen whose read fails on the second cell. `count`
* is in every cooler by construction, so the intersection is never empty.
*/
scoreColumns(): string[] {
const [first, ...rest] = this.cells
if (first === undefined) return [DEFAULT_COUNT_COLUMN]
let columns = first.store.scoreColumns()
for (const cell of rest) {
const theirs = new Set(cell.store.scoreColumns())
columns = columns.filter((name) => theirs.has(name))
}
return columns.length > 0 ? columns : [DEFAULT_COUNT_COLUMN]
}
/**
* Cumulative stored-pixel counts over the whole set.
*
* The sum of the cells' own `indexes/bin1_offset`, element by element, which
* is the cumulative count of the summed pixel set because a cumulative sum is
* linear. Memoised: it is `nBins + 1` doubles and every band plan reads it.
*
* This is the one place the set costs something up front — one small index
* read per cell — and it is why a selection is gated on its cell count.
*/
binOffsets(): Float64Array {
if (this.offsets !== null) return this.offsets
const total = new Float64Array(this.nBins + 1)
for (const cell of this.cells) {
const offsets = cell.store.binOffsets()
const n = Math.min(total.length, offsets.length)
for (let i = 0; i < n; i++) total[i] += offsets[i]
}
this.offsets = total
return total
}
/**
* Every selected cell's pixels for this rectangle, into one callback.
*
* The plan passes through unchanged, and that it can is the second thing the
* shared bin space buys: `ReadPlan.span` is a *bin* range, resolved against
* each cell's own `bin1_offset` when it is read (`forEachPixelSpanAsync`), so
* one plan addresses every cell correctly. It carries the storage mode's
* decisions — transpose, mirror, skip the diagonal — which an scool's cells
* share as they share their bins.
*/
async forEachWindow(
rows: Extent,
cols: Extent,
plan: ReadPlan,
onChunk: (chunk: PixelChunk) => void,
options: SpanScanOptions = {},
): Promise<boolean> {
let done = 0
for (const cell of this.cells) {
if (!(await cell.store.forEachWindow(rows, cols, plan, onChunk, options))) return false
// Between cells rather than only inside them: a set of four thousand
// sparse cells is four thousand small reads, and a cancel arriving as a
// message needs the queue to run to be seen (`../async`).
await yieldToEventLoop()
this.onCells?.(++done, this.cells.length)
}
return true
}
async forEachBands(
bands: readonly Extent[],
onChunk: (chunk: PixelChunk) => void,
options: SpanScanOptions = {},
): Promise<boolean> {
let done = 0
for (const cell of this.cells) {
if (!(await cell.store.forEachBands(bands, onChunk, options))) return false
await yieldToEventLoop()
this.onCells?.(++done, this.cells.length)
}
return true
}
/**
* Every pixel of every selected cell.
*
* Progress is reported against the set's total so a scan over four thousand
* cells advances a single bar rather than restarting it per cell.
*/
async scanAll(options: ScanAllOptions): Promise<void> {
const total = this.nnz
let done = 0
let stopped = false
for (const cell of this.cells) {
if (stopped) return
const before = done
await cell.store.scanAll({
...(options.chunkSize === undefined ? {} : { chunkSize: options.chunkSize }),
...(options.storedRows === undefined ? {} : { storedRows: options.storedRows }),
onChunk: (chunk) => {
const keepGoing = options.onChunk(chunk)
if (keepGoing === false) stopped = true
return keepGoing
},
onProgress: (cellDone) => {
done = before + cellDone
options.onProgress?.(done, total)
},
})
done = before + cell.store.nnz
options.onProgress?.(Math.min(done, total), total)
}
}
/**
* The same walk, with each chunk told which cell it came out of.
*
* {@link scanAll} exists to make a set of cells behave as one cooler, and this
* is the one question that needs it not to: *what is each cell on its own*.
* Nothing else in this application asks it, because nothing else in this
* application draws a cell — the Cells table ranks them by the attributes
* their writer recorded, which is a statement about the sequencing rather than
* about the contacts.
*
* **It costs exactly what a pseudo-bulk costs**, which is the whole argument
* for having it: the reads are the same reads in the same order, and the only
* difference is that the accumulator on the other end is one per cell instead
* of one shared. A reader who has already agreed to sum four thousand cells
* has already agreed to this.
*/
async scanPerCell(options: PerCellScanOptions): Promise<void> {
const total = this.nnz
let done = 0
for (let index = 0; index < this.cells.length; index++) {
const cell = this.cells[index]
options.onCell?.(index, cell.name)
const before = done
let stopped = false
await cell.store.scanAll({
...(options.chunkSize === undefined ? {} : { chunkSize: options.chunkSize }),
onChunk: (chunk) => {
const keepGoing = options.onChunk(index, chunk)
if (keepGoing === false) stopped = true
return keepGoing
},
onProgress: (cellDone) => {
done = before + cellDone
options.onProgress?.(done, total)
},
})
if (stopped) return
done = before + cell.store.nnz
options.onProgress?.(Math.min(done, total), total)
// Between cells, for `forEachWindow`'s reason: a cancel arrives as a
// message and a loop that never returns to the event loop never sees one.
await yieldToEventLoop()
}
}
}
/** How a caller watches a per-cell walk ({@link ScoolStore.scanPerCell}). */
export interface PerCellScanOptions {
chunkSize?: number
/** Called before each cell's own pixels, with its index into the selection. */
onCell?(index: number, name: string): void
/** Returning `false` abandons the whole walk where it stands. */
onChunk(cell: number, chunk: PixelChunk): boolean | void
onProgress?(done: number, total: number): void
}
/**
* Resolve the cells a view names against the cells a file has.
*
* `null` means *all of them*, which is what a file addressed with no cell set
* means: an scool with nothing ticked is the pseudo-bulk of everything in it,
* and that is the picture somebody who dropped one and pressed on wanted.
*
* A name the file does not carry is refused rather than skipped. Silently
* dropping it would answer a request for eight cells with the sum of seven and
* no sign that anything was missing — which is the same argument the region
* parser makes about a contig name, and the same answer.
*/
export function resolveCells(available: readonly string[], wanted: readonly string[] | null): string[] {
if (wanted === null || wanted.length === 0) return [...available]
const have = new Set(available)
const missing = wanted.filter((name) => !have.has(name))
if (missing.length > 0) {
throw new ViewerError(
'cell-not-found',
`This file has no ${missing.length === 1 ? 'cell' : 'cells'} ${pyList(missing.slice(0, 8))}. ` +
`It carries ${available.length} ${available.length === 1 ? 'cell' : 'cells'}, ` +
`starting ${pyList(available.slice(0, 3))}.`,
{ missing, available: available.length },
)
}
// Deduplicated, and in the file's own order rather than the order they were
// named: the same cell twice would be counted twice, and a set read forwards
// is one pass over the file rather than a seek per name.
const set = new Set(wanted)
return available.filter((name) => set.has(name))
}
|