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 | 2x 2x 2x 2x 18x 18x 2x 2x 2x 8x 8x 8x 8x 8x 8x 8x 8x 3x 3x 3x 3x 3x 3x 3x 8x 8x 3x 3x 3x 3x 3x 3x 3x 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 6x 6x 6x 8x 8x 8x 8x 2x 2x 2x 2x 2x 2x 2x 2x 1x 1x 2x 2x 6x 8x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 2x 2x 2x 2x 1x 1x 1x 1x 8x 1x 1x 1x 1x 1x 1x 1x 8x 8x 8x 8x 8x 8x 8x 8x 8x 8x 8x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x | /**
* The virtual 4C profile, as a *task* — which is a different object from the
* lane of the same name, and the panel is what makes that obvious (§18.4).
*
* The lane on the map is a reading of the block on screen. Its own docstring is
* careful about it: pan, and it is recomputed over the new block, because it is
* a summary of the pixels above it rather than a genome-wide track. That is the
* right thing for a lane and the wrong thing for a file — a bigWig that only
* covers whatever happened to be framed when it was written is a bigWig nobody
* can use.
*
* So this reads for itself, over a scope the reader chooses, and produces a
* track. And it is **cis by default**, because that choice is a factor of
* several thousand in what it costs:
*
* - **cis** is the anchor's own sequence, which is one bounded dense read — the
* same four slices any square costs, and free when the map is already showing
* it;
* - **with trans** is not a bounded read at all and the read planner says so. Set
* rows to the anchor band and columns to everything and `planRead` resolves to
* `symmetric` over the whole pixel table, because an upper-triangular store
* keeps every column *before* the anchor in other rows (§5.3, §6.2). Done
* properly it is one streaming pass that can stop at `bin1 >= anchor end` —
* on average half the file — accumulating into a band and materialising
* nothing large.
*
* One task with a toggle rather than two tasks, because they answer the same
* question and the toggle is the honest way to show what the second one costs.
*
* **On a metagenome there is no toggle**, and the run is always the second one.
* A MAG's contacts with the rest of the community are the reading a metagenomic
* map exists for, and the cheap half — one assembly against itself — is a
* picture of a genome pretending it was sequenced alone. The cost is stated the
* same way; what is gone is the offer to pay less for an answer to a different
* question. See {@link PARAMS}.
*/
import { contigNamed, framingCovers } from '../context'
import type {
AnalysisContext,
AnalysisTask,
CostEstimate,
ParamSpec,
Params,
} from '../types'
import { NEEDS_MAP, READY, waitingFor } from '../types'
import { BALANCE_PARAM } from './shared'
/** The anchor as typed, trimmed. */
export function anchorOf(params: Params): string {
return typeof params.anchor === 'string' ? params.anchor.trim() : ''
}
export function cisOnly(params: Params): boolean {
return params.cis !== false
}
/**
* What an anchor names, as bins, whether it is a sequence or a **genome**.
*
* The runner has always taken both — `resolveAnchor` reads a genome's extent
* and makes *cis* mean the whole MAG, deliberately, because on a metagenome the
* question "what does this organism touch" is the one being asked and answering
* for whichever contig the first bin fell in would answer a different one. This
* side did not: the estimate resolved the name against `/contigs` alone, so
* naming a MAG — which is what a click on a genome × genome map means, and what
* the map's own axis is labelled with — was refused before it could run.
*
* A genome is checked *after* a sequence, so a file where one name is both
* behaves as it always did.
*/
function anchorSpan(
context: AnalysisContext,
params: Params,
): { label: string; bins: number; extent: [number, number]; kind: 'contig' | 'genome' } | null {
const name = anchorOf(params).split(':')[0].trim()
if (name === '') return null
const contig = contigNamed(context, name)
if (contig !== undefined) {
return {
label: contig.name,
bins: contig.bins,
extent: [contig.startBin, contig.startBin + contig.bins],
kind: 'contig',
}
}
const genome = context.genomes?.find((it) => it.name === name)
if (genome === undefined) return null
return {
label: genome.name,
bins: genome.extent[1] - genome.extent[0],
extent: [genome.extent[0], genome.extent[1]],
kind: 'genome',
}
}
/**
* The anchor, whether to balance, and **two controls a metagenome does not get.**
*
* *Cis only* and *Normalisation* are real questions on a cooler. A chromosome is
* a self-contained thing, most of what a locus touches is on it, and the P(s)
* decay of that chromosome is the expectation to divide by; cis is then both the
* commonest reading and the cheap one, and the toggle is the honest way to show
* what leaving it costs.
*
* On a metacool neither question survives contact with the file. An anchor there
* is one MAG among hundreds sharing a tube, and the reading that is *wanted* —
* what this organism touches — is the trans one; a profile restricted to the
* anchor's own assembly answers a question about a genome in isolation, which is
* not why anybody opened a metagenomic map. And a decay curve accumulated from
* one MAG's contigs is not an expectation for contact *between* MAGs, so the
* ratio it would produce is a number with two meanings in one track. So both are
* absent there and the runner settles them: always trans, never O/E.
*
* `needsNoGenomes` rather than two descriptors, and rather than a `when` over
* some hidden parameter: the file decides, the rule is written where the control
* is, and the arithmetic behind both is untouched and still goldened.
*/
const PARAMS: readonly ParamSpec[] = [
{
kind: 'region',
key: 'anchor',
label: 'anchor',
pickable: true,
},
{
kind: 'choice',
key: 'normalise',
label: 'normalisation',
default: 'none',
needsNoGenomes: true,
settled: 'none',
/**
* A 4C profile is a slice of the map, and a slice of a Hi-C map is mostly
* P(s): contact falls by orders of magnitude with distance, so the profile's
* shape is dominated by *how far away* each bin is rather than by whether
* anything is happening there. On a log axis that is a clean triangle with
* the anchor at its apex, and the loops a reader ran the task to find are
* bumps on its flanks — visible near the anchor, invisible a megabase out.
*
* Dividing by the decay is the standard way to see past it, and the curve is
* taken from **the anchor's own sequence**, accumulated from the square this
* task is already reading. Not the P(s) task's genome-wide curve: a small
* chromosome divided by the average behaviour of the large ones comes out
* with a slope that is a property of the *set* rather than of this locus.
*
* `none` stays the default because the raw profile is what a bigWig of a 4C
* experiment holds and what a reader comparing against one expects — and
* because a ratio is a derived quantity that should be asked for.
*
* ## Coverage is the other bias, and it is a different one
*
* The decay is a property of *separation*; coverage is a property of the
* **target bin**. A bin in a repeat, a bin with few restriction sites, a bin
* the mapper half-lost: each is seen less by everything, so it is a dip in
* every profile drawn against it, whatever the anchor. On a raw map those
* dips are most of the small-scale structure a reader would otherwise take
* for real, and they are in exactly the places — repeats, segmental
* duplications — where a spurious interaction is most tempting to believe.
*
* Dividing by each bin's own total contact is the standard correction, and
* the divisor is **centred on one** so that it is a bias correction and not
* a change of units: the profile keeps the scale it had, and an O/E ratio
* corrected this way is still a ratio centred on one rather than one divided
* by a depth.
*
* It is offered *with* the decay rather than instead of it because the two
* are independent and a reader wanting a clean profile wants both — the
* decay out of the shape along the axis, the depth out of each bin across
* it. A balanced map has had something very like the second done to it
* already (that is what balancing is), so it is most useful on a raw one;
* the control does not hide itself there, because a reader balancing on a
* partly-balanced file is entitled to see what the second pass does.
*/
options: [
{
value: 'none',
label: 'none',
hint: 'Contacts as they were counted.',
},
{
value: 'coverage',
label: 'by coverage',
hint: 'Each bin divided by how much it is seen at all, so a dead bin is not a dip.',
},
{
value: 'ps',
label: 'by P(s)',
hint: 'Observed over expected at that separation, from this sequence’s own decay.',
},
{
value: 'ps+coverage',
label: 'by P(s) and coverage',
hint: 'Both divisions: the decay out of the shape, the depth out of each bin.',
},
],
},
{
kind: 'toggle',
key: 'cis',
label: 'cis only',
needsNoGenomes: true,
// What a metagenome's run *is*, so that is what its provenance says.
settled: false,
hint:
'On, this reads one sequence and costs a handful of requests; off, it walks the pixel ' +
'table for the contacts leaving that sequence.',
default: true,
},
/**
* The anchor against itself is the diagonal, and the bin beside it its
* shoulder — the brightest cells of the map by orders of magnitude, so in a
* track they are the spike that sets the scale every other bin is read on.
* Offered on every file: on a metagenome the anchor is a whole MAG, and what
* it touches *other* than itself is the reading this leaves.
*/
{
kind: 'toggle',
key: 'excludeAnchor',
label: 'exclude anchor region',
hint: 'Leave the anchor and one bin either side empty, so its contact with itself does not set the scale.',
default: false,
},
BALANCE_PARAM,
]
/**
* Whether this run stays on the anchor's own sequence.
*
* The file has the last word, and it is the same word the runner says: a
* metagenome is always read whole, whatever a stored session or a form on
* another file left in `cis`. See {@link PARAMS}.
*/
function cisRun(context: AnalysisContext, params: Params): boolean {
return !context.hasGenomes && cisOnly(params)
}
function estimate(context: AnalysisContext, params: Params): CostEstimate {
const anchor = anchorSpan(context, params)
if (anchorOf(params) === '') {
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Nothing to read.',
refusal: 'Name a locus to anchor the profile on — type one, or pick it on the map.',
}
}
if (anchor === null) {
const named = anchorOf(params).split(':')[0].trim()
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Nothing to read.',
refusal:
context.genomes === null
? `This file has no sequence called ${named}.`
: `This file has no sequence or genome called ${named}.`,
}
}
const balance = params.balance === true
// The metagenome's run is priced here too, and it is the same price: the walk
// is the walk whether a toggle asked for it or the file settled it. What must
// not happen is quoting four range requests for a pass the runner is going to
// make over half the file, which is what pricing `cis` from the params alone
// would do on a metacool — the key is still there, still `true` by default,
// and no longer the answer.
if (!cisRun(context, params)) {
// The coverage correction is a statement about every bin, so it is a
// statement about every pixel: the walk cannot stop at the anchor, and the
// number a reader decides on has to say that rather than quote the half-file
// pass it would otherwise have been (§18.6).
const wholeFile =
!context.hasGenomes && String(params.normalise ?? 'none').includes('coverage')
return {
class: 'scan',
requests: null,
bytes: null,
note: wholeFile
? 'One pass over the whole pixel table: the anchor reaches the contacts that leave this ' +
'sequence, and the coverage to divide by is a total over every bin.'
: 'One pass over the pixel table as far as the anchor — about half the file — to reach ' +
`the contacts that leave this ${context.hasGenomes ? 'genome' : 'sequence'}.`,
}
}
const extent = anchor.extent
const framing = context.framing
if (
context.blockLoaded &&
framing !== null &&
framing.binsPerCell === 1 &&
framing.balance === balance &&
framingCovers(framing, extent)
) {
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Computed from the block already on screen. No reads.',
}
}
const requests = balance ? 5 : 4
return {
class: 'bounded',
requests,
bytes: null,
note:
`About ${requests} range requests: one dense square over ` +
`${anchor.kind === 'genome' ? 'the whole of ' : ''}${anchor.label} ` +
`(${anchor.bins.toLocaleString()} bins).`,
}
}
export const FOURC_TASK: AnalysisTask = {
id: 'fourc',
title: '4C profile from an anchor',
blurb: 'What one locus touches, as a track you can export.',
keywords: ['virtual 4C', 'v4C', 'viewpoint', 'anchor'],
group: '2d',
venue: 'analysis',
version: 1,
availability: (context) => (context.hasMap ? READY : waitingFor(NEEDS_MAP)),
params: PARAMS,
outputs: ['profile'],
estimate,
}
|