Chem: SAR Matrix implementation

[talagayev]

Addition of Structure-Activity relationship matrix

[drizhina]

this is quite an overcomplicated way to link fragments. look for example in how its done in markush enumerator. you can just add a cyclization points to the corresponding molecules and then canonicalize smiles. will be much faster and result in same thing

[drizhina]

why two different operations? options can be passed from the addViewer call in second param

[drizhina]

make sure to check prs for such stuff)))

[drizhina]

A. The fallback path is uncapped, and it's the path big clusters take. assembleSinglePositionMatrix applies none of MAX_MATRIX_ROWS/MAX_COLUMNS_PER_POSITION: rows = every series in the cluster, columns = every distinct substituent. And decomposeCluster returns null for molIdx.length > 300 — so exactly the large clusters get routed to the unbounded builder. A 1000-molecule cluster with ~400 cores × ~600 substituents = 240 000 cell objects, and then:

computeMatrixConfidence → looPredictSlice copies the entire cells grid for every observed cell . That's O(observed × R × C) — hundreds of millions of object spreads.
linkVirtualStructures issues one RDKit parse per virtual cell — potentially 10⁵–10⁶.
The viewer then builds a with a canvas RDKit render for every cell.
This is a hard browser hang, not a slowdown. The caps belong in the fallback path too (or the fallback should refuse oversized clusters).

B. activityDirection is ignored by the initial ranking. SarMatrixParams has no direction field; runSarMatrix ends with rankMatrices(..., params.scaling === SCALING_METHODS.MINUS_LG) . The viewer's higherIsBetter getter — the one that honors "Higher is better"/"Lower is better" — is only consulted when the user touches the Rank by dropdown. So for a pre-computed pIC50 column on scaling: none + "Higher is better", the matrix list is sorted backwards on first render, while cardScore renders the correct direction — the card text and the card order disagree. The uncommitted working-tree change flipping the default scaling from -lg to none makes this the common case.

C. H is invisible. referenceValue and topValues both skip empty values (if (v)), and selectActivePositions counts only non-empty distinct values — yet isCleanSubstituent explicitly documents '' as "H at this position". Consequences: unsubstituted (H) can never be a column, never be the reference, and a position that varies H↔Me reads as having 1 distinct value and is dropped as inactive. Since the reference gates every cell (r.values[other] === refValues[other]), if H is the majority state at R2 the entire H-at-R2 population is excluded from the R1 group. This is the most likely explanation for "why is my matrix so empty" in practice, and H is the single most important reference point in real SAR.

D. Uncontrolled positions leak into the fit. heldAtRef only pins active positions , but MAX_POSITIONS = 2 means a molecule with 3+ varying R-groups contributes cells while its 3rd position varies freely. The virtual cell in the same row, meanwhile, is assembled with all non-active positions at their global reference (linkVirtualCellStructures uses Object.keys(refValues)). So the observed cells and the predicted structures aren't the same chemical series, and the Free-Wilson row/column effects are fit across compounds that differ at an unmodelled position. Either pin every decomposed position, or drop records that vary at an inactive one.

E. linkVirtualStructures runs on multi-position matrices too. In runSarMatrix it targets any virtual cell with smiles === null — including multi-position cells whose in-assembly link failed. It calls mmpLinkFragments, which does ${core}.${frag}.replaceAll('[:1]', '9') — a multi-attachment core leaves [:2] untouched, so you get either '' or a structure with a stray dummy atom, which then flows into the make-list export as a real proposal. Guard it on matrix.positions.length === 1.

F. Make-list append is schema-fragile. addAnalogToMakeList finds a table by name and does a positional rows.addNew(...) . The activity column name embeds this.scaling (Predicted activity (-lg)), so changing scaling between adds gives a mismatched schema; and any unrelated user table named SAR virtual analogs gets appended to.

G. Confidence mutates on the failure path. computeMatrixConfidence stamps cell.fit on cells before the pairs.length < MIN_CV_POINTS early return . So the chip says "Too few observations to cross-validate" while individual cells still render a red "non-additive" flag derived from a 1–3-point fit.

H. Clustering isn't permutation-stable. keyOf uses scaffolds[s.members[0].molIdx] on the premise that "all members of a single-cut series share the same ring system". That's false when a substituent contains a ring — a Murcko scaffold covers all ring systems, so Me vs. phenyl at the same site gives different scaffolds. members[0] is Map-insertion order, i.e. row order, so reordering the input table changes the clustering.

Methodologically questionable
I. The Free-Wilson estimator is the balanced-design closed form applied to a deliberately sparse design. rowMean + colMean − grandMean is the correct two-way additive estimate only for a complete/balanced matrix. A SAR matrix is by definition sparse and unbalanced, where this is biased — a row whose observations happen to sit in high-potency columns gets an inflated row effect. Since "predict the unmade analogs" is the whole product, this is the algorithm's soft centre. A ridge-regularized least-squares fit of row + column effects is a few dozen parameters, trivially solvable, strictly better, and would give real per-cell standard errors instead of the support = min(rowN, colN) heuristic. Worth noting a support = 1 prediction is literally oneRowValue + oneColValue − grandMean.

J. "SAR transfer" direction is an artifact. Pearson r is symmetric; a is simply whichever series gatherSeries emitted first (matrix index → position → row order). The UI states "a change learned on A should carry to B" and groups nav cards by transferSourceKey(t.a) — all of which is iteration order, not chemistry.

K. No multiple-testing control on a pairwise scan. computeAllTransfers tests every pair of (matrix, position, row) — easily thousands of hypotheses — with a floor of r ≥ 0.7 on as few as 3 points, where r = 0.7 has p ≈ 0.3. This will manufacture confident-looking cross-series transfers from noise. Raise MIN_COMMON, or gate on a permutation/FDR threshold.

L. foldMatch depends on arbitrary column order. It walks consecutive shared substituents, and columns are in frequency order — so the statistic changes if a substituent happens to appear one more time. And it's labelled "Fold-change match" while operating on raw differences: that's only a fold change on a lg/-lg scale, i.e. not under the new default scaling: none.

M. Ranking compares absolute potency, not SAR interest. discontinuityScore (max within-row range) and preferredScore (best column mean) are unnormalized activity magnitudes, so a series of uniformly potent compounds outranks a series with a genuine 100× cliff. Normalizing the spread by the dataset's activity SD would rank what the mode claims to rank.

• this.alignTemplate as an implicit parameter. Mutable viewer state reassigned per row in buildMatrixTable and in appendTransferSide, then read inside buildMatrixCell. Any reordering of the render silently misaligns structures. Should be a parameter.
• realCount means two different things — distinct molecules in the multi-position path (deduped via realMols), real cells in the single-position path. It drives the "N cpd" chip and the realCount > 0 filter.
• columns[ci].count++ on overwrite, so count can exceed observed cells — and SarMatrixColumn.count is never actually read by the UI (columnMeanPotency recomputes). Dead field.
• BitBIRCH gets MACCS with nFeatures left at the 2048 Morgan default . Numerically fine (the 235 pad bytes are zero on both sides of the Tanimoto), but it's a 256-byte stride for a 21-byte fingerprint — ~12× the memory and popcount work.
• coreFragId: i in the multi-position rows is a loop index, not a fragment id — the type lies.
• Murcko is a server round-trip of the entire molecule list (preceded by a full-column notation conversion), and failure degrades the analysis materially with only a logger.warning — the user sees no signal. It's also called as grok.functions.call('Chem:MurckoScaffolds', {smiles: 'smiles'}) passing a string where the generated murckoScaffolds() wrapper in package-api.ts expects a DG.Column.
• onPropertyChanged recomputes everything for every property, including rankScheme and activityDirection, which only affect sort/display. Changing "Activity direction" in the property panel re-runs fragmentation, Murcko, BitBIRCH, RGD and all the linking.
• Virtual-cell click sets dataFrame.currentRowIdx = -1 — a table-wide side effect that clears the current row for every other viewer.
• alignCache keys are templateMolblock + molStr — full multi-line molblocks as Map keys, module-global (shared across viewer instances and datasets), cleared wholesale at 4000.
• unusedTableName is an unbounded for (let i = 2; ; i++).
• seriesLabel(26) → "Series 27" right after "Series Z" (index 25); there is no "Series 26".
• Magic numbers throughout, none configurable: MIN_CV_POINTS, MIN_COMMON, MIN_CORRELATION, MAX_TRANSFERS, MAX_POSITIONS, MAX_COLUMNS_PER_POSITION, MAX_MATRIX_ROWS, MATCH_FRACTION, MAX_SAR_CLUSTER_SIZE. MATCH_FRACTION = 0.25 × matrix range in particular means a wide-range matrix essentially never flags a non-additive cell.
• Test coverage inverts the risk. All 18 tests hit pure functions with hand-built fragments; the end-to-end test is a skipReason stub. assembleMultiPositionMatrix, decomposeCluster, clusterRelatedCores and the reference-value filter — the riskiest code and the source of defects C, D, H — have zero coverage.
• Scope creep: the branch deletes initChemAutostart from package.ts, unrelated to SAR. WHY?