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 | 2x 2x 2x 35x 35x 35x 35x 6x 6x 6x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 3x 3x 3x 3x 3x 3x 4x 4x 4x 2x 6x 6x 6x 6x 6x 1x 1x 1x 1x 1x 1x 1x 1x 2x 6x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 6x 6x 6x 6x 6x 6x 6x 6x 2x 2x 2x 2x 2x 2x 2x 2x 4x 4x 4x 3x 3x 3x 3x 6x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x | /**
* Contact probability against genomic separation, as a *descriptor*.
*
* This is the first task in the catalogue (§18.8) and the first one written, and
* both of those are the same fact: P(s) is the one analysis that applies to
* every contact map there is. A cooler has a diagonal, a metacool has a
* diagonal, a file with no genomes and no tracks and no annotations loaded still
* has a diagonal — so {@link EXPECTED_TASK.availability} is the only predicate in
* the registry that is very nearly `true`, and the panel is never empty. It also
* has a second job that decides how carefully the rest of this has to be
* written: three later tasks consume the **expected** curve it produces, so a
* mistake here is a mistake in four results and it looks plausible in all four.
*
* The half of the task that lives in this file is what it *is* — its parameters,
* the files it applies to, and what it will cost before anything is read. What
* it *does* is `analysis/run/expected.ts`, and `analysis/job.ts` explains why the
* two are different modules: this one is imported by the panel, the form and the
* registry on the main thread, and it must therefore never reach `MetaCooler`,
* the analysis worker or anything under `analysis/run/`.
*/
import { extentOf, framingCovers, tickedContigs } from '../context'
import type {
AnalysisContext,
AnalysisTask,
CostEstimate,
ParamSpec,
Params,
} from '../types'
import { NEEDS_MAP, READY, waitingFor } from '../types'
import { BALANCE_PARAM } from './shared'
// ---------------------------------------------------------------------------
// Reading the form
// ---------------------------------------------------------------------------
/**
* The scope, normalised.
*
* `strParam` in `analysis/job.ts` is these same three lines, and it is
* deliberately not reached for: that module is the runner's half of the split
* and names `MetaCooler`, so a descriptor importing it is a descriptor that has
* stopped being safe to load on the main thread the first time somebody adds a
* value import to it.
*
* An unrecognised value falls back to `contig` rather than being trusted. A
* stored session written by a later build can carry a scope this one has never
* heard of, and the honest reading of that is the narrowest scope rather than a
* whole-file scan nobody asked for.
*/
function scopeOf(params: Params): 'contig' | 'genome' | 'file' {
const value = params.scope
return value === 'genome' || value === 'file' ? value : 'contig'
}
function balanceOf(params: Params): boolean {
return params.balance === true
}
// ---------------------------------------------------------------------------
// Parameters
// ---------------------------------------------------------------------------
const PARAMS: readonly ParamSpec[] = [
{
kind: 'choice',
key: 'scope',
label: 'scope',
default: 'contig',
/**
* Three on a metagenome, two on a cooler, and the words differ.
*
* This used to offer all three on every file, on the argument that a
* `ParamSpec` is static data shared by every open file while "does this file
* have genomes" is a property of the file — so filtering per file would mean
* building the descriptor per file, which §18.2's *a task is a descriptor,
* never a screen* refuses.
*
* The premise was right and the conclusion was not. The descriptor stays
* static; what it gained is a *rule* — `needsGenomes` — that the form reads.
* A cooler is offered "Chromosome" and "Whole genome" and nothing else,
* which is the honest list: it has no MAGs, and a reader shown "Genome"
* beside "Whole file" on a yeast cooler reasonably reads the first as the
* second. The case the old argument was really about — a stored session
* naming `genome` and reopened against a cooler — is still handled by
* {@link estimate} refusing it, and the form keeps showing it so the reader
* can see what the stored run was.
*/
options: [
{
value: 'contig',
label: 'contig',
coolLabel: 'chromosome',
hint: 'One curve per sequence.',
},
{
value: 'genome',
label: 'genome',
needsGenomes: true,
hint: 'One curve per MAG, over its own contigs.',
},
{
value: 'file',
label: 'whole file',
coolLabel: 'whole genome',
hint: 'Every sequence pooled into one curve.',
},
],
},
{
kind: 'sequences',
key: 'sequences',
label: 'sequences',
// No hint: a column of tick boxes that opens with every one ticked says
// what any sentence here would say.
// Meaningless under the other two scopes, where the sequences are decided by
// the genome table or by the file itself.
when: (params) => scopeOf(params) === 'contig',
},
{
kind: 'choice',
key: 'quantity',
label: 'quantity',
default: 'density',
/**
* Two things are called P(s) and the difference is one denominator.
*
* Both divide the contacts at a separation by something. The **density**
* divides by the contacts in the whole map and then by the separation the
* bin covers — the share of all cis contacts that fall at this distance, per
* base pair, which is `HiContacts::distanceLaw()`'s `norm_p`. The
* **expected** divides by the number of bin pairs that could have held a
* contact — `cooltools expected-cis`, the probability that two loci that far
* apart touch, and the reading a scaling exponent is compared against −1 or
* −1.5 as.
*
* So `density / expected = pairs(s) / constant`, and the pair count is flat
* while the separation is short against the sequences and falls to zero at
* the longest one. On the yeast demo the two are within 1 % below 10 kb,
* 13 % apart at 100 kb and a factor of fifty apart at 1.2 Mb — and on a
* mammalian map, where the chromosomes are two orders of magnitude longer,
* they agree over everything anybody fits. Which is why this is a choice
* rather than a correction: at the separations most figures are read at they
* are the same picture, and past that they answer different questions.
*/
options: [
{
value: 'density',
label: 'probability density',
hint: "Share of all cis contacts, per base pair. HiContacts' norm_p.",
},
{
value: 'expected',
label: 'expected (cooltools)',
hint: 'Mean contact per pair of bins that far apart. What cooltools expected-cis reports.',
},
],
},
{
kind: 'number',
key: 'factor',
label: 'log bins per step',
hint: 'Each separation bin is this many times wider than the last.',
default: 1.1,
min: 1.01,
max: 2,
step: 0.01,
// A bounded range whose only sensible use is sweeping it: the reader is
// trading the smoothness of the tail against the resolution of the shoulder,
// and that is a judgement made by watching the curve move.
slider: true,
},
{
kind: 'number',
key: 'minSeparation',
/**
* In **bins**, not base pairs.
*
* The near-diagonal is what this drops, and the near-diagonal is a fact
* about the *grid*: the first bin of separation is the one that holds every
* self-ligation and every unligated fragment, whatever resolution the file
* is at. A reader who wants "skip the first two bins" had to work out what
* two bins were in base pairs, type that, and get it back rounded — and the
* same typed number meant three bins at one resolution and one at the next.
* Bins are also what the curve's own x axis is quantised to, so this is now
* the same unit as the thing it is trimming.
*/
label: 'minimum separation',
hint: 'Drops the near-diagonal, where unligated fragments dominate.',
// **Two by default**: the diagonal and the one beside it, which hold what
// ligation could not separate — self-ligation and unligated fragments — and
// which cooltools' `expected_cis` leaves out by default for that reason
// (`ignore_diags=2`). Under the density it also floors the *total* every
// point is a share of, since the curve is a share of the contacts it was
// built from, so 0 keeps the diagonal in the first break, where HiContacts
// puts it. Under the expected the floor is one offset whatever this says:
// the diagonal is not a contact between two loci and has no place on a log
// axis. The four orientations are not floored by it (`run/expected.ts`).
default: 2,
min: 0,
integer: true,
unit: 'bins',
},
{
kind: 'toggle',
key: 'derivative',
label: 'log-log slope',
hint: 'The shoulder and the plateau are visible there and not on the curve.',
default: true,
},
// **The four orientations, retired from this form** (§18.8). They are about
// records rather than about the matrix, and they moved to Pairs QC, which
// reads the `.pairs` and draws them beside the rest of the library's
// statistics — so neither control is drawn here any more. Both stay
// parameters, because stored sessions, shared links and `chromosomix expected
// --stranded --pairsUri` lines carry them, and each of those must run as it
// did: the curve is still drawn when they ask for it, with a note saying
// where it lives now (`run/expected.ts`).
{
kind: 'toggle',
key: 'stranded',
label: 'four orientations',
default: false,
retired: true,
},
{
// The key keeps its old name, because stored sessions and `chromosomix
// --pairsUri` lines carry it. Still settled where the map's own import
// counted the four (§3.7), so a stored run on such a map is recorded as
// what it reads: the import's curves, and no file.
kind: 'file',
key: 'pairsUri',
label: '.pairs file',
accept: ['.pairs', '.pairs.gz'],
prompt: 'Drop a .pairs here, or click to choose one',
askedOf: (context) => context.binnedFrom?.orientation !== true,
settled: '',
retired: true,
},
BALANCE_PARAM,
]
// ---------------------------------------------------------------------------
// Cost
// ---------------------------------------------------------------------------
/**
* The bin range a scope covers, for the *is the block enough* question.
*
* Whole sequences, always: P(s) is a property of a sequence, and a curve
* computed from whichever window happened to be on screen but labelled `chrIV`
* is a different quantity wearing the same name. The tick boxes say which
* sequences, never how much of one.
*/
function scopeExtent(context: AnalysisContext, params: Params): [number, number] {
const whole: [number, number] = [0, context.nBins]
if (scopeOf(params) !== 'contig') return whole
const chosen = tickedContigs(context, params, 'sequences')
return chosen.length === 0 ? whole : extentOf(chosen)
}
/**
* Range requests for one contig's dense square.
*
* `MetaCooler.matrix` is "four slices and no scan" — the `bin1_offset` window,
* `bin1_id`, `bin2_id` and the value column — plus one more for the weight
* slice when balancing is on. On a local file these are `slice()` calls and cost
* nothing, which is what `requests` being a count rather than a duration is for.
*/
function boundedRequests(balance: boolean): number {
return balance ? 5 : 4
}
/**
* Ticked sequences read as their own squares, before a scan is the cheaper story.
*
* Four is where the two shapes cross in the only way that matters to a reader:
* a handful of bounded reads is a run they can start without a warning, and a
* list longer than that is on its way to being every sequence in the file —
* which the runner reads in one streaming pass whatever this says. The number
* decides which sentence is printed, never what is read.
*/
const SQUARES_WORTH_READING = 4
/**
* The run's price: the matrix's curve, and then wherever the four orientations
* are coming from.
*
* A refusal about the matrix wins, and it used to be skipped: a named `.pairs`
* was priced first and returned, so a form with no sequence ticked and a file
* named offered a Run that could only produce an empty figure.
*/
function estimate(context: AnalysisContext, params: Params): CostEstimate {
const matrix = matrixEstimate(context, params)
if (params.stranded !== true || matrix.refusal !== undefined) return matrix
// **A named `.pairs` dominates every other cost here**, so the whole run is
// priced as its pass. Reading a contact list is a pass over gigabytes of text;
// the curve it produces is a few hundred points.
const named = typeof params.pairsUri === 'string' ? params.pairsUri.trim() : ''
if (named !== '') {
return {
class: 'scan',
requests: null,
bytes: null,
note:
`${matrix.note} Then one pass over the whole of ${named} for the four orientations — ` +
'a contact list has no index to reach into.',
}
}
// The import counted them in its own pass (§3.7), so they cost nothing now —
// and the sentence says where they come from, since the form no longer asks.
const binned = context.binnedFrom
if (binned !== undefined && binned.orientation) {
return {
...matrix,
note: `${matrix.note} The four orientations were counted as ${binned.name} was binned.`,
}
}
// Reached only from a stored run or a command line now, since the form no
// longer asks (`retired`): so the sentence names where the question went
// rather than a control that is not there.
return {
...matrix,
note:
binned === undefined
? `${matrix.note} The four orientations need a .pairs, which Pairs QC asks for.`
: `${matrix.note} ${binned.name} has no strand columns, so the four orientations need ` +
'another .pairs, which Pairs QC asks for.',
}
}
/** What the matrix's own curve costs, which every run pays. */
function matrixEstimate(context: AnalysisContext, params: Params): CostEstimate {
const scope = scopeOf(params)
const balance = balanceOf(params)
const genomes = context.hasGenomes ? context.genomes : null
if (scope === 'genome' && (genomes === null || genomes.length === 0)) {
return {
class: 'free',
requests: 0,
bytes: 0,
// Costed as free because nothing is going to be read; the refusal below is
// what the drawer actually shows, and **Run** is disabled beside it.
note: 'Nothing to read.',
refusal: context.hasGenomes
? 'Every genome in this file is unassigned, and __unbinned__ is not a genome, ' +
'so there is nothing for a per-genome P(s) to be per. Use the contig or whole-file ' +
'scope.'
: 'This file is a plain cooler: it has no /genomes table, so there is nothing for a ' +
'per-genome P(s) to be per. Use the contig scope, which is the same curve over the ' +
'sequences this file does have, or the whole-file scope.',
}
}
// **Before the block is offered as free**, because nothing ticked is nothing to
// read even when the block is already on screen: a "no reads" price with no
// refusal beside it is a Run button that produces an empty figure.
if (scope === 'contig' && tickedContigs(context, params, 'sequences').length === 0) {
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Nothing to read.',
refusal: 'Tick at least one sequence.',
}
}
// The view's block is only an honest input when it covers the whole of what
// was asked for, was read at the file's own resolution, and was balanced the
// way this run wants. A pooled block is coarser than the file and a P(s) taken
// from it would be plotted against separations it cannot resolve — the framing
// records `binsPerCell` precisely so a task can decline to pretend otherwise.
const framing = context.framing
const free =
context.blockLoaded &&
framing !== null &&
framing.binsPerCell === 1 &&
framing.balance === balance &&
framingCovers(framing, scopeExtent(context, params))
if (free) {
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Computed from the block already on screen. No reads.',
}
}
if (scope === 'contig') {
const chosen = tickedContigs(context, params, 'sequences')
// **One square each, up to a handful.** A square per sequence is four
// requests per sequence, which is the cheap shape while the reader has
// ticked a few of them and is 154 k requests on a metagenome with every box
// ticked — so past a handful the streaming pass below is priced instead,
// which is the shape the runner takes there anyway.
if (chosen.length <= SQUARES_WORTH_READING) {
const requests = boundedRequests(balance) * chosen.length
const names = chosen.map((contig) => contig.name).join(', ')
const bins = chosen.reduce((sum, contig) => sum + contig.bins, 0)
return {
class: 'bounded',
requests,
bytes: null,
note:
`About ${requests} range requests: one dense square over ${names} ` +
`(${bins.toLocaleString()} bins).`,
}
}
}
// Everything left is more than one square, and reading it a square at a time
// would be one request per contig — 154 k of them on a real metagenome. One
// streaming pass over the pixel table is the cheaper shape as well as the
// honest one, and §18.6 says what that costs the reader in exchange: an
// explicit action, a progress bar and a cancel.
return {
class: 'scan',
requests: null,
bytes: null,
note:
'One pass over the whole pixel table, accumulating every cis contact by its separation.',
}
}
// ---------------------------------------------------------------------------
// The task
// ---------------------------------------------------------------------------
export const EXPECTED_TASK: AnalysisTask = {
id: 'expected',
title: 'Contact probability P(s)',
blurb: 'How contact frequency falls with genomic separation.',
keywords: ['P(s)', 'Ps', 'decay', 'distance', 'scaling', 'expected', 'slope'],
group: '2d',
venue: 'analysis',
version: 2,
// Every contact map has a diagonal, so this is the one row that is offered on
// everything — a cooler, a metacool, a file with nothing else loaded beside
// it. The guard is against a mount with no bins at all, which is a file the
// panel should be empty for.
availability: (context) => (context.hasMap ? READY : waitingFor(NEEDS_MAP)),
params: PARAMS,
outputs: ['curve'],
estimate,
}
|