Convert SMILES format to Amino Acid masses using RDKit - Cyclopeptides

Question

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.

Answer 1

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')
am = np.array(Chem.GetAdjacencyMatrix(m))

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
    aa_atom_idx.add(peptide_bond[0]) #adding a single nitrogen back in

    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]