""" pdbfixer.py: Fixes problems in PDB files This is part of the OpenMM molecular simulation toolkit originating from Simbios, the NIH National Center for Physics-Based Simulation of Biological Structures at Stanford, funded under the NIH Roadmap for Medical Research, grant U54 GM072970. See https://simtk.org. Portions copyright (c) 2013-2023 Stanford University and the Authors. Authors: Peter Eastman Contributors: Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ from __future__ import absolute_import __author__ = "Peter Eastman" __version__ = "1.7" import openmm as mm import openmm.app as app import openmm.unit as unit from openmm.app.internal.pdbstructure import PdbStructure from openmm.app.internal.pdbx.reader.PdbxReader import PdbxReader from openmm.app.element import hydrogen, oxygen from openmm.app.forcefield import NonbondedGenerator # Support Cythonized functions in OpenMM 7.3 # and also implementations in older versions. try: from openmm.app.internal import compiled matchResidue = compiled.matchResidueToTemplate except ImportError: matchResidue = app.forcefield._matchResidue import numpy as np import numpy.linalg as lin import sys import os import os.path import math if sys.version_info >= (3,0): from urllib.request import urlopen from io import StringIO else: from urllib2 import urlopen from cStringIO import StringIO substitutions = { '2AS':'ASP', '3AH':'HIS', '5HP':'GLU', '5OW':'LYS', 'ACL':'ARG', 'AGM':'ARG', 'AIB':'ALA', 'ALM':'ALA', 'ALO':'THR', 'ALY':'LYS', 'ARM':'ARG', 'ASA':'ASP', 'ASB':'ASP', 'ASK':'ASP', 'ASL':'ASP', 'ASQ':'ASP', 'AYA':'ALA', 'BCS':'CYS', 'BHD':'ASP', 'BMT':'THR', 'BNN':'ALA', 'BUC':'CYS', 'BUG':'LEU', 'C5C':'CYS', 'C6C':'CYS', 'CAS':'CYS', 'CCS':'CYS', 'CEA':'CYS', 'CGU':'GLU', 'CHG':'ALA', 'CLE':'LEU', 'CME':'CYS', 'CSD':'ALA', 'CSO':'CYS', 'CSP':'CYS', 'CSS':'CYS', 'CSW':'CYS', 'CSX':'CYS', 'CXM':'MET', 'CY1':'CYS', 'CY3':'CYS', 'CYG':'CYS', 'CYM':'CYS', 'CYQ':'CYS', 'DAH':'PHE', 'DAL':'ALA', 'DAR':'ARG', 'DAS':'ASP', 'DCY':'CYS', 'DGL':'GLU', 'DGN':'GLN', 'DHA':'ALA', 'DHI':'HIS', 'DIL':'ILE', 'DIV':'VAL', 'DLE':'LEU', 'DLY':'LYS', 'DNP':'ALA', 'DPN':'PHE', 'DPR':'PRO', 'DSN':'SER', 'DSP':'ASP', 'DTH':'THR', 'DTR':'TRP', 'DTY':'TYR', 'DVA':'VAL', 'EFC':'CYS', 'FLA':'ALA', 'FME':'MET', 'GGL':'GLU', 'GL3':'GLY', 'GLZ':'GLY', 'GMA':'GLU', 'GSC':'GLY', 'HAC':'ALA', 'HAR':'ARG', 'HIC':'HIS', 'HIP':'HIS', 'HMR':'ARG', 'HPQ':'PHE', 'HTR':'TRP', 'HYP':'PRO', 'IAS':'ASP', 'IIL':'ILE', 'IYR':'TYR', 'KCX':'LYS', 'LLP':'LYS', 'LLY':'LYS', 'LTR':'TRP', 'LYM':'LYS', 'LYZ':'LYS', 'MAA':'ALA', 'MEN':'ASN', 'MHS':'HIS', 'MIS':'SER', 'MK8':'LEU', 'MLE':'LEU', 'MPQ':'GLY', 'MSA':'GLY', 'MSE':'MET', 'MVA':'VAL', 'NEM':'HIS', 'NEP':'HIS', 'NLE':'LEU', 'NLN':'LEU', 'NLP':'LEU', 'NMC':'GLY', 'OAS':'SER', 'OCS':'CYS', 'OMT':'MET', 'PAQ':'TYR', 'PCA':'GLU', 'PEC':'CYS', 'PHI':'PHE', 'PHL':'PHE', 'PR3':'CYS', 'PRR':'ALA', 'PTR':'TYR', 'PYX':'CYS', 'SAC':'SER', 'SAR':'GLY', 'SCH':'CYS', 'SCS':'CYS', 'SCY':'CYS', 'SEL':'SER', 'SEP':'SER', 'SET':'SER', 'SHC':'CYS', 'SHR':'LYS', 'SMC':'CYS', 'SOC':'CYS', 'STY':'TYR', 'SVA':'SER', 'TIH':'ALA', 'TPL':'TRP', 'TPO':'THR', 'TPQ':'ALA', 'TRG':'LYS', 'TRO':'TRP', 'TYB':'TYR', 'TYI':'TYR', 'TYQ':'TYR', 'TYS':'TYR', 'TYY':'TYR' } proteinResidues = ['ALA', 'ASN', 'CYS', 'GLU', 'HIS', 'LEU', 'MET', 'PRO', 'THR', 'TYR', 'ARG', 'ASP', 'GLN', 'GLY', 'ILE', 'LYS', 'PHE', 'SER', 'TRP', 'VAL'] rnaResidues = ['A', 'G', 'C', 'U', 'I'] dnaResidues = ['DA', 'DG', 'DC', 'DT', 'DI'] class Sequence(object): """Sequence holds the sequence of a chain, as specified by SEQRES records.""" def __init__(self, chainId, residues): self.chainId = chainId self.residues = residues class ModifiedResidue(object): """ModifiedResidue holds information about a modified residue, as specified by a MODRES record.""" def __init__(self, chainId, number, residueName, standardName): self.chainId = chainId self.number = number self.residueName = residueName self.standardName = standardName def _guessFileFormat(file, filename): """Guess whether a file is PDB or PDBx/mmCIF based on its filename and contents.""" filename = filename.lower() if '.pdbx' in filename or '.cif' in filename: return 'pdbx' if '.pdb' in filename: return 'pdb' for line in file: if line.startswith('data_') or line.startswith('loop_'): file.seek(0) return 'pdbx' if line.startswith('HEADER') or line.startswith('REMARK') or line.startswith('TITLE '): file.seek(0) return 'pdb' # It's certainly not a valid PDBx/mmCIF. Guess that it's a PDB. file.seek(0) return 'pdb' def _overlayPoints(points1, points2): """Given two sets of points, determine the translation and rotation that matches them as closely as possible. Parameters ---------- points1 (numpy array of openmm.unit.Quantity with units compatible with distance) - reference set of coordinates points2 (numpy array of openmm.unit.Quantity with units compatible with distance) - set of coordinates to be rotated Returns ------- translate2 - vector to translate points2 by in order to center it rotate - rotation matrix to apply to centered points2 to map it on to points1 center1 - center of points1 Notes ----- This is based on W. Kabsch, Acta Cryst., A34, pp. 828-829 (1978). """ if len(points1) == 0: return (mm.Vec3(0, 0, 0), np.identity(3), mm.Vec3(0, 0, 0)) if len(points1) == 1: return (points1[0], np.identity(3), -1*points2[0]) # Compute centroids. center1 = unit.sum(points1)/float(len(points1)) center2 = unit.sum(points2)/float(len(points2)) # Compute R matrix. R = np.zeros((3, 3)) for p1, p2 in zip(points1, points2): x = p1-center1 y = p2-center2 for i in range(3): for j in range(3): R[i][j] += y[i]*x[j] # Use an SVD to compute the rotation matrix. (u, s, v) = lin.svd(R) return (-1*center2, np.dot(u, v).transpose(), center1) def _findUnoccupiedDirection(point, positions): """Given a point in space and a list of atom positions, find the direction in which the local density of atoms is lowest.""" point = point.value_in_unit(unit.nanometers) direction = mm.Vec3(0, 0, 0) for pos in positions.value_in_unit(unit.nanometers): delta = pos-point distance = unit.norm(delta) if distance > 0.1: distance2 = distance*distance direction -= delta/(distance2*distance2) direction /= unit.norm(direction) return direction class PDBFixer(object): """PDBFixer implements many tools for fixing problems in PDB and PDBx/mmCIF files. """ def __init__(self, filename=None, pdbfile=None, pdbxfile=None, url=None, pdbid=None): """Create a new PDBFixer instance to fix problems in a PDB or PDBx/mmCIF file. Parameters ---------- filename : str, optional, default=None The name of the file to read. The format is determined automatically based on the filename extension, or if that is ambiguous, by looking at the file content. pdbfile : file, optional, default=None A file-like object from which the PDB file is to be read. The file is not closed after reading. pdbxfile : file, optional, default=None A file-like object from which the PDBx/mmCIF file is to be read. The file is not closed after reading. url : str, optional, default=None A URL specifying the internet location from which the file contents should be retrieved. The format is determined automatically by looking for a filename extension in the URL, or if that is ambiguous, by looking at the file content. pdbid : str, optional, default=None A four-letter PDB code specifying the structure to be retrieved from the RCSB. Notes ----- Only one of structure, filename, pdbfile, pdbxfile, url, or pdbid may be specified or an exception will be thrown. Examples -------- Start from a filename. >>> filename = 'tests/data/test.pdb' >>> fixer = PDBFixer(filename=filename) Start from a file object. >>> with open(filename) as f: ... fixer = PDBFixer(pdbfile=f) Start from a URL. >>> fixer = PDBFixer(url='http://www.rcsb.org/pdb/files/1VII.pdb') Start from a PDB code. >>> fixer = PDBFixer(pdbid='1VII') """ # Check to make sure only one option has been specified. if bool(filename) + bool(pdbfile) + bool(pdbxfile) + bool(url) + bool(pdbid) != 1: raise Exception("Exactly one option [filename, pdbfile, pdbxfile, url, pdbid] must be specified.") self.source = None if pdbid: # A PDB id has been specified. url = 'http://www.rcsb.org/pdb/files/%s.pdb' % pdbid if filename: # A local file has been specified. self.source = filename file = open(filename, 'r') if _guessFileFormat(file, filename) == 'pdbx': self._initializeFromPDBx(file) else: self._initializeFromPDB(file) file.close() elif pdbfile: # A file-like object has been specified. self._initializeFromPDB(pdbfile) elif pdbxfile: # A file-like object has been specified. self._initializeFromPDBx(pdbxfile) elif url: # A URL has been specified. self.source = url file = urlopen(url) contents = file.read().decode('utf-8') file.close() file = StringIO(contents) if _guessFileFormat(file, url) == 'pdbx': self._initializeFromPDBx(contents) else: self._initializeFromPDB(StringIO(contents)) # Check the structure has some atoms in it. atoms = list(self.topology.atoms()) if len(atoms) == 0: raise Exception("Structure contains no atoms.") # Load the templates. self.templates = {} templatesPath = os.path.join(os.path.dirname(__file__), 'templates') for file in os.listdir(templatesPath): templatePdb = app.PDBFile(os.path.join(templatesPath, file)) name = next(templatePdb.topology.residues()).name self.templates[name] = templatePdb def _initializeFromPDB(self, file): """Initialize this object by reading a PDB file.""" structure = PdbStructure(file) pdb = app.PDBFile(structure) self.topology = pdb.topology self.positions = pdb.positions self.sequences = [Sequence(s.chain_id, s.residues) for s in structure.sequences] self.modifiedResidues = [ModifiedResidue(r.chain_id, r.number, r.residue_name, r.standard_name) for r in structure.modified_residues] def _initializeFromPDBx(self, file): """Initialize this object by reading a PDBx/mmCIF file.""" pdbx = app.PDBxFile(file) self.topology = pdbx.topology self.positions = pdbx.positions # PDBxFile doesn't record the information about sequence or modified residues, so we need to read them separately. file.seek(0) reader = PdbxReader(file) data = [] reader.read(data) block = data[0] # Load the sequence data. sequenceData = block.getObj('entity_poly_seq') sequences = {} if sequenceData is not None: entityIdCol = sequenceData.getAttributeIndex('entity_id') residueCol = sequenceData.getAttributeIndex('mon_id') for row in sequenceData.getRowList(): entityId = row[entityIdCol] residue = row[residueCol] if entityId not in sequences: sequences[entityId] = [] sequences[entityId].append(residue) # Sequences are stored by "entity". There could be multiple chains that are all the same entity, so we need to # convert from entities to chains. asymData = block.getObj('struct_asym') self.sequences = [] if asymData is not None: asymIdCol = asymData.getAttributeIndex('id') entityIdCol = asymData.getAttributeIndex('entity_id') for row in asymData.getRowList(): asymId = row[asymIdCol] entityId = row[entityIdCol] if entityId in sequences: self.sequences.append(Sequence(asymId, sequences[entityId])) # Load the modified residues. modData = block.getObj('pdbx_struct_mod_residue') self.modifiedResidues = [] if modData is not None: asymIdCol = modData.getAttributeIndex('label_asym_id') resNameCol = modData.getAttributeIndex('label_comp_id') resNumCol = modData.getAttributeIndex('auth_seq_id') standardResCol = modData.getAttributeIndex('parent_comp_id') if -1 not in (asymIdCol, resNameCol, resNumCol, standardResCol): for row in modData.getRowList(): self.modifiedResidues.append(ModifiedResidue(row[asymIdCol], int(row[resNumCol]), row[resNameCol], row[standardResCol])) def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules): """Create a new Topology in which missing atoms have been added. Parameters ---------- heavyAtomsOnly : bool If True, only heavy atoms will be added to the topology. omitUnknownMolecules : bool If True, unknown molecules will be omitted from the topology. Returns ------- newTopology : openmm.app.Topology A new Topology object containing atoms from the old. newPositions : list of openmm.unit.Quantity with units compatible with nanometers Atom positions for the new Topology object. newAtoms : openmm.app.Topology.Atom New atom objects. existingAtomMap : dict Mapping from old atoms to new atoms. """ newTopology = app.Topology() newPositions = []*unit.nanometer newAtoms = [] existingAtomMap = {} addedAtomMap = {} addedOXT = [] residueCenters = [self._computeResidueCenter(res).value_in_unit(unit.nanometers) for res in self.topology.residues()]*unit.nanometers for chain in self.topology.chains(): if omitUnknownMolecules and not any(residue.name in self.templates for residue in chain.residues()): continue chainResidues = list(chain.residues()) newChain = newTopology.addChain(chain.id) for indexInChain, residue in enumerate(chain.residues()): # Insert missing residues here. if (chain.index, indexInChain) in self.missingResidues: insertHere = self.missingResidues[(chain.index, indexInChain)] endPosition = self._computeResidueCenter(residue) if indexInChain > 0: startPosition = self._computeResidueCenter(chainResidues[indexInChain-1]) loopDirection = _findUnoccupiedDirection((startPosition+endPosition)/2, residueCenters) else: outward = _findUnoccupiedDirection(endPosition, residueCenters)*unit.nanometers norm = unit.norm(outward) if norm > 0*unit.nanometer: outward *= len(insertHere)*0.5*unit.nanometer/norm startPosition = endPosition+outward loopDirection = None firstIndex = int(residue.id)-len(insertHere) self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, loopDirection, residue, newAtoms, newPositions, firstIndex) # Create the new residue and add existing heavy atoms. newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode) for atom in residue.atoms(): if not heavyAtomsOnly or (atom.element is not None and atom.element != hydrogen): if atom.name == 'OXT' and (chain.index, indexInChain+1) in self.missingResidues: continue # Remove terminal oxygen, since we'll add more residues after this one newAtom = newTopology.addAtom(atom.name, atom.element, newResidue) existingAtomMap[atom] = newAtom newPositions.append(self.positions[atom.index]) if residue in self.missingAtoms: # Find corresponding atoms in the residue and the template. template = self.templates[residue.name] atomPositions = dict((atom.name, self.positions[atom.index]) for atom in residue.atoms()) points1 = [] points2 = [] for atom in template.topology.atoms(): if atom.name in atomPositions: points1.append(atomPositions[atom.name].value_in_unit(unit.nanometer)) points2.append(template.positions[atom.index].value_in_unit(unit.nanometer)) # Compute the optimal transform to overlay them. (translate2, rotate, translate1) = _overlayPoints(points1, points2) # Add the missing atoms. addedAtomMap[residue] = {} for atom in self.missingAtoms[residue]: newAtom = newTopology.addAtom(atom.name, atom.element, newResidue) newAtoms.append(newAtom) addedAtomMap[residue][atom] = newAtom templatePosition = template.positions[atom.index].value_in_unit(unit.nanometer) newPositions.append((mm.Vec3(*np.dot(rotate, templatePosition+translate2))+translate1)*unit.nanometer) if residue in self.missingTerminals: terminalsToAdd = self.missingTerminals[residue] else: terminalsToAdd = None # If this is the end of the chain, add any missing residues that come after it. if residue == chainResidues[-1] and (chain.index, indexInChain+1) in self.missingResidues: insertHere = self.missingResidues[(chain.index, indexInChain+1)] if len(insertHere) > 0: startPosition = self._computeResidueCenter(residue) outward = _findUnoccupiedDirection(startPosition, residueCenters)*unit.nanometers norm = unit.norm(outward) if norm > 0*unit.nanometer: outward *= len(insertHere)*0.5*unit.nanometer/norm endPosition = startPosition+outward firstIndex = int(residue.id)+1 self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, None, residue, newAtoms, newPositions, firstIndex) newResidue = list(newChain.residues())[-1] if newResidue.name in proteinResidues: terminalsToAdd = ['OXT'] else: terminalsToAdd = None # If a terminal OXT is missing, add it. if terminalsToAdd is not None: atomPositions = dict((atom.name, newPositions[atom.index].value_in_unit(unit.nanometer)) for atom in newResidue.atoms()) if 'OXT' in terminalsToAdd: newAtom = newTopology.addAtom('OXT', oxygen, newResidue) newAtoms.append(newAtom) addedOXT.append(newAtom) d_ca_o = atomPositions['O']-atomPositions['CA'] d_ca_c = atomPositions['C']-atomPositions['CA'] d_ca_c /= unit.sqrt(unit.dot(d_ca_c, d_ca_c)) v = d_ca_o - d_ca_c*unit.dot(d_ca_c, d_ca_o) newPositions.append((atomPositions['O']+2*v)*unit.nanometer) newTopology.setUnitCellDimensions(self.topology.getUnitCellDimensions()) newTopology.createStandardBonds() newTopology.createDisulfideBonds(newPositions) # Add the bonds between atoms in heterogens. for a1,a2 in self.topology.bonds(): if a1 in existingAtomMap and a2 in existingAtomMap and (a1.residue.name not in app.Topology._standardBonds or a2.residue.name not in app.Topology._standardBonds): newTopology.addBond(existingAtomMap[a1], existingAtomMap[a2]) # Return the results. return (newTopology, newPositions, newAtoms, existingAtomMap) def _computeResidueCenter(self, residue): """Compute the centroid of a residue.""" return unit.sum([self.positions[atom.index] for atom in residue.atoms()])/len(list(residue.atoms())) def _addMissingResiduesToChain(self, chain, residueNames, startPosition, endPosition, loopDirection, orientTo, newAtoms, newPositions, firstIndex): """Add a series of residues to a chain.""" orientToPositions = dict((atom.name, self.positions[atom.index]) for atom in orientTo.atoms()) if loopDirection is None: loopDirection = mm.Vec3(0, 0, 0) # We'll add the residues in an arc connecting the endpoints. Figure out the height of that arc. length = unit.norm(endPosition-startPosition) numResidues = len(residueNames) if length > numResidues*0.3*unit.nanometers: loopHeight = 0*unit.nanometers else: loopHeight = (numResidues*0.3*unit.nanometers-length)/2 # Add the residues. for i, residueName in enumerate(residueNames): template = self.templates[residueName] # Find a translation that best matches the adjacent residue. points1 = [] points2 = [] for atom in template.topology.atoms(): if atom.name in orientToPositions: points1.append(orientToPositions[atom.name].value_in_unit(unit.nanometer)) points2.append(template.positions[atom.index].value_in_unit(unit.nanometer)) (translate2, rotate, translate1) = _overlayPoints(points1, points2) # Create the new residue. newResidue = chain.topology.addResidue(residueName, chain, "%d" % ((firstIndex+i)%10000)) fraction = (i+1.0)/(numResidues+1.0) translate = startPosition + (endPosition-startPosition)*fraction + loopHeight*math.sin(fraction*math.pi)*loopDirection templateAtoms = list(template.topology.atoms()) if newResidue == next(chain.residues()): templateAtoms = [atom for atom in templateAtoms if atom.name not in ('P', 'OP1', 'OP2')] for atom in templateAtoms: newAtom = chain.topology.addAtom(atom.name, atom.element, newResidue) newAtoms.append(newAtom) templatePosition = template.positions[atom.index].value_in_unit(unit.nanometer) newPositions.append(mm.Vec3(*np.dot(rotate, templatePosition))*unit.nanometer+translate) def removeChains(self, chainIndices=None, chainIds=None): """Remove a set of chains from the structure. Parameters ---------- chainIndices : list of int, optional, default=None List of indices of chains to remove. chainIds : list of str, optional, default=None List of chain ids of chains to remove. Examples -------- Load a PDB file with two chains and eliminate the second chain. >>> fixer = PDBFixer(pdbid='4J7F') >>> fixer.removeChains(chainIndices=[1]) Load a PDB file with two chains and eliminate chains named 'B' and 'D'. >>> fixer = PDBFixer(pdbid='4J7F') >>> fixer.removeChains(chainIds=['B','D']) """ modeller = app.Modeller(self.topology, self.positions) allChains = list(self.topology.chains()) if chainIndices == None: chainIndices = list() if chainIds != None: # Add all chains that match the selection to the list. for (chainNumber, chain) in enumerate(allChains): if chain.id in chainIds: chainIndices.append(chainNumber) # Ensure only unique entries remain. chainIndices = list(set(chainIndices)) # Do nothing if no chains will be deleted. if len(chainIndices) == 0: return modeller.delete(allChains[i] for i in chainIndices) self.topology = modeller.topology self.positions = modeller.positions return def findMissingResidues(self): """Find residues that are missing from the structure. The results are stored into the missingResidues field, which is a dict. Each key is a tuple consisting of the index of a chain, and the residue index within that chain at which new residues should be inserted. The corresponding value is a list of the names of residues to insert there. Examples -------- >>> fixer = PDBFixer(pdbid='1VII') >>> fixer.findMissingResidues() >>> missing_residues = fixer.missingResidues """ chains = [c for c in self.topology.chains() if len(list(c.residues())) > 0] chainWithGaps = {} # Find the sequence of each chain, with gaps for missing residues. for chain in chains: residues = list(chain.residues()) ids = [int(r.id) for r in residues] for i, res in enumerate(residues): if res.insertionCode not in ('', ' '): for j in range(i, len(residues)): ids[j] += 1 minResidue = min(ids) maxResidue = max(ids) chainWithGaps[chain] = [None]*(maxResidue-minResidue+1) for r, id in zip(residues, ids): chainWithGaps[chain][id-minResidue] = r.name # Try to find the chain that matches each sequence. chainSequence = {} chainOffset = {} for sequence in self.sequences: for chain in chains: if chain.id != sequence.chainId: continue if chain in chainSequence: continue for offset in range(len(sequence.residues)-len(chainWithGaps[chain])+1): if all(a == b or b == None for a,b in zip(sequence.residues[offset:], chainWithGaps[chain])): chainSequence[chain] = sequence chainOffset[chain] = offset break if chain in chainSequence: break # Now build the list of residues to add. self.missingResidues = {} for chain in self.topology.chains(): if chain in chainSequence: offset = chainOffset[chain] sequence = chainSequence[chain].residues gappedSequence = chainWithGaps[chain] index = 0 for i in range(len(sequence)): if i < offset or i >= len(gappedSequence)+offset or gappedSequence[i-offset] is None: key = (chain.index, index) if key not in self.missingResidues: self.missingResidues[key] = [] residueName = sequence[i] if residueName in substitutions: residueName = substitutions[sequence[i]] self.missingResidues[key].append(residueName) else: index += 1 def findNonstandardResidues(self): """Identify non-standard residues found in the structure, and select standard residues to replace them with. The results are stored into the nonstandardResidues field, which is a map of Residue objects to the names of suggested replacement residues. Examples -------- Find nonstandard residues. >>> fixer = PDBFixer(pdbid='1YRI') >>> fixer.findNonstandardResidues() >>> nonstandard_residues = fixer.nonstandardResidues """ # First find residues based on our table of standard substitutions. nonstandard = dict((r, substitutions[r.name]) for r in self.topology.residues() if r.name in substitutions) # Now add ones based on MODRES records. modres = dict(((m.chainId, str(m.number), m.residueName), m.standardName) for m in self.modifiedResidues) for chain in self.topology.chains(): for residue in chain.residues(): key = (chain.id, residue.id, residue.name) if key in modres: replacement = modres[key] if replacement == 'DU': replacement = 'DT' if replacement in self.templates: nonstandard[residue] = replacement self.nonstandardResidues = [(r, nonstandard[r]) for r in sorted(nonstandard, key=lambda r: r.index)] def replaceNonstandardResidues(self): """Replace every residue listed in the nonstandardResidues field with the specified standard residue. Notes ----- You must have first called findNonstandardResidues() to identify nonstandard residues. Examples -------- Find and replace nonstandard residues using replacement templates stored in the 'templates' field of PDBFixer object. >>> fixer = PDBFixer(pdbid='1YRI') >>> fixer.findNonstandardResidues() >>> fixer.replaceNonstandardResidues() """ if len(self.nonstandardResidues) > 0: deleteAtoms = [] # Find atoms that should be deleted. for residue, replaceWith in self.nonstandardResidues: residue.name = replaceWith template = self.templates[replaceWith] standardAtoms = set(atom.name for atom in template.topology.atoms()) for atom in residue.atoms(): if atom.element in (None, hydrogen) or atom.name not in standardAtoms: deleteAtoms.append(atom) # Delete them. modeller = app.Modeller(self.topology, self.positions) modeller.delete(deleteAtoms) self.topology = modeller.topology self.positions = modeller.positions def applyMutations(self, mutations, chain_id): """Apply a list of amino acid substitutions to make a mutant protein. Parameters ---------- mutations : list of strings Each string must include the resName (original), index, and resName (target). For example, ALA-133-GLY will mutate alanine 133 to glycine. chain_id : str String based chain ID of the single chain you wish to mutate. Notes ----- We require three letter codes to avoid possible ambiguitities. We can't guarantee that the resulting model is a good one; for significant changes in sequence, you should probably be using a standalone homology modelling tool. Examples -------- Find nonstandard residues. >>> fixer = PDBFixer(pdbid='1VII') >>> fixer.applyMutations(["ALA-57-GLY"], "A") >>> fixer.findMissingResidues() >>> fixer.findMissingAtoms() >>> fixer.addMissingAtoms() >>> fixer.addMissingHydrogens(7.0) """ # Retrieve all residues that match the specified chain_id. # NOTE: Multiple chains may have the same chainid, but must have unique resSeq entries. resSeq_to_residue = dict() # resSeq_to_residue[resid] is the residue in the requested chain corresponding to residue identifier 'resid' for chain in self.topology.chains(): if chain.id == chain_id: for residue in chain.residues(): resSeq_to_residue[int(residue.id)] = residue # Make a map of residues to mutate based on requested mutation list. residue_map = dict() # residue_map[residue] is the name of the new residue to mutate to, if a mutation is desired for mut_str in mutations: old_name, resSeq, new_name = mut_str.split("-") resSeq = int(resSeq) if resSeq not in resSeq_to_residue: raise(KeyError("Cannot find chain %s residue %d in system!" % (chain_id, resSeq))) residue = resSeq_to_residue[resSeq] # retrieve the requested residue if residue.name != old_name: raise(ValueError("You asked to mutate chain %s residue %d name %s, but that residue is actually %s!" % (chain_id, resSeq, old_name, residue.name))) try: template = self.templates[new_name] except KeyError: raise(KeyError("Cannot find residue %s in template library!" % new_name)) # Store mutation residue_map[residue] = new_name # If there are mutations to be made, make them. if len(residue_map) > 0: deleteAtoms = [] # list of atoms to delete # Find atoms that should be deleted. for residue in residue_map.keys(): replaceWith = residue_map[residue] residue.name = replaceWith template = self.templates[replaceWith] standardAtoms = set(atom.name for atom in template.topology.atoms()) for atom in residue.atoms(): if atom.element in (None, hydrogen) or atom.name not in standardAtoms: deleteAtoms.append(atom) # Delete atoms queued to be deleted. modeller = app.Modeller(self.topology, self.positions) modeller.delete(deleteAtoms) self.topology = modeller.topology self.positions = modeller.positions def findMissingAtoms(self): """Find heavy atoms that are missing from the structure. The results are stored into two fields: missingAtoms and missingTerminals. Each of these is a dict whose keys are Residue objects and whose values are lists of atom names. missingAtoms contains standard atoms that should be present in any residue of that type. missingTerminals contains terminal atoms that should be present at the start or end of a chain. Notes ----- You must have first called findMissingResidues(). Examples -------- Find missing heavy atoms in Abl kinase structure. >>> fixer = PDBFixer(pdbid='2F4J') >>> fixer.findMissingResidues() >>> fixer.findMissingAtoms() >>> # Retrieve missing atoms. >>> missingAtoms = fixer.missingAtoms >>> # Retrieve missing terminal atoms. >>> missingTerminals = fixer.missingTerminals """ missingAtoms = {} missingTerminals = {} # Loop over residues. for chain in self.topology.chains(): chainResidues = list(chain.residues()) for residue in chain.residues(): if residue.name in self.templates: template = self.templates[residue.name] atomNames = set(atom.name for atom in residue.atoms()) templateAtoms = list(template.topology.atoms()) if residue == chainResidues[0] and (chain.index, 0) not in self.missingResidues: templateAtoms = [atom for atom in templateAtoms if atom.name not in ('P', 'OP1', 'OP2')] # Add atoms from the template that are missing. missing = [] for atom in templateAtoms: if atom.name not in atomNames: missing.append(atom) if len(missing) > 0: missingAtoms[residue] = missing # Add missing terminal atoms. terminals = [] if residue == chainResidues[-1] and (chain.index, len(chainResidues)) not in self.missingResidues: templateNames = set(atom.name for atom in template.topology.atoms()) if 'OXT' not in atomNames and all(name in templateNames for name in ['C', 'O', 'CA']): terminals.append('OXT') if len(terminals) > 0: missingTerminals[residue] = terminals self.missingAtoms = missingAtoms self.missingTerminals = missingTerminals def addMissingAtoms(self, seed=None): """Add all missing heavy atoms, as specified by the missingAtoms, missingTerminals, and missingResidues fields. Parameters ---------- seed : int Integer to set the random seed number of the integrator used in the minimization of the coordinates of the newly-added atoms. Notes ----- You must already have called findMissingAtoms() to have identified atoms to be added. Examples -------- Find missing heavy atoms in Abl kinase structure. >>> fixer = PDBFixer(pdbid='2F4J') >>> fixer.findMissingResidues() >>> fixer.findMissingAtoms() >>> fixer.addMissingAtoms() """ # Create a Topology that 1) adds missing atoms, 2) removes all hydrogens, and 3) removes unknown molecules. (newTopology, newPositions, newAtoms, existingAtomMap) = self._addAtomsToTopology(True, True) if len(newAtoms) == 0: # No atoms were added, but new bonds might have been created. newBonds = set(newTopology.bonds()) for atom1, atom2 in self.topology.bonds(): if atom1 in existingAtomMap and atom2 in existingAtomMap: a1 = existingAtomMap[atom1] a2 = existingAtomMap[atom2] if (a1, a2) in newBonds: newBonds.remove((a1, a2)) elif (a2, a1) in newBonds: newBonds.remove((a2, a1)) # Add the new bonds to the original Topology. inverseAtomMap = dict((y,x) for (x,y) in existingAtomMap.items()) for atom1, atom2 in newBonds: self.topology.addBond(inverseAtomMap[atom1], inverseAtomMap[atom2]) else: # Create a System for energy minimizing it. forcefield = self._createForceField(newTopology, False) system = forcefield.createSystem(newTopology) # Set any previously existing atoms to be massless, they so won't move. for atom in existingAtomMap.values(): system.setParticleMass(atom.index, 0.0) # If any heavy atoms were omitted, add them back to avoid steric clashes. nonbonded = [f for f in system.getForces() if isinstance(f, mm.CustomNonbondedForce)][0] for atom in self.topology.atoms(): if atom.element not in (None, hydrogen) and atom not in existingAtomMap: system.addParticle(0.0) nonbonded.addParticle([]) newPositions.append(self.positions[atom.index]) # For efficiency, only compute interactions that involve a new atom. nonbonded.addInteractionGroup([atom.index for atom in newAtoms], range(system.getNumParticles())) # Do an energy minimization. integrator = mm.LangevinIntegrator(300*unit.kelvin, 10/unit.picosecond, 5*unit.femtosecond) if seed is not None: integrator.setRandomNumberSeed(seed) context = mm.Context(system, integrator) context.setPositions(newPositions) mm.LocalEnergyMinimizer.minimize(context) state = context.getState(getPositions=True) if newTopology.getNumResidues() > 1: # When looking for pairs of atoms that are too close to each other, exclude pairs that # are in the same residue or are directly bonded to each other. exclusions = dict((atom, {a.index for a in atom.residue.atoms()}) for atom in newAtoms) for a1, a2 in newTopology.bonds(): if a1 in exclusions: exclusions[a1].add(a2.index) if a2 in exclusions: exclusions[a2].add(a1.index) cutoff = 0.13 nearest = self._findNearestDistance(context, newAtoms, cutoff, exclusions) if nearest < cutoff: # Some atoms are very close together. Run some dynamics while slowly increasing the strength of the # repulsive interaction to try to improve the result. for i in range(10): context.setParameter('C', 0.15*(i+1)) integrator.step(200) d = self._findNearestDistance(context, newAtoms, cutoff, exclusions) if d > nearest: nearest = d state = context.getState(getPositions=True) if nearest >= cutoff: break context.setState(state) context.setParameter('C', 1.0) mm.LocalEnergyMinimizer.minimize(context) state = context.getState(getPositions=True) # Now create a new Topology, including all atoms from the original one and adding the missing atoms. (newTopology2, newPositions2, newAtoms2, existingAtomMap2) = self._addAtomsToTopology(False, False) # Copy over the minimized positions for the new atoms. for a1, a2 in zip(newAtoms, newAtoms2): newPositions2[a2.index] = state.getPositions()[a1.index] self.topology = newTopology2 self.positions = newPositions2 def removeHeterogens(self, keepWater=True): """Remove all heterogens from the structure. Parameters ---------- keepWater : bool, optional, default=True If True, water molecules will not be removed. Examples -------- Remove heterogens in Abl structure complexed with imatinib. >>> fixer = PDBFixer(pdbid='2F4J') >>> fixer.removeHeterogens(keepWater=False) """ keep = set(proteinResidues).union(dnaResidues).union(rnaResidues) keep.add('N') keep.add('UNK') if keepWater: keep.add('HOH') toDelete = [] for residue in self.topology.residues(): if residue.name not in keep: toDelete.append(residue) modeller = app.Modeller(self.topology, self.positions) modeller.delete(toDelete) self.topology = modeller.topology self.positions = modeller.positions def addMissingHydrogens(self, pH=7.0, forcefield=None): """Add missing hydrogen atoms to the structure. Parameters ---------- pH : float, optional, default=7.0 The pH based on which to select hydrogens. forcefield : ForceField, optional, default=None The forcefield used when adding and minimizing hydrogens. If None, a default forcefield is used. Notes ----- No extensive electrostatic analysis is performed; only default residue pKas are used. Examples -------- Examples -------- Add missing hydrogens appropriate for pH 8. >>> fixer = PDBFixer(pdbid='1VII') >>> fixer.addMissingHydrogens(pH=8.0) """ modeller = app.Modeller(self.topology, self.positions) modeller.addHydrogens(pH=pH, forcefield=forcefield) self.topology = modeller.topology self.positions = modeller.positions def addSolvent(self, boxSize=None, padding=None, boxVectors=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*unit.molar, boxShape='cube'): """Add a solvent box surrounding the structure. Parameters ---------- boxSize : openmm.Vec3, optional, default=None The size of the box to fill with water. If specified, padding and boxVectors must not be specified. padding : openmm.unit.Quantity compatible with nanometers, optional, default=None Padding around macromolecule for filling box with water. If specified, boxSize and boxVectors must not be specified. boxVectors : 3-tuple of openmm.Vec3, optional, default=None Three vectors specifying the geometry of the box. If specified, padding and boxSize must not be specified. positiveIon : str, optional, default='Na+' The type of positive ion to add. Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+'. negativeIon : str, optional, default='Cl-' The type of negative ion to add. Allowed values are 'Cl-', 'Br-', 'F-', and 'I-'. ionicStrength : openmm.unit.Quantity with units compatible with molar, optional, default=0*molar The total concentration of ions (both positive and negative) to add. This does not include ions that are added to neutralize the system. boxShape: str='cube' the box shape to use. Allowed values are 'cube', 'dodecahedron', and 'octahedron'. If padding is None, this is ignored. Examples -------- Add missing residues, heavy atoms, and hydrogens, and then solvate with 10 A padding. >>> fixer = PDBFixer(pdbid='1VII') >>> fixer.findMissingResidues() >>> fixer.findMissingAtoms() >>> fixer.addMissingAtoms() >>> fixer.addMissingHydrogens(pH=8.0) >>> fixer.addSolvent(padding=10*unit.angstrom, ionicStrength=0.050*unit.molar) """ modeller = app.Modeller(self.topology, self.positions) forcefield = self._createForceField(self.topology, True) modeller.addSolvent(forcefield, padding=padding, boxSize=boxSize, boxVectors=boxVectors, boxShape=boxShape, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength) chains = list(modeller.topology.chains()) if len(chains) == 1: chains[0].id = 'A' else: chains[-1].id = chr(ord(chains[-2].id)+1) self.topology = modeller.topology self.positions = modeller.positions def addMembrane(self, lipidType='POPC', membraneCenterZ=0*unit.nanometer, minimumPadding=1*unit.nanometer, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*unit.molar): """Add a lipid membrane to the structure. This method adds both lipids and water, so you should call either addSolvent() or addMembrane(), but not both. See Modeller.addMembrane() for more details. Parameters ---------- lipidType : string='POPC' the type of lipid to use. Supported values are 'POPC', 'POPE', 'DLPC', 'DLPE', 'DMPC', 'DOPC', and 'DPPC'. membraneCenterZ: distance=0*nanometer the position along the Z axis of the center of the membrane minimumPadding : distance=1*nanometer the padding distance to use positiveIon : str, optional, default='Na+' The type of positive ion to add. Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+'. negativeIon : str, optional, default='Cl-' The type of negative ion to add. Allowed values are 'Cl-', 'Br-', 'F-', and 'I-'. ionicStrength : openmm.unit.Quantity with units compatible with molar, optional, default=0*molar The total concentration of ions (both positive and negative) to add. This does not include ions that are added to neutralize the system. """ modeller = app.Modeller(self.topology, self.positions) forcefield = self._createForceField(self.topology, True) modeller.addMembrane(forcefield, lipidType=lipidType, minimumPadding=minimumPadding, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength) chains = list(modeller.topology.chains()) if len(chains) == 1: chains[0].id = 'A' else: chains[-1].id = chr(ord(chains[-2].id)+1) self.topology = modeller.topology self.positions = modeller.positions def _createForceField(self, newTopology, water): """Create a force field to use for optimizing the positions of newly added atoms.""" if water: forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3p.xml') nonbonded = [f for f in forcefield._forces if isinstance(f, NonbondedGenerator)][0] radii = {'H':0.198, 'Li':0.203, 'C':0.340, 'N':0.325, 'O':0.299, 'F':0.312, 'Na':0.333, 'Mg':0.141, 'P':0.374, 'S':0.356, 'Cl':0.347, 'K':0.474, 'Br':0.396, 'Rb':0.527, 'I':0.419, 'Cs':0.605} else: forcefield = app.ForceField(os.path.join(os.path.dirname(__file__), 'soft.xml')) # The Topology may contain residues for which the ForceField does not have a template. # If so, we need to create new templates for them. atomTypes = {} bondedToAtom = [] for atom in newTopology.atoms(): bondedToAtom.append(set()) for atom1, atom2 in newTopology.bonds(): bondedToAtom[atom1.index].add(atom2.index) bondedToAtom[atom2.index].add(atom1.index) for residue in newTopology.residues(): # Make sure the ForceField has a template for this residue. signature = app.forcefield._createResidueSignature([atom.element for atom in residue.atoms()]) if signature in forcefield._templateSignatures: if any(matchResidue(residue, t, bondedToAtom) is not None for t in forcefield._templateSignatures[signature]): continue # Create a new template. resName = "extra_"+residue.name template = app.ForceField._TemplateData(resName) forcefield._templates[resName] = template indexInResidue = {} for atom in residue.atoms(): element = atom.element typeName = 'extra_'+element.symbol if element not in atomTypes: atomTypes[element] = app.ForceField._AtomType(typeName, '', 0.0, element) forcefield._atomTypes[typeName] = atomTypes[element] if water: # Select a reasonable vdW radius for this atom type. if element.symbol in radii: sigma = radii[element.symbol] else: sigma = 0.5 nonbonded.registerAtom({'type':typeName, 'charge':'0', 'sigma':str(sigma), 'epsilon':'0'}) indexInResidue[atom.index] = len(template.atoms) template.atoms.append(app.ForceField._TemplateAtomData(atom.name, typeName, element)) for atom in residue.atoms(): for bondedTo in bondedToAtom[atom.index]: if bondedTo in indexInResidue: b = (indexInResidue[atom.index], indexInResidue[bondedTo]) if b[0] < b[1]: template.bonds.append(b) template.atoms[b[0]].bondedTo.append(b[1]) template.atoms[b[1]].bondedTo.append(b[0]) else: b = indexInResidue[atom.index] template.externalBonds.append(b) template.atoms[b].externalBonds += 1 if signature in forcefield._templateSignatures: forcefield._templateSignatures[signature].append(template) else: forcefield._templateSignatures[signature] = [template] return forcefield def _findNearestDistance(self, context, newAtoms, cutoff, exclusions): """Given a set of newly added atoms, find the closest distance between one of those atoms and another atom.""" positions = context.getState(getPositions=True).getPositions(asNumpy=True).value_in_unit(unit.nanometer) boxSize = np.max(positions, axis=0)-np.min(positions, axis=0) boxVectors = [(boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2])] cells = app.modeller._CellList(positions, cutoff, boxVectors, False) nearest_squared = sys.float_info.max for atom in newAtoms: excluded = exclusions[atom] for i in cells.neighbors(positions[atom.index]): if i not in excluded: p = positions[atom.index]-positions[i] dist_squared = np.dot(p, p) if dist_squared < nearest_squared: nearest_squared = dist_squared return np.sqrt(nearest_squared) def main(): if len(sys.argv) < 2: # Display the UI. from . import ui ui.launchUI() else: # Run in command line mode. from optparse import OptionParser parser = OptionParser(usage="Usage: %prog\n %prog filename [options] \n\nWhen run with no arguments, it launches the user interface. If any arguments are specified, it runs in command line mode.") parser.add_option('--pdbid', default=None, dest='pdbid', metavar='PDBID', help='PDB id to retrieve from RCSB [default: None]') parser.add_option('--url', default=None, dest='url', metavar='URL', help='URL to retrieve PDB from [default: None]') parser.add_option('--output', default='output.pdb', dest='output', metavar='FILENAME', help='output pdb file [default: output.pdb]') parser.add_option('--add-atoms', default='all', dest='atoms', choices=('all', 'heavy', 'hydrogen', 'none'), help='which missing atoms to add: all, heavy, hydrogen, or none [default: all]') parser.add_option('--keep-heterogens', default='all', dest='heterogens', choices=('all', 'water', 'none'), metavar='OPTION', help='which heterogens to keep: all, water, or none [default: all]') parser.add_option('--replace-nonstandard', action='store_true', default=False, dest='nonstandard', help='replace nonstandard residues with standard equivalents') parser.add_option('--add-residues', action='store_true', default=False, dest='residues', help='add missing residues') parser.add_option('--water-box', dest='box', type='float', nargs=3, metavar='X Y Z', help='add a water box. The value is the box dimensions in nm [example: --water-box=2.5 2.4 3.0]') parser.add_option('--ph', type='float', default=7.0, dest='ph', help='the pH to use for adding missing hydrogens [default: 7.0]') parser.add_option('--positive-ion', default='Na+', dest='positiveIon', choices=('Cs+', 'K+', 'Li+', 'Na+', 'Rb+'), metavar='ION', help='positive ion to include in the water box: Cs+, K+, Li+, Na+, or Rb+ [default: Na+]') parser.add_option('--negative-ion', default='Cl-', dest='negativeIon', choices=('Cl-', 'Br-', 'F-', 'I-'), metavar='ION', help='negative ion to include in the water box: Cl-, Br-, F-, or I- [default: Cl-]') parser.add_option('--ionic-strength', type='float', default=0.0, dest='ionic', metavar='STRENGTH', help='molar concentration of ions to add to the water box [default: 0.0]') parser.add_option('--verbose', default=False, action='store_true', dest='verbose', metavar='VERBOSE', help='Print verbose output') (options, args) = parser.parse_args() if (len(args) == 0) and (options.pdbid==None) and (options.url==None): parser.error('No filename specified') if len(args) > 1: parser.error('Must specify a single filename or --pdbid or --url') if options.pdbid != None: if options.verbose: print('Retrieving PDB "' + options.pdbid + '" from RCSB...') fixer = PDBFixer(pdbid=options.pdbid) elif options.url != None: if options.verbose: print('Retrieving PDB from URL "' + options.url + '"...') fixer = PDBFixer(url=options.url) else: fixer = PDBFixer(filename=sys.argv[1]) if options.residues: if options.verbose: print('Finding missing residues...') fixer.findMissingResidues() else: fixer.missingResidues = {} if options.nonstandard: if options.verbose: print('Finding nonstandard residues...') fixer.findNonstandardResidues() if options.verbose: print('Replacing nonstandard residues...') fixer.replaceNonstandardResidues() if options.heterogens == 'none': fixer.removeHeterogens(False) elif options.heterogens == 'water': fixer.removeHeterogens(True) if options.verbose: print('Finding missing atoms...') fixer.findMissingAtoms() if options.atoms not in ('all', 'heavy'): fixer.missingAtoms = {} fixer.missingTerminals = {} if options.verbose: print('Adding missing atoms...') fixer.addMissingAtoms() if options.atoms in ('all', 'hydrogen'): if options.verbose: print('Adding missing hydrogens...') fixer.addMissingHydrogens(options.ph) if options.box is not None: if options.verbose: print('Adding solvent...') fixer.addSolvent(boxSize=options.box*unit.nanometer, positiveIon=options.positiveIon, negativeIon=options.negativeIon, ionicStrength=options.ionic*unit.molar) with open(options.output, 'w') as f: if options.verbose: print('Writing output...') if fixer.source is not None: f.write("REMARK 1 PDBFIXER FROM: %s\n" % fixer.source) app.PDBFile.writeFile(fixer.topology, fixer.positions, f, True) if options.verbose: print('Done.') if __name__ == '__main__': main()