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 | 2x 2x 2x 2x 2x 2x 6x 6x 6x 6x 6x 1x 1x 1x 5x 6x 6x 6x 6x 6x 6x 6x 6x 6x 6x 6x 6x 5x 5x 5x 5x 6x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 4x 5x 1x 1x 1x 1x 1x 1x 5x 5x 5x 4x 4x 4x 4x 4x 4x 4x 6x 6x 6x 4x 4x 4x 4x 4x 4x 6x 6x 6x 6x 6x 6x 6x 6x 6x 6x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 6x 6x 6x 14x 10x 10x 10x 1630x 1630x 1630x 1534x 1534x 1534x 10x 10x 10x 10x 10x 10x 10x 10x 10x 14x 6x 4x 4x 4x 4x 9x 9x 6x 6x 6x 6x 6x 9x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 9x 6x 6x 6x 6x 6x 6x 9x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 9x 6x 6x 6x 3006x 3006x 3006x 1368x 1638x 3006x 6x 6x 6x 9x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 20x 20x 20x 20x 8x 8x 8x 8x | /**
* Reproducibility between the maps on screen, run.
*
* The arithmetic is `analysis/hicrep.ts`. What is here is the read, and the read
* is the whole of why this is affordable: a stratified correlation needs the
* cells within `maxDistance` of the diagonal and **nothing else**, which is a
* band rather than a square — the same shape insulation and the pattern detector
* already read, and the reason a mammalian chromosome is a few dozen bounded
* reads instead of a refusal on size.
*
* ## The order of the loops, and what it costs
*
* Sequences outermost, maps innermost. Every map's strata for **one** sequence
* are held at once, because an N×N coefficient matrix needs all of them
* together; nothing holds a second sequence. So the memory is
* `maps × strata × bins` and not `maps × strata × genome`, which is the
* difference between forty megabytes and a gigabyte on a human assembly.
*
* The peers are opened once for the whole run rather than once per sequence —
* mounting a remote cooler is a metadata walk §12 budgets at twenty requests,
* and paying that per chromosome would be most of the run. They are closed in a
* `finally`, because a cancelled job must not leave three files mounted in a
* worker that goes on living.
*/
import {
contigIdFor,
numParam,
boolParam,
tickedSequences,
type JobContext,
type TaskRunner,
} from '../job'
import { emptyStrata, smoothTileInto, stratumAdjusted, type Strata } from '../hicrep'
import { strataOf } from '../tasks/hicrep'
import { DEFAULT_COUNT_COLUMN } from '../../metacooler/matrix'
import type { MetaCooler } from '../../metacooler/MetaCooler'
import type { CurveOutput, MatrixOutput, TableOutput, TaskOutput } from '../types'
import { forEachDiagonalTile, missingSequencesNote, usableBlock, type Span } from './blocks'
/** One map in the comparison, with the name the figure calls it. */
interface Compared {
label: string
mc: MetaCooler
/** Whether this one is the mounted map, which is never closed here. */
lead: boolean
close?: () => void
}
export const runHicrep: TaskRunner = async (job) => {
const peers = job.inputs.peers ?? []
if (peers.length === 0) {
job.note(
'No second map is open. Press + along the edge of the map to open a cell beside it, and ' +
'drop another file in.',
)
return []
}
if (job.peer === undefined) {
throw new Error(
'Reproducibility reads a second file, and this venue cannot open one. Run it in the app, ' +
'or with `chromosomix hicrep <file> --against <other>`.',
)
}
// Unanswered is every sequence and an empty list is nothing (§18.2): a reader
// who pressed **none** is told so rather than handed every sequence of every
// map. The estimate refuses that run before it is sent; one that arrives
// anyway — a stored session, a hand-written line — hears it before a second
// file is mounted for it.
const ticked = tickedSequences(job.params, 'sequences')
if (ticked !== null && ticked.length === 0) {
job.note('No sequence is ticked, so there is nothing to compare over.')
return []
}
const balance = boolParam(job.params, 'balance', false)
const value = job.framing.value === '' ? DEFAULT_COUNT_COLUMN : job.framing.value
const h = Math.max(0, Math.round(numParam(job.params, 'smooth', 1)))
const maxK = strataOf(job.params)
// Named as the reader names it (§16.10), like the peers beside it: a batch
// runs this once per map, and a table that called each run's own map *This
// map* named a different map by the same words in every run.
const maps: Compared[] = [{ label: job.inputs.mapLabel ?? 'This map', mc: job.mc, lead: true }]
try {
for (let i = 0; i < peers.length; i++) {
job.progress(i, peers.length, `opening ${peers[i].label}`)
const opened = await job.peer(peers[i])
// **Refused by name, not compared anyway.** A stratum is a diagonal offset
// in bins, so the k-th diagonal of a 1 kb map and a 5 kb map are two
// different separations — and §16's lock is in base pairs precisely so a
// deck may hold both. Comparing them would report a low coefficient that
// reads as a biological difference.
if (opened.mc.binsize !== job.mc.binsize) {
job.note(
`${peers[i].label} is at ${opened.mc.binsize?.toLocaleString() ?? '?'} bp and this map ` +
`is at ${job.mc.binsize.toLocaleString()} bp, so their strata are different ` +
'separations and they are left out. Open them at one resolution.',
)
opened.close()
continue
}
maps.push({ label: peers[i].label, mc: opened.mc, lead: false, close: opened.close })
}
if (maps.length < 2) {
job.note('No map beside this one could be compared with it, so there is nothing to report.')
return []
}
return await compare(job, maps, { balance, value, h, maxK })
} finally {
for (const map of maps) map.close?.()
}
}
async function compare(
job: JobContext,
maps: Compared[],
options: { balance: boolean; value: string; h: number; maxK: number },
): Promise<TaskOutput[]> {
const { balance, value, h, maxK } = options
// Every sequence where the run was never asked, and the ticked ones this
// file has where it was; an empty list has already been turned away in
// `runHicrep`. A ticked name this file has not got is named and left out
// (§5.9) rather than dropped in silence — and a list of nothing but those is
// said as that, not as nothing ticked.
const wanted = tickedSequences(job.params, 'sequences')
const chosen: number[] = []
const absent: string[] = []
if (wanted === null) {
for (let cid = 0; cid < job.mc.nContigs; cid++) chosen.push(cid)
} else {
for (const name of wanted) {
const cid = contigIdFor(job.mc, name)
if (cid >= 0) chosen.push(cid)
else absent.push(name)
}
}
const said = missingSequencesNote(job.mc, absent, chosen.length)
if (said !== null) job.note(said)
if (chosen.length === 0) return []
const n = maps.length
// Pooled across sequences by the same weights the strata are pooled by, which
// is what makes the genome-wide number a weighted mean of the per-sequence
// ones rather than an unweighted average over chromosomes of wildly different
// length. A mean of per-chromosome coefficients gives a 200 kb contig the same
// say as chromosome 1, which on a metagenome is every contig.
const pooledWeighted = new Float64Array(n * n)
const pooledWeights = new Float64Array(n * n)
const perSequence: Array<Record<string, string | number>> = []
// Summed across sequences so the curve is one line per pair rather than one
// per pair per chromosome.
const strataWeighted: Float64Array[] = []
const strataWeights: Float64Array[] = []
for (let p = 0; p < (n * (n - 1)) / 2; p++) {
strataWeighted.push(new Float64Array(maxK + 1))
strataWeights.push(new Float64Array(maxK + 1))
}
// The view's block, where it honestly covers what is being read — and for the
// **lead only**, because the block on screen is one panel's pixels and the
// peers are other files (`AnalysisRunner`'s own note about a deck sharing one
// block makes the same point one layer up).
const leadBlock = usableBlock(
job,
chosen.map((cid) => spanOf(job.mc, cid)),
balance,
)
for (let at = 0; at < chosen.length; at++) {
await job.breathe()
const cid = chosen[at]
const name = job.mc.chromnames[cid]
job.progress(at, chosen.length, name)
const span = spanOf(job.mc, cid)
const bins = span.to - span.from
if (bins < 3) continue
// One sequence, every map. Held together because the coefficient matrix
// needs them together, and freed before the next sequence is read.
const strata: Strata[] = []
let missing = false
for (const map of maps) {
const theirs = map.lead ? span : spanIn(map.mc, name)
if (theirs === null) {
missing = true
break
}
if (theirs.to - theirs.from !== bins) {
// Same name, different length: two assemblies that share a chromosome's
// name and not its sequence. Left out with a sentence rather than
// compared over whichever is shorter.
job.note(
`${name} is ${bins.toLocaleString()} bins here and ` +
`${(theirs.to - theirs.from).toLocaleString()} in ${map.label}, so it is left out.`,
)
missing = true
break
}
strata.push(
await stratifyOne(job, map, theirs, {
balance,
value,
h,
maxK: Math.min(maxK, bins - 1),
bins,
block: map.lead ? leadBlock : null,
}),
)
}
if (missing) continue
let pair = 0
for (let a = 0; a < n; a++) {
for (let b = a + 1; b < n; b++) {
const result = stratumAdjusted(strata[a], strata[b])
const key = a * n + b
for (const stratum of result.strata) {
// Pooled as each sequence's own number is: a stratum with no
// correlation counts as zero, with its weight.
pooledWeighted[key] += stratum.weight * (Number.isFinite(stratum.r) ? stratum.r : 0)
pooledWeights[key] += stratum.weight
// The curve is the correlation *at* each separation, so a stratum that
// has none is left out of it rather than drawn as a zero.
if (!Number.isFinite(stratum.r)) continue
strataWeighted[pair][stratum.k] += stratum.weight * stratum.r
strataWeights[pair][stratum.k] += stratum.weight
}
perSequence.push({
sequence: name,
a: maps[a].label,
b: maps[b].label,
scc: result.scc,
strata: result.strata.length,
})
pair++
}
}
}
job.progress(chosen.length, chosen.length, 'done')
job.throwIfCancelled()
const values = new Float32Array(n * n)
for (let a = 0; a < n; a++) {
values[a * n + a] = 1
for (let b = a + 1; b < n; b++) {
const key = a * n + b
const scc = pooledWeights[key] > 0 ? pooledWeighted[key] / pooledWeights[key] : Number.NaN
values[a * n + b] = scc
values[b * n + a] = scc
}
}
const labels = maps.map((map) => map.label)
const outputs: TaskOutput[] = [
{
kind: 'matrix',
id: 'hicrep',
label: 'Stratum-adjusted correlation',
nrows: n,
ncols: n,
values,
xLabel: '',
yLabel: '',
xTicks: labels,
yTicks: labels,
colormap: 'RdBu_r',
// A correlation is read about no-agreement, and `[0, 1]` is fixed rather
// than fitted so that two runs on two pairs of files are the same picture
// where they are the same result.
range: [0, 1],
stats: [{ label: 'sequences', value: chosen.length, digits: 0 }],
// Every map against every other, whichever map this run was about — so
// the all-maps tab shows it once rather than once per map (§18.7).
acrossMaps: true,
} satisfies MatrixOutput,
]
const pooledRows: Array<Record<string, string | number>> = []
for (let a = 0; a < n; a++) {
for (let b = a + 1; b < n; b++) {
pooledRows.push({
a: labels[a],
b: labels[b],
scc: values[a * n + b],
})
}
}
outputs.push({
kind: 'table',
id: 'hicrep-pairs',
label: 'Pooled over every sequence',
columns: [
{ key: 'a', label: 'Map', kind: 'text' },
{ key: 'b', label: 'Map', kind: 'text' },
{
key: 'scc',
label: 'SCC',
kind: 'number',
digits: 4,
hint: 'Every stratum of every sequence, each weighted by the positions it compares.',
},
],
rows: pooledRows.sort((x, y) => Number(y.scc) - Number(x.scc)),
acrossMaps: true,
} satisfies TableOutput)
if (perSequence.length > 0) {
outputs.push({
kind: 'table',
id: 'hicrep-sequences',
label: 'Per sequence',
columns: [
{ key: 'sequence', label: 'Sequence', kind: 'text' },
{ key: 'a', label: 'Map', kind: 'text' },
{ key: 'b', label: 'Map', kind: 'text' },
{ key: 'scc', label: 'SCC', kind: 'number', digits: 4 },
{ key: 'strata', label: 'Strata', kind: 'number', digits: 0 },
],
rows: perSequence,
acrossMaps: true,
} satisfies TableOutput)
}
// The figure that says more than the scalar does: *where* two maps stop
// agreeing. Short-range agreement with long-range disagreement is a depth
// difference; disagreement everywhere is two different experiments.
const series = []
let pair = 0
for (let a = 0; a < n; a++) {
for (let b = a + 1; b < n; b++) {
const x = new Float64Array(maxK + 1)
const y = new Float64Array(maxK + 1)
for (let k = 0; k <= maxK; k++) {
x[k] = job.toBp(k)
y[k] =
strataWeights[pair][k] > 0
? strataWeighted[pair][k] / strataWeights[pair][k]
: Number.NaN
}
series.push({ name: `${labels[a]} vs ${labels[b]}`, x, y })
pair++
}
}
if (series.length > 0) {
outputs.push({
kind: 'curve',
id: 'hicrep-strata',
label: 'Correlation by separation',
series,
xLabel: 'Genomic separation (bp)',
yLabel: 'Pearson r on that stratum',
xScale: 'log',
yScale: 'linear',
yDomain: [-0.2, 1],
acrossMaps: true,
} satisfies CurveOutput)
}
return outputs
}
/** One map's smoothed strata for one sequence. */
async function stratifyOne(
job: JobContext,
map: Compared,
span: Span,
options: {
balance: boolean
value: string
h: number
maxK: number
bins: number
block: Parameters<typeof forEachDiagonalTile>[2]['block']
},
): Promise<Strata> {
const into = emptyStrata(options.bins, options.maxK)
// **The kernel's whole reach**, which is the band plus the filter: a cell on
// the `maxK`-th diagonal is smoothed from cells `h` further out still, so a
// halo of `maxK` alone would smooth the outermost stratum from a window that
// is half outside the tile.
const halo = options.maxK + options.h
// A shim rather than a second reader: `forEachDiagonalTile` needs a job only
// for its `mc` and its cancellation, and handing it one pointing at the peer
// is what lets the two maps be read by one piece of code.
const shim: JobContext = map.lead ? job : { ...job, mc: map.mc }
await forEachDiagonalTile(
shim,
span,
{ halo, block: options.block, balance: options.balance, value: options.value },
(tile) => {
job.throwIfCancelled()
smoothTileInto(tile, options.h, options.maxK, into)
},
)
return into
}
function spanOf(mc: MetaCooler, cid: number): Span {
const [from, to] = mc.contigExtent(cid)
return { from, to }
}
/** The same sequence in another file, resolved by name (§5.9), or `null`. */
function spanIn(mc: MetaCooler, name: string): Span | null {
const cid = contigIdFor(mc, name)
return cid < 0 ? null : spanOf(mc, cid)
}
|