Convert SMILES format to Amino Acid masses using RDKit - Cyclopeptides

This is a variation of a question asked on October 30, 2021 here, but it is different enough that I think it warrants a new post. I'm also trying to use RDKit.

Is it possible to convert a SMILES notation of a cyclopeptide into a list of masses, rather than a list of amino acid names? For example, Tyrocidine B1 has the following isomeric SMILES format:

CC(C)C[C@H]1C(=O)N[C@@H](C(=O)N2CCC[C@H]2C(=O)N[C@H](C(=O)N[C@@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N1)CCCCN)C(C)C)CC3=CC=C(C=C3)O)CCC(=O)N)CC(=O)N)CC4=CC=CC=C4)CC5=CNC6=CC=CC=C65)CC7=CC=CC=C7

The above question would try to write out the specific amino acids. Tyrocidine B1 would look something like the following:

cyclo[Asn-Gln-Tyr-Val-Lys-Leu-D-Phe-Pro-Trp-D-Phe]

I believe the cyclic nature of the peptide invalidates the answers provided in the above question. Furthermore, rather than listing the amino acids by name, I was hoping to obtain a list of amino acid masses. For Tyrocidine B1 it would look something like the following in python:

[114.04, 128.06, 163.06, 99.07, 128.09, 113.08, 147.07, 97.05, 186.08, 147.07]

I'm thinking that, because I'm focusing on the masses specifically and may deal with unknown amino acids in the future, it may be more straightforward to identify the peptide bonds and add up the mass of atoms in between them to achieve this goal. For example, if I found peptide bond 1 and peptide bond 2, I would add up all carbon, hydrogen, and nitrogen atoms (and possibly others) between those bonds to obtain the mass of that amino acid. I've been researching how best to do so, and feel like I've come close with the below code. However, I'm struggling with understanding the output of running a Mol object through GetSubstructMatches, and don't know how to use that to obtain amino acid masses.

from rdkit import Chem

tyroB1 = 'CC(C)C[C@H]1C(=O)N[C@@H](C(=O)N2CCC[C@H]2C(=O)N[C@H](C(=O)N[C@@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N1)CCCN)C(C)C)CC3=CC=C(C=C3)O)CCC(=O)N)CC(=O)N)CC4=CC=CC=C4)CC5=CNC6=CC=CC=C65)CC7=CC=CC=C7'
m = Chem.MolFromSmiles(tyroB1)

# identify peptide bonds
substructure = Chem.MolFromSmarts('NCC(=O)N')
print(m.GetSubstructMatches(substructure))


output:

((7, 8, 9, 10, 11), (11, 15, 16, 17, 18), (18, 19, 20, 21, 22), (22, 23, 24, 25, 26), (26, 27, 28, 29, 30), (30, 31, 32, 33, 34), (34, 35, 36, 37, 38), (38, 39, 40, 41, 42), (42, 43, 44, 45, 46), (46, 4, 5, 6, 7))


I'm optimistic in that the above code clearly identifies the correct number of peptide bonds (and does so with the different cyclopeptides I've tried), but I've stuggled moving from there to assembling a list of amino acid masses.

Would you have recommendations or ideas on how best to move forward, either in the direction I'm going or somewhere new? I would be obtaining information like SMILES formatting from PubChem, so anything that appears on cyclopeptide pages (like Tyrocidine B1 above) is fair game for use.

• Welcome to the site! It may help if you clarify what you are intending to use the masses for. While I don't know if the code from the other answer will convert a cycle to the amino acid names, if it can, you can just store a small dictionary of masses (and any other properties you care about) for each amino acid.
– Tyberius
Mar 25 at 13:42
• Thank you Tyberius! This is for compiling a dataset to test algorithms designed to sequence cyclopeptides from mass spec data (see bioinformaticsalgorithms.org/bioinformatics-chapter-4 for details). I found some cyclic peptide files on GNPS, but need a consistent way to obtain expected masses for testing. I haven't gotten the other answer to work on cyclics yet.Even if I get that working, It would be a middle man that needs to be updated with new amino acids, because cyclics are weird that way. Better, if possible, to just sum atoms, though mapping can work if there's no other way. Mar 26 at 15:15

I'm not sure if it's improper to answer my own question, but I think I've cracked it.

For each peptide bond match identified by GetSubstructMatches the index [1] is the atom index of the alpha carbon. If you find all atoms connected to that carbon, excluding the nitrogens from the peptide bond, you get all of the atoms contained in the amino acid. Throw in one of the excluded nitrogens and you can calculate the mass using the rdkit.Chem.Descriptors.ExactMolWt function. It accurately determined the sequences of Tyrocidine B1, Surugamide A and Surugamide B.

The below code hasn't been optimized, and it does NOT account for the possibility of connected side chains. I can't think of a reliable method to account for connected side chains just by comparing atoms and bonds, a prior knowledge of the specific amino acids in question would likely be required. Fortunately, that is not required for my current project, though additional insights would be welcome.

from rdkit import Chem
from rdkit.Chem.Descriptors import ExactMolWt
import numpy as np
from itertools import combinations

tyroB1 = 'CC(C)C[C@H]1C(=O)N[C@@H](C(=O)N2CCC[C@H]2C(=O)N[C@H](C(=O)N[C@@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N[C@H](C(=O)N1)CCCCN)C(C)C)CC3=CC=C(C=C3)O)CCC(=O)N)CC(=O)N)CC4=CC=CC=C4)CC5=CNC6=CC=CC=C65)CC7=CC=CC=C7'

m = Chem.MolFromSmiles(tyroB1)

peptide_bond_representation = Chem.MolFromSmarts('NCC(=O)N')

masses = []
for peptide_bond in m.GetSubstructMatches(peptide_bond_representation):

alpha = peptide_bond[1]
nitrogens = set([peptide_bond[0],peptide_bond[-1]])

aa_atom_idx = set([alpha])
set2 = set()
while aa_atom_idx != set2: # while loop repeats until all atoms in the amino acid are accounted for
set2 = aa_atom_idx.copy()
temp_am = am[:, list(aa_atom_idx)]
aa_atom_idx = set(np.where(temp_am==1)[0]) | aa_atom_idx
aa_atom_idx -= nitrogens # nitrogens from the peptide bond excluded to isolate atoms of this amino acid

bonds = [] # Chem.PathToSubmol requires a list of bonds, not atoms
for i,j in combinations(aa_atom_idx, 2):
b = m.GetBondBetweenAtoms(int(i),int(j))
if b: bonds.append(b.GetIdx())

tempMol = Chem.PathToSubmol(m, bonds)
mass = ExactMolWt(tempMol)
masses.append(mass)
print(masses)


Output: [147.068413908, 97.052763844, 186.07931294, 147.068413908, 114.042927432, 128.058577496, 163.063328528, 99.068413908, 128.094963004, 113.084063972]

• Self answers are fine, glad to see you were able to figure this out. I imagine you are right about the issue of connected side chains, as the only reason you can reliably divide up the cycle is that you know every atom between peptide bonds must belong to one amino acid. With connected side chains, you would have to decide arbitrarily where to split them or actually be able to identify the individual groups.
– Tyberius
Mar 27 at 13:19