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 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 | 2x 2x 2x 2x 4x 4x 2x 2x 10x 10x 10x 10x 2x 2x 2x 2x 27x 27x 27x 27x 10995x 10995x 10995x 12539017x 12539017x 12539016x 12539016x 12539016x 10995x 10995x 27x 27x 2x 15x 15x 15x 15x 15x 15x 7257x 7257x 7257x 8354887x 8354887x 7257x 15x 15x 30x 30x 30x 30x 14514x 14212x 14212x 14212x 30x 30x 30x 30x 14514x 14514x 14514x 30x 30x 2x 10x 10x 10x 8x 8x 10x 10x 10x 10x 10x 10x 10x 10x 10x 10x 300x 300x 300x 187753x 187753x 183639x 183639x 187753x 218614308x 218614308x 213510529x 213510529x 213510529x 213510529x 213510529x 213510529x 183639x 183639x 300x 267x 267x 267x 267x 184773x 180659x 184773x 180659x 180659x 180659x 180659x 267x 267x 259x 259x 33x 33x 300x 33x 300x 2980x 2980x 2980x 2980x 2980x 2980x 300x 2147x 2116x 2147x 2116x 2116x 2116x 300x 33x 300x 10x 10x 3734x 3734x 3734x 4184120x 4184120x 4184120x 3734x 10x 10x 333x 333x 333x 333x 189900x 185755x 185755x 185755x 333x 333x 2x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 5466x 5466x 5466x 5466x 5639x 5639x 5639x 173x 974x 974x 952x 952x 952x 173x 173x 5466x 5639x 1863282x 1863282x 1863248x 1863248x 1863248x 1863248x 5639x 1869014x 1869014x 1868962x 1868962x 1868962x 1868962x 5466x 5466x 66x 66x 66x 66x 2x 22x 22x 22x 22x 22x 22x 22x 22x 22x 22x 22x 4730x 4730x 4730x 4730x 3683118x 5x 5x 5x 3683113x 3683118x 3682985x 3682985x 3683118x 128x 128x 3683118x 3683118x 4730x 22x 22x 2x 36x 36x 36x 36x 36x 36x 36x 36x 36x 36x 36x 36x 36x 36x 36x 665x 665x 665x 665x 754x 754x 754x 89x 89x 89x 665x 754x 754x 665x 665x 36x 36x 49175x 49175x 49175x 102x 102x 102x 102x 102x 102x 102x 204x 204x 187x 187x 95x 95x 95x 187x 187x 187x 204x 204x 102x 102x 102x 102x 102x 102x 102x 11x 11x 11x 11x 39888x 39888x 39888x 11x 11x 18x 18x 18x 18x 18x 18x 18x 18x 3703251x 3703251x 3703251x 3703251x 3703251x 3703251x 18x 18x 2x 146x 146x 146x 146x 146x 146x 146x 146x 146x 5x 5x 5x 146x 4x 4x 4x 146x 5x 5x 5x 5x 5x 5x 146x 132x 146x 146x 146x 146x 146x 146x 146x 122x 122x 119x 14x 14x 14x 14x 6x 8x 14x 14x 14x 14x 10x 27x 2x 2x 2x 2x 2x 2x 2x 8x 8x 4x 146x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 146x 146x 146x 146x 11x 11x 11x 11x 146x 37x 37x 5x 5x 5x 5x 5x 5x 5x 5x 37x 19x 19x 30x 13x 13x 13x 13x 37x 135x 98x 98x 98x 146x 146x 1x 1x 146x 146x 146x 146x 146x 146x 146x 146x 2x 9x 7x 9x 2421458x 2421458x 2421458x 7x 7x 2x 93x 93x 93x 66x 66x 66x 66x 66x 66x 66x 66x 66x 66x 62x 62x 44x 44x 44x 30x 30x 30x 44x 44x | /**
* Values a file does not contain: coverage normalisation, iterative balancing,
* and the distance-decay baseline every Hi-C map is read against.
*
* A cooler stores one number per pixel and, if somebody ran `cooler balance`,
* one weight per bin. Everything else a reader wants — normalised counts on a
* file nobody balanced, the ratio of a cell to what its genomic distance would
* predict — is arithmetic over numbers that are already on screen. `ui/scores.ts`
* offers the file's own columns because a file that went to the trouble of
* storing an expected value had no way of saying so; this is the other half of
* that argument, for the far more common file that stored nothing but counts.
*
* **Why this is honest here and would not be in HiGlass.** Every function below
* measures a normalisation from the block it is handed. That is a genome-wide
* statement only if the block is the genome — and spec §2 is precisely the claim
* that it is: a MAG is 2–6 Mb, its matrix is loaded in one bounded read, and pan
* and zoom are canvas operations over a block that never changes. A tile server
* computing coverage from the 256×256 tile under the pointer would produce a
* different normalisation at every zoom level and every pan; here there is one
* block, so there is one answer. Where that stops being true — an `__unbinned__`
* genome, a windowed region of a whole cooler chromosome — the number is a
* statement about the window, and {@link DerivedValues.scope} is what says so
* rather than leaving the caption to imply otherwise.
*
* Nothing here draws, and nothing here reads. Each function takes a dense block
* and returns a new array of the same shape, which is what keeps the renderer,
* the pooling layer and the reader unaware that derived values exist.
*/
// The one import from `analysis/`, and a pure one — arithmetic over arrays, no
// worker and no wire format. It is shared rather than copied for the reason
// `correlationField` is: the O/E a reader looks at and the one the compartment
// chain correlates must be one implementation, or they agree until they do not.
import { armExpected, type ArmFrame } from '../analysis/arms'
import {
adaptiveSmoothValues,
parseSmooth,
smoothingCaption,
smoothValues,
SMOOTH_INLINE_CELLS,
type Smoothing,
} from './smooth'
/** The part of a read this module needs: shape, values, and where they sit. */
export interface BlockValues {
nrows: number
ncols: number
/** Row-major, at least `nrows * ncols` long. */
values: Float32Array
rowExtent: readonly [number, number]
colExtent: readonly [number, number]
}
/**
* A stretch of one axis belonging to one sequence.
*
* Structurally a `worker/protocol.ts` `AxisGroup` minus everything this module
* has no use for, and declared here rather than imported so that `render/` keeps
* knowing nothing about the wire format — the same separation `render/heatmap.ts`
* observes by taking a `DenseMatrix` rather than a `MatrixResult`.
*/
export interface AxisRun {
label: string
/** First row/column of the run within the block. */
start: number
/** One past the last row/column of the run within the block. */
end: number
}
/**
* How the cells were put on a comparable footing, before anything is read off
* them.
*
* `none` is the file as written and is the default everywhere: a viewer that
* silently normalised would be answering a question nobody asked, and the count
* in a cell is the one number the file actually contains.
*
* The other three are the standard ladder, cheapest first, and they differ in
* what they assume rather than in quality:
*
* - **`vc`** — vanilla coverage. Divide each cell by its two bins' coverage.
* One pass, no iteration, and it assumes the only bias is how often each bin
* was seen at all.
* - **`vcsqrt`** — the same correction at half strength, `√` of each marginal.
* The convention exists because full VC over-corrects on sparse data, where a
* bin's marginal is as much sampling noise as it is bias.
* - **`ice`** — the fixed point rather than the single step: rescale, remeasure,
* repeat until every bin contributes equally. This is what `cooler balance`
* writes into `/bins/weight`, so on a balanced file `ice` and `balanced`
* should agree to the tolerance below — and where they do not, the difference
* is the filtering a real balancing run does and this one deliberately does
* not (see {@link balanceBlock}).
*
* Juicer's `KR` is not a fourth entry because it is not a fourth *answer*.
* Knight–Ruiz and iterative correction solve the same problem — the unique
* diagonal scaling that makes a symmetric non-negative matrix doubly stochastic
* — and converge to the same weights; KR is a Newton solver that gets there in
* fewer passes. At §2's scale the whole matrix is at most low millions of cells
* and the fixed point converges in milliseconds, so the faster solver would buy
* a name in a menu and nothing a reader could see. A file balanced with `KR`
* elsewhere and read here under `ice` shows the same picture.
*/
export type Normalisation = 'none' | 'vc' | 'vcsqrt' | 'ice'
/** Every normalisation, in the order they should be offered. */
export const NORMALISATIONS: readonly Normalisation[] = ['none', 'vc', 'vcsqrt', 'ice']
export function isNormalisation(value: string): value is Normalisation {
return (NORMALISATIONS as readonly string[]).includes(value)
}
/**
* Rounds of {@link balanceBlock} before it gives up and returns what it has.
*
* Iterative correction on a well-connected matrix converges geometrically and is
* done in well under twenty rounds; the cap is for the matrices that do not
* converge at all, which are the ones made almost entirely of zeros. Returning a
* partly-converged scaling is right for those: it is still a better picture than
* raw counts, and refusing would mean a reader with a sparse map gets an error
* where every other normalisation gives them something.
*/
export const BALANCE_MAX_ROUNDS = 200
/**
* Cell visits one balancing pass is allowed in total, across all its rounds.
*
* The round cap alone is a cap on *iterations*, and iterations are not what
* costs: one round is one pass over the block, so 200 rounds over a 4,096²
* whole-genome overview is 3.3 billion cell visits and a freeze measured in
* seconds — on the main thread, at the moment a reader picks an entry from a
* dropdown. The blocks this is really for are far smaller (§2: a MAG at 5 kb is
* ~1,000 bins, so ~10⁶ cells, and the full 200 rounds cost less than the frame
* they land in), and the budget is what keeps the pathological case from being
* paid for by the ordinary one.
*
* Spent rather than refused, and a floor of a few rounds regardless: iterative
* correction converges geometrically, so a dozen rounds over a large block
* removes most of the bias the fixed point would, and stopping early is much
* closer to right than falling back to raw counts would be. What it is not is
* *converged*, which is why {@link balanceBlock} reports it.
*/
export const BALANCE_CELL_BUDGET = 2.5e8
/** Rounds this block can afford, between a usable floor and the cap. */
function roundsFor(cells: number): number {
if (cells <= 0) return BALANCE_MAX_ROUNDS
return Math.max(4, Math.min(BALANCE_MAX_ROUNDS, Math.floor(BALANCE_CELL_BUDGET / cells)))
}
/**
* When the marginals are flat enough to call the block balanced.
*
* The largest relative departure of any bin's marginal from their mean. 1e-5 is
* far below what a 256-entry colormap can show — the picture stops changing
* several rounds earlier — so a block that gets this far is drawn and captioned
* as balanced, without the qualification {@link BalanceResult.converged} puts on
* one the rounds ran out on. It is not where the loop stops while it can afford
* to go on: that is {@link BALANCE_PRECISION}.
*/
export const BALANCE_TOLERANCE = 1e-5
/**
* When the marginals are flat enough to stop, rounds permitting.
*
* The cells are handed on as Float32, and someone comparing them with a
* reference — cooltools' iterative correction, `cooler balance` — is comparing
* to that precision: the weights left 1e-9 from the fixed point put each cell
* within the half unit in the last place its own rounding costs
* (`test/unit/normalise-cooler.test.ts`). Stopping at {@link BALANCE_TOLERANCE}
* left them 1e-5 apart. Iterative correction converges geometrically, so the
* four decades between the two are some thirty rounds, and the cell budget still
* caps them: a block that runs out between the two is balanced to what it
* reached and still called converged.
*/
export const BALANCE_PRECISION = 1e-9
/**
* What a derived array is, beside the numbers.
*
* The values alone cannot say what colour scale they want or what the colourbar
* should be headed, and every caller needs both. `neutral` is the load-bearing
* one: a ratio is not a quantity with a bright end and a dark end but one with a
* *middle*, and a picture of one drawn on a sequential ramp anchored at zero
* hides the only thing it is for. See {@link isDiverging}.
*/
export interface DerivedValues {
values: Float32Array
/**
* The value that means "nothing to report", or `null` where the result is an
* ordinary one-ended quantity.
*
* `1` for a ratio, `0` for a difference. The colour scale is made symmetric
* about it and the default colormap becomes a diverging one; both follow from
* this single number rather than from a flag per derivation.
*/
neutral: number | null
/** What the colourbar is headed and what the readout calls the number. */
label: string
/**
* What the number was measured over, where that is narrower than the label
* implies. `null` when the block is the whole of what the statement is about.
*/
scope: string | null
/**
* A smoothing this derivation wanted and did not have.
*
* Present only when an `adaptive` setting is chosen on a block big enough to be
* worth moving and the caller said it could move it. The field above is the
* honest unsmoothed picture; a caller that can run the filter starts one and
* derives again. See {@link SmoothingRequest}.
*/
pending?: SmoothingRequest
}
/** Whether a derived quantity is read against a middle rather than against zero. */
export function isDiverging(derived: Pick<DerivedValues, 'neutral'>): boolean {
return derived.neutral !== null
}
// ---------------------------------------------------------------------------
// Coverage and balancing
// ---------------------------------------------------------------------------
/**
* The two marginals of a block, with non-finite cells contributing nothing.
*
* A `NaN` is a bin the balancer threw out (§6.4) and it has to be skipped rather
* than propagated: one such bin would otherwise turn the marginal of every bin
* it touches into `NaN`, and a single discarded row would blank the whole map.
*/
export function marginalsOf(block: BlockValues): { rows: Float64Array; cols: Float64Array } {
const { nrows, ncols, values } = block
const rows = new Float64Array(nrows)
const cols = new Float64Array(ncols)
for (let r = 0; r < nrows; r++) {
const base = r * ncols
let sum = 0
for (let c = 0; c < ncols; c++) {
const value = values[base + c]
if (!Number.isFinite(value)) continue
sum += value
cols[c] += value
}
rows[r] = sum
}
return { rows, cols }
}
/**
* Coverage normalisation: each cell divided by the coverage of its two bins.
*
* The marginals are put on a mean of 1 before they divide anything, which is
* what keeps the result in the units it started in. Dividing by raw marginals
* would leave a matrix whose values are around `1 / total`, so every colour
* scale, every stored clip and every number in the readout would change meaning
* with the depth of the library — and a reader flipping between `counts` and
* `VC` would see the picture blank rather than flatten.
*
* `power` is `1` for VC and `0.5` for VC √; nothing else is offered, because the
* exponent is a convention rather than a parameter and two conventions exist.
*
* A bin with no coverage at all gets `NaN` rather than a division by zero. That
* is the same value a bin the balancer discarded already carries, so the raster
* leaves both transparent and the readout says `NaN (unbalanced)` for both —
* which is the truth in both cases: there is no normalised value for a bin
* nothing was observed in.
*/
export function coverageNormalise(block: BlockValues, power: 0.5 | 1): Float32Array {
const { nrows, ncols, values } = block
const { rows, cols } = marginalsOf(block)
const rowScale = scaleToUnitMean(rows, power)
const colScale = scaleToUnitMean(cols, power)
const out = new Float32Array(nrows * ncols)
for (let r = 0; r < nrows; r++) {
const base = r * ncols
const rw = rowScale[r]
for (let c = 0; c < ncols; c++) {
out[base + c] = rw > 0 && colScale[c] > 0 ? values[base + c] / (rw * colScale[c]) : Number.NaN
}
}
return out
}
/**
* A marginal rescaled to a mean of one over the bins that have any, raised to
* `power`.
*
* The mean is taken over the *non-empty* bins only. A block half of which is a
* contig nothing mapped to would otherwise have its mean halved by the empty
* half, and every value in the half that does carry data would be scaled by two
* for a reason that is a fact about the other contig.
*/
function scaleToUnitMean(marginal: Float64Array, power: 0.5 | 1): Float64Array {
let total = 0
let seen = 0
for (const value of marginal) {
if (!(value > 0)) continue
total += value
seen++
}
const out = new Float64Array(marginal.length)
if (seen === 0) return out
const mean = total / seen
for (let i = 0; i < marginal.length; i++) {
const value = marginal[i]
out[i] = value > 0 ? (power === 1 ? value / mean : Math.sqrt(value / mean)) : 0
}
return out
}
/**
* Iterative correction: the diagonal scaling under which every bin contributes
* equally.
*
* The assumption is the one every balancing method makes and it is worth stating
* because it is the reason a balanced map is read differently from a raw one:
* *every bin should see the same total number of contacts*. Whatever it does not
* see is bias — mappability, restriction-site density, GC, copy number — and the
* weights are whatever makes the marginals flat. Where that assumption is wrong,
* balancing is wrong in exactly the same way, and the classic case is real: a
* region genuinely present at twice the copy number is corrected until it is not.
*
* Two updates, chosen by the block rather than offered:
*
* - **Symmetric**, for a block whose two axes cover the same bins. One weight
* vector applied to both sides, stepped by `√(target / marginal)`, so the
* result stays symmetric to the last bit. Alternating row and column passes
* over a symmetric matrix does *not* — each half-step breaks the symmetry the
* other restores — and a contact map that is not symmetric draws a visible
* seam along its own diagonal, which is the one artefact a reader is
* guaranteed to notice and to misread as structure.
* - **Alternating** (Sinkhorn–Knopp), for a rectangular trans block, where there
* is no shared bin to give one weight to and the two axes are scaled
* independently.
*
* Rows that are entirely empty are held at zero weight throughout and emitted as
* `NaN`; they are excluded from the target, or a matrix that is mostly empty
* would drive the target towards zero and take every real row with it.
*
* `converged` is `false` where the rounds ran out — see
* {@link BALANCE_CELL_BUDGET}. The values are still the best scaling found and
* are still worth drawing; what a caller must not do is describe them as
* balanced without qualification, which is what {@link DerivedValues.scope}
* carries to the caption.
*/
export interface BalanceResult {
values: Float32Array
/**
* Whether the marginals reached {@link BALANCE_TOLERANCE} before the rounds
* ran out — whether or not they went on to {@link BALANCE_PRECISION}.
*/
converged: boolean
}
export function balanceBlock(block: BlockValues): BalanceResult {
const { nrows, ncols, values } = block
const symmetric =
nrows === ncols &&
block.rowExtent[0] === block.colExtent[0] &&
block.rowExtent[1] === block.colExtent[1]
const rowWeight = new Float64Array(nrows)
const colWeight = new Float64Array(ncols)
const { rows: rowMass, cols: colMass } = marginalsOf(block)
for (let r = 0; r < nrows; r++) rowWeight[r] = rowMass[r] > 0 ? 1 : 0
for (let c = 0; c < ncols; c++) colWeight[c] = colMass[c] > 0 ? 1 : 0
const rowSum = new Float64Array(nrows)
const colSum = new Float64Array(ncols)
const rounds = roundsFor(nrows * ncols)
let converged = false
for (let round = 0; round < rounds; round++) {
rowSum.fill(0)
colSum.fill(0)
for (let r = 0; r < nrows; r++) {
const rw = rowWeight[r]
if (rw === 0) continue
const base = r * ncols
let sum = 0
for (let c = 0; c < ncols; c++) {
const cw = colWeight[c]
if (cw === 0) continue
const value = values[base + c]
if (!Number.isFinite(value)) continue
const scaled = rw * value * cw
sum += scaled
colSum[c] += scaled
}
rowSum[r] = sum
}
if (symmetric) {
const target = meanOfActive(rowSum, rowWeight)
// Nothing left to balance — an empty block, or one every bin of which was
// dropped. Converged in the only sense available: another round would
// change nothing.
if (!(target > 0)) {
converged = true
break
}
let drift = 0
for (let r = 0; r < nrows; r++) {
if (rowWeight[r] === 0) continue
const mass = rowSum[r]
if (!(mass > 0)) {
rowWeight[r] = 0
colWeight[r] = 0
continue
}
drift = Math.max(drift, Math.abs(mass / target - 1))
// √ of the ratio, because the weight is applied to *both* ends of every
// cell: a full step would overshoot by exactly the square.
rowWeight[r] *= Math.sqrt(target / mass)
colWeight[r] = rowWeight[r]
}
if (drift < BALANCE_TOLERANCE) converged = true
if (drift < BALANCE_PRECISION) break
continue
}
const rowTarget = meanOfActive(rowSum, rowWeight)
const colTarget = meanOfActive(colSum, colWeight)
if (!(rowTarget > 0) || !(colTarget > 0)) {
converged = true
break
}
let drift = 0
for (let r = 0; r < nrows; r++) {
if (rowWeight[r] === 0) continue
const mass = rowSum[r]
if (!(mass > 0)) {
rowWeight[r] = 0
continue
}
drift = Math.max(drift, Math.abs(mass / rowTarget - 1))
rowWeight[r] *= rowTarget / mass
}
for (let c = 0; c < ncols; c++) {
if (colWeight[c] === 0) continue
const mass = colSum[c]
if (!(mass > 0)) {
colWeight[c] = 0
continue
}
drift = Math.max(drift, Math.abs(mass / colTarget - 1))
colWeight[c] *= colTarget / mass
}
if (drift < BALANCE_TOLERANCE) converged = true
if (drift < BALANCE_PRECISION) break
}
const out = new Float32Array(nrows * ncols)
for (let r = 0; r < nrows; r++) {
const base = r * ncols
const rw = rowWeight[r]
for (let c = 0; c < ncols; c++) {
const cw = colWeight[c]
out[base + c] = rw > 0 && cw > 0 ? rw * values[base + c] * cw : Number.NaN
}
}
return { values: out, converged }
}
/** The mean of the entries whose weight is still live. */
function meanOfActive(sums: Float64Array, weights: Float64Array): number {
let total = 0
let seen = 0
for (let i = 0; i < sums.length; i++) {
if (weights[i] === 0) continue
total += sums[i]
seen++
}
return seen === 0 ? 0 : total / seen
}
// ---------------------------------------------------------------------------
// Expected, and observed over it
// ---------------------------------------------------------------------------
/**
* The distance-decay baseline: what a cell would hold if the only thing
* governing it were how far apart its two ends are.
*
* This is the single strongest signal in any Hi-C map and it is not the one
* anybody is looking for. Contact frequency falls by two to three orders of
* magnitude from the diagonal outwards, so a raw map is a bright diagonal
* fading to nothing and every feature — a loop, a domain boundary, a plasmid
* that belongs to a different genome — is a small deviation riding on it.
* Dividing it out is what makes those deviations the picture (§7.1's colour
* scale is doing the same job far more crudely, by clipping).
*
* **One curve, pooled over sequences.** Cells whose two ends are in the same
* sequence give a mean per bin separation; cells whose ends are in different
* sequences give a single scalar between them all. Per-sequence curves would be
* better on an assembly whose contigs are genuinely different things, and worse
* on the case this viewer is for: a MAG's contigs are pieces of *one* genome
* broken by assembly rather than by biology, and a 40 kb contig on its own
* yields eight points to fit a decay curve to. Pooling spends the whole genome's
* cells on one curve, which is the estimate a reader of a MAG wants.
*
* Separation is measured in **absolute bins**, so a rectangular block covering
* two distant windows of one contig is placed on the curve where it belongs
* rather than at the separations its row and column indices suggest — the same
* distinction `render/compare.ts` draws when it decides which side of the
* diagonal a cell is on.
*/
export interface ExpectedModel {
/**
* Mean value per bin separation, for cells within one sequence, indexed from
* {@link base} rather than from zero.
*
* The offset is what bounds the allocation. Indexed by absolute separation, a
* 100 × 100 block covering two windows of one chromosome three million bins
* apart would need a three-million-entry curve to hold ten thousand cells —
* 24 MB per array, for a block that is 40 KB. Separations within a block span
* at most `nrows + ncols`, whatever their absolute size, so storing them
* relative to the smallest one the block can produce makes the curve the size
* of the block instead of the size of the genome.
*/
cis: Float64Array
/** The separation `cis[0]` describes. */
base: number
/** Mean value between any two different sequences, or `NaN` where there are none. */
trans: number
}
export function expectedModel(
block: BlockValues,
runs: { rows: readonly AxisRun[]; cols: readonly AxisRun[] },
): ExpectedModel {
const { nrows, ncols, values } = block
const { rowIds, colIds, colSpans } = sequenceIds(runs, nrows, ncols)
const rowStart = block.rowExtent[0]
const colStart = block.colExtent[0]
const base = minSeparation(rowStart, nrows, colStart, ncols)
const span = maxSeparation(rowStart, nrows, colStart, ncols) - base
const sums = new Float64Array(span + 1)
const counts = new Float64Array(span + 1)
let transSum = 0
let transCount = 0
for (let r = 0; r < nrows; r++) {
const row = r * ncols
// The column this row's own bin is at, where the separation is zero.
const diagonal = rowStart + r - colStart
const mine = rowIds[r]
for (let s = 0; s + 1 < colSpans.length; s++) {
const c0 = colSpans[s]
const c1 = colSpans[s + 1]
if (colIds[c0] !== mine) {
for (let c = c0; c < c1; c++) {
const value = values[row + c]
if (!Number.isFinite(value)) continue
transSum += value
transCount++
}
continue
}
// Either side of the diagonal, so the separation is a subtraction.
const split = Math.min(c1, Math.max(c0, diagonal))
for (let c = c0; c < split; c++) {
const value = values[row + c]
if (!Number.isFinite(value)) continue
const d = diagonal - c - base
sums[d] += value
counts[d]++
}
for (let c = split; c < c1; c++) {
const value = values[row + c]
if (!Number.isFinite(value)) continue
const d = c - diagonal - base
sums[d] += value
counts[d]++
}
}
}
const cis = new Float64Array(span + 1)
for (let d = 0; d <= span; d++) cis[d] = counts[d] > 0 ? sums[d] / counts[d] : Number.NaN
return { cis, base, trans: transCount > 0 ? transSum / transCount : Number.NaN }
}
/**
* The baseline itself, cell by cell: what {@link observedOverExpected} divides by.
*
* Every cell is replaced by its own expected value — the mean of the block at
* that bin separation within a sequence, the mean between the two sequences
* across one — so the picture is the distance decay with nothing else in it: a
* bright diagonal fading outward, banded, with a flat plate over every off-
* diagonal pair of sequences. It is not a reading of the file. It is the model
* the file's own cells were averaged into, and the reason to look at it is that
* an O/E is only as good as this: a block whose baseline is estimated from a
* handful of cells at long separation produces a ratio that swings, and the
* swing is in the denominator rather than in the biology.
*
* **Masked exactly where the ratio is.** A cell the block has no value for stays
* absent — the baseline exists there arithmetically, but painting a model over a
* balancer-dropped row would put colour where every other map this viewer draws
* shows paper. A separation the model has no mean for, and one whose mean is not
* positive, both go the same way, and those are precisely the two cases
* {@link observedOverExpected} refuses to divide by. So every cell drawn here is
* a cell that has an O/E, which is what makes this the denominator rather than a
* picture that resembles it.
*
* The `0` and the infinity that function blanks are *not* mirrored: those are
* properties of the ratio — nothing observed, or a baseline too small to divide
* by — and neither says the baseline is missing.
*/
export function expectedValues(
block: BlockValues,
runs: { rows: readonly AxisRun[]; cols: readonly AxisRun[] },
model: ExpectedModel = expectedModel(block, runs),
): Float32Array {
const { nrows, ncols, values } = block
const rowLabel = labelsByIndex(runs.rows, nrows)
const colLabel = labelsByIndex(runs.cols, ncols)
const rowStart = block.rowExtent[0]
const colStart = block.colExtent[0]
const out = new Float32Array(nrows * ncols)
for (let r = 0; r < nrows; r++) {
const row = r * ncols
const rowBin = rowStart + r
const mine = rowLabel[r]
for (let c = 0; c < ncols; c++) {
if (!Number.isFinite(values[row + c])) {
out[row + c] = Number.NaN
continue
}
let expected: number
if (mine === colLabel[c]) {
const d = Math.abs(rowBin - (colStart + c)) - model.base
expected = d >= 0 && d < model.cis.length ? model.cis[d] : Number.NaN
} else {
expected = model.trans
}
out[row + c] = expected > 0 ? expected : Number.NaN
}
}
return out
}
/**
* Observed over expected, cell by cell.
*
* The result is a **ratio**, neutral at `1`, and deliberately not a logarithm:
* `types.ts`'s `Transform` already owns the question of what scale a picture is
* drawn on, and `log2` is one of its answers. Returning `log2(O/E)` here would
* put the same decision in two places and make the readout report a number no
* reader asked for — "2.3× expected" is the sentence, and the colour scale is
* what turns it into a position on a ramp.
*
* A separation nothing was observed at, and a cell whose expected value is zero,
* both give `NaN`: there is no ratio to a baseline that does not exist, and the
* raster already leaves `NaN` transparent.
*
* **So do the two ends of the ratio itself**, and that is the one place this
* departs from the arithmetic. A cell where nothing was observed is `0`, which
* on a scale centred at `1` is the strongest depletion the ramp can draw — so an
* O/E map of a sparse block came out as a wall of "maximally depleted" over
* every empty pixel, which is the loudest possible statement of the one thing
* the file does not say. Zero is *nothing seen here*, which is what paper means
* on every other map this viewer draws, so it is drawn as paper. An infinite
* ratio is the same claim from the other end — a baseline too small to divide by
* rather than a measurement — and it goes the same way. Neither is a number a
* colour could honestly stand for, and both would drag the auto range with them.
*/
export function observedOverExpected(
block: BlockValues,
runs: { rows: readonly AxisRun[]; cols: readonly AxisRun[] },
model: ExpectedModel = expectedModel(block, runs),
/** Keep a zero ratio as zero — for the correlation only ({@link Derivation.keepZeros}). */
keepZeros = false,
): Float32Array {
const { nrows, ncols, values } = block
const { rowIds, colIds, colSpans } = sequenceIds(runs, nrows, ncols)
const rowStart = block.rowExtent[0]
const colStart = block.colExtent[0]
const out = new Float32Array(nrows * ncols)
// What each separation is divided by, `NaN` where there is nothing to divide
// by — no mean measured there, or one that is not positive. A value divided by
// `NaN` is `NaN`, so the cell needs no test of its own for either, and none
// for its own absence: `NaN` and the infinities stay what they are divided.
const cis = new Float64Array(model.cis.length)
for (let d = 0; d < cis.length; d++) cis[d] = model.cis[d] > 0 ? model.cis[d] : Number.NaN
const trans = model.trans > 0 ? model.trans : Number.NaN
const base = model.base
for (let r = 0; r < nrows; r++) {
const row = r * ncols
const diagonal = rowStart + r - colStart
const mine = rowIds[r]
for (let s = 0; s + 1 < colSpans.length; s++) {
const c0 = colSpans[s]
const c1 = colSpans[s + 1]
if (colIds[c0] !== mine) {
for (let c = c0; c < c1; c++) out[row + c] = paper(values[row + c] / trans, keepZeros)
continue
}
// A separation outside the curve reads `undefined`, which divides to `NaN`
// as a missing mean does.
const split = Math.min(c1, Math.max(c0, diagonal))
for (let c = c0; c < split; c++) out[row + c] = paper(values[row + c] / cis[diagonal - c - base], keepZeros)
for (let c = split; c < c1; c++) out[row + c] = paper(values[row + c] / cis[c - diagonal - base], keepZeros)
}
}
return out
}
/**
* A ratio as the map draws it: itself, or paper at either end.
*
* Not `> 0`: an extra `/pixels` column may legitimately hold negative numbers,
* and blanking those would be this rule deciding what a file is allowed to
* contain. Exactly the two ends go — zero (and `-0`, which `!== 0` is false
* for) and either infinity — with `isFinite` carrying a `NaN` through
* unchanged. The zero stays where the ratio is for the correlation
* ({@link Derivation.keepZeros}).
*/
function paper(ratio: number, keepZeros: boolean): number {
return (keepZeros || ratio !== 0) && Number.isFinite(ratio) ? ratio : Number.NaN
}
/**
* The sequence of every row and column as a number both axes share, and where
* the columns change sequence.
*
* Numbers rather than {@link labelsByIndex}'s names because the O/E asks the
* question of every cell — sixteen million at the cell cap, twice — and the
* spans because it need not: within a span of columns the answer is the same.
* An axis with no runs at all — a plain cooler read as one span — is sequence
* `0` on both sides, so it is one sequence as far as the O/E is concerned.
*/
function sequenceIds(
runs: { rows: readonly AxisRun[]; cols: readonly AxisRun[] },
nrows: number,
ncols: number,
): { rowIds: Int32Array; colIds: Int32Array; colSpans: Int32Array } {
const ids = new Map<string, number>()
const fill = (axis: readonly AxisRun[], length: number): Int32Array => {
const out = new Int32Array(length)
for (const run of axis) {
let id = ids.get(run.label)
if (id === undefined) {
id = ids.size + 1
ids.set(run.label, id)
}
const end = Math.min(run.end, length)
for (let i = Math.max(0, run.start); i < end; i++) out[i] = id
}
return out
}
const rowIds = fill(runs.rows, nrows)
const colIds = fill(runs.cols, ncols)
const spans = [0]
for (let c = 1; c < ncols; c++) if (colIds[c] !== colIds[c - 1]) spans.push(c)
if (ncols > 0) spans.push(ncols)
return { rowIds, colIds, colSpans: Int32Array.from(spans) }
}
// ---------------------------------------------------------------------------
// The two of them together
// ---------------------------------------------------------------------------
/** What a caller asks for: a normalisation, a baseline, either, or neither. */
export interface Derivation {
normalise?: Normalisation
/** Divide by what the genomic distance predicts. */
expected?: boolean
/**
* Keep a cell where nothing was observed as the ratio it is — zero — rather
* than as paper (§7.5).
*
* For the map's correlation, which is the one reader of the ratio that is not
* a colour ramp. Paper is right for the eye and wrong for the statistic: an
* empty pixel is an observation of no contacts, and the compartment call
* correlates it as one (§18.8, as cooltools does). Blanked, every empty pixel
* left the correlation — a different matrix from the one the task solves,
* over the cells each pair happened to share — and on a sparse map made every
* pair's arithmetic an exception to correct for.
*/
keepZeros?: boolean
/**
* Show what that division would be *by*, instead of doing it.
*
* Overrides {@link expected} rather than combining with it, and there is
* nothing to combine: a baseline over its own baseline is 1 in every cell that
* has one. The two are one control in the interface for the same reason — the
* score picker offers this and the O/E toggle is inert while it is chosen
* (`ui/scores.ts`).
*/
baseline?: boolean
/**
* How much despeckling, as the token the picker and the link carry — a fixed
* radius (`1` is 3 × 3) or an adaptive target (`a200`); `0` or absent for none
* (`render/smooth.ts`).
*/
smooth?: string
/**
* A smoothing already computed for this exact block and setting.
*
* Absent on the first derivation of a block, which is what produces
* {@link DerivedValues.pending}; supplied on the second, once the answer is
* back. Nothing here checks that it belongs to this block — the caller holds
* the key that says so, exactly as it does for the correlation.
*/
smoothed?: SmoothedField
/**
* Whether this caller can run the adaptive filter off this thread.
*
* Absent for a caller that cannot — a test, a panel with no analysis thread
* behind it, or **either file of a compared pair**, where two blocks would need
* two jobs and a half-finished pair would put a smoothed map over an unsmoothed
* one and call the result a ratio. Those compute it inline at whatever it
* costs, which is right: the alternative is not doing it at all.
*/
deferSmoothing?: boolean
/**
* Take the baseline against **both** of a bacterial chromosome's diagonals
* rather than the main one alone (`analysis/arms.ts`, §8.15).
*
* Replaces the distance-decay model wherever {@link expected} or
* {@link baseline} would use it, and nowhere else. `rowBins` and `colBins`
* say which bin of the sequence each row and column of the block is — which
* is what lets a block rotated to put the origin in the middle, or a window
* of one, be read in the one frame.
*/
arms?: {
frame: ArmFrame
rowBins: Int32Array
colBins: Int32Array
/** Half-width along the arms, in bins; `0` for the whole sequence. */
window: number
}
}
/**
* The arm-frame expectation of every cell, masked exactly as
* {@link expectedValues} masks the distance decay: a cell with no value stays
* absent, and so does one whose expectation is not positive.
*/
function armValues(block: BlockValues, arms: NonNullable<Derivation['arms']>): Float32Array {
const expected = armExpected(block.values, arms.rowBins, arms.colBins, arms.frame, arms.window)
const out = new Float32Array(block.nrows * block.ncols)
for (let i = 0; i < out.length; i++) {
const e = expected[i]
out[i] = Number.isFinite(block.values[i]) && e > 0 ? e : Number.NaN
}
return out
}
/**
* Apply a normalisation, then a baseline, then the despeckle, in that order,
* and describe the result.
*
* The order is not a preference. The first two are estimates *from the block*,
* and the distance-decay curve is an average over cells — so measuring it before
* the coverage bias has been taken out means measuring it partly on the bias,
* and every cell of a badly-covered bin is then compared against a baseline its
* own under-coverage helped set. Normalise first, and the curve is a statement
* about distance alone. It is also the order cooltools uses, for the same
* reason.
*
* The despeckle goes **last**, and by the same argument turned around: the other
* two are estimators and this is not. A coverage marginal and a distance curve
* are sums over many cells and are already as smooth as their sample size makes
* them, so feeding them pre-averaged cells would not steady them — it would only
* make each cell's own neighbours count twice toward the baseline it is then
* compared against. Smoothing is a statement about the picture, so it is applied
* to the picture.
*
* It is applied here, per file, rather than after two files are combined, which
* is the whole reason it can be applied to a ratio at all — see
* `render/smooth.ts`.
*
* The identity case returns the block's own array rather than a copy. Nothing
* downstream mutates it, and the alternative is a needless pass over up to 16 M
* cells on the single commonest path there is — a reader who asked for none of
* this.
*/
/** What {@link deriveBlock}'s despeckle step produces, with its baseline beside it. */
interface Smoothed {
values: Float32Array
/** The matched-window baseline, when one was asked for. */
over: Float32Array | null
suffix: string
scope: string | null
/** Set when the work was handed back rather than done — see {@link SmoothingRequest}. */
pending?: SmoothingRequest
}
/**
* Everything the adaptive filter needs, for a caller that can run it elsewhere.
*
* It reads the block twice — once to measure every cell against its neighbourhood and
* once to average the capped field — and at the cell cap that is seconds rather
* than the pass or two every other step here costs. So above
* {@link SMOOTH_INLINE_CELLS} it is not done here. This module runs the pipeline
* to the point the smoothing would happen, hands the caller exactly the arrays it
* would have used, and carries on with the **unsmoothed** field so there is a
* real picture on screen meanwhile. When the answer arrives the caller derives
* again with {@link Derivation.smoothed} supplied and the step is a lookup.
*
* That is the correlation's arrangement (`ui/MapDerivation.ts`) and it is here for
* the correlation's reason: every other step in `deriveBlock` is synchronous, and
* a function that returned a promise for the sake of one of them would make every
* caller wait on arithmetic that was already finished.
*/
export interface SmoothingRequest {
/** The field to average, normalised and before any division. */
values: Float32Array
/** The block as read, which sets the threshold and drives the criterion. */
counts: Float32Array
/** The baseline to average over the same windows, when there is one. */
over: Float32Array | null
nrows: number
ncols: number
spec: Extract<Smoothing, { kind: 'adaptive' }>
}
/** A smoothing that has come back, ready to be substituted for the step. */
export interface SmoothedField {
values: Float32Array
over: Float32Array | null
/** What was *done* — the fixed floor where the filter refused the block. */
applied: Smoothing
scope: string | null
}
/**
* A ratio field under the same two rules {@link observedOverExpected} applies to
* its own: a cell where nothing was seen and a baseline too small to divide by
* are both paper rather than the loudest colour on the ramp (§7.5).
*/
function ratioOf(
numerator: Float32Array,
denominator: Float32Array | null,
length: number,
keepZeros: boolean,
): Float32Array {
const out = new Float32Array(length)
for (let i = 0; i < length; i++) {
const n = numerator[i]
const d = denominator === null ? Number.NaN : denominator[i]
if (!Number.isFinite(n)) {
out[i] = Number.NaN
continue
}
const ratio = d > 0 ? n / d : Number.NaN
out[i] = (keepZeros || ratio !== 0) && Number.isFinite(ratio) ? ratio : Number.NaN
}
return out
}
export function deriveBlock(
block: BlockValues,
derivation: Derivation,
runs: { rows: readonly AxisRun[]; cols: readonly AxisRun[] },
): DerivedValues {
const parts: string[] = []
let scope: string | null = null
let values = block.values
switch (derivation.normalise) {
case 'vc':
values = coverageNormalise(block, 1)
parts.push('VC')
break
case 'vcsqrt':
values = coverageNormalise(block, 0.5)
parts.push('VC √')
break
case 'ice': {
const balanced = balanceBlock(block)
values = balanced.values
parts.push('ICE')
// Said rather than hidden: a partly-converged scaling is a better picture
// than raw counts and a worse one than the label promises, and the reader
// is the one who can decide whether that matters.
if (!balanced.converged) scope = 'stopped before convergence — block too large'
break
}
default:
break
}
// How much despeckling, and what it did — resolved before the baseline below
// because the O/E has to be able to run *after* it. Nothing is applied yet.
const smoothing = parseSmooth(derivation.smooth)
const despeckle = (
source: Float32Array,
/** A baseline the caller is about to divide by, averaged over the same windows. */
over?: Float32Array,
): Smoothed => {
if (smoothing === null) {
return { values: source, over: over ?? null, suffix: '', scope: null }
}
if (smoothing.kind === 'fixed') {
return {
values: smoothValues({ ...block, values: source }, smoothing.radius),
over:
over === undefined
? null
: smoothValues({ ...block, values: over }, smoothing.radius),
suffix: ` · smoothed ${smoothingCaption(smoothing)}`,
scope: null,
}
}
// Already computed elsewhere, for this exact block and setting.
const held = derivation.smoothed
if (held !== undefined) {
return {
values: held.values,
over: held.over,
suffix: ` · smoothed ${smoothingCaption(held.applied)}`,
scope: held.scope,
}
}
// Big enough to be worth moving, and a caller that can move it. The
// unsmoothed field goes out with **no suffix**, because the caption has to
// describe the picture actually drawn: a map labelled as smoothed that is not
// smoothed is the one claim an exported figure must never carry.
if (
derivation.deferSmoothing === true &&
block.nrows * block.ncols > SMOOTH_INLINE_CELLS
) {
return {
values: source,
over: over ?? null,
suffix: '',
scope: null,
pending: {
values: source,
counts: block.values,
over: over ?? null,
nrows: block.nrows,
ncols: block.ncols,
spec: smoothing,
},
}
}
// The adaptive width measures on the block's **as-read** values whatever
// stage this is called from, because the thing it compensates for — counts
// thinning out with distance — is a property of the measurement rather than
// of the arithmetic done to it since.
const adapted = adaptiveSmoothValues(
{ ...block, values: source },
block.values,
smoothing,
over,
)
return {
values: adapted.values,
over: adapted.over,
suffix: ` · smoothed ${smoothingCaption(adapted.applied)}`,
scope: adapted.scope,
}
}
// The baseline first, and exclusive: it is the same estimate the ratio would
// have used, shown instead of applied. `expected` is not consulted at all
// under it, so a reader who leaves the O/E toggle on and picks the baseline
// gets the baseline rather than a field of ones.
let smoothed: Smoothed
const arms = derivation.arms
if (derivation.baseline === true) {
values = arms === undefined ? expectedValues({ ...block, values }, runs) : armValues({ ...block, values }, arms)
parts.push(arms === undefined ? 'expected' : 'expected along the arms')
// Last, and it barely matters: a distance curve is a mean over many cells
// and is already smooth. Kept in the order for consistency rather than for
// effect — there is no speckle in a baseline to remove.
smoothed = despeckle(values)
values = smoothed.values
} else if (derivation.expected === true) {
// **The despeckle runs before the division, and this is the whole reason a
// smoothed O/E map is worth looking at.**
//
// `observedOverExpected` blanks a cell where nothing was observed, because
// zero over a positive baseline is the strongest depletion a ramp centred on
// 1 can draw and that is the loudest possible statement of the one thing the
// file does not say (§7.5). Away from the diagonal *most* cells are that
// cell — 98% of them at 200 bins out on a yeast map at 500 bp — so an O/E
// computed first hands the filter a field that is almost entirely `NaN`, and
// `NaN` is the one thing no despeckle will touch. Every setting then drew
// the same picture, because the blanks that make up the picture were not
// reachable by any of them.
//
// Smoothed first, a zero cell with contacts around it is no longer zero, so
// it is never blanked and the map fills in: the same 98% falls to 1% under
// the adaptive settings and to 39% at 7 × 7. This is exactly the argument
// `render/smooth.ts` already makes for putting the filter ahead of a *ratio
// between two files* — make the awkward cells finite before the division
// rather than trying to average them afterwards — applied to the division
// this module does itself.
//
// **The curve is still estimated from the unsmoothed block.** That is the
// half of "the despeckle goes last" worth keeping: a decay curve is an
// average over cells, and fitting it to pre-averaged ones would let each
// cell's neighbours count twice toward the baseline it is then compared
// against. Estimating before and dividing after keeps the estimator honest
// and still lets the filter do its work.
const stepped = { ...block, values }
if (arms !== undefined) {
// Both diagonals divided out: the same division as below, against the
// arm-frame model, which is estimated from the unsmoothed block and — once
// there is a despeckle — averaged over the same windows as the map, for
// exactly the reason the decay curve is.
const model = armValues(stepped, arms)
smoothed = smoothing === null ? despeckle(values) : despeckle(values, model)
values = ratioOf(
smoothing === null ? stepped.values : smoothed.values,
smoothing === null ? model : smoothed.over,
block.nrows * block.ncols,
derivation.keepZeros === true,
)
} else if (smoothing === null) {
smoothed = despeckle(values)
values = observedOverExpected(stepped, runs, undefined, derivation.keepZeros === true)
} else {
// **And the baseline is averaged over the same windows as the map.**
// Dividing a smoothed map by an unsmoothed curve is a different operation
// and a biased one: near the diagonal the decay is steep, so a window that
// spans it averages cells whose baselines differ several-fold and the
// ratio comes out systematically high — 24-32% at radii 3 to 7, an
// enrichment produced by the filter and by nothing in the file. Averaged
// over the same cells with the same divisor the bias is in both halves and
// cancels, and the mean sits back at 1 where an O/E belongs.
const model = expectedModel(stepped, runs)
smoothed = despeckle(values, expectedValues(stepped, runs, model))
values = ratioOf(smoothed.values, smoothed.over, block.nrows * block.ncols, derivation.keepZeros === true)
}
parts.push(arms === undefined ? 'O/E' : 'O/E along the arms')
} else {
smoothed = despeckle(values)
values = smoothed.values
}
// Named in the label like the steps above it, because it is the one of the
// three that a reader cannot see they asked for: a normalisation changes the
// numbers and a baseline changes the ramp, while a despeckled map looks like a
// map. A figure exported without the word on it is a claim about single-cell
// features that were averaged away.
//
// Appended *after* the name is settled rather than pushed onto `parts`, and
// that is the difference between an operation and an annotation: the steps
// above change what the quantity **is**, so any of them makes the word
// `counts` wrong; this one leaves it exactly what it was and says something
// about how it was measured.
const suffix = smoothed.suffix
if (smoothed.scope !== null) {
scope = scope === null ? smoothed.scope : `${scope} · ${smoothed.scope}`
}
return {
values,
// Only the baseline produces a middle. A normalisation rescales a count and
// leaves it a count — still a quantity with a floor at zero — which is why
// `VC` alone stays on a sequential ramp and `O/E` does not. The despeckle is
// not in this either: a mean of neighbours is in the units of its
// neighbours, so it moves no floor and creates no middle.
// The baseline is not in this either, and it is the near miss worth naming:
// it comes *out* of the ratio's denominator, but it is a mean of counts and
// so is a count-like quantity with a floor at zero and no middle. A neutral
// of 1 would put it on a diverging ramp centred where nothing in it lives.
neutral: derivation.expected === true && derivation.baseline !== true ? 1 : null,
label: (parts.length === 0 ? 'counts' : parts.join(' · ')) + suffix,
scope,
...(smoothed.pending === undefined ? {} : { pending: smoothed.pending }),
}
}
/**
* A field derived with {@link Derivation.keepZeros}, as the same derivation
* without it would have drawn it: every zero of the ratio made paper.
*
* That is the whole of the flag's effect, and only where the derivation divides
* — so the map's correlation derives its block once, sends it, and draws this
* under the veil, rather than deriving the same block twice (`ui/MapDerivation.ts`).
* Elsewhere a zero is a count and is left alone, and the array is returned as it
* came.
*/
export function paperZeros(values: Float32Array, derivation: Derivation): Float32Array {
if (derivation.expected !== true || derivation.baseline === true) return values
const out = new Float32Array(values.length)
for (let i = 0; i < values.length; i++) {
const value = values[i]
// `-0` too, which `=== 0` is true for, exactly as the division's own test.
out[i] = value === 0 ? Number.NaN : value
}
return out
}
/**
* The transform a derived quantity should be read on, where it has an opinion.
*
* A ratio wants `log2`, so that twice and half are the same distance either side
* of neutral; a difference is already symmetric and wants none. Offered as a
* *default* the panel adopts on the way in and gives back on the way out —
* nothing here forbids reading a ratio on a linear scale, it only refuses to
* make that the thing a reader is shown without asking.
*/
export function preferredTransform(neutral: number | null): 'log2' | 'none' | null {
if (neutral === null) return null
return neutral === 1 ? 'log2' : 'none'
}
/** The largest bin separation any cell of the block can have. */
function maxSeparation(rowStart: number, nrows: number, colStart: number, ncols: number): number {
const a = Math.abs(rowStart - (colStart + Math.max(0, ncols - 1)))
const b = Math.abs(rowStart + Math.max(0, nrows - 1) - colStart)
return Math.max(a, b)
}
/**
* The smallest bin separation any cell of the block can have.
*
* Zero wherever the two axes overlap, which is every on-diagonal block; the gap
* between them otherwise. Subtracting it is what keeps {@link ExpectedModel.cis}
* the size of the block rather than the size of the distance between its two
* windows.
*/
function minSeparation(rowStart: number, nrows: number, colStart: number, ncols: number): number {
const rowEnd = rowStart + Math.max(0, nrows - 1)
const colEnd = colStart + Math.max(0, ncols - 1)
if (rowStart > colEnd) return rowStart - colEnd
if (colStart > rowEnd) return colStart - rowEnd
return 0
}
/**
* The sequence each row or column belongs to, as an index into the run array.
*
* Compared by *label* rather than by run index, because the two axes of a
* rectangular block have two different run arrays and run 0 of one is not run 0
* of the other. `null` is "no run covers this index", which compares equal to
* itself and so keeps an unlabelled block on the cis curve — the right answer
* for a plain cooler read as one span, where every cell is within one sequence.
*/
function labelsByIndex(runs: readonly AxisRun[], length: number): Array<string | null> {
const out = new Array<string | null>(length).fill(null)
for (const run of runs) {
const end = Math.min(run.end, length)
for (let i = Math.max(0, run.start); i < end; i++) out[i] = run.label
}
return out
}
|