"""Peptides and nucleic acids from a sequence or HELM, as ordinary structures (#98).
A one-letter sequence or a HELM string becomes the full skeletal structure,
which is then an ordinary species: it is laid out, exported and calculated
like any other. RDKit builds both, with ``MolFromSequence`` and
``MolFromHELM``; this module checks the letters first, so a refusal names
the residue and its position instead of RDKit's silent failure.
HELM is Pistoia Alliance's Hierarchical Editing Language for Macromolecules:
Zhang et al., J. Chem. Inf. Model. 2012, 52, 2796 (doi:10.1021/ci3001925).
"""
from __future__ import annotations
import re
from dataclasses import dataclass, field, replace
from typing import Any, Literal
from vibesketch.document.model import Document, Species
from vibesketch.document.monomers import MODIFIED, MonomerRefused, curated_peptide
__all__ = [
"ALPHABETS",
"MAX_RESIDUES",
"SEQUENCE_FORMATS",
"Printed",
"SequenceRefused",
"helm_of",
"helm_molecules",
"plain_sequence",
"printed",
"sequence_export",
"sequence_molecules",
"with_sequences",
]
#: The import formats of this module.
SEQUENCE_FORMATS: tuple[str, ...] = ("peptide", "dna", "rna", "helm")
Kind = Literal["peptide", "dna", "rna"]
#: The one-letter codes each kind takes: the 20 standard amino acids, L, and
#: the four bases of DNA or RNA.
ALPHABETS: dict[str, str] = {
"peptide": "ACDEFGHIKLMNPQRSTVWY",
"dna": "ACGT",
"rna": "ACGU",
}
#: RDKit's ``MolFromSequence`` flavours: L-protein, uncapped RNA and DNA, so
#: peptides have free N and C termini and nucleic acids free 5'- and 3'-OH.
_FLAVOURS = {"peptide": 0, "rna": 2, "dna": 6}
_NAMES = {"peptide": "peptide", "dna": "DNA", "rna": "RNA"}
#: The longest sequence read; longer ones exceed the layout limit anyway.
MAX_RESIDUES = 2_500
[docs]
class SequenceRefused(ValueError):
"""The text is not a sequence of the stated kind."""
def _records(text: str) -> list[tuple[str, str]]:
"""``(where, sequence)`` per sequence: FASTA records, or else one per line."""
lines = [line.strip() for line in text.splitlines() if line.strip()]
if not any(line.startswith(">") for line in lines):
return [(f"line {i}", line) for i, line in enumerate(lines, start=1)]
records: list[tuple[str, str]] = []
for line in lines:
if line.startswith(">"):
records.append((f"record {line[1:].strip() or len(records) + 1!r}", ""))
elif not records:
raise SequenceRefused("a FASTA file starts with a '>' header line")
else:
where, sequence = records[-1]
records[-1] = (where, sequence + line)
return records
[docs]
def sequence_molecules(text: str, kind: Kind) -> list[Any]:
"""One molecule per sequence in ``text``: one per line, or FASTA records.
Letters may be in either case and spaced or numbered as in GenBank, so
``acgt acgt`` and ``1 ACGT`` both read; anything else is refused.
"""
from rdkit import Chem
alphabet = ALPHABETS[kind]
records = _records(text)
if not records:
raise SequenceRefused(f"no {_NAMES[kind]} sequence given")
molecules = []
for where, raw in records:
sequence = "".join(c for c in raw.upper() if not c.isspace() and not c.isdigit())
if not sequence:
raise SequenceRefused(f"{where} holds no {_NAMES[kind]} sequence")
for position, letter in enumerate(sequence, start=1):
if letter not in alphabet:
raise SequenceRefused(
f"{where}: {letter!r} at position {position} is not one of the "
f"{_NAMES[kind]} letters {alphabet}"
)
if len(sequence) > MAX_RESIDUES:
raise SequenceRefused(
f"{where} has {len(sequence):,} residues, more than the "
f"{MAX_RESIDUES:,} vibe-sketch reads"
)
mol = Chem.MolFromSequence(sequence, flavor=_FLAVOURS[kind])
if mol is None: # pragma: no cover - every checked letter is one RDKit knows
raise SequenceRefused(f"RDKit could not build the {_NAMES[kind]} at {where}")
molecules.append(mol)
return molecules
#: A HELM polymer: its type and number, and its units.
_POLYMER = re.compile(r"(PEPTIDE|RNA|CHEM|BLOB)(\d+)\{([^{}]*)\}")
def _check_every_residue(helm: str, mol: Any) -> None:
"""Refuse ``mol`` when RDKit left a unit of ``helm`` out of it (#115).
RDKit reads some monomers of its library, such as ``[Nal]`` in 2026.09,
without building them and joins the rest around the gap, so the
structure is silently another one. RDKit numbers each atom's residue
within its chain (A, B, ... in the order of the polymers), as in PDB
files: an N-terminal cap such as ``[ac]`` is residue -1 and the units
after it count from 1. A chain with fewer numbered residues than units
has lost one. Chains of other types than peptides and nucleic acids are
not numbered so and are not checked.
"""
polymers = _POLYMER.findall(helm.split("$", 1)[0])
if any(kind not in ("PEPTIDE", "RNA") for kind, _, _ in polymers):
return
found: dict[str, set[int]] = {}
for atom in mol.GetAtoms():
info = atom.GetMonomerInfo()
if info is None:
return # pragma: no cover - RDKit labels every atom it builds from HELM
found.setdefault(info.GetChainId(), set()).add(info.GetResidueNumber())
for index, (kind, number, text) in enumerate(polymers):
units = text.split(".")
present = found.get(chr(ord("A") + index), set())
if len(present) >= len(units):
continue
# Name the first unit whose number is missing, past a cap numbered -1.
shift = 1 if -1 in present else 0
missing = [
(position, unit)
for position, unit in enumerate(units, start=1)
if (-1 if shift and position == 1 else position - shift) not in present
]
position, unit = missing[0] if missing else (None, None)
where = f"{unit} at position {position}" if unit else "a unit"
raise SequenceRefused(
f"RDKit's HELM reader leaves out {where} of {kind}{number} in {helm!r}, so the "
"structure would be wrong; that monomer has no usable structure here"
)
[docs]
def helm_molecules(text: str) -> list[Any]:
"""One molecule per HELM string in ``text``, one per line."""
from rdkit import Chem
lines = [line.strip() for line in text.splitlines() if line.strip()]
if not lines:
raise SequenceRefused("no HELM given")
molecules = []
for line in lines:
if len(line) > 8 * MAX_RESIDUES:
raise SequenceRefused("the HELM string is longer than vibe-sketch reads")
mol = Chem.MolFromHELM(line)
if mol is not None:
_check_every_residue(line, mol)
if mol is None:
try:
mol = curated_peptide(line)
except MonomerRefused as error:
raise SequenceRefused(f"the HELM {line!r}: {error}") from None
if mol is None:
raise SequenceRefused(
f"RDKit could not read the HELM {line!r}: its syntax, or a monomer "
"neither RDKit's HELM library nor vibe-sketch's modifications have"
)
mol.SetProp(_GIVEN, line)
molecules.append(mol)
return molecules
# ── The recorded monomer form (#98) ──────────────────────────────────────
#: IUPAC three-letter symbols of the standard amino acids: "Nomenclature and
#: symbolism for amino acids and peptides", IUPAC-IUB JCBN Recommendations
#: 1983, Pure Appl. Chem. 56, 595 (doi:10.1351/pac198456050595).
_THREE_LETTER = {
"A": "Ala", "C": "Cys", "D": "Asp", "E": "Glu", "F": "Phe", "G": "Gly", "H": "His",
"I": "Ile", "K": "Lys", "L": "Leu", "M": "Met", "N": "Asn", "P": "Pro", "Q": "Gln",
"R": "Arg", "S": "Ser", "T": "Thr", "V": "Val", "W": "Trp", "Y": "Tyr",
} # fmt: skip
#: HELM monomers beyond the standard ones that figures print, by their IUPAC
#: symbols: each D-amino acid as D- and its L-form's symbol, and the less
#: common amino acids of the same recommendations that RDKit's HELM library
#: has (ornithine, sarcosine, norleucine, norvaline, 2-aminobutanoic acid).
_PRINTED_MONOMERS = {
**{f"d{letter}": f"D-{symbol}" for letter, symbol in _THREE_LETTER.items() if letter != "G"},
"Orn": "Orn", "Sar": "Sar", "Nle": "Nle", "Nva": "Nva", "Abu": "Abu",
# The curated modifications (vibesketch.document.monomers).
**{symbol: change.printed for symbol, change in MODIFIED.items()},
} # fmt: skip
#: Peptide chains, their links in the second section: the HELM a figure can
#: print as sequences when the links are disulfides within or between them.
_PEPTIDES = re.compile(r"((?:PEPTIDE\d+\{[^{}$]+\}\|?)+)\$([^$]*)\$\$\$(?:V2\.0)?")
_CHAIN = re.compile(r"PEPTIDE(\d+)\{([^{}$]+)\}")
#: A disulfide: residue to residue through their side chains (R3), in one chain or two.
_DISULFIDE = re.compile(r"PEPTIDE(\d+),PEPTIDE(\d+),(\d+):R3-(\d+):R3")
#: The monomers a disulfide joins: L- and D-cysteine.
_CYSTEINE = ("C", "[dC]")
[docs]
@dataclass(frozen=True)
class Printed:
"""A sequence as a figure prints it, with what it draws above the text.
A peptide of several chains is its first chain here, the others in
``chains``, each printed on a line of its own below the one before.
"""
text: str
#: Per residue of a peptide, the start and end of its symbol in ``text``.
spans: tuple[tuple[int, int], ...] = ()
#: Disulfide bridges, as pairs of residue indices into ``spans``, first < second.
bridges: tuple[tuple[int, int], ...] = ()
#: The chains after the first, in the order of the HELM.
chains: tuple[Printed, ...] = ()
#: Disulfides between chains, as (chain, residue) pairs, the earlier chain first;
#: chain 0 is this one and chain ``n`` is ``chains[n - 1]``.
links: tuple[tuple[tuple[int, int], tuple[int, int]], ...] = ()
#: What a peptide's ends print as: the free termini, or a curated cap.
_Termini = tuple[str, str]
@dataclass
class _Chain:
"""One peptide chain while its HELM is read."""
symbols: list[str]
termini: _Termini
#: How many positions the HELM counts before the first printed residue (an acetyl cap).
shift: int
monomers: list[str]
bridges: list[tuple[int, int]] = field(default_factory=list)
def _chain(text: str) -> _Chain | None:
"""One chain's printed residues and termini, or ``None`` for a monomer with no symbol here.
The curated caps (acetyl first, amide last) print as Ac- and -NH2 in
place of H- and -OH.
"""
from vibesketch.document.monomers import CAPS
monomers = text.split(".")
acetyl = monomers[0] == "[ac]"
amide = len(monomers) > 1 and monomers[-1] == "[am]"
termini = (CAPS["ac"] if acetyl else "H", CAPS["am"] if amide else "OH")
# Link positions count the caps; the residues between them are what is printed.
shift = int(acetyl)
monomers = monomers[shift : len(monomers) - int(amide)]
if not monomers:
return None
symbols = []
for monomer in monomers:
if monomer in _THREE_LETTER:
symbols.append(_THREE_LETTER[monomer])
elif (
monomer.startswith("[") and monomer.endswith("]") and monomer[1:-1] in _PRINTED_MONOMERS
):
symbols.append(_PRINTED_MONOMERS[monomer[1:-1]])
else:
return None
return _Chain(symbols, termini, shift, monomers)
def _residues(
helm: str,
) -> tuple[list[_Chain], list[tuple[tuple[int, int], tuple[int, int]]]] | None:
"""The peptide chains of ``helm`` with their disulfides, and the disulfides between them.
``None`` when a chain has a monomer without an IUPAC symbol here, or a
link is anything but a disulfide between two cysteines, each cysteine in
one at most, joining one chain or two printed one under the other.
Chains are in the HELM's order, numbered PEPTIDE1, 2, … .
"""
match = _PEPTIDES.fullmatch(helm)
if match is None:
return None
found = _CHAIN.findall(match.group(1))
if [int(number) for number, _ in found] != list(range(1, len(found) + 1)):
return None
chains = []
for _, text in found:
chain = _chain(text)
if chain is None:
return None
chains.append(chain)
links: list[tuple[tuple[int, int], tuple[int, int]]] = []
used: set[tuple[int, int]] = set()
for link in filter(None, match.group(2).split("|")):
bridge = _DISULFIDE.fullmatch(link)
if bridge is None:
return None
ends = []
for number, position in ((bridge[1], bridge[3]), (bridge[2], bridge[4])):
index = int(number) - 1
if not 0 <= index < len(chains):
return None
residue = int(position) - 1 - chains[index].shift
if not 0 <= residue < len(chains[index].monomers):
return None
if chains[index].monomers[residue] not in _CYSTEINE:
return None
ends.append((index, residue))
first, second = sorted(ends)
if first == second or {first, second} & used:
return None
used |= {first, second}
if first[0] == second[0]:
chains[first[0]].bridges.append((first[1], second[1]))
elif second[0] == first[0] + 1:
links.append((first, second))
else: # its line would cross the chains printed between them
return None
return chains, links
def _line(chain: _Chain) -> Printed:
"""One chain printed with its termini, H-Cys-…-OH."""
start, end = chain.termini
text, spans = start, []
for symbol in chain.symbols:
text += "-"
spans.append((len(text), len(text) + len(symbol)))
text += symbol
return Printed(f"{text}-{end}", tuple(spans), tuple(chain.bridges))
#: The molecule property that keeps the HELM a molecule was read from.
_GIVEN = "vibesketch_helm"
[docs]
def helm_of(mol: Any) -> str:
"""The HELM of a molecule RDKit built from a sequence or HELM.
RDKit's own HELM where it writes one; RDKit reads some monomers it does
not write (Sar, Abu, D-Thr, D-Ile in 2026.09), and those keep the HELM
they were read from.
"""
from rdkit import Chem
helm = str(Chem.MolToHELM(mol))
if not helm and mol.HasProp(_GIVEN):
helm = str(mol.GetProp(_GIVEN))
if not helm: # pragma: no cover - every sequence RDKit builds has a HELM it writes
raise SequenceRefused("RDKit could not write the structure's HELM")
return helm
[docs]
def plain_sequence(helm: str) -> tuple[Kind, str] | None:
"""The kind and one-letter sequence of ``helm``, or ``None`` if it says more.
A HELM is plain when it is one chain of standard residues without links:
exactly what that one-letter sequence builds again. D-amino acids,
other monomers, links such as disulfides and several chains are not.
"""
from rdkit import Chem
mol = Chem.MolFromHELM(helm)
if mol is None:
return None
letters = str(Chem.MolToSequence(mol))
kind: Kind
if helm.startswith("PEPTIDE1{"):
kind = "peptide"
elif helm.startswith("RNA1{"):
kind = "dna" if "[dR]" in helm else "rna"
else:
return None
if not letters or any(c not in ALPHABETS[kind] for c in letters):
return None
again = Chem.MolFromSequence(letters, flavor=_FLAVOURS[kind])
if again is None or str(Chem.MolToHELM(again)) != helm:
return None
return kind, letters
[docs]
def printed(helm: str) -> Printed:
"""``helm`` as a figure prints it: H-Tyr-D-Ala-Gly-OH, 5'-ACGU-3' or 5'-d(ACGT)-3'.
Peptides take IUPAC's three-letter symbols with their termini, D-amino
acids prefixed D- and the less common ones by their own symbols, and a
disulfide as a bridge above the two cysteines, as "Nomenclature and
symbolism for amino acids and peptides" (IUPAC-IUB JCBN 1983) draws
cystine. A peptide of several chains, such as insulin, prints each on a
line of its own, with a line joining the cysteines of each disulfide
between two chains. Nucleic acids take their one-letter bases, 5' to 3', DNA marked
d( ) as in IUPAC-IUB CBN's "Abbreviations and symbols for nucleic acids,
polynucleotides and their constituents" (1970).
"""
plain = plain_sequence(helm)
if plain is not None and plain[0] != "peptide":
kind, letters = plain
return Printed(f"5'-d({letters})-3'" if kind == "dna" else f"5'-{letters}-3'")
found = _residues(helm)
if found is None:
raise SequenceRefused(
f"{helm!r} is not one chain of residues, or several printed one under the other, "
"with IUPAC symbols linked at most by disulfides, so it cannot be printed as a "
"sequence; other links and other monomers need the structure"
)
chains, links = found
first, *others = (_line(chain) for chain in chains)
return replace(first, chains=tuple(others), links=tuple(links))
[docs]
def with_sequences(document: Document, molecules: list[Any]) -> Document:
"""``document`` with each species built from one of ``molecules`` recording its HELM."""
from vibesketch.chem.species import identity
from vibesketch.document.model import Sequence
found: dict[str, str] = {}
for mol in molecules:
key = identity(mol).inchikey
if key is not None:
found[key] = helm_of(mol)
species = [
s.model_copy(update={"sequence": Sequence(helm=found[s.inchikey], inchikey=s.inchikey)})
if s.inchikey in found
else s
for s in document.species
]
return document.model_copy(update={"species": species})
def _recorded(document: Document) -> list[tuple[Species, str]]:
"""Every species with its HELM; refuses when one has none."""
missing = [s.label or s.id for s in document.species if s.sequence is None]
if not document.species:
raise SequenceRefused("the document has no species")
if missing:
raise SequenceRefused(
f"{', '.join(missing)} {'has' if len(missing) == 1 else 'have'} no sequence: "
"only species imported from a sequence or HELM can be written as one"
)
return [(s, s.sequence.helm) for s in document.species if s.sequence is not None]
[docs]
def sequence_export(document: Document, fmt: Literal["helm", "fasta"]) -> str:
"""Every species as HELM, one per line, or as FASTA records headed by its label or id."""
recorded = _recorded(document)
if fmt == "helm":
return "".join(f"{helm}\n" for _, helm in recorded)
lines = []
for species, helm in recorded:
plain = plain_sequence(helm)
if plain is None:
raise SequenceRefused(
f"{species.label or species.id} is not one chain of standard residues, "
"which FASTA cannot hold; export it as HELM"
)
lines += [f">{species.label or species.id}", plain[1]]
return "\n".join(lines) + "\n"