from copy import * import math import random maxatoms = 29997 class CharmmPDB: def __init__(self, filename, aareplace=[], chargeexeptions=[]): """ Reads and parses a pdb file """ self.ssbonds = [] self.atoms = [] pdbfile = open(filename, "r") lines = pdbfile.readlines() pdbfile.close() self._parsepdb(lines) def _parsepdb(self, lines): for line in lines: if line[0:4] == "ATOM" or line[0:6] == "HETATM": self.atoms.append(Atom(line)) elif line[0:6] == "SSBOND": self.ssbonds.append(SSbond(line)); def write(self, filename): """ Write the structure to file. Since charmm has an arbitrary limit of the length of a structure, we cut it in multiple files below this limit. This is controlled by the global variable maxatoms. """ blocks = len(self.atoms) / maxatoms for i in range(0, blocks): subset = self.atoms[i*maxatoms:(i+1)*maxatoms] self._renumber_residues(subset) self._write(filename + str(i), subset) rest = self.atoms[blocks*maxatoms:len(self.atoms)] self._renumber_residues(rest) if len(rest) > 0 and blocks > 0: self._write(filename + str(blocks), rest) elif len(rest) > 0: self._write(filename, rest) def _write(self, filename, atoms): """ Write a file of atom records in a pdb-like format to file """ print "Writing file " + filename pdbout = open(filename, "w") for atom in atoms: pdbout.write(atom.writeatom() +"\n") pdbout.write("TER\n") pdbout.close() def atomfilter(self, predicate): """ Delete all atoms where the predicate return false """ self.atoms = filter(predicate, self.atoms) atommap = {} for atom in self.atoms: atommap[atom.resnum.strip()] = True def copy(self): """ Make a deep copy of this structure """ new_self = copy(self) new_self.atoms = deepcopy(self.atoms) new_self.ssbonds = deepcopy(self.ssbonds) return new_self def renumber_residues(self): return self._renumber_residues(self.atoms) def _renumber_residues(self,atoms): """ Renumber all residues and atoms to start from 1 """ oldnewmap = {} newresnum = 0 atomnum = 1 if len(atoms) == 0: return oldresnum = 'xxxxxxx' for atom in atoms: if atom.resnum != oldresnum: newresnum = newresnum + 1 oldresnum = atom.resnum oldnewmap[oldresnum.strip()] = str(newresnum) atom.resnum=str(newresnum) atom.atomnum=atomnum atomnum = atomnum + 1 atom.renumbered = True return oldnewmap class Atom: renumbered = False def __init__(self, line): if line[0:6] == "HETATM": self.hetatm = True else: self.hetatm = False if line[6] == "*": self.atomnum=99999 else: self.atomnum = int(line[6:11]) self.rawatomname = line[12:16] self.atomname = self.rawatomname.strip() self.multiocc = line[16].strip() self.restype = line[17:21].strip() self.chain = line[21].strip() self.resnum = line[22:27] self.x = float(line[30:38]) self.y = float(line[38:46]) self.z = float(line[46:54]) self.occ = float(line[56:60]) self.b = float(line[60:66]) self.segment = line[72:76].strip() if self.segment == "": self.segment = "SEGM" spaces = " "*80 def _wf(self, field, length, aln="r"): if(len(field) > length): print "Warning field truncated: " + str(field) if aln == 'r': return (self.spaces + field)[-length:] else: return (field + self.spaces)[:length] def writeatom(self): spaces = " "*80 line = self._wf("ATOM", 6, aln='l') line += self._wf("%d" % self.atomnum, 5) line += " " line += self._wf(self.atomname, 4) line += self._wf(self.multiocc, 1) line += self._wf(self.restype, 4, aln='l') line += self._wf(self.chain,1) line += self._wf(str(self.resnum), 5) line += " "*3 line += self._wf("%8.3F" % self.x, 8) line += self._wf("%8.3F" % self.y, 8) line += self._wf("%8.3F" % self.z, 8) line += " "*2 line += self._wf("%4.2F" % self.occ, 4) line += self._wf("%6.2F" % self.b, 6) line += " "*7 line += self._wf(self.segment, 4, aln='l') return line class SSbond: def __init__(self, line): res1 = line[11:23].strip() res2 = line[23:].strip() self.bondnum = line[6:10].strip() self.res1 = self._parseres(res1) self.res2 = self._parseres(res2) def _parseres(self, res): reslist = filter(lambda(x): x != '' , res.split(" ")) return reslist def lipid_replace(structure, lipids, aasubst, lipid_type): for lipid in lipids: for atom in structure.atoms: append=True #print "atom" + str(atom.resnum) + " " + atom.atomname if atom.resnum == lipid: atom.restype=lipid_type for subst in aasubst: if atom.atomname == subst[0]: #print "delete atom" + str(atom.resnum) + " " + atom.atomname atom.hetatm=True structure.atomfilter(is_peptide) for atom in structure.atoms: if atom.hetatm == True: print "atom" + str(atom.resnum) + " " + atom.atomname return structure headgroup_fixes = [["C15" ],["N"],[ "H15A" ],[ "H15B" ],[ "H15C" ],[ "C14" ], [ "H14A" ],[ "H14B" ],[ "H14C" ],[ "C13" ],[ "H13A" ],[ "H13B" ],[ "H13C" ], [ "C12" ],[ "H12A" ],[ "H12B" ],[ "C11" ],[ "H11A" ],[ "H11B"]] #Functions for replacing residues to the ones charmm expects def res_replace(structure, aasubst, chain=""): atoms = structure.atoms for subst in aasubst: if subst[3] == chain or chain == "": for atom in atoms: if atom.restype == subst[0] and atom.resnum.strip() == str(subst[1]): atom.restype = subst[2] def fix_charges(structure,aasubst): atoms = structure.atoms charge = 0 for atom in atoms: for subst in aasubst: if atom.restype == subst[0] and atom.atomname == subst[1]: charge+=1 return charge charge_fixes = [["ASP","HD2","1"], ["GLU","HE2","1"]] def fix_residues(structure, fixes): for fix in fixes: for atom in structure.atoms: if atom.restype == fix[0] and (atom.atomname == fix[1] or fix[1] == "*"): atom.restype = fix[2] if fix[3] != "*": atom.atomname = fix[3] atom.rawatomname = (" " + fix[3])[-5:] atom_fixes = [["ILE", "CD1","ILE", "CD"], ["HOH", "*", "TIP3","OH2"], ["NA" , "*", "SOD", "SOD"], ["CL" , "*", "CLA", "CLA"], ["CA" , "*", "CAL", "CAL"], ["K" , "*", "POT", "POT"], ["HIS", "*", "HSD", "*"]] #Finding the charge balance residues_charges = {"ASP" : -1, "GLU" : -1, "ARG" : 1, "LYS" : 1, "HSP" : 1, "DMPS": -1, "DMPG": -1 } ion_charges = {"NA" : 1, "SOD" : 1, "MG" : 2, "CL" : -1, "CLA" : -1, "K" : 1, "POT" : 1, "CA" : 2, "CAL" : 2, "FE" : 2, "ZN" : 2} def charge_balance(structure, chargehash): totalcharge = 0 reshash = {} segments = {} for atom in structure.atoms: segments[atom.segment] = atom.segment print segments for segname in segments.keys(): segfilter = make_segment_filter(segname) segment = structure.copy() segment.atomfilter(segfilter) segment.renumber_residues() totalcharge+=fix_charges(segment,charge_fixes) for atom in segment.atoms: reshash[atom.resnum.strip()] = atom.restype for resnum in reshash.keys(): restype = reshash.get(resnum) if chargehash.has_key(restype): totalcharge += chargehash[restype] return totalcharge #Find all the chains, and return each of them as a single peptide def make_chain_filter(chain): def chain_filter(atom): return atom.chain == chain return chain_filter def split_peptides(peptides): chainnames = {} chains = [] for atom in peptides.atoms: chainnames[atom.chain] = atom.chain for chainname in chainnames.keys(): chain_filter = make_chain_filter(chainname) chain = peptides.copy() chain.atomfilter(chain_filter) chains.append(chain) return chains #Functions for splitting in peptide, ions, crystalwater def is_water(atom): return atom.restype == "HOH" or atom.restype == "TIP3" def is_ion(atom): return ion_charges.has_key(atom.restype) def is_cofactor(atom): rt = atom.restype return atom.hetatm and not (rt=="HOH" or rt=="TIP3" or ion_charges.has_key(rt)) def is_peptide(atom): return not atom.hetatm def prepare_build(pdbpath, aasubst=[], chargediffs=[], custom_filters=[]): """ Prepare a PDB structure for charmm: Some atoms/residues has different names in charmm and pdb, change to the charmm name This is the case for ions and waters HIS does not exists, should be replaced by HSD, HDE or HSP depending on hydrogen bonds and protonation. The aasubst list is for this Find charged amino acids Find ion charges and compute total charge balance The CD1 atom of ILE residues is called CD in charmm -> rename Suggest disulfide patch sentences Suggest ion patch sentences split after peptide, cofactor, ions and crystal water """ pdb = CharmmPDB(pdbpath) allchains = [] peptidechains = {} for filter in custom_filters: pdb.atomfilter(filter) #write the peptides - one chain in one file. peptides = pdb.copy() peptides.atomfilter(is_peptide) peptide_list = split_peptides(peptides) for peptide in peptide_list: chainname = peptide.atoms[0].chain pepname = "peptide_" + peptide.atoms[0].chain.lower() + ".pdb" allchains.append([peptide, chainname, pepname]) peptidechains[chainname] = chainname for atom in peptide.atoms: atom.segment = "PEP" + chainname #filter out the ions ions = pdb.copy() ions.atomfilter(is_ion) if len(ions.atoms) > 0: chainname = ions.atoms[0].chain allchains.append([ions, chainname, "ions.pdb"]) #filter out the cofactor cofactors = pdb.copy() cofactors.atomfilter(is_cofactor) if len(cofactors.atoms) > 0: chainname = cofactors.atoms[0].chain allchains.append([cofactors, chainname, "cofactors.pdb"]) #filter out the crystal waters crystwater = pdb.copy() crystwater.atomfilter(is_water) chainname = crystwater.atoms[0].chain allchains.append([crystwater, chainname, "crystalwater.pdb"]) #Renumber the residues and print out all chains renumbered = {} for chainent in allchains: chain = chainent[0] chainname = chainent[1] res_replace(chain, aasubst, chainname) fix_residues(chain, atom_fixes) if chain.atoms[0].segment[0:3] == "PEP": renum = chain.renumber_residues() renumbered[chainname] = renum chain.write(chainent[2]) print "Add the following disulfide patches to the charmm build script:\n" for ssbond in pdb.ssbonds: chain1 = ssbond.res1[1] chain2 = ssbond.res2[1] oldresnum1 = ssbond.res1[2] oldresnum2 = ssbond.res2[2] renum1 = renumbered[chain1] renum2 = renumbered[chain2] str1 = " PEP" + chain1 + " " + renum1[oldresnum1] str2 = " PEP" + chain2 + " " + renum2[oldresnum2] print "patch disu" + str1 + str2 #Functions for solvation preparations def is_water_oxygen(atom): return is_solvbox_water(atom) and atom.atomname == "OH2" def is_solvbox_water(atom): return atom.segment[0:2] == "WT" def is_crystal_water(atom): return atom.segment[0:3] == "SOL" def is_protein2(atom): return atom.segment[0:3] == "PRO" def is_protpept(atom): return atom.segment[0:3] == "PRO" or atom.segment[0:3] == "PEP" def is_peptide2(atom): return atom.segment[0:3] == "PEP" def is_membrane(atom): return atom.segment[0:3]=="MEM" def is_memb_resid(atom, resid): return atom.segment[0:3]=="MEM" and atom.resnum == resid def make_segment_filter(segname): def is_segment(atom): return atom.segment == segname return is_segment #Suggestion of counter ions def atom_distance(atom1, atom2): xdist = atom2.x - atom1.x ydist = atom2.y - atom1.y zdist = atom2.z - atom1.z return math.sqrt(xdist*xdist + ydist*ydist + zdist*zdist) def suggest_negative_lipids(membrane,peptides, lipid_type, total_number,mindist=5.0): #pick a random lipid # check if upper or lower membrane leaflet # check if far enough away from protein nolipid=total_number/2 lipids= _suggest_lipids(membrane, peptides, nolipid, False,mindist) lipids += _suggest_lipids(membrane, peptides, nolipid, True,mindist) return lipid_replace(membrane, lipids, headgroup_fixes, lipid_type) def _suggest_lipids(membrane,peptides, num_lipids, above,mindist=5.0): nolipids=num_lipids replaced_lipids = [] while nolipids>0: new=True selected_atom = random.choice(membrane.atoms) selected_resid = selected_atom.resnum print "random resid " + str(selected_resid) #if lipid already in replacing for replaced in replaced_lipids: print "replaced " + replaced if replaced == selected_resid: new=False print "already chosen" if new==True: too_near=False position=True for latom in membrane.atoms: if latom.resnum==selected_resid: if latom.z < 0 and latom.atomname=="P": position=False print "below Z" #print "look for atoms in peptide" + str(len(peptides.atoms)) for atom in peptides.atoms: #print str(atom.atomname) if atom_distance(atom, latom) < mindist: too_near = True #print "Atom too close: " + str(atom.atomnum) else: too_near=True if too_near == False and above!=position: replaced_lipids.append(selected_resid) nolipids = nolipids - 1 print "nolipids " + str(nolipids) return replaced_lipids def suggest_counter_ions(peptide, solvboxwaters, total_charge, mindist=5.0): #pick a random water #check if it is too near the peptide or prev selected counterions #if not add to counterions #continue till neutral abscharge = abs(total_charge) replaced_waters = [] while abscharge > 0: selected_atom = random.choice(solvboxwaters) too_near = False for atom in peptide: if atom_distance(atom, selected_atom) < mindist: too_near = True for atom in replaced_waters: if atom_distance(atom, selected_atom) < mindist: too_near = True if too_near == False: replaced_waters.append(selected_atom) abscharge = abscharge - 1 else: print "Atom too close: " + str(atom.atomnum) resnums = [] for atom in replaced_waters: resnums.append(atom.resnum) return resnums def prepare_mixedmemb(pdbpath, lipid_type, num_lipids): pdb = CharmmPDB(pdbpath) peptides = pdb.copy() peptides.atomfilter(is_protpept) #handle membrane memb = pdb.copy() memb.atomfilter(is_membrane) newmemb = suggest_negative_lipids(memb, peptides, lipid_type, num_lipids) newmemb.write("mixed_memb.pdb") peptides=pdb.copy() peptides.atomfilter(is_peptide2) proteins=pdb.copy() proteins.atomfilter(is_protein2) crystalwater=pdb.copy() crystalwater.atomfilter(is_water) crystalwater.renumber_residues() ions=pdb.copy() ions.atomfilter(is_ion) if len(ions.atoms) > 0: ions.renumber_residues() ions.write("ions2.pdb") crystalwater.write("crystalwater_allh2.pdb") #filter out the cofactor cofactors = pdb.copy() cofactors.atomfilter(is_cofactor) if len(cofactors.atoms) > 0: chainname = cofactors.atoms[0].chain allchains.append([cofactors, chainname, "cofactors.pdb"]) proteins = pdb.copy() proteins.atomfilter(is_protein2) protsegments = {} for atom in proteins.atoms: protsegments[atom.segment] = atom.segment print protsegments for segname in protsegments.keys(): segfilter = make_segment_filter(segname) segment = proteins.copy() segment.atomfilter(segfilter) segment.renumber_residues() segment.write("protein_" + segname[3].lower() + "_allh2.pdb") peptides = pdb.copy() peptides.atomfilter(is_peptide2) pepsegments = {} for atom in peptides.atoms: pepsegments[atom.segment] = atom.segment print pepsegments for segname in pepsegments.keys(): segfilter = make_segment_filter(segname) segment = peptides.copy() segment.atomfilter(segfilter) segment.renumber_residues() segment.write("peptide_" + segname[3].lower() + "_allh2.pdb") def prepare_solvation(pdbpath, chargediff=[]): """ Find the ion balance - to find how many counter ions to place, and of which type (Na or Cl) Then find candidate solvent waters to be replaced by counter ions not too near each others Split in peptide, cofactor, ions, crystal water and solvent water if more than 30000 solvent atoms, split in several files. """ pdb = CharmmPDB(pdbpath) solvent = pdb.copy() solvent.atomfilter(is_solvbox_water) solvent.renumber_residues() solvboxox = solvent.copy() solvboxox.atomfilter(is_water_oxygen) print "Number of solvbox oxygens: %d" % len(solvboxox.atoms) proteins = pdb.copy() proteins.atomfilter(is_protein2) protsegments = {} for atom in proteins.atoms: protsegments[atom.segment] = atom.segment print protsegments for segname in protsegments.keys(): segfilter = make_segment_filter(segname) segment = proteins.copy() segment.atomfilter(segfilter) segment.renumber_residues() segment.write("protein_" + segname[3].lower() + "_allh2.pdb") peptides = pdb.copy() peptides.atomfilter(is_peptide2) pepsegments = {} for atom in peptides.atoms: pepsegments[atom.segment] = atom.segment print pepsegments for segname in pepsegments.keys(): segfilter = make_segment_filter(segname) segment = peptides.copy() segment.atomfilter(segfilter) segment.renumber_residues() segment.write("peptide_" + segname[3].lower() + "_allh2.pdb") ions = pdb.copy() ions.atomfilter(is_ion) memb = pdb.copy() memb.atomfilter(is_membrane) memb.renumber_residues() memb.write("membrane_allh2.pdb") #find charge balance protein_charge = charge_balance(proteins, residues_charges) peptide_charge = charge_balance(peptides, residues_charges) ion_charge = charge_balance(ions, ion_charges) membrane_charge=charge_balance(memb, residues_charges) print "protein charge is %d" %protein_charge print "peptide charge is %d" % peptide_charge print "ION charge is %d" % ion_charge print "MEMBRANE charge is %d" % membrane_charge charge_bal = protein_charge + peptide_charge + ion_charge + membrane_charge print "total charge is %d" % charge_bal if abs(charge_bal) > 0: suggestions = suggest_counter_ions(peptides.atoms, solvboxox.atoms, charge_bal) print "****** Suggested counter ion patches *******\n" waterpatch = "define watrepl select ( " count = 0 for wat in suggestions[1:]: count+=1 if count > 50: count=0 id = int(int(wat) *3. / maxatoms ) + 1 wat = int(wat) % int(maxatoms / 3) waterpatch += "(segid wat%s .and. " % id waterpatch += " resi %s) " % wat waterpatch += ") end \n\n" waterpatch += "define watrepl select watrepl .or. ( " else: id = int(int(wat) *3. / maxatoms ) + 1 wat = int(wat) % int(maxatoms / 3) waterpatch += "(segid wat%s .and. " % id waterpatch += " resi %s) .or. - \n " % wat id = int(int(suggestions[0]) *3. / maxatoms ) + 1 wat = int(suggestions[0]) % int(maxatoms / 3) waterpatch += "(segid wat%s .and. " % id waterpatch += " resi %s) " % wat waterpatch += ") end \n\n" #print waterpatch if charge_bal < 0: waterpatch += "set ion SOD \n" else: waterpatch += "set ion CLA \n" waterpatch += "stream ../lib/add-ions.str\n\n" createions_out = open("createions.str", "w") createions_out.write(waterpatch) createions_out.close() else: print "Charge balance 0, no counter ions needed\n" print waterbox_dimensions(solvboxox.atoms) crystalwater = pdb.copy() crystalwater.atomfilter(is_crystal_water) crystalwater.renumber_residues() if len(ions.atoms) > 0: ions.renumber_residues() ions.write("ions2.pdb") crystalwater.write("crystalwater_allh2.pdb") solvent.write("solvboxwater_allh2.pdb") #The dimension of the water box def waterbox_dimensions(solvboxwaters): minx = maxx = miny = maxy = minz = maxz = 0.0 for atom in solvboxwaters: if atom.x < minx: minx = atom.x if atom.y < miny: miny = atom.y if atom.z < minz: minz = atom.z if atom.x > maxx: maxx = atom.x if atom.y > maxy: maxy = atom.y if atom.z > maxz: maxz = atom.z return [minx, maxx, miny, maxy, minz, maxz] def pme_suggestions(proteinpath): pdb = CharmmPDB(proteinpath) dim = waterbox_dimensions(pdb.atoms) xdim = abs(dim[1] - dim[0]) ydim = abs(dim[3] - dim[2]) zdim = abs(dim[5] - dim[4]) print "Please suggest to insert the following lines into the NAMD setup file:\n" print "cellbasisvector1 %5.1F 0 0 " % xdim print "cellbasisvector2 0 %5.1F 0 " % ydim print "cellbasisvector3 0 0 %5.1F " % zdim print "" print "PMEGridSizeX %d" % int(round(xdim)) print "PMEGridSizeY %d" % int(round(ydim)) print "PMEGridSizeZ %d" % int(round(zdim))