PDB Reader

Background

The Protein Data Bank (PDB) format allows for a relatively straightforward process of describing biomolecular structures. The RCSB protein data bank (RCSB PDB) has information on thousands of structures. Unfortunately, simulating these structures is not always straightforward, as their structures do not always model all relevant info explicitly. Modifying the structures can be especially difficult, even if you know what you are doing.

Most MolCube projects start with a PDB Reader project. PDB Reader parses a .pdb or .cif file, determines what structures already exist in the file, and allows manipulating the structure to add or remove features like Disulfide bonds, Mutations, Glycosylations, etc.

A finished PDB Reader project can then be used for more complex operations like Solvation, embedding in a Membrane, and calculating Free Energy.

Overview

The general procedure for using PDB Reader in MolCube-API Client goes like this:

Authenticate server connection.
Create project.
Select chains.
Manipulate structure (optional).
Finalize model.
Download project (optional).

The example below shows how this works in the simplest case by using only the default settings that you’d see on the MolCube Apps site. This is equivalent to entering a PDB ID and just clicking “Next” until you get to the final page, then clicking “Download Project”.

import molcube as mc
from pprint import pprint

api_access_key = 'your-api-key'
molcube = mc.API('api.molcube.com', 443)
molcube.authenticate(api_token=api_access_key)

#
# Initialize project by downloading structure from RCSB
#
pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='test-defaults', ff='charmmff', pdbId='2hac')
pdbreader.set_defaults()

#
# if modifying chain selection, do so here
#
assert pdbreader.confirm_chains()

#
# if modifying manipulation options, do so here
#
assert pdbreader.model_pdb()

pdbreader.download_project('myproject.tgz')

The assert keyword is prepended to each command that submits a step. This prevents the script from proceeding if a step fails, though it is not required.

Setup PDB Reader project

Creating a PDB Reader project requires setting a force field and project title, and either providing a RCSB pdbId or uploading a customPdb. This is done by calling the create_project() method.

The following arguments are available for the method:

correct_topo (bool): Correct chains and bonds information using distance between each atom. (default: False)
rename_dupl_atoms (bool): Rename hetero atoms if there are duplicate atom names. (default: True)
calc_pka (bool): Calculate pKa of protein residues to apply system pH. (default: False)

Alternatively, if your project already exists, e.g. because you created it interactively in MolCube Apps, you can use resume_project() as shown below.

Create PDB Reader Project

Let’s walk through how to create a PDB Reader project.

Two force field options are available: charmmff and amberff. Although MolCube supports martiniff and drudeff, these options are not yet supported in the Python API client.

Most code examples assume you are using charmmff. For amberff, you can select different force field options by passing an amberOptions dict to create_project() like so:

pdbreader = molcube.create_pdb_reader_project()
pdbreader.create_project(title='test', ff='amberff', pdbId="2hac", amberOptions={
    "protein": "FF19SB",
    "dna": "OL15",
    "rna": "OL3",
    "glycan": "GLYCAM_06j",
    "lipid": "Lipid21",
    "water": "OPC"
})

Here are all available choices for amberOptions:

Protein: [FF19SB, FF14SB, FF14SBonlysc]
DNA: [OL15, BSC1]
RNA: [OL3, YIL, Shaw]
Glycan: [GLYCAM_06j]
Lipid: [Lipid21, Lipid17]
Water: [OPC, TIP3P, TIP4PEW, TIP4PD]

You can also find these options in the molcube.pdbreader.enums module:

from molcube.pdbreader import enums

print(f"Protein options: [{', '.join(enums.Protein)}]")
print(f"DNA options: [{', '.join(enums.DNA)}]")
print(f"RNA options: [{', '.join(enums.RNA)}]")
print(f"Glycan options: [{', '.join(enums.Glycan)}]")
print(f"Lipid options: [{', '.join(enums.Lipid)}]")
print(f"Water options: [{', '.join(enums.Water)}]")

Fetch PDB from RCSB

Using the pdbId keyword argument will attempt to obtain the PDB from RCSB automatically. create_project() returns True on success:

# Create a PDB Reader project and fetch PDB from RCSB using PDB ID
pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='test', ff='charmmff', pdbId="2hac")

Upload a custom PDB file

If you already have a local copy of your structure, you can pass the path to your structure with the customPdb keyword argument:

pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='test', ff='charmmff', customPdb="files/2hac.cif")

MolCube recognizes structures in PDB (.pdb), PDBx/mmCIF (.cif), and GROMACS (.gro) formats.

Resume an existing project

An existing project can be resumed by passing the project ID to resume_project():

pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.resume_project(project_id='b33384ed-7e4e-48cb-9afd-b2f0fe6456a2')

Search project list

If you don’t already know the ID, you can use search_projects().

Acceptable args:

page (int): page to return (default: 1)
perPage (int): number of results per page (default: 10)
keyword (str): limit results to those containing a keyword
searchKey (title|pk): restrict keyword search to either the title string or project ID (pk).

For the complete list of accepted arguments, see search_projects() in the API Reference for more detail.

Example:

>>> search_results = molcube.search_projects()
>>> search_results
{'projects': [{'pk': '1eebb792-f267-4c01-9f9f-1b179819c3f3',
   'createdAt': '2026-04-02 13:51:33',
   'forcefieldType': 'charmmff',
   'projectType': 'PDB Reader',
   'title': 'My Test Project',
   'step': 2,
   'status': 'Success',
   'fileName': '2hac.cif',
   'sideChainOriented': False,
   'tag': None,
   'user': 'Your User Name',
   'team': None,
   'teamId': None,
   'workspace': 'Personal'},
  {'pk': '205dffb4-2b15-43ce-9228-305e6ae510a6', ...},
  ...,
 ],
 'totalPages': 2,
 'currentPage': 1,
 'totalCount': 11,
 'hasNext': True,
 'hasPrevious': False}

E.g., to resume the most recent project, use the first pk from the returned object:

my_projects = search_results['projects']
project_id = my_projects[0]['pk']

pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.resume_project(project_id=project_id)

Check available info about PDB

The get_chains() method returns a list of chains for each chain type and (where applicable) the available terminal caps:

>>> pdbreader.get_chains()
{'protein': [{'chainIndex': 'PROT_A',
   'terminal': {'nter': ['NTER', 'NNEU', 'ACE', 'NONE'],
    'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE']},
   'nsdTerminal': {'nter': ['ACE'], 'cter': ['NONE']},
   'chainId': 'A'},
  {'chainIndex': 'PROT_B',
   'terminal': {'nter': ['NTER', 'NNEU', 'ACE', 'NONE'],
    'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE']},
   'nsdTerminal': {'nter': ['ACE'], 'cter': ['NONE']},
   'chainId': 'B'}],
 'nucleicAcid': [],
 'standaloneLigand': [],
 'heme': [],
 'ion': [],
 'glycan': [],
 'water': []}

The get_pdb_info() method returns a large dict with all info the MolCube server was able to parse from the structure:

>>> pdb_info = pdbreader.get_pdb_info()
>>> pdb_info.keys()

dict_keys(['ph', 'pdbId', 'source',
   'forceFieldType', 'models', 'availResnames',
   'resnames', 'titrableResidues',
   'ptmResidues', 'protonationStates',
   'ssbondResidues', 'phosphorylatableResidues',
   'phosphorylationStates', 'staplingPatches',
   'missingResidues', 'ssbonds',
   'glycosylations', 'hemes', 'staplings',
   'covalentLigands', 'nonStandards',
   'terminalCappings', 'acidsOptions',
   'surfaceProteinResidues', 'calcPka',
   'invalidCovalentLigands', 'selectedChains',
   'ffGeneration', 'ffGenAtomType'])

Using default settings

Initially, a PdbReaderProject contains no modifications, even those that would automatically be applied by MolCube Apps. To see the current settings applied to a PdbReaderProject, you can simply print the object. If you are in an interactive console, then entering just the object as a statement will accomplish the same thing:

>>> pdbreader.create_project(title='defaults-example', ff='charmmff', pdbId='2hac')
[ status update message, which should eventually say: ]
task finished...: 100%|███████████████████████████████| 5/5
>>> pdbreader
{'projectPk': '',
 'ph': 7.0,
 'chain': {},
 'glycosylation': []}

As you can see, before applying any settings, there isn’t much to show.

To apply the default modifications, you can use set_defaults():

>>> pdbreader.set_defaults()
>>> pdbreader
<PdbReaderProject with settings: {'chain': {'calcPka': False,
           'ffGeneration': None,
           'glycan': [],
           'glycosylation': [],
           'heme': [],
           'ion': [],
           'nucleicAcid': [],
           'ph': 7.0,
           'projectPk': 'your-project-id',
           'protein': [{'chainIndex': 'PROT_A',
                        'missing': [],
                        'selected': True,
                        'terminal': {'cter': 'CTER', 'nter': 'NTER'}},
                       {'chainIndex': 'PROT_B',
                        'missing': [],
                        'selected': True,
                        'terminal': {'cter': 'CTER', 'nter': 'NTER'}}],
           'ssbond': [{'residue1': {'chainIndex': 'PROT_A', 'resid': '2'},
                       'residue2': {'chainIndex': 'PROT_B', 'resid': '2'}}],
           'standaloneLigand': [],
           'water': []},
 'ffGeneration': None,
 'glycosylation': [],
 'heme': [],
 'ph': 7.0,
 'projectPk': 'your-project-id',
 'ssbond': [{'residue1': {'chainIndex': 'PROT_A', 'resid': '2'},
             'residue2': {'chainIndex': 'PROT_B', 'resid': '2'}}]}>

If, for some reason, you just want to see what settings would be applied, you can instead use the get_defaults() method.

Confirm chain selection (required)

After using set_defaults() and (optionally) a toggle function, your settings are still local to your machine. To tell the MolCube server to apply your chain selection, you must use the confirm_chains() method.

Args it accepts:

ph (float): pH to use (default: 7.0)
model (int): PDB model to use (default: 1st model)

This is equivalent to pressing “Submit” on the chain selection page of PDB Reader:

assert pdbreader.confirm_chains()

Model manipulation (required)

The sections below demonstrate model manipulation. Each of them is optional.

Use model_pdb() to confirm model manipulation. This must be performed after confirm_chains().

In the simplest case where you want to use default chain selection and default manipulations, then this is all you need to do:

import molcube as mc
from pprint import pprint

api_access_key = 'your-api-key'
molcube = mc.API('api.molcube.com', 443)
molcube.authenticate(api_token=api_access_key)

# simplest possible case: use defaults for everything
pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='test-defaults', ff='charmmff', pdbId='2hac')
pdbreader.set_defaults()

#
# if modifying chain selection, do so here
#

assert pdbreader.confirm_chains()

#
# if modifying manipulation options, do so here
#

assert pdbreader.model_pdb()

Supported manipulations

A quick summary of supported manipulations: The following sections explain basic usage of the manipulations supported by MolCube-API Client. Here is a summary of relevant methods with links to the API Reference:

Chain selection: toggle_chain(), toggle_chains_by_type()
Terminal patching: set_terminal_patch(), get_terminal_residues()
Mutations: add_mutation(), remove_mutation()
Phosphorylation: add_phosphorylation(), remove_phosphorylation()
Protonation: add_protonation(), remove_protonation()
Disulfide bonds: add_ssbond(), remove_ssbond()
Peptide stapling: add_staple(), remove_staple(), get_valid_staples()
Missing residue modeling: get_missing_residues(), add_missing_residues(), remove_missing_residues(), get_valid_missing_terminals()
Side chain orientation: orient_side_chains()
Glycosylation: See examples in Glycosylation.

Selecting different chains

The chain selection functions are toggle_chain() and toggle_chains_by_type(). The default chain selection when using set_defaults() is to select all chains except for water. You only need to call one of the these methods if deviating from this default.

Both functions take up to two arguments:

toggle_chain() args:
    enable (str | list[str]): chain ID or list of chain IDs to enable
    disable (str | list[str]): chain ID or list of chain IDs to disable

toggle_chains_by_type() args:
    enable (str | list[str]): a category or list or categories to enable
    disable (str | list[str]): same as above

Chain types:

protein
nucleicAcid
standaloneLigand
ion
water
glycan

See molcube.pdbreader.enums.CHAIN_TYPE():

>>> from molcube.pdbreader import enums
>>> print(f"Chain types: [{', '.join(enums.CHAIN_TYPE)}]")
Chain types: [protein, nucleicAcid,
   standaloneLigand, heme, ion, water, glycan]

As a more complex example, let’s consider 3PQR. This PDB has several chains, as shown below:

>>> import molcube as mc
>>> from pprint import pprint
>>> molcube = mc.API('api.molcube.com', 443)
>>> molcube.authenticate(api_token='your-api-key')
>>> pdbreader = molcube.create_pdb_reader_project()
>>> pdbreader.create_project(title='test-chains', ff='charmmff', pdbId='3pqr')
[ progress updates ]
>>> pdbreader.set_defaults()
Applied Chain Indices: GLYC_A, GLYC_B, GLYC_C, GLYC_D, GLYC_E
>>> pdbreader  # example output below is truncated for clarity
<PdbReaderProject with settings: {'chain': {'calcPka': False,
           'glycosylation': [ ... ],
           'ssbond': [ ... ],
           'heme': [],
           'ion': [],
           'nucleicAcid': [],
           'ph': 7.0,
           'projectPk': 'your-project-id',
           'glycan': [{'chainIndex': 'GLYC_A', 'selected': True},
                      {'chainIndex': 'GLYC_B', 'selected': True},
                      {'chainIndex': 'GLYC_C', 'selected': True},
                      {'chainIndex': 'GLYC_D', 'selected': True},
                      {'chainIndex': 'GLYC_E', 'selected': True}],
           'protein': [{'chainIndex': 'PROT_A',
                        'missing': [],
                        'selected': True,
                        'terminal': {'cter': 'CTER', 'nter': 'NTER'}},
                       {'chainIndex': 'PROT_B',
                        'missing': [],
                        'selected': True,
                        'terminal': {'cter': 'CTER', 'nter': 'NTER'}}],
           'standaloneLigand': [
               {'chainIndex': 'HETE_C', 'resname': 'SO4', 'selected': True},
               {'chainIndex': 'HETE_D', 'resname': 'ACT', 'selected': True},
               {'chainIndex': 'HETE_E', 'resname': 'ACT', 'selected': True}],
           'water': [{'chainIndex': 'WATE_A', 'selected': False},
                     {'chainIndex': 'WATE_B', 'selected': False}]},
 'ffGeneration': None, 'glycosylation': [ ... ], 'heme': [], 'ph': 7.0,
 'projectPk': 'your-project-id',
 'ssbond': [ ... ]}>

Chains can be enabled/disabled individually or by category:

# enable or disable a single chain
pdbreader.toggle_chain(enable='PROT_A')
pdbreader.toggle_chain(disable='GLYC_A')
# same as above
pdbreader.toggle_chain(enable='PROT_A', disable='GLYC_B')

# enable or disable multiple chains
pdbreader.toggle_chain(enable=['PROT_A', 'GLYC_C'],
                       disable='PROT_B')

# disable everything except protein
pdbreader.toggle_chains_by_type(enable='protein',
   disable=['glycan', 'water', 'ion', 'standaloneLigand'])

Terminal patching

Each protein chain returned by get_chains() has a 'terminal' key with a list of valid N-/C-terminals. The default terminal patch is the first one in the list. E.g., defaults for PROT_A below are CTER and NTER:

>>> from pprint import pprint
>>> chains = pdbreader.get_chains()
>>> pprint(chains['protein'])
{'chainId': 'A',
  'chainIndex': 'PROT_A',
  'nsdTerminal': {'cter': ['CT3'], 'nter': ['ACE']},
  'terminal': {'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE'],
               'nter': ['NTER', 'NNEU', 'ACE', 'NONE']}},
 {'chainId': 'B',
  'chainIndex': 'PROT_B',
  'nsdTerminal': {'cter': ['NONE'], 'nter': ['ACE']},
  'terminal': {'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE'],
               'nter': ['NTER', 'NNEU', 'ACE', 'NONE']}}]

>>> pprint({chain: pdbreader._option_by_chain[chain] for chain in ('PROT_A', 'PROT_B')})

{'PROT_A': {'chainIndex': 'PROT_A',
            'missing': [],
            'selected': True,
            'terminal': {'cter': 'CTER', 'nter': 'NTER'}},
 'PROT_B': {'chainIndex': 'PROT_B',
            'missing': [],
            'selected': True,
            'terminal': {'cter': 'CTER', 'nter': 'NTER'}}}

To set a different terminal patch, use set_terminal_patch().

Expected arguments:

chain_id (str, required): chain to set
nter (str, optional): use this patch, if given; else use default patch
cter (str, optional): use this patch, if given; else use default patch

Mutations

Point mutations are added with add_mutation() and removed with remove_mutation().

Expected arguments:

add_mutation()

chain_id (str): chain containing residue to mutate
resid (str): residue ID to mutate
new_resname (str): name of residue to mutate to

remove_mutation()

chain_id (str): chain containing mutated residue
resid (str): residue ID to restore

Example usage:

import molcube as mc

api_access_key = 'your-api-key'
molcube = mc.API('api.molcube.com', 443)
molcube.authenticate(api_token=api_access_key)

# simplest possible case: use defaults for everything
pdbreader = molcube.create_pdb_reader_project()
pdbreader.create_project(title='test-mutation', ff='charmmff', pdbId='2klu')

pdbreader.set_defaults()
assert pdbreader.confirm_chains()

pdbreader.add_mutation(chain_id='PROT_A', resid='364', new_resname='ASN')  # GLY 364 -> ASN
pdbreader.add_mutation(chain_id='PROT_A', resid='365', new_resname='ALA')  # PRO 365 -> ALA

assert pdbreader.model_pdb()

Phosphorylation and protonation

add_phosphorylation() and add_protonation() take the same arguments and differ only by what patch residues are considered valid.

Args:

chain_id (str): chain index, e.g. PROT_A
resid (str): resid to protonate
patch (str): name of phosphorylation/titration patch to apply

You can find valid options in the dict returned by get_pdb_info():

>>> from pprint import pprint
>>> pdb_info = pdbreader.get_pdb_info()
>>> print('protonations:')
>>> pprint( pdb_info['titrableResidues'] )
protonations:
{'ARG': ['RN1', 'RN2', 'RN3'],
 'ASP': ['ASPP'],
 'CYS': ['CYM'],
 'GLU': ['GLUP'],
 'HIE': ['HSP', 'HSD', 'HSE'],
 'HIP': ['HSP', 'HSD', 'HSE'],
 'HIS': ['HSP', 'HSD', 'HSE'],
 'HSD': ['HSP', 'HSE'],
 'HSE': ['HSP', 'HSD'],
 'HSP': ['HSD', 'HSE'],
 'LYS': ['LSN']}

>>> print('phosphorylations:')
>>> pprint( pdb_info['phosphorylatableResidues'] )
phosphorylations:
{'SER': ['SP2', 'SP1'], 'THR': ['THPB',
 'THP1'], 'TYR': ['TP2', 'TP1']}

Example usage:

import molcube as mc

api_access_key = 'your-api-key'
molcube = mc.API('api.molcube.com', 443)
molcube.authenticate(api_token=api_access_key)

pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='test-phosphorylation-2klu', ff='charmmff', pdbId='2klu')

pdbreader.set_defaults()
assert pdbreader.confirm_chains()

pdbreader.add_phosphorylation(chain_id='PROT_A', resid='394', patch='SP1')
pdbreader.add_phosphorylation(chain_id='PROT_A', resid='415', patch='SP1')
pdbreader.add_phosphorylation(chain_id='PROT_A', resid='431', patch='SP1')
assert pdbreader.model_pdb()

Disulfide bonds

Disulfide bonds present in the PDB are shown in pdb_info:

>>> pdb_info = pdbreader.get_pdb_info()
>>> pdb_info['ssbonds']
{'residue1': {'chainIndex': 'PROT_A', 'resid': '2'},
  'residue2': {'chainIndex': 'PROT_B', 'resid': '2'}}]

Disulfide bonds are added with add_ssbond() and removed with remove_ssbond().

Both require the same arguments:

residue1 (str|dict): first ssbond residue

residue2 (str|dict): second ssbond residue

The two acceptable formats are shown below:

# string format: "chain_id residue_id"
pdbreader_project.add_ssbond('PROT_A 50', 'PROT_A 62')

# if passing dict, it must be structed like below
pdbreader_project.add_ssbond(
    residue1={'chainIndex': 'PROT_A', 'resid': '50'},
    residue2={'chainIndex': 'PROT_A', 'resid': '62'})

Peptide Stapling

Staples are added with add_staple() and removed with remove_staple(). Usage is almost exactly like with disulfide bonds, except for an additional argument when adding a staple: the staple_type.

To see the valid staple types, use get_valid_staples():

>>> pdbreader.get_valid_staples()
['META3', 'META4', 'META5', 'META6', 'META7',
 'RMETA3', 'RMETA4', 'RMETA5', 'RMETA6',
 'RMETA7', 'DIBM', 'DIBP', 'CR12', 'CR21']

Example usage:

import molcube as mc

api_access_key = 'your-api-key'
molcube = mc.API('api.molcube.com', 443)
molcube.authenticate(api_token=api_access_key)

pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='test-staple', ff='charmmff', pdbId='1ubq')

pdbreader.set_defaults()
assert pdbreader.confirm_chains()

pdbreader.add_staple('RMETA3', 'PROT_A 1', 'PROT_A 3')
assert pdbreader.model_pdb()

Missing residues

RCSB entries frequently include proteins with unmodeled protein regions, e.g., because those regions could not be resolved crystallographically. These regions are indicated as ‘missing’ in the PDB. There are two ways of handling them:

Leave the section unmodeled. This is the default for missing residues located at either the C-terminus or N-terminus (we call these regions “terminal”).
Attempt to model the missing residues. This is the default for all other missing residues (we call these regions “non-terminal”).

Which approach you choose depends on how much you think the unmodeled region influences the dynamics that you are interested in. We explore each case below.

Case 1: Model the missing residues

Terminal residues: Modeling terminal residues usually means that the residue located at the N- or C-terminus will be a different amino acid type—with possibly different terminal patch options.

To see the valid terminal patches for missing residues, use get_valid_missing_terminals(). When called without arguments, it returns a dict of options for each missing residue range. Example (using 2ZFF):

>>> from pprint import pprint
>>> pdbreader.set_default()
>>> pprint(pdbreader.get_valid_missing_terminals())
{'PROT_A': {   '1H_1D': {'nter': ['NTER', 'NNEU', 'ACE', 'NONE']},
             '14L_14N': {'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE']}},
 'PROT_B': {'148_149E': {'nter': ['GLYP', 'NGNE', 'ACE', 'NONE'],
                         'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE']},
             '247_247': {'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE']}},
 'PROT_C': {   '64_64': {'cter': ['CTER', 'CNEU', 'CT1', 'CT2', 'CT3', 'NONE']}}}

The function can be called with chain_id to give only the results for a given chain, or both chain_id and residue_range to get just a single range’s options. The sub-dicts that contain only a nter or a cter key correspond to residue ranges located at the corresponding terminus, whereas those with both keys are non-terminal.

As shown above, the residue range consists of two residue IDs joined by _. These resids correspond to the start and end of the missing region. In the case of ranges with only a single missing residue, the start and end will be the same residue (e.g., 247_247 and 64_64 above).

To model a terminal residue range, (1) add it to the list of residues to model with add_missing_residues() and (2) give a valid option for set_terminal_patch. For example, to model the range 1H~1D on PROT_A with the NNEU patch and 14L_14N with CNEU:

>>> pdbreader.add_missing_residues(chain_id='PROT_A', residue_range='1H_1D')
>>> pdbreader.add_missing_residues(chain_id='PROT_A', residue_range='14L_14N')
>>> pdbreader.set_terminal_patch(chain_id='PROT_A', nter='NNEU', cter='CNEU')

As of this release (0.4.2), if both missing termini on the same protein chain will be modeled, then each terminus should be added with add_missing_residues() before using set_terminal_patch(), and both nter and cter args to set_terminal_patch() should be explicitly set. Currently, set_terminal_patch() otherwise uses the previous (non-missing) terminal residue’s default terminal patch.

Non-terminal residues: The procedure below is only necessary if you do not call set_defaults() or to reverse the effect of remove_missing_residues().

Modeling non-terminal residues is simpler than terminal residues, as there are no terminal patches to set. E.g., to model the non-terminal missing residues of 2ZFF, it’s just one command:

>>> pdbreader.add_missing_residues(chain_id='PROT_B', residue_range='148_149E')

Case 2: Don’t model the missing residues

Terminal residues: The procedure below is only necessary to reverse the effect of add_missing_residues().

To undo modeling a terminal missing residue range, use remove_missing_residues(). For example:

>>> pdbreader.remove_missing_residues(chain_id='PROT_A', residue_range='14L_14N')

If you previously used set_terminal_patch(), then refer to Terminal patching to find and set a valid terminal patch.

Non-terminal residues: Since non-terminal missing residues occur (by definition) somewhere in the middle of the amino acid sequence, choosing not to model them implies splitting the original chain into two disconnected chains—thus creating two more terminal sites. Unlike with terminal patches, the nter and cter keys returned by get_valid_missing_terminals() are only relevant for non-terminal patches when the given range is disabled.

To disable modeling of a non-terminal missing residue range, remove it from the list of residues to model and set the terminal patches for the new disconnected region with remove_missing_residues(). For example, to omit 2ZFF’s non-terminal missing residues and model the new disconnected terminals with GLYP (N-terminus) and CTER (C-terminus), do this:

>>> pdbreader.remove_missing_residues(chain_id='PROT_B',
...     residue_range='148_149E', nter='GLYP', cter='CTER')

Handling ligands and nonstandard residues

There are currently two supported ways to handle ligands using MolCube-API Client:

Use the default SDF structure that MolCube obtains from RCSB.
Upload a custom SDF structure

Although MolCube Apps also supports uploading custom force fields, this is not yet available in MolCube-API Client.

For nonstandard protein residues in particular, only the first option is supported. However, nonstandard amino acids can be substituted with standard amino acids, as shown in Substitute nonstandard amino acids.

Using default SDF

To demonstrate default ligand handling, we’ll use the RCSB ID 2ZFF. This structure contains a modified amino acid (TYS) and a ligand (53U).

Default SDF for ligand

Standalone (non-covalent) ligands are handled with validate_standalone_sdf(). Since the 53U ligand is available on RCSB, we can tell MolCube to use that. This is the default action when only the residue name is provided. For example:

pdbreader.validate_standalone_sdf(resname='53U')
pdbreader.confirm_chains()

Standalone ligands are handled in the chain selection step, so confirm_chains() should be called before handling nonstandard residues.

Default SDF for nonstandard residues

Engineered or modified amino acids are considered “nonstandard” residues. They need special handling due to the covalent connections to other amino acids. After doing confirm_chains(), we can parameterize TYS like so:

pdbreader.get_nonstandard_sdf()
pdbreader.model_pdb()

When called without arguments, get_nonstandard_sdf() uses the default FF generation method for all nonstandard residues and uses energy minimization to refine their conformation.

If your PDB contains multiple engineered residues that must each be handled differently, then you should specify each one separately using the corresponding resname. E.g.:

pdbreader.get_nonstandard_sdf(resname='ABC', ff_type='XFF')
pdbreader.get_nonstandard_sdf(resname='DEF', ff_type='fastGAFF')
pdbreader.model_pdb()

Full example: 2ZFF with defaults

import molcube as mc

api_access_key = 'your-api-key'
molcube = mc.API('api.molcube.com', 443)
molcube.authenticate(api_token=api_access_key)

#
# Initialize project by downloading structure from RCSB
#
pdbreader = molcube.create_pdb_reader_project()
assert pdbreader.create_project(title='ligand-example', ff='charmmff', pdbId='2zff')
pdbreader.set_defaults()

# handle ligand: 53U
assert pdbreader.validate_standalone_sdf(resname='53U')
assert pdbreader.confirm_chains()

# handle nonstandard residues: TYS
assert pdbreader.get_nonstandard_sdf()
assert pdbreader.model_pdb()

pdbreader.download_project('myproject.tgz')

Using custom SDF file

To use a custom SDF file (standalone ligands, only), use the sdf argument. If you already have the contents of the SDF stored in a string variable, you can pass that:

# option 1: sdf is contained in string
pdbreader.get_nonstandard_sdf(sdf='the SDF file contents...')

Otherwise, pass an open file handle:

# option 2: file handle
with open('path/to/my.sdf') as sdf_file:
   pdbreader.get_nonstandard_sdf(sdf=sdf_file)

Substitute nonstandard amino acids

There are some residues that MolCube cannot generate FFs for. In the web UI, you might see an explanation like, “This modified residue has modification made in backbone atoms. It cannot generate FF at this moment.” Instead, you can swap the residue with an amino acid or engineered residue that already exists in the FF. E.g., suppose you have a residue ASH with unsupported modifications and you want to replace it with aspartate (ASP). This can be done like so:

pdbreader.get_nonstandard_sdf(resname='ASH', replace='ASP')
pdbreader.model_pdb()

Side chain orientation

The orientation of amino acid side chains can be modified by rotating the dihedral angle of the bond between the side chain and protein backbone. To do this, you must first pass sideChainOrient=True to model_pdb():

>>>
# [...]
# ... setup pdbreader project ...
# [...]
#
# set PDB options
>>> pdbreader.set_defaults()
>>> assert pdbreader.confirm_chains()
>>> assert pdbreader.model_pdb(sideChainOrient=True)
>>>
# [...]
# ... do side chain settings ...
# [...]

To view the current orientation settings, You can view the current orientation settings with get_available_dihedrals():

>>> from pprint import pprint
>>> dihe = pdbreader.get_available_dihedrals()
>>> pp(dihe)
{'dihedrals': [['PROT_B', '-3', -62.721],
               ['PROT_B', '19', -57.573],
               ['PROT_B', '16', -137.886],
               ...
               ... output truncated for this example
               ...
               ['PROT_A', '25', -66.206],
               ['PROT_B', '17', -51.048],
               ['PROT_B', '3', -171.18]],
 'availableResidues': [['PROT_A', '30', 'ASP'],
                       ['PROT_A', '10', 'PHE'],
                       ['PROT_B', '23', 'VAL'],
                       ...
                       ... output truncated for this example
                       ...
                       ['PROT_B', '-2', 'SER'],
                       ['PROT_B', '27', 'ARG'],
                       ['PROT_B', '8', 'ILE']],
 'sideChainOrient': True}

The server’s response is a dict with the following keys:

dihedrals: a list of [chain_id, resid, angle]
- shows the current dihedral angles in degrees.
- if called as get_available_dihedrals(as_string=True), the response format is a list of strings: chain_id resid angle
availableResidues: a list of [chain_id, resid, resname]
- lists all protein residues whose side chains can be oriented
- if called as get_available_dihedrals(as_string=True), the response format is a list of strings: chain_id resid resname
sideChainOrient: always True

Use the orient_side_chains() method to set the new rotations. This function accepts a single argument, which can either be a string, a list of strings, or a list of dicts. For example, the following all do the same thing:

pdbreader.orient_side_chains(['PROT_A 1 30.0'])

# in dict format, 'angle' should be a float
pdbreader.orient_side_chains({'chainIndex': 'PROT_A', 'resid': '1', 'angle': 30.0})

# this only works because we are only changing one angle; for multiple angles,
# use one of the above formats
pdbreader.orient_side_chains('PROT_A 1 30.0')

Calling this method runs a task to set the new side chains. When it finishes, you can download your PDB/project files as usual:

pdbreader.download_pdb("oriented.pdb")