[lucas-diedrich]

Reviewing #457, I noticed that the logic in the same smiles.py module got reworked for the ~third time within a few releases. I believe that this is probably indicative that the current approach has some limitations

• It relies on a custom placeholder convention that we have to maintain
• Using placeholder atoms prevents us from doing a proper validation of the SMILES strings in our database (as the utilized atoms would actually follow a very different chemistry)
• To ensure the correct bond order etc. relies on a lot of custom logic. Ideally, we would like to "mark" the correct connection point explicitly (which we know in advance).

I wanted to quickly explore some alternatives.

This PR is fully vibed and is not meant to be merged but to start the discussion whether it makes sense to refactor the module. I think as a proof of principle it revealed some nice properties:

• RDKit can actually "mark" connection points with an Atom-map number, which makes joining PTMs with peptides much easier
• A lot of the handling logic for H-atoms etc. can be done by RDKit under the hood

Overall, this makes the merging of molecular subgraphs generalizable, i.e. one function (merge_at_anchors) can handle the connection of Modification-N terminus, Modification-C terminus, or residue-residue; and each modification in modification.tsv is a valid molecule on its own.

On this PR

The next comment summarizes my conversation with Claude and gives more technical details on the implementation of the old and new approach.

For reviewers, I would suggest to look at:

• smiles.py: Major changes to logic, move to RDKit logic. Slight reduction in LOC
• peptide.py: Highest impact, reduces the required custom logic the most
• adding_smiles.py: See the diff to get an intuition for the changes to the SMILES strings

The rest is just to have an end-to-end solution

[lucas-diedrich]

Summary of a design discussion and prototype on branch refactor/smiles. The prototype is meant to start a discussion and is not intended to be merged; it stays a draft PR.

1. Purpose of the modules

Module / file	Purpose
alphabase/smiles/smiles.py	Core chemistry. AminoAcidModifier loads amino-acid and modification SMILES from the constants tables and applies N-/C-terminal modifications to a residue (modify_amino_acid).
alphabase/smiles/peptide.py	PeptideSmilesEncoder turns a sequence plus alphabase modification strings (e.g. "Acetyl@Any_N-term;Oxidation@M", sites "0;2") into a peptide SMILES or RDKit molecule, or into one SMILES per residue.
alphabase/smiles/adding_smiles.py	Curation script. Holds dictionaries of new modification SMILES, derives each one's elemental composition from the SMILES, and writes it into modification.tsv.
constants/const_files/amino_acid.tsv, modification.tsv	The data. The smiles column holds residue and modification SMILES.
tests/test_smiles.py, tests/test_peptide_smiles.py	Unit tests for both modules.
docs/tutorials/adding_smiles.rst, docs/nbs/tutorial_smiles.ipynb, docs/nbs/peptide_smiles.ipynb	User guide and tutorials for the SMILES format.
scripts/convert_legacy_smiles.py (new in the prototype)	One-off utility that converts the tables from the old placeholder format to the new anchor format. It isn't imported by the library.

2. Current approach

Logical flow

Format: superheavy elements as placeholders for join points.

• Residues, including PTM residues, are written as the whole amino acid. Each amine hydrogen becomes [Fl] and the carboxylic OH becomes [Ts], e.g. alanine is N([Fl])([Fl])[C@@]([H])(C)C(=O)[Ts]. Proline has a single [Fl].
• N-terminal modifications follow one of two conventions:
  • With [Lv], which stands for the amine nitrogen itself: CC(=O)[Lv], C[Lv]C.
  • Without any placeholder: C(=O)CC. The first atom of the SMILES is implicitly bonded to the nitrogen.
• C-terminal modifications: the first atom replaces the OH, e.g. OC (methyl ester), N (amidation).

Applying an N-terminal modification (_apply_n_terminal_modification):

• [Lv] branch:
  1. Bond every neighbour of [Lv] to the residue N, copying the bond order.
  2. Copy [Lv]'s charge onto the N.
  3. Delete as many [Fl] as the total bond order (a double bond uses up two).
  4. Delete [Lv].
• No-placeholder branch: ReplaceSubstructs swaps the first [Fl] for the whole modification, attaching it through its atom 0.
• In both branches, any leftover [Fl] becomes an explicit [H].

Applying a C-terminal modification: ReplaceSubstructs replaces [Ts] with the modification, or with O if there is none.

Building a peptide (_build_peptide, _connect_amino_acids):

1. Residue 1 gets the N-terminal modification.
2. For each further residue:
   • find [Ts] in the chain and the first [Fl] in the next residue;
   • remove both and bond the carbonyl C to the N;
   • turn the remaining [Fl] into [H].
3. Apply the C-terminal modification last.

Issues

• Three placeholder conventions, each with its own code path and its own hydrogen bookkeeping: counting [Fl], bond-order arithmetic, transferring charge, cleaning up placeholders.
• The attachment atom is implicit for N-terminal mods without [Lv] and for all C-terminal mods (atom 0 of the SMILES). Wrong attachment points therefore go unnoticed. Because the elemental formula doesn't change, mass-based tests never catch them. Entries found that way:
  • Benzyl@Any_N-term and Benzyl-TA* attach through an aromatic ring carbon (N-tolyl, not N-benzyl).
  • SATA-/SATP-Succinamide attach through the thioacetyl methyl.
  • The *-TA-Succinamide glycan mods form an N–O bond to the carboxylic acid oxygen.
• Other data issues found:
  • Succinyl@K has the succinyl group on the α-amine instead of the ε-amine.
  • 12 residue entries lack α-carbon stereo: Oxidation@M, Phospho@T/Y, Carbamidomethyl@C/M, Pyridylethyl@C, Propionyl/Succinyl/Butyryl/Crotonyl/hydroxyisobutyryl@K, Pro->3,4-Dehydro-L-proline@P.
  • Carbamidomethyl@M is written as a sulfur ylide rather than a sulfonium.
• Latent stereo bug: a stereocentre directly bonded to [Lv] can invert depending on how the SMILES is written. No current entry is affected.
• Output artefacts: leftover [Fl] become explicit hydrogens, giving [H]N([H])….
• Table entries aren't real molecules, so they can't be checked against ChEBI or PubChem, and the composition isn't directly the formula.
• adding_smiles.py had drifted from the table. Three isotope-labelled Dimethyl entries held single-methyl SMILES that would overwrite the correct table values if the script were run.

3. New approach

RDKit concepts that simplify it

The rework uses hardly any new RDKit API; the old code already used CombineMols, RWMol editing and SanitizeMol. It works by delegating to RDKit's chemical model what the old code emulated with placeholder atoms.

RDKit concept	Effect
Implicit hydrogens: sanitizing recomputes hydrogen counts from valence.	Open valences are ordinary implicit hydrogens, reset after a merge (SetNumExplicitHs(0), SetNoImplicit(False)). All [Fl] bookkeeping, the bond-order arithmetic and the [H] artefacts disappear. This is the biggest simplification.
Atom-map numbers ([NH2:1]): a non-chemical annotation kept through edits and written to SMILES.	They mark join points on real atoms, so every table entry is a valid molecule. The map number also tells which site the SMILES belongs to.
Atoms can be edited in place (SetAtomicNum, SetFormalCharge, SetIsotope).	The anchor takes on the fragment anchor's element, charge and isotope. One rule covers amidation (O→N), imines, aromatic ring nitrogens and ¹⁵N labels.
Aromaticity: Kekulize, then re-detection while sanitizing.	Aromatic bonds become explicit single/double bonds before editing, and aromaticity is re-detected afterwards.
How chirality is stored: a chiral tag is relative to the order of an atom's bonds.	Moving a bond can change that order, so the merge inverts the tag when needed (InvertChirality). This fixes the latent stereo bug.
Sanitizing as validation (MolSanitizeException).	Impossible chemistry, such as dimethyl on proline, raises a clear ValueError.

Rejected alternatives:

• Chem.molzip: each pair of placeholder atoms makes exactly one bond, so it can't express two bonds to one N (dimethyl) or transfer a charge.
• Reaction SMARTS: these handle stereo internally, but the pattern matches any N, so the α-amine has to be targeted with extra code, and the templates are harder to curate.

Reactions remain an option if the team prefers not to own the stereo logic. The table format can stay as it is, because each mod SMILES can be turned into a reaction automatically.

New format

Every entry is a complete, valid molecule. Its join points carry atom-map numbers: :1 marks the α-amine nitrogen, and :2 marks the oxygen of the carboxylic OH.

Kind	Anchors	Example	Old (placeholder format)	New (anchor format)
Residue	:1 and :2	Alanine	N([Fl])([Fl])[C@@]([H])(C)C(=O)[Ts]	[NH2:1][C@@H](C)C(=O)[OH:2]
Residue	:1 and :2	Proline	N1([Fl])[C@@]([H])(CCC1)C(=O)[Ts]	[NH:1]1CCC[C@H]1C(=O)[OH:2]
PTM residue	:1 and :2	Phospho@S	O=P(O)(O)OC[C@@H](C(=O)[Ts])N([Fl])([Fl])	[NH2:1][C@@H](COP(=O)(O)O)C(=O)[OH:2]
N-terminal modification	:1	Acetyl@Any_N-term	CC(=O)[Lv]	[NH2:1]C(C)=O
N-terminal modification	:1	Dimethyl@Any_N-term	C[Lv]C	[NH:1](C)C
N-terminal modification	:1	Propionyl@Any_N-term (implicit attachment: atom 0)	C(=O)CC	[NH2:1]C(=O)CC
C-terminal modification	:2	Methyl@Any_C-term	OC	[OH:2]C
C-terminal modification	:2	Amidated@Any_C-term	N	[NH3:2]

New flow

One operation, merge_at_anchors(target, a, fragment, b), fuses the fragment's anchor atom into the target's anchor atom:

1. Combine the two molecules and convert aromatic bonds to explicit single/double bonds.
2. Move every bond of the fragment's anchor onto the target's anchor, keeping its bond order. Correct the chirality of any affected stereocentre.
3. The target anchor takes the fragment anchor's element, charge and isotope. Hydrogens are reset so RDKit recomputes them.
4. The fused atom keeps the label if a == b; otherwise it gets label 0.
5. Delete the fragment's anchor and sanitize. A valence violation raises a ValueError.

Everything is built from this one operation:

• N-terminal modification: merge(residue, :1, mod, :1).
• C-terminal modification: merge(residue, :2, mod, :2).
• Peptide bond: merge(chain, :2, next_residue, :1). The chain's OH oxygen becomes the next residue's nitrogen; recomputing hydrogens accounts for the lost water. The chain always keeps exactly one :1 and one :2.
• modify_amino_acid: parse → N-terminal merge → C-terminal merge → SMILES without atom maps.
• PeptideSmilesEncoder: fold the residues with peptide-bond merges, then apply the terminal mods. The per-residue variant processes each residue separately.

Example, acetyl-alanine: [NH2:1][C@@H](C)C(=O)[OH:2] + [NH2:1]C(C)=O → CC(=O)N[C@@H](C)C(=O)O. The old approach produced [H]N(C(C)=O)[C@@H](C)C(=O)O.

Prototype state (draft PR on refactor/smiles)

• Self-contained module: smiles.py has no knowledge of the old format and no backwards compatibility. Old-format SMILES fail the normal anchor validation.
• Public API unchanged: AminoAcidModifier and PeptideSmilesEncoder keep their signatures.
• Tables: converted with scripts/convert_legacy_smiles.py. Only the smiles column changed, and running the script on the original tables reproduces the converted ones byte for byte.
• Equivalence with the old implementation:
  • 2,992 modify_amino_acid calls (every amino acid, PTM residue, and terminal mod on every amino acid): identical molecules, stereo included. The 11 chemically impossible combinations raise in both versions.
  • 8 peptides, whole and per residue: identical.
• Tests: 1427 passed, 12 xfailed. New tests cover:
  • stereo next to an anchor, over shuffled atom orders (8 of them fail if the correction is disabled);
  • anchor and valence errors;
  • checks across every table entry, including a strict expected-failure (xfail) test for the 12 residues with missing α-carbon stereo.
• Data issues are deliberately left for a follow-up PR.

[lucas-diedrich]

Prompts

1. We try to represent post-translational modifications in Proteomics measurements as SMILES strings in constants/modifications.tsv to make them accessible for chemo-informatic processing. The logic to parse these modifications is encoded in @alphabase/smiles/smiles.py

   We have to jump through multiple hoops to get this to work and e.g. encode protein N/C-termini (via a placeholder atom), amines, etc.

   I am wondering whether there is a more elegant way to do this, e.g. by writing SMILES with placeholders and using the rdkit-intrinsic replacement logic to replace them with a suitable chemical strucuture.

   If you were to burn down the current module and reimplement it completely from scratch, what would you do? Look at the tests that handle the current logic and various edgecases.

2. Move this to a PR-ready change on the current branch

3. Why are the legacy functions in the smiles.py module still necessary? Please make the module self-contained (i.e. no backwards compatibility). If conversion etc. is needed, add this as a utility script to scripts/.

   This PR is just for educational purposes to spark a discussion how the module could get refactored. We will not merge it, it will remain a draft.

4. Now that you've done the rework, please explain

   1. The purpose of all modules that you touched
   2. The old logic
   3. The new logic

5. On what kind of RDKit-logic does this rely? Are there version constraints?

6. So assuming that the previous version of the modules was not fully aware of the full RDKit API - which core new logic does this rework rely on?

7. Summarize our conversation into a markdown file.

   1. Purpose of the modules
   2. Logical flow of the current approach, Issues with the current approach
   3. RDKit concepts that simplify the approach, new suggested flow

   Then add all prompts that I added to the bottom of the document