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 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 | 2x 2x 841x 841x 2x 936823x 936823x 1522x 1522x 1522x 1522x 936823x 2x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 936823x 2x 209x 209x 2x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 833x 497x 336x 336x 336x 336x 336x 336x 833x 833x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 936824x 28322x 28322x 28322x 28322x 29768x 29768x 29768x 29768x 28322x 936824x 14227x 85122x 85122x 85122x 85122x 14227x 936824x 145605x 145605x 145605x 145605x 145605x 936824x 936824x 936824x 936824x 936824x 802x 1521x 1521x 1521x 833x 833x 38x 38x 38x 38x 38x 38x 38x 833x 336x 336x 336x 336x 336x 336x 833x 833x 2x 2x 2x 199x 199x 199x 199x 199x 199x 198x 198x 199x 198x 198x 198x 198x 198x 396x 396x 396x 396x 794x 794x 355x 794x 396x 396x 2x 31x 31x 2x 63x 63x 63x 63x 63x 63x 63x 63x 63x 63x 63x 63x 63x 63x 63x 2x 2x 36x 36x 36x 2x 17050x 17050x 103x 103x 103x 103x 17050x 2x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 2x 21x 21x 2x 35x 35x 35x 35x 35x 35x 35x 35x 35x 35x 35x 35x 35x 35x 35x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 8525x 35x 35x 2x | /**
* Building and querying a {@link FeatureSet}.
*
* Two jobs, kept together because they are two halves of one contract: the
* builder establishes the sort and the per-contig index, and the query relies on
* exactly those. Splitting them would let one drift from the other, and the way
* that failure presents — a long gene silently missing from a window it overlaps
* — is invisible rather than loud.
*/
import type { BarChartBars, FeatureSet, PairSet } from './types'
/** One feature as the parsers hand it over, before anything is sorted. */
interface Raw {
contig: number
start: number
end: number
strand: number
score: number
name: string
id: string
detail: string
exons: number[] | null
cds: [number, number] | null
/** A barChart's per-category values, or `null` on every other format (§8.17). */
bars: readonly number[] | null
/** A structure's bead position, or `null` on every other format (§25). */
xyz: readonly [number, number, number] | null
/** Which haplotype the bead belongs to, as an index; `0` everywhere else. */
hap: number
}
/**
* What a structure knows about itself before any bead is placed (§25).
*
* The counterpart of {@link BarChartSchema}: the spacing and the category names
* come from a pass over the whole file and the coordinates come from every row,
* and they meet once, after the sort, where the arrays can be laid out in the
* order the beads will be drawn and searched in.
*/
export interface StructureSchema {
/** The spacing between consecutive beads, in base pairs; see `StructureModel.width`. */
width: number
/** The category names, in the order the `hap` indices count them. */
haplotypes: string[]
}
/**
* What a barChart knows about its categories before any row is read (§8.17).
*
* Everything in {@link BarChartBars} except the numbers themselves, because the
* numbers are what {@link FeatureBuilder.finish} is for: the labels, the
* colours and the units come from one `track` line and the values come from
* every row, and the two meet exactly once — after the sort, where the values
* can be laid out in the order they will be drawn in.
*/
export type BarChartSchema = Omit<BarChartBars, 'values' | 'max'> & {
/**
* `maxLimit` from the file, or `NaN` to fit the bars to the file's own
* largest value. Never fitted to a window; see {@link BarChartBars.max}.
*/
limit: number
}
/**
* Accumulates features in whatever order a file lists them, then sorts once.
*
* A file is not required to be sorted and most are only sorted by accident:
* this Ensembl GTF opens on chromosome IV. Sorting here — once, in the worker,
* at parse time — is what lets every later query be a binary search.
*/
export class FeatureBuilder {
private readonly rows: Raw[] = []
private readonly contigIds = new Map<string, number>()
readonly contigs: string[] = []
/** The id for a contig name, minting one on first sight. */
contigId(name: string): number {
const existing = this.contigIds.get(name)
if (existing !== undefined) return existing
const id = this.contigs.length
this.contigIds.set(name, id)
this.contigs.push(name)
return id
}
add(
contig: string,
start: number,
end: number,
extra: {
strand?: number
score?: number
name?: string
/** The accession behind the name, where the format has one; see `FeatureSet.id`. */
id?: string
detail?: string
/** Flattened `[start, end, start, end, …]`, in file order. */
exons?: number[]
cds?: [number, number]
/**
* One value per category, for a barChart (§8.17).
*
* Carried through the builder rather than collected beside it, because
* {@link finish} **sorts**: a parallel array written in file order and a
* set read in genomic order are two orders, and the bug that makes is a
* chart drawn on the wrong gene — which looks like data.
*/
bars?: readonly number[]
/**
* Where this bead is, for a structure (§25). Carried through the builder
* for the reason `bars` is: {@link finish} sorts, and a coordinate kept
* in file order beside a set in genomic order is a bead drawn at the
* wrong locus — which looks exactly like a model.
*/
xyz?: readonly [number, number, number]
/** The bead's haplotype, as an index into the schema's list. */
hap?: number
} = {},
): void {
// A zero-length feature is legal in BED and unpaintable: give it one base so
// it has somewhere to be drawn rather than dropping a row the file asserts.
const lo = Math.min(start, end)
const hi = Math.max(lo + 1, Math.max(start, end))
this.rows.push({
contig: this.contigId(contig),
start: lo,
end: hi,
strand: extra.strand ?? 0,
score: extra.score ?? Number.NaN,
name: extra.name ?? '',
id: extra.id ?? '',
detail: extra.detail ?? '',
exons: extra.exons ?? null,
cds: extra.cds ?? null,
bars: extra.bars ?? null,
xyz: extra.xyz ?? null,
hap: extra.hap ?? 0,
})
}
get length(): number {
return this.rows.length
}
/**
* Sort, index, and flatten into typed arrays.
*
* The exons of one feature are sorted too. A GTF lists a minus-strand
* transcript's exons in transcript order, which is descending in genomic
* coordinates, and a renderer that joined consecutive pairs would draw the
* introns as overlapping backwards spans.
*
* `bars` is a barChart's categories (§8.17). Passed to `finish` rather than
* set on the result afterwards so that there is one owner of the whole
* record: the stride is `labels.length` and the values are laid out against
* the sort this method performs, and a caller that attached them later would
* have to reproduce that sort to do it.
*/
finish(bars: BarChartSchema | null = null, model: StructureSchema | null = null): FeatureSet {
const rows = this.rows
// The haplotype last, as a tiebreak and nothing more: a phased structure
// has two rows at every `(contig, start, end)`, and without it their order
// would be whatever the sort happened to leave, which is a model whose
// homologues swap places between two loads of the same file.
rows.sort(
(a, b) => (a.contig - b.contig) || (a.start - b.start) || (a.end - b.end) || (a.hap - b.hap),
)
const n = rows.length
const set: FeatureSet = {
contigs: this.contigs,
length: n,
contig: new Int32Array(n),
start: new Int32Array(n),
end: new Int32Array(n),
strand: new Int8Array(n),
score: new Float32Array(n),
name: new Array<string>(n),
id: new Array<string>(n),
detail: new Array<string>(n),
exonAt: new Int32Array(n + 1),
exonStart: new Int32Array(0),
exonEnd: new Int32Array(0),
cdsStart: new Int32Array(n),
cdsEnd: new Int32Array(n),
contigAt: new Int32Array(this.contigs.length + 1),
maxSpan: new Int32Array(this.contigs.length),
}
let exonTotal = 0
for (const row of rows) exonTotal += row.exons === null ? 0 : row.exons.length >> 1
set.exonStart = new Int32Array(exonTotal)
set.exonEnd = new Int32Array(exonTotal)
// One flat `Float32Array` with a stride, filled in the sorted order below.
// `NaN` throughout to start with, so a row that carried no vector — or one
// whose `expCount` disagreed with the file's — is *no value* rather than a
// run of zeroes, which would draw as a real measurement of nothing.
const stride = bars === null ? 0 : bars.labels.length
const values = new Float32Array(stride * n).fill(Number.NaN)
let widest = 0
// Three columns for a structure, filled in the sorted order below — the
// whole reason the coordinates ride through the builder (§25).
const placed =
model === null
? null
: {
x: new Float32Array(n),
y: new Float32Array(n),
z: new Float32Array(n),
hap: new Int8Array(n),
}
let exonCursor = 0
for (let i = 0; i < n; i++) {
const row = rows[i]
set.contig[i] = row.contig
set.start[i] = row.start
set.end[i] = row.end
set.strand[i] = row.strand
set.score[i] = row.score
set.name[i] = row.name
set.id[i] = row.id
set.detail[i] = row.detail
set.cdsStart[i] = row.cds === null ? -1 : row.cds[0]
set.cdsEnd[i] = row.cds === null ? -1 : row.cds[1]
set.exonAt[i] = exonCursor
if (row.exons !== null) {
const pairs: Array<[number, number]> = []
for (let k = 0; k + 1 < row.exons.length; k += 2) pairs.push([row.exons[k], row.exons[k + 1]])
pairs.sort((a, b) => a[0] - b[0])
for (const [lo, hi] of pairs) {
set.exonStart[exonCursor] = lo
set.exonEnd[exonCursor] = hi
exonCursor++
}
}
if (stride > 0 && row.bars !== null && row.bars.length === stride) {
for (let c = 0; c < stride; c++) {
const value = row.bars[c]
values[i * stride + c] = value
if (value > widest) widest = value
}
}
if (placed !== null) {
// A row with no coordinate in a structure is `NaN`, which every reader of
// the model treats as *not placed* — the same answer a bar chart gives a
// row whose vector was missing, for the same reason.
placed.x[i] = row.xyz === null ? Number.NaN : row.xyz[0]
placed.y[i] = row.xyz === null ? Number.NaN : row.xyz[1]
placed.z[i] = row.xyz === null ? Number.NaN : row.xyz[2]
placed.hap[i] = row.hap
}
const span = row.end - row.start
if (span > set.maxSpan[row.contig]) set.maxSpan[row.contig] = span
}
set.exonAt[n] = exonCursor
// Half-open run per contig. Written by counting forward rather than by
// searching, because the rows are already grouped by the sort above.
let cursor = 0
for (let c = 0; c < this.contigs.length; c++) {
set.contigAt[c] = cursor
while (cursor < n && set.contig[cursor] === c) cursor++
}
set.contigAt[this.contigs.length] = cursor
if (bars !== null) {
const { limit, ...rest } = bars
set.bars = {
...rest,
values,
// The file's own ceiling where it declares one, and the tallest bar in
// the file where it does not. Zero for a file of zeroes and empties,
// which would otherwise divide every bar by nothing.
max: Number.isFinite(limit) && limit > 0 ? limit : widest,
}
}
if (placed !== null && model !== null) {
set.model = {
...placed,
haplotypes: [...model.haplotypes],
width: model.width,
}
}
return set
}
}
/**
* How one feature names itself in prose: `LipL1 (PF3D7_1314600)`.
*
* The symbol *and* the accession, wherever the file gave two different strings,
* because a tooltip that said only `VAR` over a *Plasmodium* map would be
* answering "which gene is this" with the name of sixty of them. Just the one
* string where there is only one to have: a BED row has no accession, and a gene
* with no symbol is drawn as its accession already, so `name === id` there
* (`annotations/types.ts`).
*
* The glyph on the map is *not* this — it is the name alone. A label drawn
* inside a box that is often eight pixels wide has room for one thing, and the
* one thing a reader scanning a track wants is the name.
*/
export function featureLabel(set: FeatureSet, at: number): string {
const name = set.name[at]
const id = set.id[at]
if (id === '' || id === name) return name
return name === '' ? id : `${name} (${id})`
}
/**
* The half-open index range worth *examining* for features overlapping
* `[lo, hi)` on `contig`.
*
* Not the answer — the answer is whichever of those rows actually overlaps, and
* the caller tests that as it draws. The rows are sorted by `start`, so a
* feature can begin before `lo` and still end inside the window; backing the
* lower bound off by the contig's longest feature is what makes the range a
* superset rather than a plausible-looking undercount.
*/
export function featureRange(
set: FeatureSet,
contig: number,
lo: number,
hi: number,
): [number, number] {
if (contig < 0 || contig >= set.contigs.length) return [0, 0]
const from = set.contigAt[contig]
const to = set.contigAt[contig + 1]
if (to <= from) return [0, 0]
const back = set.maxSpan[contig]
const first = lowerBound(set.start, from, to, lo - back)
const last = lowerBound(set.start, from, to, hi)
return [first, last]
}
/** First index in `[from, to)` whose value is `>= target`, else `to`. */
function lowerBound(values: Int32Array, from: number, to: number, target: number): number {
let low = from
let high = to
while (low < high) {
const mid = (low + high) >> 1
if (values[mid] >= target) high = mid
else low = mid + 1
}
return low
}
/** An empty set, for a file that parsed cleanly and held nothing. */
export function emptyFeatureSet(): FeatureSet {
return new FeatureBuilder().finish()
}
/**
* The names in a set that a map has never heard of, and the other way round.
*
* Spec §8.2: an overlap of zero has to be *said*, because it looks exactly like
* a broken feature. Both directions are reported — a GTF using `chrI` against a
* map using `I` is a naming mismatch, a GTF covering only chromosome I against
* a whole-genome map is a partial annotation, and the two want different fixes.
*/
export function overlapWith(
fileContigs: readonly string[],
mapContigs: readonly string[],
sample = 4,
): { matched: number; total: number; missing: string[]; extra: string[] } {
const inFile = new Set(fileContigs)
const inMap = new Set(mapContigs)
const missing = mapContigs.filter((name) => !inFile.has(name))
const extra = fileContigs.filter((name) => !inMap.has(name))
return {
matched: mapContigs.length - missing.length,
total: mapContigs.length,
missing: missing.slice(0, sample),
extra: extra.slice(0, sample),
}
}
/** Accumulates pairs; the bedpe counterpart of {@link FeatureBuilder}. */
export class PairBuilder {
private readonly rows: Array<{
c1: number
s1: number
e1: number
c2: number
s2: number
e2: number
score: number
name: string
detail: string
}> = []
private readonly contigIds = new Map<string, number>()
readonly contigs: string[] = []
contigId(name: string): number {
const existing = this.contigIds.get(name)
if (existing !== undefined) return existing
const id = this.contigs.length
this.contigIds.set(name, id)
this.contigs.push(name)
return id
}
add(
contig1: string,
start1: number,
end1: number,
contig2: string,
start2: number,
end2: number,
extra: { score?: number; name?: string; detail?: string } = {},
): void {
this.rows.push({
c1: this.contigId(contig1),
s1: start1,
e1: Math.max(start1 + 1, end1),
c2: this.contigId(contig2),
s2: start2,
e2: Math.max(start2 + 1, end2),
score: extra.score ?? Number.NaN,
name: extra.name ?? '',
detail: extra.detail ?? '',
})
}
get length(): number {
return this.rows.length
}
finish(): PairSet {
const n = this.rows.length
const set: PairSet = {
contigs: this.contigs,
length: n,
contig1: new Int32Array(n),
start1: new Int32Array(n),
end1: new Int32Array(n),
contig2: new Int32Array(n),
start2: new Int32Array(n),
end2: new Int32Array(n),
score: new Float32Array(n),
name: new Array<string>(n),
detail: new Array<string>(n),
}
for (let i = 0; i < n; i++) {
const row = this.rows[i]
set.contig1[i] = row.c1
set.start1[i] = row.s1
set.end1[i] = row.e1
set.contig2[i] = row.c2
set.start2[i] = row.s2
set.end2[i] = row.e2
set.score[i] = row.score
set.name[i] = row.name
set.detail[i] = row.detail
}
return set
}
}
|