SequenceUtils

SequenceUtils is a package for the Datagrok platform that provides sequence alignment utilities for HELM-notation peptides and other macromolecules.

HELM MSA

Multiple Sequence Alignment engine for HELM-notation peptide sequences. Supports linear, cyclic, branched, and multi-chain structures with any number of non-canonical amino acids.

Features

• Aligns any HELM sequences: linear, cyclic (head-to-tail, lariat, stapled), branched, multi-chain
• Handles non-canonical amino acids as opaque tokens (no substitution matrix needed)
• Chains are paired across molecules by what they contain, not by polymer numbering
• Chains are emitted in the order the molecule reads, N-terminus first
• Automatic cyclic rotation detection and normalization
• CHEM-bridged macrocycle support (ring-closing metathesis, PEG linkers, etc.)
• Affine gap penalties with separate terminal gap control
• Gap chains (CHEM1{*}) inserted for missing branches and bonded into the place they stand for
• Connection positions automatically remapped after gap insertion and rotation
• Empty/null values in input columns are gracefully skipped
• Processes 2,800 sequences in under 1 second

Platform Integration

The HELM MSA function is registered with meta.role: sequenceMSA and is automatically discovered by the MSA dialog in the Bio package (Bio | Analyze | MSA...). When a user opens the MSA dialog on a HELM or non-canonical sequence column, HELM MSA appears as an available engine with its parameters editable in the dialog.

How It Works

                    HELM strings (input)
                          |
                [1] Parse each HELM
                          |
                [2] Match chains across molecules
                     by what they are made of
                          |
                [3] Order chains along the backbone
                     (N-terminus first)
                          |
                [4] Detect macrocycles
                     (head-to-tail, lariat, CHEM-bridge)
                          |
                [5] Rotate cyclic sequences
                     to best match reference
                          |
                [6] MSA: UPGMA (n<=50) or
                     center-star (n>50)
                          |
                [7] Reconstruct HELM
                     (gaps, connections, virtual chains)
                          |
                    HELM strings (output)

Step 1 - Parse: Each HELM string is parsed into polymers, connections, and metadata. The parser handles both bracketed ([meI]) and unbracketed (meI) monomer notation, splitting by dot separators.

Step 2 - Match Chains: Every polymer type gets as many alignment slots as the busiest molecule has chains of that type, and each molecule puts each of its chains into a distinct slot. Which chain goes where is decided by what the chain contains, because HELM numbering carries no meaning across molecules — see Chain Correspondence.

Step 3 - Order Chains: Slots are ordered the way the molecules read rather than by polymer type and number, so a conjugate bonded to a peptide's N-terminus takes its columns in front of the peptide and one bonded to the C-terminus takes them behind it. The order comes from the bonds: a bond with exactly one R1 end puts the R1 chain downstream. Every molecule votes and the majority direction wins.

Step 4 - Detect Macrocycles: The aligner checks for macrocyclic topology:

| Topology | Detection Rule | Rotatable? | |---|---|---| | Head-to-tail (R2->R1 self-connection) | Always | Yes | | Lariat (R1->R3 or R2->R3 self) | If span >= 70% of chain | Yes | | Backbone CHEM-bridge (R1/R2 through CHEM) | If ring span >= 50% | Yes | | Side-chain staple (R3->R3 through CHEM) | Never | No |

Step 5 - Cyclic Rotation: Cyclic peptides can be written starting at any ring position. The aligner tries all rotations and picks the one that best aligns to a reference.

Step 6 - MSA: Uses UPGMA-guided progressive alignment for small sets (n <= 50) and center-star alignment for large sets (n > 50). Both use Needleman-Wunsch with affine gap penalties, integer monomer encoding, and frequency-vector profile scoring for efficiency.

Step 7 - Reconstruct: Gap monomers (*) are inserted, connection positions remapped (terminal caps point to first/last positions, side-chains follow their monomers), and virtual gap chains added for missing polymers.

Chain Correspondence

Nothing in HELM says that the PEPTIDE1 of one molecule plays the same role as the PEPTIDE1 of the next. Numbering is whatever the tool that wrote the string happened to emit. A payload conjugate is the case where reading it as meaningful goes badly wrong:

PEPTIDE1{[PEGAA].Y.A.A.C.Y.Y.K.C.C}$$$$V2.0
CHEM1{[MCC]}|PEPTIDE1{E}|PEPTIDE2{[PEGAA].Y.A.A.C.Y.Y.K.C.C}
$CHEM1,PEPTIDE1,1:R2-1:R1|PEPTIDE1,PEPTIDE2,1:R3-1:R1$$$V2.0

The second molecule's PEPTIDE1 is a one-residue stub the linker bonds to; its side chain carries the peptide, written as PEPTIDE2. Pairing the two PEPTIDE1 chains asks the aligner to line one residue up against ten, which it can only do by opening a peptide-length gap in both rows:

*     PEGAA Y A A C Y Y K C C  *  *  *  *  *  *  *  *  *  *
MCC   *     * * * * * * * * *  E  PEGAA Y A A C Y Y K C C

Matching by content pairs the peptide with the peptide and gives the stub a column of its own:

*     *  PEGAA Y A A C Y Y K C C
MCC   E  PEGAA Y A A C Y Y K C C

The pairing is chosen to minimise the number of monomers left unpaired, over all one-to-one assignments of a molecule's chains to slots. Names are still the default: the content pairing has to cut the name-order pairing's cost by at least half (chainMatchMargin), and by at least a few monomers, before it is taken. So sets that were already matched correctly come out byte-identical, and a permutation of one-monomer staples — worth a monomer either way — never triggers a reshuffle.

Where the bonds settle a chain's place in its molecule, that place counts too. The two halves of a dimer are joined N-to-C, so which is which is never in doubt, and they are not swapped because their composition happens to be similar. Where the bonds settle nothing — three staples branching off one peptide have no order among themselves — the numbering is left alone rather than guessed at.

Set matchChainsByContent to false for the pre-existing name-order behaviour.

Gap Placeholders

A molecule missing a chain the others have gets an all-gap placeholder for it, so every row of the alignment has the same columns. Those placeholders are bonded into the place the chain they stand for occupies, rather than left as loose fragments — the row above reads as one connected molecule:

CHEM1{*}|CHEM2{*}|PEPTIDE2{*}|PEPTIDE1{[PEGAA].Y.A.A.C.Y.Y.K.C.C}
$CHEM1,CHEM2,1:R2-1:R1|CHEM2,PEPTIDE2,1:R2-1:R1|PEPTIDE2,PEPTIDE1,1:R3-1:R1$$$V2.0

The limit is reversibility: contracting the placeholders back out of an aligned molecule has to reproduce exactly the bonds the original had. A placeholder with one real neighbour dangles and contracts away to nothing, so it is always wired in. One with two real neighbours would contract to a bond between them, which is right only where it replaces a bond that is really there — a linker between two chains — and an invention anywhere else. A staple is the clearest case: it is bonded to its peptide twice, so a placeholder wired the same way would claim a ring between two residues that were never joined. Those are left unbonded.

Scoring

Identity-based scoring with affine gap model:

| Parameter | Default | Description | |---|---|---| | matchScore | 1 | Score when two monomers are identical | | mismatchScore | 0 | Score when two monomers differ | | gapOpenPenalty | -3 | Penalty for starting a new internal gap | | gapExtendPenalty | -0.5 | Penalty for extending an internal gap | | terminalGapOpenPenalty | -3 | Penalty for starting a gap at sequence ends | | terminalGapExtendPenalty | -0.5 | Penalty for extending a terminal gap |

Opening a gap is expensive; extending it is cheap. This encourages fewer, longer gaps.

Parameters

| Parameter | Type | Default | Description | |---|---|---|---| | sequenceCol | column | (required) | Column with HELM sequences to align | | gapOpen | double | -3 | Penalty for opening a new gap (internal) | | gapExtend | double | -0.5 | Penalty for extending an existing gap | | termGapOpen | double | -3 | Penalty for terminal gap opening | | termGapExtend | double | -0.5 | Penalty for terminal gap extension | | alignAllChains | bool | true | Align all chains or only the longest | | alignConjugateRuns | bool | true | Align a backbone-linked CHEM stack as one unit | | matchChainsByContent | bool | true | Pair chains by content rather than polymer numbering | | orderSlotsAlongBackbone | bool | true | Emit chains in molecular reading order |

When alignAllChains is true (default), every polymer chain type is aligned independently. Sequences missing a chain get a virtual gap chain. When false, only the longest chain is aligned and no virtual gap chains are created.

matchChainsByContent and orderSlotsAlongBackbone each restore the previous behaviour exactly when turned off, so they double as a way to A/B a dataset against the older output.

The API also accepts chainMatchMargin (default 0.5), the share of the name-order pairing's cost the content pairing must cut away before it is used. Raise it towards 1 to trust polymer numbering more, lower it to trust chain content more.

Supported HELM Topologies

Linear peptides:

PEPTIDE1{A.G.meI.hHis.C.K}$$$$

Cyclic peptides (head-to-tail):

PEPTIDE1{ClAc.F.Y.S.W.C.NH2}$PEPTIDE1,PEPTIDE1,6:R2-1:R1$$$V2.0

Lariat / lasso peptides:

PEPTIDE1{ClAc.F.Y.S.W.Sar.N.Y.W.S.Y.Y.Sar.W.C.NH2}$PEPTIDE1,PEPTIDE1,1:R1-15:R3$$$V2.0

Side-chain stapled (through CHEM linker):

PEPTIDE1{Ac.LysN3.R.R.dC.P.L.Y.I.S.Y.D.P.V.Pen.R.R.NH2}|CHEM1{CH2}
$PEPTIDE1,CHEM1,5:R3-1:R1|PEPTIDE1,CHEM1,15:R3-1:R2$$$V2.0

CHEM-bridged macrocycles:

PEPTIDE1{A.G.C.K.L.M.N.P}|CHEM1{CH2}
$PEPTIDE1,CHEM1,1:R1-1:R1|PEPTIDE1,CHEM1,8:R2-1:R2$$$

Multi-chain with branches:

PEPTIDE1{A.G.C.K.L}|PEPTIDE2{F.W.Y}|CHEM1{PEG}
$PEPTIDE1,PEPTIDE2,3:R3-2:R3|PEPTIDE1,CHEM1,4:R3-1:R1$$$

Payload conjugates (linker, stub, peptide):

CHEM1{[MCC]}|PEPTIDE1{E}|PEPTIDE2{[PEGAA].Y.A.A.C.Y.Y.K.C.C}
$CHEM1,PEPTIDE1,1:R2-1:R1|PEPTIDE1,PEPTIDE2,1:R3-1:R1$$$V2.0

Capped peptides:

PEPTIDE1{dAsp.dPro.dAla...}|CHEM1{Ac}|CHEM2{NHMe}
$CHEM1,PEPTIDE1,1:R1-1:R1|CHEM2,PEPTIDE1,1:R1-18:R2$$$V2.0

Performance

| Optimization | Impact | |---|---| | Monomer integer encoding | All comparisons are === on Int32Array | | K-mer distance | O(L) per pair instead of O(L^2) for guide tree | | Center-star MSA (n > 50) | O(n * L^2) instead of O(n^2 * L^2) | | TypedArray DP matrices | Flat Float64Array, cache-friendly | | Frequency vector profiles | O(|alphabet|) scoring instead of O(k1 * k2) |

Benchmark: 2,800 sequences of ~18 monomers each align in ~0.7 seconds.

Source Structure

src/
  package.ts              Entry point, registers HELM MSA function
  helm-msa/
    types.ts              Type definitions, AlignmentOptions, defaults
    helm-parser.ts        Parse/serialize HELM notation (dot-separator-aware)
    encoding.ts           Map monomer strings to integers for fast comparison
    nw-align.ts           Needleman-Wunsch (full, score-only, integer-encoded, profile-profile)
    msa.ts                Progressive MSA (UPGMA + center-star), k-mer distance
    chain-matcher.ts      Match polymer chains across molecules by content
    linear-conjugates.ts  Recover the backbone reading order from the connection list
    index.ts              Orchestrator: parse -> match -> order -> rotate -> align -> reconstruct
  tests/
    helm-msa-tests.ts         Algorithm units: parser, NW, MSA, rotation, remapping
    helm-msa-csv-tests.ts     Integration against files/*.csv
    linear-conjugate-tests.ts Conjugate runs: terminal stacks, linkers, gap wiring
    chain-matching-tests.ts   Chain correspondence, numbering, slot order
    helm-test-utils.ts        Shared invariant checks

Testing

grok test --host localhost                           # Run all tests
grok test --host localhost --category "HELM MSA: Parser"  # Run a specific category

277 tests across 64 categories covering:

• HELM parsing and serialization (bracketed, unbracketed, multi-char monomers)
• Needleman-Wunsch correctness and score optimality
• Progressive MSA (UPGMA and center-star)
• Cyclic rotation detection and alignment
• CHEM-bridged cycles and side-chain staples
• Connection remapping (caps, cycles, branches)
• Terminal cap R-group correctness
• Gap penalty behavior (strict vs relaxed)
• Chain correspondence: payload conjugates, renumbered polymers, dimer halves, staple permutations
• Backbone slot order for N-terminal and C-terminal conjugates
• Edge cases (single monomer, CHEM-only, empty input)
• CSV integration tests against real datasets (rotated_helm, pseqs_random, bicyclic, multiple_conjugates, 2800+ sequences)

Please note, these test csv files are not included in the package bundle. Contact the developer for files