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 | 2x 2x 2x 2x 2x 2x 2x 2x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 66x 66x 66x 66x 66x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 2x 2x 2x 2x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 132x 132x 132x 132x 66x 66x 102x 102x 102x 102x 102x 102x 102x 66x 4x 4x 4x 102x 102x 2710x 2608x 154x 77x 102x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 32x 32x 32x 9071x 9071x 3443835x 3443835x 3443835x 9071x 32x 32x 8x 8x 8x 8x 8x 618x 8x 8x 4x 4x 4x 4x 102x 102x 102x 102x 102x 102x 102x 102x 102x 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 102x 102x 102x 102x 102x 102x 102x 102x 102x 4x 4x 4x | /**
* TAD annotation, run: Arrowhead over each ticked sequence.
*
* `analysis/arrowhead.ts` is the method and does not know where a square comes
* from; this file hands it each window Juicer would read — out of the block on
* screen where that block can answer for it, from the file where it cannot —
* and turns what comes back into the two things a domain call leaves as: a set
* drawn as corners on the map, and Juicer's own BEDPE columns as a table.
*/
import { FeatureBuilder } from '../../annotations/featureset'
import { DEFAULT_COUNT_COLUMN } from '../../metacooler/matrix'
import { ARROWHEAD_DEFAULTS, arrowheadWindows, callDomains, type ArrowheadSettings, type HighScore } from '../arrowhead'
import { boolParam, numParam, tickedSequences, type TaskRunner } from '../job'
import { matrixSizeOf } from '../tasks/arrowhead'
import type { IntervalsOutput, TableOutput } from '../types'
import { missingSequencesNote, spanOfContig, squareFor, tickedContigIds, usableBlock, validBins, type Span } from './blocks'
import { formatBp } from './labels'
interface Unit {
name: string
span: Span
length: number
}
/** The domains of one sequence, as Juicer writes them: `x1` to `x2`, in base pairs. */
interface Domain extends HighScore {
sequence: string
start: number
end: number
}
export const runArrowhead: TaskRunner = async (job) => {
const mc = job.mc
const balance = boolParam(job.params, 'balance', false)
const value = job.framing.value === '' ? DEFAULT_COUNT_COLUMN : job.framing.value
const settings: ArrowheadSettings = {
matrixSize: matrixSizeOf(job.params),
varThreshold: numParam(job.params, 'varThreshold', ARROWHEAD_DEFAULTS.varThreshold),
highSignThreshold: numParam(job.params, 'highSignThreshold', ARROWHEAD_DEFAULTS.highSignThreshold),
maxLowSignThreshold: numParam(job.params, 'maxLowSignThreshold', ARROWHEAD_DEFAULTS.maxLowSignThreshold),
minLowSignThreshold: numParam(job.params, 'minLowSignThreshold', ARROWHEAD_DEFAULTS.minLowSignThreshold),
decrementLowSignThreshold: numParam(
job.params,
'decrementLowSignThreshold',
ARROWHEAD_DEFAULTS.decrementLowSignThreshold,
),
minBlockSize: Math.max(0, Math.round(numParam(job.params, 'minBlockSize', ARROWHEAD_DEFAULTS.minBlockSize))),
}
const ticked = tickedSequences(job.params, 'sequences')
if (ticked !== null && ticked.length === 0) {
job.note('No sequence is ticked, so there is nothing to call domains on.')
return []
}
const { ids, missing } = tickedContigIds(mc, ticked)
const lengths = mc.chromsizes
const units: Unit[] = []
for (const cid of ids) {
const span = spanOfContig(mc, cid)
if (span === null) continue
const name = mc.chromnames[cid]
units.push({ name, span, length: lengths.get(name) ?? (span.to - span.from) * mc.binsize })
}
const said = missingSequencesNote(mc, missing, units.length)
if (said !== null) job.note(said)
if (units.length === 0) return []
// What Juicer was tuned at, said rather than enforced: it picks 5 or 10 kb
// for itself, and a reader who opened a yeast map at 1 kb asked for 1 kb.
if (mc.binsize < 5_000 || mc.binsize > 25_000) {
job.note(
`Arrowhead was tuned on maps at 5 and 10 kb, and Juicer picks one of the two when it is not told; ` +
`this map is at ${formatBp(mc.binsize)}, which is what it was run at. Its window of ` +
`${settings.matrixSize.toLocaleString()} bins is ${formatBp(settings.matrixSize * mc.binsize)} here.`,
)
}
if (!balance) {
job.note(
'Read without balancing weights, which is Juicer’s -k NONE. Juicer itself reads a normalised map ' +
'unless told otherwise.',
)
}
// The counts, and the weights applied here rather than by the reader, so a
// cell is `NaN` only where it holds a contact and a bin of it was dropped —
// the window Juicer reads. A balanced read is `NaN` across every dropped
// row and column, empty cells included, and there the transform differs.
const block = usableBlock(
job,
units.map((unit) => unit.span),
false,
)
const weights = validBins(mc, balance)
const total = units.reduce((sum, unit) => sum + arrowheadWindows(unit.span.to - unit.span.from, settings.matrixSize).length, 0)
let done = 0
const domains: Domain[] = []
for (const unit of units) {
await job.breathe()
const bins = unit.span.to - unit.span.from
const before = done
const called = await callDomains(
bins,
mc.binsize,
async (from, to) => {
const first = unit.span.from + from
const square = await squareFor(
job,
{ from: first, to: unit.span.from + to },
{ block, balance: false, value },
)
return weights === null ? square.values : weighted(square.values, to - from, weights, first)
},
settings,
async (window, windows) => {
await job.breathe()
done = before + window
job.progress(done, total, `${unit.name}, window ${Math.min(window + 1, windows)} of ${windows}`)
},
)
for (const score of called) {
domains.push({
...score,
sequence: unit.name,
start: score.i,
end: Math.min(unit.length, score.j),
})
}
}
job.progress(total, total)
if (domains.length === 0) {
job.note(
'No domain was called. Arrowhead needs a deep map: Juicer warns that one with fewer than ' +
'about 300 million contacts may give few domains.',
)
return []
}
// How many sit inside another on their sequence: the one thing a tiling of
// boundaries cannot say, so worth a number.
const nested = domains.filter((inner) =>
domains.some(
(outer) =>
outer !== inner &&
outer.sequence === inner.sequence &&
outer.start <= inner.start &&
outer.end >= inner.end &&
(outer.start < inner.start || outer.end > inner.end),
),
).length
const count = domains.length
job.note(
`${count.toLocaleString()} ${count === 1 ? 'domain' : 'domains'} on ` +
`${new Set(domains.map((domain) => domain.sequence)).size.toLocaleString()} of ` +
`${units.length.toLocaleString()} ${units.length === 1 ? 'sequence' : 'sequences'}` +
(nested === 0 ? ', none inside another.' : `, ${nested.toLocaleString()} of them inside another.`),
)
return [domainSet(domains), domainTable(domains, mc.binsize)]
}
/**
* A square of counts balanced as the reader balances one — `count × w[i] × w[j]`,
* rounded to float32 — but `0` where there is no contact, whatever its weights.
*/
function weighted(raw: Float32Array, n: number, weights: Float64Array, first: number): Float32Array {
const out = new Float32Array(n * n)
for (let r = 0; r < n; r++) {
const wr = weights[first + r]
for (let c = 0; c < n; c++) {
const count = raw[r * n + c]
out[r * n + c] = count === 0 ? 0 : count * wr * weights[first + c]
}
}
return out
}
/** In genomic order: the sequence as the file has it, then where a domain starts, then where it ends. */
function genomic(domains: readonly Domain[]): Domain[] {
const rank = new Map<string, number>()
for (const domain of domains) if (!rank.has(domain.sequence)) rank.set(domain.sequence, rank.size)
return [...domains].sort(
(a, b) =>
(rank.get(a.sequence) ?? 0) - (rank.get(b.sequence) ?? 0) || a.start - b.start || a.end - b.end,
)
}
/** The domains as a set drawn on the map, one corner each. */
function domainSet(domains: readonly Domain[]): IntervalsOutput {
const builder = new FeatureBuilder()
let named = 0
for (const domain of genomic(domains)) {
if (domain.end <= domain.start) continue
builder.add(domain.sequence, domain.start, domain.end, {
// The corner score, which is what Juicer ranks a domain's evidence by.
score: domain.score,
name: `domain_${++named}`,
detail:
`corner score ${domain.score.toFixed(3)} · variance ${domain.uVarScore.toFixed(3)} up, ` +
`${domain.lVarScore.toFixed(3)} down · sign ${domain.upSign.toFixed(3)} up, ${domain.loSign.toFixed(3)} down`,
})
}
return {
kind: 'intervals',
id: 'arrowhead-domains',
label: 'Arrowhead domains',
payload: { kind: 'features', data: builder.finish() },
noun: { one: 'domain', many: 'domains' },
// On the map, as the corner of each square — the picture the corner score
// is read from, and where a nested domain is seen to be nested.
display: 'overlay',
look: 'corner',
}
}
/** Juicer's `_blocks.bedpe`, as a table: its columns, under its names. */
function domainTable(domains: readonly Domain[], binsize: number): TableOutput {
return {
kind: 'table',
id: 'arrowhead-table',
label: 'Arrowhead domains, as Juicer writes them',
columns: [
{ key: 'sequence', label: 'Sequence', kind: 'text' },
{ key: 'start', label: 'Start', kind: 'number', digits: 0, coordinate: 'start' },
{ key: 'end', label: 'End', kind: 'number', digits: 0, coordinate: 'end' },
{ key: 'width', label: 'Width (bins)', kind: 'number', digits: 0 },
{
key: 'score',
label: 'Corner score',
kind: 'number',
digits: 6,
hint: 'The arrowhead difference and sign difference of the two triangles, less their variance.',
},
{ key: 'uVarScore', label: 'uVarScore', kind: 'number', digits: 6 },
{ key: 'lVarScore', label: 'lVarScore', kind: 'number', digits: 6 },
{ key: 'upSign', label: 'upSign', kind: 'number', digits: 6 },
{ key: 'loSign', label: 'loSign', kind: 'number', digits: 6 },
],
rows: genomic(domains).map((domain) => ({
sequence: domain.sequence,
start: domain.start,
end: domain.end,
width: Math.round((domain.j - domain.i) / binsize),
score: domain.score,
uVarScore: domain.uVarScore,
lVarScore: domain.lVarScore,
upSign: domain.upSign,
loSign: domain.loSign,
})),
}
}
|