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 | 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 4x 4x 4x 4x 4x 4x 4x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x | /**
* The HiC score, run.
*
* The same plan as the directionality index beside it — a band per sequence, one
* kernel, one lane — and the shallowest read of the three: the score is the
* bin's own line across the near-diagonal, so the band is read exactly as deep
* as the line reaches, `diagonals − 1`. At the bacterial default that is two
* bins, which on most files is one strip and no tiling at all.
*
* What this file owns beyond that is the arithmetic between the reader's number
* and the kernel's: a count of diagonals becomes a reach in bins, and the notes
* say what that came to in base pairs on the file actually open and how many
* cells each score is a mean of. Both are one line and both are things a reader
* would otherwise have to work out from the label.
*/
import { DEFAULT_COUNT_COLUMN } from '../../metacooler/matrix'
import { hicScore, lineCells } from '../hicscore'
import {
boolParam,
numParam,
readsAsCircular,
tickedSequences,
type JobContext,
type TaskRunner,
} from '../job'
import { diagonalsOf, reachOf } from '../tasks/hicscore'
import type { ProfileOutput, TaskOutput } from '../types'
import {
bandFor,
missingSequencesNote,
spanOfContig,
tickedContigIds,
usableBlock,
type Span,
} from './blocks'
import { formatBp } from './labels'
interface Unit {
name: string
span: Span
circular: boolean
}
/** The sequences to score: the ticked ones, or every one where none was. */
function planUnits(job: JobContext, names: readonly string[] | null): {
units: Unit[]
missing: string[]
} {
const mc = job.mc
const { ids, missing } = tickedContigIds(mc, names)
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, circular: readsAsCircular(job, name) })
}
return { units, missing }
}
export const runHicScore: TaskRunner = async (job) => {
const balance = boolParam(job.params, 'balance', false)
const value = job.framing.value === '' ? DEFAULT_COUNT_COLUMN : job.framing.value
const minValid = Math.max(0, Math.min(1, numParam(job.params, 'minValid', 0.5)))
const ignoreDiags = Math.max(0, Math.round(numParam(job.params, 'ignoreDiags', 0)))
const diagonals = diagonalsOf(job.params)
const depth = reachOf(job.params)
const reachBp = depth * job.mc.binsize
const cells = lineCells(diagonals, ignoreDiags)
// What the score reaches in base pairs on *this* file, and how many cells are
// behind each point of the track. The form asks in diagonals so the number a
// reader types is the number the method is described with (see
// `tasks/hicscore.ts`), and the cost of that is exactly this line.
const included = ignoreDiags > 0 ? 'less the diagonals cut below' : 'the main diagonal included'
job.note(
`${diagonals.toLocaleString()} ${diagonals === 1 ? 'diagonal reaches' : 'diagonals reach'} ` +
`${formatBp(reachBp)} either way at this file’s ` +
`${job.mc.binsize.toLocaleString()} bp resolution, and each score is the mean of ` +
`${cells.toLocaleString()} ${cells === 1 ? 'cell' : 'cells'} — the bin’s own line ` +
`across them, ${included}.`,
)
if (cells === 0) {
job.note(
`Ignoring ${ignoreDiags.toLocaleString()} diagonals leaves nothing of a ` +
`${diagonals.toLocaleString()}-diagonal line, so there is no cell left to average. ` +
'Count more diagonals, or ignore fewer.',
)
return []
}
const ticked = tickedSequences(job.params, 'sequences')
if (ticked !== null && ticked.length === 0) {
job.note('No sequence is ticked, so there is nothing to score.')
return []
}
const { units: planned, missing } = planUnits(job, ticked)
const said = missingSequencesNote(job.mc, missing, planned.length)
if (said !== null) job.note(said)
if (planned.length === 0) return []
// The line spans `2 × depth + 1` bins, so a sequence shorter than that has no
// bin with a whole line through it — on a circle as on a line, because past
// half the circle the short way round is the other way.
const units = planned.filter((unit) => unit.span.to - unit.span.from > 2 * depth)
if (units.length === 0) {
job.note(
`No sequence in this file is longer than the ${(2 * depth + 1).toLocaleString()} bins this ` +
'line spans, so no bin on any of them has a whole line through it. Use fewer diagonals.',
)
return []
}
const nBins = job.mc.nBins
const scores = new Float32Array(nBins).fill(Number.NaN)
const block = usableBlock(
job,
units.map((unit) => unit.span),
balance,
)
let circularUnits = 0
for (let u = 0; u < units.length; u++) {
const unit = units[u]
await job.breathe()
job.progress(u, units.length, unit.name)
if (unit.circular) circularUnits++
const band = await bandFor(job, unit.span, {
hi: depth,
block,
balance,
value,
circular: unit.circular,
onTile: async (done, total) => {
await job.breathe()
job.progress(
u,
units.length,
total === 1 ? unit.name : `${unit.name} · strip ${Math.min(done + 1, total)} of ${total}`,
)
},
})
const local = hicScore(band, { diagonals, minValid, ignoreDiags })
for (let i = 0; i < local.length; i++) scores[unit.span.from + i] = local[i]
}
job.progress(units.length, units.length, 'done')
if (circularUnits > 0) {
job.note(
`${circularUnits.toLocaleString()} of ${units.length.toLocaleString()} ` +
`${units.length === 1 ? 'sequence was' : 'sequences were'} read as circular, so the ` +
'line wraps past the end and every bin is scored.',
)
}
if (circularUnits < units.length && depth > 0) {
// The other half of the same fact, and the one a reader is more likely to be
// puzzled by: a linear sequence has no whole line through its first and last
// few bins, so the track begins and ends with a gap rather than a fade.
job.note(
`The first and last ${depth.toLocaleString()} ${depth === 1 ? 'bin' : 'bins'} of a linear ` +
'sequence have no whole line through them and are left as gaps, rather than scored on ' +
'the part of it that fits — which would fade both ends towards zero.',
)
}
const profile: ProfileOutput = {
kind: 'profile',
id: `hicscore-d${diagonals}`,
label: `HiC score · ${diagonals} diag`,
values: scores,
extent: [0, nBins],
binsize: job.mc.binsize,
// A mean of cells, so it carries the same unit one cell does — counts on a
// raw map, and a balanced map's own dimensionless numbers where the run was
// balanced. `contacts` is what every other lane of cell values in the
// catalogue says and what the axis should be read as; that this one is an
// average rather than a total is in the notes and in the label's window.
unit: 'contacts',
}
const outputs: TaskOutput[] = [profile]
return outputs
}
|