"""Fabricate a peptide from alpha- or beta-amino acid residues."""
import warnings
try:
from pymol import cmd
import chempy.models
import pymol.plugins
except ImportError:
warnings.warn(
'Cannot import PyMOL: functionality will suffer (you can ignore this if you are just building the documentation).')
cmd = None
import re
import os
from .utils import tempObjectName, fetchModel, set_dihedral
from typing import Optional, List, Dict, Union
import logging
from .secstructdb import SecondaryStructureDB
import itertools
import random
logger = logging.getLogger()
logger.setLevel(logging.INFO)
class BetaPeptideException(Exception):
pass
class BetaPeptideConsistencyError(BetaPeptideException):
pass
class BetaPeptide:
"""A class for constructing alpha/beta peptides.
"""
GREEKLETTERS = 'ABGDEZHTIKLMNP'
SIDECHAINS = {
'A': 'ala',
'C': 'cys',
'CM': 'cys', # this will be treated differently in the code, since PyMOL does not have a ready-made fragment
'D': 'asp',
'DH': 'asph',
'E': 'glu',
'EH': 'gluh',
'F': 'phe',
'G': 'gly',
'HD': 'hid',
'HE': 'hie',
'HH': 'hip',
'H': 'his',
'I': 'ile',
'K': 'lys',
'O': 'orp', # this will be treated differently in the code, since PyMOL does not have a ready-made fragment
'KN': 'lysn',
'ON': 'orn', # this will be treated differently in the code, since PyMOL does not have a ready-made fragment
'L': 'leu',
'M': 'met',
'N': 'asn',
'P': 'pro', # proline, only available for alpha-peptides!
'Q': 'gln',
'R': 'arg',
'S': 'ser',
'T': 'thr',
'V': 'val',
'W': 'trp',
'Y': 'tyr',
}
def __init__(self, modelname: Optional[str], isbeta: bool = True, alphasidechain: str = 'G', alphastereo: str = 'S',
betasidechain: str = 'G', betastereo: str = 'S'):
"""Create a single amino-acid. Peptides can be built by *adding* them
:param modelname: the name under which this peptide will be known by PyMOL
:type modelname: str
:param isbeta: True if this is to be a beta-amino acid, False if an alpha
:type isbeta: bool
:param alphasidechain: the sidechain on the C-alpha
:type alphasidechain: str
:param alphastereo: chirality of the C-alpha
:type alphastereo: str, "R" or "S" (or "D" or "L" if `isbeta` is False, i.e. alpha-amino acids)
:param betasidechain: the sidechain on the C-beta
:type betasidechain: str
:param betastereo: chirality of the C-beta
:type betastereo: str, "R" or "S"
"""
self._modelname = modelname
if modelname is None:
return
if isbeta:
self.initializeBetaResidue(alphasidechain, alphastereo, betasidechain, betastereo)
else:
self.initializeAlphaResidue(alphasidechain, alphastereo)
def copy(self, newname: str) -> "BetaPeptide":
"""Make a deep copy of ourselves
:param newname: the new object name
:type newname: str
:return: the new BetaPeptide object
:rtype: BetaPeptide
"""
model = cmd.get_model('model {}'.format(self._modelname))
cmd.delete('model {}'.format(newname))
cmd.load_model(model, newname)
obj = BetaPeptide(None)
obj._modelname = newname
return obj
def initializeAlphaResidue(self, sidechain: str = 'G', stereo: str = 'S'):
"""Initialize this residue to be an alpha-amino acid
:param sidechain: designation of the sidechain
:type sidechain: str
:param stereo: absolute chirality of the sidechain ('R' or 'S')
:type stereo: str
:return: None
:rtype: NoneType
"""
cmd.delete('model {}'.format(self._modelname))
try:
cmd.fragment(self.SIDECHAINS[sidechain], self._modelname)
except pymol.CmdException as exc:
# could not load fragment: this might be a special one. Try to load it from our resource folder
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
cmd.load(os.path.join(resourcedir, '{}.pkl'.format(self.SIDECHAINS[sidechain])), self._modelname)
if sidechain == 'CM':
cmd.remove('model {} and name HG'.format(self._modelname))
cmd.alter('model {} and name SG'.format(self._modelname), 'formal_charge=-1')
cmd.alter('model {} and name SG'.format(self._modelname), 'partial_charge=-1')
cmd.alter('model {}'.format(self._modelname), 'resn="CYM"')
cmd.alter('model {}'.format(self._modelname), 'resv=1')
cmd.alter('model {}'.format(self._modelname), 'chain="A"')
# alpha-residues are by default all have L chirality. Now let's fix the chirality. The case of cisteine is
# complicated...
if sidechain != 'G' and ( # skip the achiral glycine
(stereo == 'D') or # all amino-acids if D is requested
(stereo == 'R' and sidechain not in ['C',
'CM']) or # for most amino-acids, L is S, thus for R, we have to flip
(stereo == 'S' and sidechain in ['C', 'CM']) # for Cys, L is R, thus for S, we have to flip
):
cmd.edit('model {} and resi {} and name CA'.format(self._modelname, 1),
'model {} and resi {} and name HA'.format(self._modelname, 1),
'model {} and resi {} and name C'.format(self._modelname, 1))
cmd.invert()
cmd.unpick()
# now rename the atoms
with fetchModel(self._modelname) as model:
for a in model.atom:
if a.name in ['N', 'H', 'C', 'O']:
a.name = a.name + '_sc0'
else:
a.name = a.name + '_sc2'
self.renameAtomsInResidue(1)
def checkPeptide(self):
"""Check this peptide for validity
The following tests are done:
1. each non-N-terminal peptide nitrogen must have one neighbour named "H" or "HN", except alpha-proline
2. each peptide nitrogen must have only one neighbour named "CB" or "CB1"
3. each peptide carbon must have only one neighbour named "CA"
4. each non-C-terminal peptide carbon must have exactly one neighbour named "N" with resi=resi+1
5. each non-C-terminal peptide carbon must have exactly one neighbour named "O"
6. each non-N-terminal peptide nitrogen must have exactly one neighbour named "C" with resi=resi-1
7. residue numbering must be continuous and start from one.
Capping groups are exempt from the above rules:
. ACE from rule 3.
. NME from rule 2.
:raises BetaPeptideConsistencyError: when one of the above tests fails
"""
with fetchModel(self._modelname) as model:
minresidue = min([a.resi_number for a in model.atom])
maxresidue = max([a.resi_number for a in model.atom])
if minresidue != 1:
raise BetaPeptideConsistencyError(
'Lowest residue number is {} instead of the expected 1.'.format(minresidue))
for i in range(1, maxresidue + 1):
resname = cmd.get_model('model {} and resi {}'.format(self._modelname, i)).atom[0].resn
logger.debug('Residue {}: {}'.format(i, resname))
if (i > 1) and resname not in ['PRO', 'DPRO']: # not the N-terminal and not proline
# count amide hydrogens on the nitrogen: must be one on a non-N-terminal residue
hcount = cmd.count_atoms('model {} and resi {} and (neighbor name N ) and (name H+HN)'.format(
self._modelname, i))
if hcount != 1:
raise BetaPeptideConsistencyError('Amide N in residue {} has {} hydrogens.'.format(i, hcount))
# count amide carbons bonded to the nitrogen: must be one on a non-N-terminal residue. This also ensures
# continuity of residue numbering (partly)
ccount = cmd.count_atoms('(neighbor (model {} and resi {} and name N)) and (resi {} and name C)'.format(
self._modelname, i, i - 1))
if ccount != 1:
raise BetaPeptideConsistencyError(
'Amide N in residue {} has {} amide carbons from the previous residue.'.format(i, ccount))
if i < maxresidue:
# not the C-terminal
# count amide oxygens on the carbon: must be one on a non-C-terminal residue
ocount = cmd.count_atoms('model {} and resi {} and (neighbor name C) and name O'.format(
self._modelname, i))
if ocount != 1:
raise BetaPeptideConsistencyError('Amide C in residue {} has {} oxygens.'.format(i, ocount))
# count amide nitrogens bonded to the carbon: must be one on a non-C-terminal residue. Together with
# the previous similar check, this ensures the continuity of residue numbering.
ncount = cmd.count_atoms(
'(neighbor (model {} and resi {} and name C)) and (resi {} and name N)'.format(
self._modelname, i, i + 1))
if ncount != 1:
raise BetaPeptideConsistencyError(
'Amide C in residue {} has {} amide nitrogens from the next residue'.format(i, ncount)
)
# all residues, including terminals
cbcount = cmd.count_atoms('model {} and resi {} and (neighbor name N) and (name CB+CB1+CA)'.format(
self._modelname, i))
if (cbcount != 1) and (resname != 'NME') and (resname != 'ACE') and (resname != 'BUT'):
raise BetaPeptideConsistencyError('Amide N in residue {} has {} CA or CB neighbours.'.format(
i, cbcount))
cacount = cmd.count_atoms('model {} and resi {} and (neighbor name C) and (name CA)'.format(
self._modelname, i))
if (cacount != 1) and (resname != 'ACE') and (resname != 'NME'):
raise BetaPeptideConsistencyError('Amide C in residue {} has {} CA neighbours.'.format(i, cacount))
@property
def maxResidue(self) -> int:
"""Get the highest residue number
:return: the highest residue number
:rtype: int
"""
with fetchModel(self._modelname) as model:
return max([a.resi_number for a in model.atom])
def isNProtected(self) -> bool:
"""Check if the N-terminal is protected
:return: True if the N-terminal is protected
:rtype: bool
"""
return cmd.count_atoms(
'(neighbor (model {} and name N and resi 1)) and not name CB+CB1+CA+H+HN+CH3'.format(self._modelname)) > 0
def isCProtected(self) -> bool:
"""Check if the C-terminal is protected
:return: True if the C-terminal is protected
:rtype: bool
"""
maxresidue = self.maxResidue
return cmd.count_atoms(
'(neighbor (model {} and name C and resi {})) and (not name O+CA+CH3)'.format(self._modelname,
maxresidue)) > 0
def startsWithProline(self) -> bool:
"""Check if the N-terminus is a proline residue
:return: True if the N-terminus is an alpha-proline
:rtype: bool
"""
with fetchModel(self._modelname) as m:
minresi = min({a.resi_number for a in m.atom})
firstresidue = {a.resn for a in m.atom if a.resi_number == minresi}
if len(firstresidue) > 1:
raise BetaPeptideConsistencyError('Atoms in the same residue must have the same residue names.')
else:
assert len(firstresidue) == 1 # len()==0 cannot happen.
return firstresidue.pop() in ['PRO', 'DPRO']
def __iadd__(self, right: "BetaPeptide") -> "BetaPeptide":
"""Append another peptide to our C-terminal
:param right: the C-terminal part
:type right: BetaPeptide
:return: the concatenated beta-peptide (the present one is overwritten).
:rtype: BetaPeptide
"""
left = self # the code below will be more logical this way.
if not isinstance(left, BetaPeptide) or not isinstance(right, BetaPeptide):
raise TypeError('Only peptides can be concatenated!')
if left.isCProtected():
raise ValueError('The N-terminal part is C-protected!')
if right.isNProtected() and (not right.startsWithProline()):
raise ValueError('The C-terminal part is N-protected!')
# do some validation
left.checkPeptide()
right.checkPeptide()
with tempObjectName() as mouldname, tempObjectName() as lefttemp, tempObjectName() as righttemp:
self.loadPeptideBond(mouldname)
mould = BetaPeptide(None)
mould._modelname = mouldname
left = left.copy(lefttemp)
right = right.copy(righttemp)
# increment the residue numbers in right
maxleft = left.maxResidue
with fetchModel(righttemp) as model:
for a in model.atom:
a.resi_number += maxleft
# mould is a peptide bond, with the following layout:
#
# O Cnext
# \\ /
# C --- N
# / \
# Cprev H
#
# align the following:
# C of the last residue of "left" to C of mould
# O of the last residue of "left" to O of mould
# CA or CH3 of the last residue of "left" to Cprev of mould
# in case of alpha-proline:
#
# O CA
# \\ / \
# C --- N CB
# / \ /
# Cprev CD - CG
logger.debug('Fitting left side')
logger.debug('Models: \n' + '\n'.join(cmd.get_object_list()))
logger.debug('Mould name: {}'.format(mouldname))
logger.debug('Lefttemp: {}'.format(lefttemp))
logger.debug('Righttemp: {}'.format(righttemp))
self.safe_pairfit('model {}'.format(lefttemp),
'model {}'.format(mouldname),
maxleft,
1,
['CA+CH3', 'C', 'O'],
['Cprev', 'C', 'O'],
0.1)
# find the neighbour of N
cbname = cmd.get_model(
'(neighbor (model {} and resi {} and name N)) and (symbol C) and (not name C)'.format(righttemp,
1 + maxleft)).atom[
0].name
# N of the first residue of "right" to N of mould
# H of the first residue of "right" to H of mould
# CB, CB1 or CA or CH3 of the first residue of "right" to Cnext of mould
logger.debug('Fitting right side')
self.safe_pairfit('model {}'.format(righttemp),
'model {}'.format(mouldname),
1 + maxleft,
1,
['N', 'H' if not right.startsWithProline() else 'CD', cbname],
['N', 'H', 'Cnext'],
0.1 if not right.startsWithProline() else 0.3)
cmd.fuse("model {} and resi {} and name C".format(lefttemp, maxleft),
"model {} and resi {} and name N".format(righttemp, 1 + maxleft),
mode=3)
cmd.bond("model {} and resi {} and name C".format(righttemp, maxleft),
"model {} and resi {} and name N".format(righttemp, 1 + maxleft), quiet=True)
obj = right.copy(self._modelname)
return obj
@staticmethod
def safe_pairfit(model1: str, model2: str, resi1: int, resi2: int, atoms1: List[str], atoms2: List[str],
rmstolerance: float = 0.1, ntries: int = 10):
"""Align model1 and model2 by the given atoms. Try to find an optimal tolerance.
The problem with the original pair_fit of PyMOL was that sometimes it found a local minimum, not
the desired one. If it converged into a local minimum and moved out by a random rotation, it
can find the global minimum.
:param model1: the name of the model to be moved
:type model1: str
:param model2: the name of the stationary model
:type model2: str
:param resi1: residue in model 1
:type resi1: int
:param resi2: residue in model 2
:type resi2: int
:param atoms1: name of atoms in the movable model
:type atoms1: list of str
:param atoms2: name of atoms in the stationary model
:type atoms2: list of str
:param rmstolerance: try to fit the models until the RMS is less than this
:type rmstolerance: float
:param ntries: the number of trials
:type ntries: int
:return: the final rms
:rtype: float
"""
selections1 = ['({}) and resi {} and name {}'.format(model1, resi1, name) for name in atoms1]
selections2 = ['({}) and resi {} and name {}'.format(model2, resi2, name) for name in atoms2]
selections = list(itertools.chain.from_iterable(zip(selections1, selections2)))
for i in range(ntries):
rms = cmd.pair_fit(*selections, quiet=True)
if rms < rmstolerance:
return rms
else:
# it may not have converged: do a random rotation and try to fit
cmd.rotate('xyz'[random.randint(0, 2)], random.randrange(0, 360), model1)
else:
raise RuntimeError('Could not fit with RMS less than {}. RMS is: {}'.format(rmstolerance, rms))
@staticmethod
def betaResidueName(sidechain2: Optional[str] = None, sidechain3: Optional[str] = None) -> str:
"""Get a residue name usable in GROMACS/CHARMM input files for a beta-amino acid
:param sidechain2: one (or two-)-letter abbreviation of the alpha-sidechain
:type sidechain2: str
:param sidechain3: one (or two-)-letter abbreviation of the beta-sidechain
:type sidechain3: str
:return: the residue name
:rtype: str
"""
if sidechain2 in [None, 'G'] and sidechain3 in [None, 'G']:
return 'BALA'
elif sidechain2 in [None, 'G']:
return 'B3' + sidechain3
elif sidechain3 in [None, 'G']:
return 'B2' + sidechain2
else:
return 'B0' + sidechain2 + sidechain3
def initializeBetaResidue(self, sidechain2: str = 'G', stereo2: str = 'S', sidechain3: str = 'G',
stereo3: str = 'S'):
"""Initialize this to be a beta-amino acid residue
:param sidechain2: designation of the side chain of the alpha carbon
:type sidechain2: str
:param stereo2: absolute conformation of the alpha side chain
:type stereo2: str, either 'R' or 'S'
:param sidechain3: designation of the side chain of the beta carbon
:type sidechain3: str
:param stereo3: absolute conformation of the beta side chain
:type stereo3: str, either 'R' or 'S'
:return: None
:rtype: NoneType
:raises ValueError: on invalid input
"""
residuename = self.betaResidueName(sidechain2, sidechain3)
if sidechain2 == 'P' or sidechain3 == 'P':
raise NotImplementedError('Beta-Proline residues not yet supported')
# delete the model
cmd.delete('model {}'.format(self._modelname))
# prime us with a homo-glycine
self.loadBetaBackbone(self._modelname)
# set the residue numbers to 1
with fetchModel(self._modelname) as model:
for a in model.atom:
a.resi_number = 1
a.resn = residuename
a.chain = 'A'
a.name = a.name + '_sc0' # add a tag identifying this atom to be part of the backbone
# attach alpha sidechain if needed
if sidechain2 != 'G':
self.attachBetaSideChain(1, 2, stereo2, sidechain2)
# rename the atoms
# attach beta sidechain if needed
if sidechain3 != 'G':
self.attachBetaSideChain(1, 3, stereo3, sidechain3)
# rename the atoms
self.renameAtomsInResidue(1)
cmd.alter('model {} and resi {}'.format(self._modelname, 1), 'chain="A"')
# Adjustment of the chirality is needed in the following special cases:
if sidechain2 in ['C', 'CM', 'S', 'T']:
# switch chirality of the alpha carbon
cmd.edit('model {} and resi {} and name CA'.format(self._modelname, 1),
'model {} and resi {} and name CB+CB1'.format(self._modelname, 1),
'model {} and resi {} and name C'.format(self._modelname, 1))
cmd.invert()
cmd.unpick()
if sidechain3 in ['C', 'CM', 'S', 'T']:
# switch chirality of the beta carbon
cmd.edit('model {} and resi {} and name CB+CB1'.format(self._modelname, 1),
'model {} and resi {} and name CA'.format(self._modelname, 1),
'model {} and resi {} and name N'.format(self._modelname, 1))
cmd.invert()
cmd.unpick()
if sidechain3 in ['M', 'D', 'DH'] and sidechain2 in ['G']:
# switch chirality of the beta carbon
cmd.edit('model {} and resi {} and name CB+CB1'.format(self._modelname, 1),
'model {} and resi {} and name CA'.format(self._modelname, 1),
'model {} and resi {} and name N'.format(self._modelname, 1))
cmd.invert()
cmd.unpick()
def attachBetaSideChain(self, residue: int, site: int, stereo: str, sidechain: str):
"""Attach a side-chain to the beta-backbone. The atoms in the side-chain (including the attach site and its other
hydrogen) will be correctly named.
:param residue: number of the residue to be modified
:type residue: int
:param site: the site to attach the sidechain: either 2 or 3
:type site: int
:param stereo: absolute conformation
:type stereo: str, either 'R' or 'S'
:param sidechain: sidechain designation
:type sidechain: str
:return: None
:rtype: NoneType
:raises ValueError: if `site` or `stereo` is incorrect
:raises KeyError: if `sidechain` is unknown
"""
with tempObjectName() as fragmentname:
# load the appropriate amino-acid fragment
try:
cmd.fragment(self.SIDECHAINS[sidechain], fragmentname)
except pymol.CmdException as exc:
# could not load fragment: this might be a special one. Try to load it from our resource folder
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
cmd.load(os.path.join(resourcedir, '{}.pkl'.format(self.SIDECHAINS[sidechain])), fragmentname)
# keep only the side-chain and Calpha
cmd.remove('model {} and (name C+O+N+H+HA)'.format(fragmentname))
# hack for residue 'CM'
if sidechain == 'CM':
cmd.remove('model {} and name HG'.format(fragmentname))
cmd.alter('model {} and name SG'.format(fragmentname), 'formal_charge=-1')
has_beta2sidechain = cmd.count_atoms(
'model {} and resi {} and name CB1'.format(self._modelname, residue)) == 1
# find the name of the backbone atoms where the sidechain will be attached
if site == 2 and stereo == 'S':
attachatoms = ('CA_sc0', 'HA1_sc0')
elif site == 2 and stereo == 'R':
attachatoms = ('CA_sc0', 'HA2_sc0')
elif site == 3 and stereo == 'S':
attachatoms = (
'CB1_sc0' if has_beta2sidechain else 'CB_sc0', 'HB11_sc0' if has_beta2sidechain else 'HB1_sc0')
elif site == 3 and stereo == 'R':
attachatoms = (
'CB1_sc0' if has_beta2sidechain else 'CB_sc0', 'HB12_sc0' if has_beta2sidechain else 'HB2_sc0')
else:
raise ValueError('Invalid site/stereo combination: {}/{}'.format(site, stereo))
# relabel the backbone atoms to avoid name clashes: simply add a "_sc0" tag at the end.
with fetchModel(self._modelname) as model:
for a in model.atom:
assert isinstance(a, chempy.Atom)
if a.resi_number != residue:
continue
if '_' not in a.name: # don't add two tags if one is already present.
a.name = a.name + '_sc0'
# relabel the sidechain atoms to avoid name clashes: add a "_sc<i>" tag at the end, where <i> is 2 or 3
with fetchModel(self._modelname) as model:
residuename = [a for a in model.atom if a.resi_number == residue][0].resn
with fetchModel(fragmentname) as model:
for a in model.atom:
assert isinstance(a, chempy.Atom)
a.name = a.name + '_sc{}'.format(site)
a.resi_number = residue
a.resn = residuename
# select the first heavy atom (i.e. not hydrogen) attached to the CA in the backbone-stripped sidechain
# model (typically CB)
with tempObjectName(nocreate=True) as tempsel:
cmd.select(tempsel,
'neighbor (model {0} and name CA_sc{1}) and not symbol H'.format(fragmentname, site))
caneighbour = cmd.get_model('%{}'.format(tempsel)).atom[0].name
# align the two models: fit CA and CB (or the first non-hydrogen attached to CA) in the backbone-stripped
# sidechain to the carbon and hydrogen on the beta-backbone chosen above for attaching.
cmd.pair_fit('(model {0} and name CA_sc{1}+{2})'.format(fragmentname, site, caneighbour),
'model {} and resi {} and name {}+{}'.format(self._modelname, residue, *attachatoms),
quiet=True)
cmd.remove('model {} and name CA_sc{}'.format(fragmentname, site))
cmd.remove('model {} and resi {} and name {}'.format(self._modelname, residue, attachatoms[1]))
cmd.fuse('model {} and name {}'.format(fragmentname, caneighbour),
'model {} and resi {} and name {}'.format(self._modelname, residue, attachatoms[0]))
# we are done with fusing, only renaming is needed.
def renameAtomsInResidue(self, residue: int) -> None:
"""Rename atoms in a residue
:param int residue: the residue number.
"""
with fetchModel(self._modelname) as model:
# first deconstruct the names of the atoms in the sidechains and the backbone into the form:
# <symbol>:<greek>:<index>:<sidechain>
#
# for example:
# CD1 in the alpha-sidechain -> C:D:1:2
# NE2 in the beta-sidechain -> N:E:2:3
# hydrogens are skipped for the time being.
for a in model.atom:
assert isinstance(a, chempy.Atom)
if a.resi_number != residue:
# skip atoms which do not belong to this residue
continue
if a.symbol == 'H':
# this is a hydrogen, skip it for now.
continue
m = re.match(
'^(?P<symbol>[BCNOFPS])(?P<greek>[{}])(?P<index>[123456789])?_sc(?P<sidechain>[023])$'.format(
self.GREEKLETTERS), a.name)
if m is None:
if a.name in ['C_sc0', 'O_sc0', 'N_sc0', 'H_sc0']:
# these are legal
continue
raise ValueError('Cannot match atom name: {}'.format(a.name))
m = m.groupdict()
if m['index'] is None:
m['index'] = 0
if m['sidechain'] == '3':
# increment the greek letters by one
m['greek'] = self.GREEKLETTERS[self.GREEKLETTERS.index(m['greek']) + 1]
a.name = '{0[symbol]}_{0[greek]}_{0[index]}_{0[sidechain]}'.format(m)
# now fix the numbering of heavy atoms. Sorting categories are (from strongest to weakest):
# - backbone, sidechain3, sidechain2
# - original index
# the names of the heavy atoms at this point are: <symbol>_<greek>_<index>_<sidechain>, where:
# - symbol is the element (e.g. C, N, O, S...)
# - greek is the greek letter (e.g A, B, G, D, E, Z, H, T ...)
# - index is the number of this atom in the original residue (e.g. the "2" from CD2)
# - sidechain is the sidechain designation, '0' being the backbone, '2' the alpha- and '3' the beta-sidechain
scsort = {'0': 0, '3': 1, '2': 2}
for g in self.GREEKLETTERS:
atoms = [a for a in model.atom if
a.symbol != 'H' and '_' in a.name and a.resi_number == residue and a.name.split('_')[1] == g]
if not atoms:
continue
elif len(atoms) == 1:
# this is the simplest case, no indexing required
a = atoms[0]
symbol, greek, index, sidechain = a.name.split('_')
a.name = symbol + greek
else:
logger.debug('Sorting heavy atoms for greek letter {}:'.format(g))
for i, a in enumerate(
sorted(atoms, key=lambda a: (scsort[a.name.split('_')[3]], a.name.split('_')[2]))):
logger.debug('{}: {}'.format(i + 1, a.name))
symbol, greek, index, sidechain = a.name.split('_')
a.name = symbol + greek + str(i + 1)
# re-fetch it from PyMOL, because we will use cmd.select for finding hydrogen neighbours of heavy atoms
with fetchModel(self._modelname) as model:
# now fix the hydrogens
for a in model.atom:
assert isinstance(a, chempy.Atom)
if a.symbol == 'H':
continue
if a.resi_number != residue:
continue
# find the atom to which this is bonded
atomindex = model.atom.index(a)
pairindices = [b.index for b in model.bond if atomindex in b.index]
neighbourindices = [x[0] if x[1] == atomindex else x[1] for x in pairindices]
hydrogens = [model.atom[idx] for idx in neighbourindices if model.atom[idx].symbol == 'H']
if not hydrogens:
continue
elif len(hydrogens) == 1:
hydrogens[0].name = 'H' + a.name[1:]
else:
# hydrogens are numbered as 1HG1, 2HG1, 3HG1 originally, i.e. the first character is the index.
# Another nomenclature is HG11, HG12, HG13. Either way, an alphabetical sort should do the trick of
# preserving the original order.
for i, h in enumerate(sorted(hydrogens, key=lambda h: h.name)):
h.name = 'H' + a.name[1:] + str(i + 1)
with fetchModel(self._modelname) as model:
# see if some atoms are still having some :tags at the end of their names. Must be the N:sc0, H:sc0, O:sc0, C:sc0
for a in model.atom:
if a.resi_number != residue:
continue
if a.name.endswith('_sc0'):
assert a.name.split('_')[0] in ['N', 'H', 'O', 'C']
a.name = a.name.split('_')[0]
assert '_' not in a.name
def isBetaResidue(self, residue: int) -> bool:
"""Check if the residue is a beta-amino acid
:param residue: residue index
:type residue: int
:return: True if it is a beta-amino acid
:rtype: bool
"""
if cmd.count_atoms('model {0} and resi {1} and name N+C+H+CA+O'.format(self._modelname, residue)) != 5:
# this might not even be an amino-acid!
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name N)) and resi {1} and name H'.format(self._modelname,
residue)) != 1:
# N-H bond missing
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name C)) and resi {1} and name CA'.format(self._modelname,
residue)) != 1:
# N-CA bond missing
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name C)) and (neighbor (model {0} and resi {1} and name N)) and (resi {1} and name CA)'.format(
self._modelname, residue)) == 1:
# N and C have a common neighbour, named CA. This is an alpha-amino acid
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name CA)) and (neighbor (model {0} and resi {1} and name N)) and (resi {1} and name CB+CB1)'.format(
self._modelname, residue)) == 1:
# N and CA have a common neighbour, named either CB or CB1. This is a beta-amino acid.
return True
return False
def isAlphaResidue(self, residue: int) -> bool:
"""Check if the residue is an alpha-amino acid
:param residue: residue index
:type residue: int
:return: True if it is an alpha-amino acid
:rtype: bool
"""
if cmd.count_atoms('model {0} and resi {1} and name N+C+H+CA+O'.format(self._modelname, residue)) != 5:
# this might not even be an amino-acid!
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name N)) and resi {1} and name H'.format(self._modelname,
residue)) != 1:
# N-H bond missing
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name C)) and resi {1} and name CA'.format(self._modelname,
residue)) != 1:
# N-CA bond missing
return False
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name C)) and (neighbor (model {0} and resi {1} and name N)) and (resi {1} and name CA)'.format(
self._modelname, residue)) == 1:
# N and C have a common neighbour, named CA. This is an alpha-amino acid
return True
if cmd.count_atoms(
'(neighbor (model {0} and resi {1} and name CA)) and (neighbor (model {0} and resi {1} and name N)) and (resi {1} and name CB+CB1)'.format(
self._modelname, residue)) == 1:
# N and CA have a common neighbour, named either CB or CB1. This is a beta-amino acid.
return False
return False
def fold(self, residue: int, dihedrals: Union[List[float], str]):
"""Fold a residue into the desired secondary structure
The secondary structure can be given in `dihedrals` in two ways:
1. a list of the torsion angles in degrees (2 or 3, depending on the type
of the amino-acid
2. the name of a secondary structure in the BETAFAB_HELIXTYPES PyMOL Plugin
variable
:param residue: the number of the residue
:type residue: int
:param dihedrals: secondary structure designation
:type dihedrals: list of floats or a str
:raises ValueError: if the residue is neither an alpha-, nor a beta-amino acid
"""
if isinstance(dihedrals, str):
dihedrals = SecondaryStructureDB.dihedrals(dihedrals)
logger.debug('Folding alpha-residue {} to {}'.format(residue, dihedrals))
if self.isAlphaResidue(residue):
logger.debug('Folding alpha-residue {} to {}'.format(residue, dihedrals))
logger.debug(
set_dihedral(self._modelname, ('C', residue - 1), ('N', residue), ('CA', residue), ('C', residue),
dihedrals[0]))
logger.debug(
set_dihedral(self._modelname, ('N', residue), ('CA', residue), ('C', residue), ('N', residue + 1),
dihedrals[1]))
elif self.isBetaResidue(residue):
logger.debug('Folding beta-residue {} to {}'.format(residue, dihedrals))
logger.debug(
set_dihedral(self._modelname, ('C', residue - 1), ('N', residue), ('CB+CB1', residue), ('CA', residue),
dihedrals[0]))
logger.debug(
set_dihedral(self._modelname, ('N', residue), ('CB+CB1', residue), ('CA', residue), ('C', residue),
dihedrals[1]))
logger.debug(
set_dihedral(self._modelname, ('CB+CB1', residue), ('CA', residue), ('C', residue), ('N', residue + 1),
dihedrals[2]))
else:
raise ValueError('Residue {} is neither an alpha-, nor a beta-amino acid'.format(residue))
@staticmethod
def loadBetaBackbone(objectname: str) -> str:
"""Create a new bare beta-backbone as a new object
:param objectname: desired name of the new object
:type objectname: str
:return: the name of the new object
:rtype: str
"""
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
cmd.load(os.path.join(resourcedir, 'betabackbone.pkl'), objectname)
return objectname
@staticmethod
def loadPeptideBond(objectname: str) -> str:
"""Create a new object: a bare peptide bond
:param objectname: desired name of the new object
:type objectname: str
:return: the name of the new object
:rtype: str
"""
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
cmd.load(os.path.join(resourcedir, 'peptidebond.pkl'), objectname)
return objectname
def _debugModel(self, model=None):
if model is None:
logger.debug('Fetching model')
with fetchModel(self._modelname) as model:
logger.debug('Fetched model.')
self._debugModel(model)
logger.debug('End of recursive call.')
else:
logger.debug('Atoms in the model:')
for a in model.atom:
logger.debug(' - name: {}, resi: {}, resn: {}'.format(a.name, a.resi_number, a.resn))
logger.debug('Bonds in the model:')
for b in model.bond:
logger.debug(' - {} - {} (order {})'.format(model.atom[b.index[0]].name, model.atom[b.index[1]].name,
b.order))
@classmethod
def But(cls, modelname=None):
obj = cls(None)
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
cmd.load(os.path.join(resourcedir, 'butyrate.pkl'), modelname)
cmd.alter('model {}'.format(modelname), 'resv=1')
cmd.alter('model {}'.format(modelname), 'chain="A"')
obj._modelname = modelname
return obj
@classmethod
def Ace(cls, modelname=None):
obj = cls(None)
cmd.fragment('ace', modelname)
for i in range(1, 4):
cmd.alter('model {} and name {}HH3'.format(modelname, i), 'name="HH3{}"'.format(i))
cmd.alter('model {}'.format(modelname), 'resv=1')
cmd.alter('model {}'.format(modelname), 'chain="A"')
obj._modelname = modelname
return obj
@classmethod
def NMe(cls, modelname=None):
obj = cls(None)
cmd.fragment('nme', modelname)
for i in range(1, 4):
cmd.alter('model {} and name {}HH3'.format(modelname, i), 'name="HH3{}"'.format(i))
cmd.alter('model {}'.format(modelname), 'resv=1')
cmd.alter('model {}'.format(modelname), 'chain="A"')
obj._modelname = modelname
return obj
@classmethod
def ACHC(cls, stereo2: str, stereo3: str, modelname: str = None):
obj = cls(None)
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
try:
cmd.load(os.path.join(resourcedir, 'ACHC_2{}3{}.pkl'.format(stereo2.upper(), stereo3.upper())), modelname)
except pymol.CmdException:
raise ValueError('Invalid stereochemistry: {} and {}'.format(stereo2, stereo3))
cmd.alter('model {}'.format(modelname), 'resv=1')
cmd.alter('model {}'.format(modelname), 'chain="A"')
cmd.alter('model {}'.format(modelname), 'resn="ACHC"')
obj._modelname = modelname
return obj
@classmethod
def ACPC(cls, stereo2: str, stereo3: str, modelname: str = None):
obj = cls(None)
resourcedir = os.path.join(os.path.split(pymol.plugins.plugins['pmlbeta'].filename)[0], 'resource')
try:
cmd.load(os.path.join(resourcedir, 'ACPC_2{}3{}.pkl'.format(stereo2.upper(), stereo3.upper())), modelname)
except pymol.CmdException:
raise ValueError('Invalid stereochemistry: {} and {}'.format(stereo2, stereo3))
cmd.alter('model {}'.format(modelname), 'resv=1')
cmd.alter('model {}'.format(modelname), 'chain="A"')
cmd.alter('model {}'.format(modelname), 'resn="ACPC"')
obj._modelname = modelname
return obj
@staticmethod
def parseBetaPeptideSequence(sequencestr: str) -> List[Dict[str, str]]:
"""Parse a beta-peptide sequence
The markup is the following (case sensitive):
1. beta-2 amino-acid (kind: B2):
({S|R})B2h<sc> e.g. (S)B2hV, (R)B2hA
2. beta-3 amino-acid (kind: B3):
({S|R})B3h<sc> e.g. (S)B3hL, (R)B3hR
3. beta-2,3 amino-acid (kind: B23):
(2{S|R}3{S,R})B23h(2<sc1>3<sc2>) e.g. (2S3R)B23h(2A3L)
4. bare beta-backbone (kind: BA):
BA or (2S3R)B23h(2G3G) (with any combination of R and S
in the first pair of parentheses)
5. alpha-amino acid (kind: A):
({S|R|L|D})A<sc> e.g. (S)AL, (R)AV, (L)AA, (D)AC
6. capping groups:
ACE: acetyl on the N-terminus (kind: ACE)
NME: N-methylamide on the C-terminus (kind: NME)
BUT: butyl on the N-terminus (kind: BUT)
Subsequent residues must be separated with comma. Whitespace around commas
and at the ends of the sequence string will be ignored.
After each residue a set of backbone torsion angles can be optionally
supplied in the form:
[ <angle1> <angle2> <angle3>]
i.e. a list enclosed in square brackets, separated by whitespace only. The angles
must be expressed in degrees, in normal notation (no exponential!) For example:
[-140.0 66 -136]
For alpha-amino acids two angles must be given, for betas three.
Employing the secondary structure database, the name of the desired fold can also be given
in curly braces, e.g.:
{H14M} or {Straight}
:param sequence: a string representation of the peptide sequence, following the above markup
:type sequence: str
:return: parsed alpha- or beta-amino acids
:rtype: a list of dicts having the following keys: 'kind', 'stereo2', 'stereo3', 'sidechain2', 'sidechain3'
:raises ValueError: if an element of a sequence cannot be interpreted
"""
float_regex = r"""([+-]\s*)? # optional sign and whitespace
( # start of the mantissa
(\d+(\.\d*)?) # one option: some digits, then optionally some decimals
| or
(\.\d+) # second option: only the decimals, led in by a decimal point
) # end of the mantissa
([eE][+-]?\d+)? # optionally, an exponent
"""
regex = re.compile(r"""^( # begin
(?P<ace>ACE) # acetyl group on the N-terminus
|
(?P<but>BUT) # butyl group on the N-terminus
|
(?P<nme>NME) # N-methylamide group on the C-terminus
|
\(2(?P<achcstereo2>[RS])3(?P<achcstereo3>[RS])\)ACHC # 2-aminocyclohexanecarboxylic acid
|
\(2(?P<acpcstereo2>[RS])3(?P<acpcstereo3>[RS])\)ACPC # 2-aminocyclopentanecarboxylic acid
|
( # begin beta23-amino acid regex
\(2(?P<stereo2>[RS])3(?P<stereo3>[RS])\) # chirality for beta2,3
B23h
\(2(?P<sc2>[A-Z]{{1,2}})3(?P<sc3>[A-Z]{{1,2}})\) # sidechain designation
) # end beta23-amino acid regex
| # or
( # begin monosubstituted beta-amino acid regex
\((?P<monostereo>[RS])\) # chirality for mono-substituted beta-amino acids
B(?P<monokind>[23])h
(?P<monosc>[A-Z]{{1,2}}) # sidechain designation
) # end monostubstituted beta-amino acid regex
| # or
(?P<bare>BA) # bare beta-amino acid (homo-glycine)
| # or
( # begin alpha-amino acid regex
\((?P<alphastereo>[RSLD])\) # chirality for alpha-amino acids
A
(?P<alphasc>[A-Z]{{1,2}}) # sidechain designation
))
\s* # allow for whitespace between the amino-acid designation and the secondary structure info
( # begin optional torsion angle declaration
(\[\s*(?P<phi>{0})\s+(?P<theta>{0})(\s+(?P<psi>{0}))?\s*\]) # two or three torsion angles in square brackets
| # or
(\{{\s*(?P<ssname>[-a-zA-Z0-9_+ ]+)\s*\}}) # name of a secondary structure database entry in curly braces
)? # end optional torsion angle declaration
$ # end alpha-amino acid regex
""".format(float_regex), re.X)
sequence = []
for aminoacid in [a.strip() for a in sequencestr.split(',')]:
m = regex.match(aminoacid)
if not m:
raise ValueError('Invalid amino-acid designation: {}'.format(aminoacid))
if m['phi'] is not None and m['theta'] is not None and m['psi'] is not None:
dihedrals = [float(m['phi']), float(m['theta']), float(m['psi'])]
elif m['phi'] is not None and m['theta'] is not None: # implies that 'psi' not in m
dihedrals = [float(m['phi']), float(m['theta'])]
elif m['ssname'] is not None:
dihedrals = SecondaryStructureDB.dihedrals(m['ssname'])
else:
dihedrals = None
if m['bare'] is not None:
if dihedrals is not None and len(dihedrals) != 3:
raise ValueError('Not enough dihedral angles specified in residue {}'.format(aminoacid))
sequence.append(
{'kind': 'BA', 'sidechain2': 'G', 'stereo2': '', 'sidechain3': 'G', 'stereo3': '',
'dihedrals': dihedrals})
elif m['alphastereo'] is not None:
if dihedrals is not None and len(dihedrals) != 2:
raise ValueError('Too many dihedral angles specified in residue {}'.format(aminoacid))
sequence.append(
{'kind': 'A', 'sidechain2': m['alphasc'], 'stereo2': m['alphastereo'],
'sidechain3': '', 'stereo3': '', 'dihedrals': dihedrals})
elif m['monostereo'] is not None:
if dihedrals is not None and len(dihedrals) != 3:
raise ValueError('Not enough dihedral angles specified in residue {}'.format(aminoacid))
if m['monokind'] == '2':
sequence.append(
{'kind': 'B2', 'sidechain2': m['monosc'], 'stereo2': m['monostereo'],
'sidechain3': 'G', 'stereo3': '', 'dihedrals': dihedrals})
else:
assert m['monokind'] == '3'
sequence.append(
{'kind': 'B3', 'sidechain2': 'G', 'stereo2': '',
'sidechain3': m['monosc'], 'stereo3': m['monostereo'], 'dihedrals': dihedrals})
elif m['ace'] is not None:
if dihedrals is not None:
raise ValueError('Residue ACE does not supports custom dihedrals.')
sequence.append({'kind': 'ACE', 'sidechain2': '', 'stereo2': '', 'sidechain3': '', 'stereo3': '',
'dihedrals': dihedrals})
elif m['but'] is not None:
if dihedrals is not None:
raise ValueError('Residue BUT does not supports custom dihedrals.')
sequence.append({'kind': 'BUT', 'sidechain2': '', 'stereo2': '', 'sidechain3': '', 'stereo3': '',
'dihedrals': dihedrals})
elif m['nme'] is not None:
if dihedrals is not None:
raise ValueError('Residue NME does not supports custom dihedrals.')
sequence.append({'kind': 'NME', 'sidechain2': '', 'stereo2': '', 'sidechain3': '', 'stereo3': '',
'dihedrals': dihedrals})
elif m['acpcstereo2'] is not None and m['acpcstereo3'] is not None:
sequence.append({'kind': 'ACPC', 'sidechain2': '', 'stereo2': m['acpcstereo2'], 'sidechain3': '',
'stereo3': m['acpcstereo3'], 'dihedrals': dihedrals})
elif m['achcstereo2'] is not None and m['achcstereo3'] is not None:
sequence.append({'kind': 'ACHC', 'sidechain2': '', 'stereo2': m['achcstereo2'], 'sidechain3': '',
'stereo3': m['achcstereo3'], 'dihedrals': dihedrals})
else:
if dihedrals is not None and len(dihedrals) != 3:
raise ValueError('Not enough dihedral angles specified in residue {}'.format(aminoacid))
assert m['stereo2'] is not None
sequence.append(
{'kind': 'B23', 'sidechain2': m['sc2'], 'stereo2': m['stereo2'],
'sidechain3': m['sc3'], 'stereo3': m['stereo3'], 'dihedrals': dihedrals})
return sequence
[docs]def betafab2(objname, *args):
""" Construct an alpha/beta peptide
Amino acids can be given in the following way:
1) beta-2 amino-acid:
({S|R})B2h<sc> e.g. (S)B2hV, (R)B2hA
2) beta-3 amino-acid:
({S|R})B3h<sc> e.g. (S)B3hL, (R)B3hR
3) beta-2,3 amino-acid:
(2{S|R}3{S,R})B23h(2<sc1>3<sc2>) e.g. (2S3R)B23h(2A3L)
4) bare beta-backbone:
BA or (2S3R)B23h(2G3G) (with any combination of R and S
in the first pair of parentheses)
5) alpha-amino acid:
({S|R|L|D})A<sc> e.g. (S)AL, (R)AV, (L)AA, (D)AC
6) capping groups:
ACE: acetyl on the N-terminus (must be first)
BUT: butyl on the N-terminus (must be first)
NME: N-methylamide on the C-terminus (must be last)
7) cyclic beta-residues:
(2{S|R}3{S|R})ACHC (2-aminocyclohexanecarboxylic acid) e.g. (2S3R)ACHC
(2{S|R}3{S|R})ACPC (2-aminocyclopentanecarboxylic acid) e.g. (2S3R)ACPC
where <sc> is the single-letter abbreviation of a natural
peptide sidechain. The following are supported:
A, C, D, E, F, G, I, K, L, M, N, O, Q, R, S, T, V, W, Y.
In addition, the following variants are also understood:
CM :
deprotonated cysteine
DH :
aspartic acid
EH :
glutamic acid
HD :
histidine protonated on the delta-nitrogen
HE :
histidine protonated on the epsilon-nitrogen
HH :
histidine protonated on both nitrogens
KN :
neutral lysine
ON :
neutral ornithine
The desired secondary structure can be given as well after
each amino-acid in two possible ways:
1. three (or two) floating point numbers in square brackets, e.g.
(S)AQ[-57 -47] or (S)B3hV[-140.3 66.5 -136.8]
2. referencing an entry of the secondary structure database in
curly braces, e.g. (S)AQ{Alpha-helix} or (S)B3hV{H14M}
:param objname: the name PyMOL will know about the resulting peptide
:type objname: str
:param args: amino-acid abbreviations
:type args: str
:return: the constructed peptide
:rtype: BetaPeptide
"""
sequence = [BetaPeptide.parseBetaPeptideSequence(a)[0] for a in args]
betapeptide = None
for ires, residue in enumerate(sequence):
with tempObjectName() as nextname:
if residue['kind'] == 'ACE':
nextresidue = BetaPeptide.Ace(nextname)
elif residue['kind'] == 'BUT':
nextresidue = BetaPeptide.But(nextname)
elif residue['kind'] == 'NME':
nextresidue = BetaPeptide.NMe(nextname)
elif residue['kind'] == 'ACHC':
nextresidue = BetaPeptide.ACHC(residue['stereo2'], residue['stereo3'], nextname)
elif residue['kind'] == 'ACPC':
nextresidue = BetaPeptide.ACPC(residue['stereo2'], residue['stereo3'], nextname)
else:
nextresidue = BetaPeptide(
nextname, residue['kind'] != 'A', residue['sidechain2'], residue['stereo2'],
residue['sidechain3'], residue['stereo3'])
if betapeptide is None:
betapeptide = nextresidue
else:
betapeptide += nextresidue
betapeptide = betapeptide.copy(objname)
# we are ready. Now we need to fold it. Folding is skipped to this point, because folding a residue needs some
# atoms from the previous and the next residues, respectively.
for ires, residue in enumerate(sequence, start=1):
if residue['dihedrals'] is not None:
betapeptide.fold(ires, residue['dihedrals'])
cmd.show_as('sticks', 'model {}'.format(objname))
return
def betafab2cmd(objname, *args, **kwargs): # kwargs is needed to swallow the _self argument given by PyMol
"""
DESCRIPTION
Build a beta peptide
USAGE
betafab2 objname [, aa1 [, aa2 [, aa3 [,...]]]]
ARGUMENTS
objname = str: target object name. Will be overwritten!
aa1, aa2 etc. = str: descriptors of the alpha/beta amino
acid residues. The following cases are understood:
1) beta-2 amino-acid:
({S|R})B2h<sc> e.g. (S)B2hV, (R)B2hA
2) beta-3 amino-acid:
({S|R})B3h<sc> e.g. (S)B3hL, (R)B3hR
3) beta-2,3 amino-acid:
(2{S|R}3{S,R})B23h(2<sc1>3<sc2>) e.g. (2S3R)B23h(2A3L)
4) bare beta-backbone:
BA or (2S3R)B23h(2G3G) (with any combination of R and S
in the first pair of parentheses)
5) alpha-amino acid:
({S|R|L|D})A<sc> e.g. (S)AL, (R)AV, (L)AA, (D)AC
6) capping groups:
ACE: acetyl on the N-terminus (must be first)
BUT: butyl on the N-terminus (must be first)
NME: N-methylamide on the C-terminus (must be last)
7) cyclic beta-residues:
(2{S|R}3{S|R})ACHC (2-aminocyclohexanecarboxylic acid) e.g. (2S3R)ACHC
(2{S|R}3{S|R})ACPC (2-aminocyclopentanecarboxylic acid) e.g. (2S3R)ACPC
where <sc> is the single-letter abbreviation of a natural
peptide sidechain. The following are supported:
A, C, D, E, F, G, I, K, L, M, N, O, Q, R, S, T, V, W, Y.
In addition, the following variants are also understood:
CM : deprotonated cysteine
DH : aspartic acid
EH : glutamic acid
HD : histidine protonated on the delta-nitrogen
HE : histidine protonated on the epsilon-nitrogen
HH : histidine protonated on both nitrogens
KN : neutral lysine
ON : neutral ornithine
The desired secondary structure can be given as well after
each amino-acid in two possible ways:
1. three (or two) floating point numbers in square brackets, e.g.
(S)AQ[-57 -47] or (S)B3hV[-140.3 66.5 -136.8]
2. referencing an entry of the secondary structure database in
curly braces, e.g.
(S)AQ{Alpha-helix} or (S)B3hV{H14M}
EXAMPLE
betafab2 test, (S)B3hV, (S)B3hA, (R)AK, (S)B3hL, (2S3S)B23h(2A3A)
betafab2 helix, (S)B3hV{H14M}, (S)B3hA{H14M}, (S)B3hL{H14M}, (2S3S)B23h(2A3A){H14M}, (S)B3hV{H14M}, (S)B3hA{H14M}, (S)B3hL{H14M}
"""
return betafab2(objname, *args)
if cmd is not None:
cmd.extend("betafab2", betafab2cmd)