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 | 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x | /**
* The pileup `chromosight quantify` draws, from the APA task.
*
* A fourth normalisation rather than a fourth task, because it *is* an aggregate
* peak analysis — the same anchors, the same averaging, the same corner
* statistic. What differs is the control, and that difference is the whole
* reason it is offered: the other three divide the pileup by something computed
* from the windows themselves (a shifted control, an expected curve), while this
* one detrends **the map** by its own distance law first and then averages
* windows that are already ratios. So the number in a cell is
* *observed / expected at that separation on that sequence*, which is what
* chromosight's figures are, and a reader can put the two side by side.
*
* Three things follow from that and none of them is optional:
*
* - **the window is the template's size**, because that is what a quantify
* window is — 17 bins for loops — and a 51-bin window detrended the same way
* is a different picture with the same name;
* - **the scan distance is the loop set's own reach**, which is `cmd_quantify`'s
* first act. The distance law is computed out to the furthest pair in the file
* rather than to the preset's two megabases, so a set of short-range loops is
* detrended against a curve fitted to short ranges;
* - **every anchor comes back**, scored or not. Quantification answers about the
* positions it was given; dropping the ones whose window is too empty would
* silently reorder somebody's loop table.
*
* And it scores as well as piles up. A pileup says what the set looks like on
* average; the per-anchor Pearson correlation with the template says which
* members of it look like that at all, which is the difference between a loop
* set and a list of coordinates.
*/
import { p2ll } from '../apa'
import {
detectableFromMarginal,
detectableFromWeights,
fdrCorrection,
pileupOf,
keepDistanceFor,
prepareSquare,
quantifyOnSquare,
symmetrise,
templatesOf,
type ScanSettings,
type Square,
} from '../chromosight'
import { CHROMOSIGHT_PRESETS, type ChromosightPreset } from '../chromosightPresets'
import { CHROMOSIGHT_SPECTRUM } from '../chromosightSpectrum'
import { boolParam, strParam, type JobContext } from '../job'
import { loopSpectrum } from '../loopSpectrum'
import { DEFAULT_COUNT_COLUMN } from '../../metacooler/matrix'
import type { CurveOutput, MatrixOutput, TableOutput, TaskOutput } from '../types'
import { bandFor, cisMarginal, usableBlock, validBins, type Span } from './blocks'
import { offsetTicks } from './labels'
/** What the APA runner hands over: an anchor pair in whole-file bins. */
export interface ChromosightAnchor {
row: number
col: number
bounds: Span
label: string
}
/** The preset an APA run is scoring against. */
export function apaPreset(job: JobContext): ChromosightPreset {
const name = strParam(job.params, 'csPattern', 'loops')
return CHROMOSIGHT_PRESETS.find((preset) => preset.name === name) ?? CHROMOSIGHT_PRESETS[0]
}
export async function chromosightPileup(
job: JobContext,
anchors: readonly ChromosightAnchor[],
): Promise<TaskOutput[]> {
const preset = apaPreset(job)
const kernel = templatesOf(preset)[0]
const balance = boolParam(job.params, 'balance', false)
const value = job.framing.value === '' ? DEFAULT_COUNT_COLUMN : job.framing.value
const binsize = job.mc.binsize
// `cmd_quantify`'s first act: the scan distance becomes the furthest pair in
// the set, so the distance law is fitted over the separations that are
// actually being asked about.
const furthest = anchors.reduce((most, anchor) => Math.max(most, anchor.col - anchor.row), 0)
const maxDist = Math.max(Math.min(furthest, job.mc.nBins), 1)
const settings: ScanSettings = {
maxDist,
pearson: preset.pearson,
missingTol: preset.maxPercUndetected / 100,
zeroTol: preset.maxPercZero / 100,
onDiagonal: preset.onDiagonal,
}
// Every anchor, with no cap: which anchors were scored is part of what the
// score means, and a subset chosen by whichever sorted first is not a subset
// anybody asked for.
const used = anchors
// One square per sequence rather than one read per window, and this is the one
// normalisation where that is required rather than merely cheaper: the
// distance law is a property of the whole sequence, and no set of windows —
// however many — can be detrended against a curve taken from itself.
const bySequence = new Map<number, ChromosightAnchor[]>()
for (const anchor of used) {
const held = bySequence.get(anchor.bounds.from)
if (held === undefined) bySequence.set(anchor.bounds.from, [anchor])
else held.push(anchor)
}
const weights = validBins(job.mc, balance)
// The one quantity a band cannot answer for itself, counted once and only
// where there are no weights. See `run/chromosight.ts`.
let detectableFor: (from: number, to: number) => Uint8Array
if (weights !== null) {
detectableFor = (from, to) => detectableFromWeights(weights, from, to)
} else {
const marginal = await cisMarginal(job)
detectableFor = (from, to) => detectableFromMarginal(marginal.subarray(from, to), 5)
}
const spans = [...bySequence.values()].map((group) => group[0].bounds)
const block = usableBlock(job, spans, balance)
const rows: Array<{ anchor: ChromosightAnchor; score: number; pvalue: number; qvalue: number }> = []
const windows: Float64Array[] = []
let done = 0
for (const [, group] of bySequence) {
await job.breathe()
job.progress(done, used.length, `${group.length.toLocaleString()} anchors`)
const span = group[0].bounds
const n = span.to - span.from
if (n <= kernel.size) {
done += group.length
continue
}
// The band the scan reaches, not the square around it — the same change
// `run/chromosight.ts` makes and for the same reason: `prepareSquare` reads
// `0 … keepDistance` and nothing else, and the square was a ceiling on the
// sequence length rather than a requirement of the arithmetic.
const contacts = await bandFor(job, span, {
hi: keepDistanceFor(n, maxDist, kernel.size),
block,
balance,
value,
onTile: async () => {
await job.breathe()
},
})
const square: Square = {
contacts,
n,
detectable: detectableFor(span.from, span.to),
}
const matrix = prepareSquare(square, { maxDist, largestKernel: kernel.size })
const scored = quantifyOnSquare(
square,
matrix,
kernel,
settings,
group.map((anchor) => ({ bin1: anchor.row - span.from, bin2: anchor.col - span.from })),
)
for (let i = 0; i < scored.length; i++) {
rows.push({
anchor: group[i],
score: scored[i].pattern.score,
pvalue: scored[i].pattern.pvalue,
qvalue: Number.NaN,
})
windows.push(scored[i].values)
}
done += group.length
}
job.progress(used.length, used.length, 'done')
job.throwIfCancelled()
if (windows.length === 0) {
job.note('No anchor landed on a sequence long enough to score against this template.')
return []
}
const qvalues = fdrCorrection(rows.map((row) => row.pvalue))
for (let i = 0; i < rows.length; i++) rows[i].qvalue = qvalues[i]
const scored = rows.filter((row) => Number.isFinite(row.score)).length
if (scored < rows.length) {
job.note(
`${(rows.length - scored).toLocaleString()} of ${rows.length.toLocaleString()} anchors could ` +
'not be scored: their window holds too many undetectable or empty pixels. They are in the ' +
'table with no score, and they contribute nothing to the pileup.',
)
}
const size = kernel.size
const raw = pileupOf(windows, size)
const values = preset.onDiagonal ? symmetrise(raw, size) : raw
const outputs: TaskOutput[] = [
pileupOutput(values, size, { preset, n: scored, binsize }),
spectrumOf(rows, preset),
]
// …and the same scores against separation, where the anchors have one. A
// pileup and a histogram both average the distance away, and a set whose near
// anchors score and whose far ones do not is the same picture as one where
// they all score — see `analysis/loopSpectrum.ts`.
if (!preset.onDiagonal) {
const spectrum = loopSpectrum(
rows.map((row) => ({
separation: (row.anchor.col - row.anchor.row) * binsize,
score: row.score,
})),
{
id: 'apa-chromosight-loop-spectrum',
label: `Loop spectrum · ${preset.name.replace(/_/g, ' ')}`,
noun: 'anchors',
yLabel: 'Pearson correlation with the template',
},
)
if (spectrum !== null) outputs.push(spectrum)
}
outputs.push(tableOf(rows, binsize))
return outputs
}
function pileupOutput(
values: Float64Array,
size: number,
about: { preset: ChromosightPreset; n: number; binsize: number },
): MatrixOutput {
return {
kind: 'matrix',
id: 'apa-chromosight',
label: `Pileup of ${about.n.toLocaleString()} ${about.preset.name.replace(/_/g, ' ')} · chromosight`,
nrows: size,
ncols: size,
values: Float32Array.from(values),
// Base pairs, not bins: the question a reader has of a pileup is how far
// from the anchor a thing lies (`run/labels.offsetTicks`).
xLabel: 'Distance from the centre',
yLabel: 'Distance from the centre',
xTicks: offsetTicks(size, about.binsize),
yTicks: offsetTicks(size, about.binsize),
// Drawn as `pileup_plot` draws it — the fixed 0–2 seismic scale — so this
// picture and the one in a chromosight PDF are the same picture.
...CHROMOSIGHT_SPECTRUM,
stats: [
{ label: 'windows', value: about.n, digits: 0, hint: 'Anchors that could be scored.' },
{
label: 'P2LL',
value: p2ll(values, size, size),
digits: 2,
hint: 'Centre against the lower-left corner. Above 1 is enrichment.',
},
],
}
}
/**
* The distribution of per-anchor scores.
*
* The half of the answer a pileup cannot carry. Two sets with the same mean
* window are the same picture, and one of them can be a hundred good loops while
* the other is four hundred coordinates of which fifty are loops — which is
* exactly the distinction anybody piling up a called set is trying to make. The
* preset's own calling threshold is marked, because that is where chromosight
* would have stopped believing the score.
*/
function spectrumOf(
rows: ReadonlyArray<{ score: number }>,
preset: ChromosightPreset,
): CurveOutput {
const scores = rows.map((row) => row.score).filter((score) => Number.isFinite(score))
const lo = scores.length > 0 ? Math.min(...scores) : 0
const hi = scores.length > 0 ? Math.max(...scores) : 1
const bins = Math.max(10, Math.min(50, Math.round(Math.sqrt(scores.length) * 2)))
const width = hi > lo ? (hi - lo) / bins : 1
const x = new Float64Array(bins)
const y = new Float64Array(bins)
for (let i = 0; i < bins; i++) x[i] = lo + width * (i + 0.5)
for (const score of scores) y[Math.min(bins - 1, Math.max(0, Math.floor((score - lo) / width)))]++
const above = scores.filter((score) => score >= preset.pearson).length
return {
kind: 'curve',
id: 'apa-chromosight-spectrum',
label:
`Score spectrum · ${above.toLocaleString()} of ${scores.length.toLocaleString()} at or above ` +
`${preset.pearson}`,
series: [{ name: 'Anchors', x, y, n: scores.length }],
xLabel: 'Pearson correlation with the template',
yLabel: 'Anchors',
xScale: 'linear',
yScale: 'linear',
}
}
function tableOf(
rows: ReadonlyArray<{ anchor: ChromosightAnchor; score: number; pvalue: number; qvalue: number }>,
binsize: number,
): TableOutput {
return {
kind: 'table',
id: 'apa-chromosight-anchors',
label: 'Per anchor',
columns: [
{ key: 'anchor', label: 'Anchor', kind: 'text' },
{ key: 'separation', label: 'Separation (bp)', kind: 'number', digits: 0 },
{ key: 'score', label: 'Score', kind: 'number', digits: 4, hint: 'Pearson r with the template.' },
{ key: 'pvalue', label: 'p', kind: 'number', digits: 6 },
{ key: 'qvalue', label: 'q', kind: 'number', digits: 6, hint: 'Benjamini-Hochberg.' },
],
rows: rows.map((row) => ({
anchor: row.anchor.label,
separation: (row.anchor.col - row.anchor.row) * binsize,
score: row.score,
pvalue: row.pvalue,
qvalue: row.qvalue,
})),
}
}
|