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 | 2x 2x 2x 2x 2x 6x 6x 6x 2x 2x 2x 2x 2x 2x 2x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 14x 24x 24x 24x 24x 24x 7x 7x 7x 7x 21x 17x 17x 17x 24x 24x 24x 24x 24x 24x 24x 24x 24x 24x 24x 24x 24x 24x 24x 14x 14x 14x 14x 3x 3x 3x 3x 14x 14x 14x 14x 53x 53x 53x 53x 53x 53x 53x 53x 53x 14x 14x 14x | /**
* Base composition, run — in the annotations worker, beside the residues.
*
* One pass per sequence counts each base per bin of the map's own grid, and
* every lane is arithmetic over those counts: the window is a sum of whole
* bins, and the four ratios are one division each. The grid is the session's,
* rebuilt from `inputs.contigLengths` exactly as the motif density is
* (`run/motif.ts`), because a lane on a private grid is a lane that lands in
* the wrong places the moment it is registered.
*
* A sequence the reference does not have is a gap on every lane rather than a
* zero, and it is named: a skew of zero is a measurement, and a sequence
* nobody could read is not one.
*/
import {
atSkewOf,
baseCounts,
cumulativeOf,
gcContentOf,
gcSkewOf,
readBaseCounts,
windowCounts,
type BaseCounts,
} from '../../annotations/composition'
import { idParam, numParam, type TrackTaskRunner } from '../job'
import { COMPOSITION_BINSIZE, compositionLanes } from '../tasks/composition'
import { coversNoMatrix } from '../cliCommand'
import type { ProfileOutput, TaskOutput } from '../types'
import { gridContigs } from './motif'
/**
* Whether this run has no map behind it — a 1D browser, or a shell with no map
* named — whose framing's bin size is a browsing grid rather than a resolution
* anybody picked.
*/
function mapless(job: Parameters<TrackTaskRunner>[0]): boolean {
return coversNoMatrix(job.framing)
}
/** What each lane is called and how it draws. */
const LANES = {
gc: { name: 'GC content', unit: '' },
gcSkew: { name: 'GC skew', unit: '' },
atSkew: { name: 'AT skew', unit: '' },
cumulative: { name: 'Cumulative GC skew', unit: '' },
} as const
export const runComposition: TrackTaskRunner = async (job) => {
const entry = idParam(job.params, 'sequence')
if (job.residues === undefined) {
job.note('This session cannot read residues, so there is no reference to read.')
return []
}
if (entry === null) {
job.note('No reference was chosen. Load a FASTA and pick it under Reference.')
return []
}
const reader = job.residues(entry)
if (reader === undefined) {
job.note('The reference chosen is no longer loaded. Drop the FASTA again, or pick another.')
return []
}
// A refusal says the residues are not here *whole*, and counting needs them
// a few megabases at a time, so a reader that can read by range reads past
// it (`SequenceReader.refusal`). It stopped here in the browser, whose
// indexed FASTA carries both: a genome dropped with its `.fai` was offered
// the task, and the run ended in the note and no tracks.
if (reader.read === undefined && reader.refusal !== undefined) {
job.note(reader.refusal)
return []
}
const lanes = compositionLanes(job.params)
if (lanes.length === 0) {
job.note('No track was ticked, so there is nothing to make.')
return []
}
// The bin size asked for; unanswered, the map's own resolution where there is
// a map, and 10 kb where there is not — which is what a shell run gets, the
// form having already answered it for a run from the drawer.
const asked = numParam(job.params, 'resolution', 0)
const binsize =
asked > 0 ? Math.round(asked) : job.framing.binsize > 0 && !mapless(job) ? job.framing.binsize : COMPOSITION_BINSIZE
const grid = gridContigs(job)
if (grid.length === 0) {
job.note('There are no sequences to lay the bins out on. Open a map, or load the FASTA itself.')
return []
}
const windowBp = Math.max(0, Math.round(numParam(job.params, 'window', 0)))
// An odd number of bins, centred: the window asked for, rounded to the grid
// and then up to the next odd count — which is the width `windowCounts`
// reads, so it is the width the note and the lane's name say.
const rounded = Math.max(1, Math.round(windowBp / binsize))
const span = rounded % 2 === 1 ? rounded : rounded + 1
const circular = new Set(job.framing.circular ?? [])
const total = grid.reduce((sum, contig) => sum + Math.max(1, Math.ceil(contig.length / binsize)), 0)
const values: Record<(typeof lanes)[number], Float32Array> = {
gc: new Float32Array(total).fill(Number.NaN),
gcSkew: new Float32Array(total).fill(Number.NaN),
atSkew: new Float32Array(total).fill(Number.NaN),
cumulative: new Float32Array(total).fill(Number.NaN),
}
// A reference read by range is read a few megabases at a time, none of them
// kept, as the motif search reads one: a sequence asked for whole is a
// chromosome held for nothing, and the pages under it cost more to find than
// its bases cost to count. A reference held is counted where it is held.
const { read } = reader
const unread: string[] = []
const resized: string[] = []
let startBin = 0
for (let i = 0; i < grid.length; i++) {
const contig = grid[i]
const bins = Math.max(1, Math.ceil(contig.length / binsize))
job.progress(i, grid.length, contig.name)
let counted: { counts: BaseCounts; length: number } | undefined
if (read !== undefined) {
counted = await readBaseCounts((start, end) => {
// A Stop is seen between two reads, as the motif search sees one: a
// mammal read by range in the browser is a minute, not a moment.
job.throwIfCancelled()
return read(contig.name, start, end)
}, binsize)
} else {
const seq = reader.bases(contig.name)
counted = seq === undefined ? undefined : { counts: baseCounts(seq, binsize), length: seq.length }
}
if (counted === undefined) {
unread.push(contig.name)
startBin += bins
continue
}
// A reference of another assembly version is read for what it has: its bins
// go where the map's are, and a length disagreement is said out loud.
if (counted.length !== contig.length) resized.push(contig.name)
const { counts } = counted
const windowed = span > 1 ? windowCounts(counts, span, circular.has(contig.name)) : counts
const per = {
gc: gcContentOf(windowed),
gcSkew: gcSkewOf(windowed),
atSkew: atSkewOf(windowed),
// Summed from each bin's own skew, not from the windowed one: a running
// sum of overlapping windows counts every base `span` times.
cumulative: cumulativeOf(gcSkewOf(counts)),
}
const kept = Math.min(bins, counts.steps)
for (const lane of lanes) values[lane].set(per[lane].subarray(0, kept), startBin)
startBin += bins
}
job.progress(grid.length, grid.length, 'done')
if (unread.length === grid.length) {
job.note(
'The reference has none of the map’s sequences under any name this viewer knows, so ' +
'nothing could be measured.',
)
return []
}
if (unread.length > 0) {
job.note(
`${unread.length.toLocaleString()} of the map’s sequences are not in this reference and ` +
`are gaps on every track: ${unread.slice(0, 6).join(', ')}${unread.length > 6 ? ', …' : ''}.`,
)
}
if (resized.length > 0) {
job.note(
`${resized.slice(0, 6).join(', ')}${resized.length > 6 ? ', …' : ''} ` +
`${resized.length === 1 ? 'is' : 'are'} not the length the map says — the reference may ` +
'be another assembly version. Bins were read from its start.',
)
}
if (span > 1) {
job.note(
`Each bin was read over ${span} bins (${(span * binsize).toLocaleString()} bp) centred on ` +
'it; the cumulative skew sums each bin’s own.',
)
}
const reference = job.inputs.labels?.[entry] ?? 'reference'
const suffix = span > 1 ? ` · ${(span * binsize).toLocaleString()} bp window` : ''
const outputs: TaskOutput[] = lanes.map(
(lane) =>
({
kind: 'profile',
id: `composition-${lane.toLowerCase()}`,
label: `${LANES[lane].name} · ${reference}${lane === 'cumulative' ? '' : suffix}`,
values: values[lane],
extent: [0, total],
binsize,
unit: LANES[lane].unit,
}) satisfies ProfileOutput,
)
return outputs
}
|