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 | 2x 2x 2x 2x 2x 2x 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 6x 6x 6x 6x 6x 1x 1x 1x 1x 1x 1x 1x 1x 5x 6x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 2x 2x 2x 2x 6x 6x 6x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 6x 6x 6x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x | /**
* The directionality index, as a descriptor.
*
* It sits directly after the diamond because it is the *other* answer to the
* question the diamond asks, and a reader who has one wants the other. Both are
* per-bin readings of a domain edge; they disagree about what an edge is, and the
* disagreement is informative rather than a nuisance. Insulation says *little
* contact crosses here* and finds the boundary whatever is on either side of it.
* The directionality index says *this bin looks one way* and finds the boundary
* from the asymmetry, which means it also marks the **inside** edges of a domain
* and gives the sign flip that Dixon's HMM segmented on. A boundary both agree
* about is a boundary.
*
* ## The window is in bins, like the diamond's
*
* Both tasks look a bounded distance up the diagonal and down it, and both are
* read side by side against the same map — so a reader comparing them is
* comparing two windows, and a form that asks for one in bins and the other in
* base pairs makes that comparison arithmetic. The diamond's argument for bins
* applies here unchanged: the window is a count of cells the statistic sums, so
* the number that decides what it measures has always been this one, and a
* reader who wants the two tasks at the same scale can now type the same number
* into both.
*
* The cost is the one bins always carry: **the default is a different genomic
* distance on every rung of an `.mcool`**, where Dixon's 2 Mb was the same
* distance everywhere. {@link DEFAULT_DI_WINDOW_BINS} is 50, which is that 2 Mb
* at the 40 kb the 2012 paper worked at, and the runner says in a note what the
* window came to in base pairs on the file actually open — because a 50-bin
* window on a 1 kb bacterial map is 50 kb and is a different reading of the
* genome, and that is worth being told once rather than inferring from a track.
*
* ## Borders and domains
*
* The same two outputs the diamond produces, called from the same kernel: see
* `analysis/di.ts` for why the index is integrated first and what the prominence
* is measured in. A border here is a **sign flip** — the bin where the track
* stops looking upstream and starts looking downstream — where the diamond's is
* a minimum, so the two sets are independent readings of the same edges and a
* border both agree on is worth more than either.
*/
import { extentOf, framingCovers, tickedContigs } from '../context'
import { bandTiles } from '../tiles'
import type {
AnalysisContext,
AnalysisTask,
CostEstimate,
ParamSpec,
Params,
} from '../types'
import { NEEDS_MAP, READY, waitingFor } from '../types'
import { BALANCE_PARAM } from './shared'
/**
* Dixon et al. 2012's 2 Mb, in bins of the 40 kb map the paper computed it on.
*
* A default has to be a number rather than a distance now that the form asks in
* bins, and this is the number that reproduces the published reading on a
* mammalian map opened at its usual resolution.
*/
export const DEFAULT_DI_WINDOW_BINS = 50
/** The window in bins, at least one. */
export function windowBinsOf(params: Params): number {
const value = typeof params.window === 'number' ? params.window : DEFAULT_DI_WINDOW_BINS
return Math.max(1, Math.round(value))
}
const PARAMS: readonly ParamSpec[] = [
{
kind: 'number',
key: 'window',
label: 'window',
hint:
'Contacts are summed this many bins each way of every bin, 50 being Dixon et al.’s 2 Mb ' +
'at the 40 kb they worked at.',
default: DEFAULT_DI_WINDOW_BINS,
min: 1,
integer: true,
unit: 'bins',
},
{
kind: 'sequences',
key: 'sequences',
label: 'sequences',
},
{
kind: 'number',
key: 'ignoreDiags',
label: 'diagonals to ignore',
/**
* `0`, where the diamond's is `2`, and `analysis/di.ts` argues the
* difference: the index differences two sums that both contain the
* near-diagonal, so it largely cancels, and Dixon's own numbers were
* computed with every cell in.
*/
hint: '0 counts every cell; the bin’s contact with itself is never in either sum.',
default: 0,
min: 0,
integer: true,
section: 'Quality',
},
{
kind: 'number',
key: 'minValid',
label: 'minimum coverage of each side',
hint: 'Below this fraction of usable cells — on either side — the bin gets no score.',
default: 0.5,
min: 0,
max: 1,
step: 0.05,
slider: true,
section: 'Quality',
},
BALANCE_PARAM,
/**
* The borders, and the domains between them.
*
* On by default, for the reason the diamond's are: the index is a reading and
* a reader almost always wants the call under it, the calls cost nothing once
* the track exists, and producing only the track left the same eight lines of
* script to be written every time. A reader who wants the curve alone unticks
* it.
*/
{
kind: 'toggle',
key: 'call',
label: 'call borders',
hint: 'Adds the borders as a set of their own, a 1D track under the index to begin with.',
default: true,
section: 'Borders',
},
{
kind: 'number',
key: 'prominence',
label: 'minimum prominence',
/**
* In bins of ordinary directional bias, which is the unit
* `integratedDirection` standardises the curve into and is not the index's
* own. A χ² grows with how deep the library was sequenced, so a
* threshold in χ² is a different threshold on every file and a reader who
* found a good one on one map would look for it again from scratch on the
* next. This one is a property of the structure: *how many bins of typical
* bias does the arm on each side of this border carry*.
*
* Ten by default — a domain arm of about ten bins at full strength, which at
* the 50-bin window above is a modest domain and calls a set a reader can
* then thin rather than one they have to go looking for.
*/
hint: 'How many bins of steady bias the track must carry on both sides of a flip.',
default: 10,
min: 0,
max: 200,
step: 1,
slider: true,
when: (params) => params.call === true,
section: 'Borders',
},
{
kind: 'number',
key: 'spacing',
label: 'minimum spacing',
hint: 'Two calls closer than this many bins keep only the stronger.',
default: 0,
min: 0,
integer: true,
unit: 'bins',
when: (params) => params.call === true,
section: 'Borders',
},
{
kind: 'toggle',
key: 'domains',
label: 'domains between them',
hint: 'A second set, the span from each border to the next, drawn as corners along the diagonal to begin with.',
default: true,
when: (params) => params.call === true,
section: 'Borders',
},
]
function estimate(context: AnalysisContext, params: Params): CostEstimate {
const balance = params.balance === true
const window = windowBinsOf(params)
const chosen = tickedContigs(context, params, 'sequences')
if (chosen.length === 0) {
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Nothing to read.',
refusal: 'Tick at least one sequence.',
}
}
// **Priced over the sequences that can carry the window, not over the ticks.**
// A bin needs a whole window each side, so a sequence no longer than twice
// the window contributes nothing at any cost — the runner's own test
// (`run/di.ts`), on a line and on a circle alike — and the run says so rather
// than reading it. With one ticked that is a refusal naming it, which is the
// same sentence the single dropdown used to produce; with several it counts
// them, since a metagenome's short contigs are the ordinary case and not a
// mistake.
const long = chosen.filter((entry) => entry.bins > window * 2)
if (long.length === 0) {
const one = chosen.length === 1 ? chosen[0] : null
return {
class: 'free',
requests: 0,
bytes: 0,
note: 'Nothing to read.',
refusal:
one !== null
? `${one.name} is ${one.bins.toLocaleString()} bins and the window is ` +
`${window.toLocaleString()}, so no bin on it has a whole window either side: it needs ` +
`more than ${(window * 2).toLocaleString()}. Use a narrower window, or a longer sequence.`
: `No ticked sequence is longer than twice a ${window.toLocaleString()}-bin window, so no ` +
'bin on any of them has a whole window either side. Use a narrower window, or tick ' +
'a longer sequence.',
}
}
const extent = extentOf(long)
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.',
}
}
// **A band, not a square**, which is what makes a 2 Mb window on a mammalian
// chromosome a sane thing to ask for at all: the index reaches `window` bins
// each way and no further, so what is read is `(window + 1) × n` cells rather
// than `n²`. `bandFor` tiles it; the count below is that walk's, and it is the
// same arithmetic `analysis/run/blocks.ts` does.
const perTile = balance ? 5 : 4
if (long.length === 1) {
const contig = long[0]
const tiles = bandTiles(contig.bins, window)
const requests = tiles * perTile
return {
class: 'bounded',
requests,
bytes: null,
note:
`About ${requests.toLocaleString()} range requests: ${contig.name} is read as a band ` +
`${window.toLocaleString()} bins deep, in ${tiles.toLocaleString()} ` +
`${tiles === 1 ? 'strip' : 'strips'}, never as the whole square.`,
}
}
const tiles = long.reduce((sum, entry) => sum + bandTiles(entry.bins, window), 0)
return {
class: 'scan',
requests: tiles * perTile,
bytes: null,
note:
`${long.length.toLocaleString()} sequences are long enough for this window, read as ` +
`${tiles.toLocaleString()} ${tiles === 1 ? 'strip' : 'strips'} along their diagonals.`,
}
}
export const DI_TASK: AnalysisTask = {
id: 'di',
title: 'Directionality index',
blurb: 'Which way each bin looks, and the domains the sign flips call.',
keywords: ['DI', 'TAD', 'domains', 'boundaries'],
group: '2d',
venue: 'analysis',
version: 1,
availability: (context) => (context.hasMap ? READY : waitingFor(NEEDS_MAP)),
params: PARAMS,
outputs: ['profile', 'intervals'],
estimate,
}
|