from __future__ import division, print_function import os, sys #This script updates the het dictionary used by Reduce to add hydrogens to het groups #It obsoletes and replaces update_het_dict.sh #Outline: #Download the lasted version of the het_dictionary from the PDB #Reformat certain columns of the het dictionary #Edit certain functional groups #Regularize endlines for Windows #Set filenames and clean up #----- Download het dict from PDB ----- import urllib import time print("Downloading most recent het dict", file=sys.stderr) #Download the lasted version of the het_dictionary from the PDB urllib.urlretrieve('ftp://ftp.wwpdb.org/pub/pdb/data/monomers/het_dictionary.txt', 'reduce_wwPDB_het_dict_download.txt') #Store time of download for versioning fetchtime = time.strftime("%Y%m%d_%H%M%S") #----- End download ----- #----- Adjust columns in het dict ----- def make_pdb3_format_atom(atomstring): #put atoms into a PDB3-like format #Note that this reformatting does not distinguish between calcium 'CA ' and # alpha carbons ' CA '. Both will be rendered as ' CA '. #No issues are known to result, but the possibility remains. atomstring = atomstring.strip() l = len(atomstring) if l == 4: return atomstring elif l == 3: return ' '+atomstring elif l == 2: return ' '+atomstring+' ' elif l == 1: return ' '+atomstring+' ' else: sys.stdout.write("ERROR: Strange atom name |"+atomstring+"|") sys.exit() infile = open("reduce_wwPDB_het_dict_download.txt") outfile = open("reduce_wwPDB_het_dict_adjusted.txt","w") print("Adjusting atom columns", file=sys.stderr) for line in infile: if line.startswith("CONECT"): x = line.split() numconnections = '%5s' % x[2] connections = x[3:] atomlist = [] for atom in connections: atomlist.append(make_pdb3_format_atom(atom)) outlist = [x[0]," ",make_pdb3_format_atom(x[1]),numconnections,' '.join(atomlist).rstrip()] print(''.join(outlist).strip(), file=outfile) elif line.startswith("RESIDUE"): x = line.split() resname = x[1].rjust(3) members = x[2].strip().rjust(7) print(x[0]+" "+resname+ members, file=outfile) elif line.startswith("END"): print(line.strip()," ", file=outfile) elif line.startswith("HETNAM") or line.startswith("HETSYN") or line.startswith("FORMUL"): outfile.write(line) else: print(line.strip(), file=outfile) infile.close() outfile.close() #----- End adjustment ----- #----- Remove H from Phosphate OH ----- #The default het dict contains some hydrogens that are not present at # physiological conditions. These hydrogens are usually present to make the # ligands charge-neutral. #The most prevalent of these is the H on phosphate OH. #Here we comment out those hydrogens so Reduce will not add them. #We also adjust SPD, SPM, and PAR nitrogens to have appropriate charge states import re class record(): def __init__(self): self.lines = [] self.has_O = False self.has_P = False def addline(self, line): self.lines.append(line) def output_record(self,outfile): for line in self.lines: outfile.write(line) outfile.write("\n") def remove_OH_on_P(self): line_index = 0 hydrogens_to_remove = [] while line_index < len(self.lines): line = self.lines[line_index] if line.startswith("CONECT"): parent_atom = line[11:15].strip() number_atoms = line[16:20].strip() connected_atoms = line[20:].split() if parent_atom in hydrogens_to_remove: self.lines[line_index] = "xx"+line #assumes hydrogen always comes after oxygen elif parent_atom.startswith("O") and number_atoms=="2": #need 2 connected atoms atom_1 = connected_atoms[0].strip() atom_2 = connected_atoms[1].strip() if atom_1.startswith("P") and atom_2.startswith("H"): hydrogens_to_remove.append(atom_2) #matches pattern for a phosphate oxygen, need to comment out the attached hydrogen line_index += 1 def fix_spd_charge(self): line_index = 0 while line_index < len(self.lines): line = self.lines[line_index] if line.startswith("CONECT N1 "): self.lines[line_index] = "CONECT N1 4 C2 HN11 HN12 HN13\n" elif line.startswith("CONECT N6 "): self.lines[line_index] = "CONECT N6 4 C5 C7 HN61 HN62\n" elif line.startswith("CONECT N10 "): self.lines[line_index] = "CONECT N10 4 C9 H101 H102 H103\n" elif line.startswith("CONECT HN12"): self.lines[line_index] = "CONECT HN12 1 N1\nCONECT NH13 1 N1\n" elif line.startswith("CONECT HN6 "): self.lines[line_index] = "CONECT HN61 1 N6\nCONECT HN62 1 N6\n" elif line.startswith("CONECT H102 1 N10"): self.lines[line_index] = "CONECT H102 1 N10\nCONECT H103 1 N10\n" line_index += 1 def fix_spm_charge(self): line_index = 0 while line_index < len(self.lines): line = self.lines[line_index] if line.startswith("CONECT N1 "): self.lines[line_index] = "CONECT N1 4 C2 HN11 HN12 HN13\n" elif line.startswith("CONECT N5 "): self.lines[line_index] = "CONECT N5 4 C4 C6 HN51 HN52\n" elif line.startswith("CONECT N10"): #HN0 -> H101 self.lines[line_index] = "CONECT N10 4 C9 C11 H101 H102\n" elif line.startswith("CONECT N14"): self.lines[line_index] = "CONECT N14 4 C13 HN41 HN42 HN43\n" elif line.startswith("CONECT HN12"): self.lines[line_index] = "CONECT HN12 1 N1\nCONECT NH13 1 N1\n" elif line.startswith("CONECT HN5"): self.lines[line_index] = "CONECT HN51 1 N5\nCONECT HN52 1 N5\n" elif line.startswith("CONECT HN0"): #HN0 -> H101 self.lines[line_index] = "CONECT H101 1 N10\nCONECT H102 1 N10\n" elif line.startswith("CONECT HN42"): self.lines[line_index] = "CONECT HN42 1 N14\nCONECT HN43 1 N14\n" line_index += 1 def fix_par_charge(self): line_index = 0 while line_index < len(self.lines): line = self.lines[line_index] if line.startswith("CONECT N21"): self.lines[line_index] = "CONECT N21 4 C21 HN21 HN22 HN23\n" elif line.startswith("CONECT N12"): self.lines[line_index] = "CONECT N12 4 C12 H121 H122 H123\n" elif line.startswith("CONECT N32"): self.lines[line_index] = "CONECT N32 4 C32 H321 H322 H323\n" elif line.startswith("CONECT N24"): self.lines[line_index] = "CONECT N24 4 C24 H241 H242 H243\n" elif line.startswith("CONECT N64"): self.lines[line_index] = "CONECT N64 4 C64 HN61 HN62 HN63\n" elif line.startswith("CONECT HN22"): self.lines[line_index] = "CONECT HN22 1 N21\nCONECT HN23 1 N21\n" elif line.startswith("CONECT H122"): self.lines[line_index] = "CONECT H122 1 N12\nCONECT H123 1 N12\n" elif line.startswith("CONECT H322"): self.lines[line_index] = "CONECT H322 1 N32\nCONECT H323 1 N32\n" elif line.startswith("CONECT H242"): self.lines[line_index] = "CONECT H242 1 N24\nCONECT H243 1 N24\n" elif line.startswith("CONECT HN62"): self.lines[line_index] = "CONECT HN62 1 N64\nCONECT HN63 1 N64\n" line_index += 1 def process(self, formul_line, outfile): phosphate_pattern = re.compile('(FORMUL)(.........)(.*H[1-9])(.*O[1-9])(.*P[1-9])') if phosphate_pattern.match(formul_line): self.remove_OH_on_P() if self.lines[0].startswith("RESIDUE SPD"): self.fix_spd_charge() elif self.lines[0].startswith("RESIDUE SPM"): self.fix_spm_charge() elif self.lines[0].startswith("RESIDUE PAR"): self.fix_par_charge() self.output_record(outfile) infile = open("reduce_wwPDB_het_dict_adjusted.txt") outfile = open("reduce_wwPDB_het_dict_no_poh.txt","w") print("Adjusting charges", file=sys.stderr) residue_line_pattern = re.compile('(RESIDUE)(...)(?P...)') formula_line_pattern = re.compile('FORMUL') for line in infile: if residue_line_pattern.match(line): #resname = residue_line_pattern.match(line).group('resname') rec = record() rec.addline(line) elif formula_line_pattern.match(line): rec.addline(line) rec.process(line, outfile) else: rec.addline(line) infile.close() outfile.close() #----- End remove H from POH ----- #----- Regularized endlines for Windows ----- #Swap newline characters for Windows-compatible ones to ensure usability on Win infile = open("reduce_wwPDB_het_dict_no_poh.txt") outfile = open("reduce_wwPDB_het_dict_endlines.txt","w") print("Switching to Windows endlines", file=sys.stderr) lines = infile.readlines() for line in lines: line=line.replace("\n","") outfile.write("%s\r\n" % line) infile.close() outfile.close() #----- End regularize endlines ----- #----- Rename and remove files ----- #The default path Reduce uses to find the het dict file can be checked by # invoking reduce -h and looking under the -DB flag. #Generally, Reduce is looking for a file named reduce_wwPDB_het_dict.txt. #Renaming the updated file is renamed to match that name allows it to be used # automatically in most circumstances. #The previous file is renamed with a datestamp. print("File cleanup", file=sys.stderr) existing_files = os.listdir(".") if "reduce_wwPDB_het_dict.txt" in existing_files: print(" The previous het dict was renamed to reduce_wwPDB_het_dict_obsoleted_"+fetchtime+".txt", file=sys.stderr) os.rename("reduce_wwPDB_het_dict.txt","reduce_wwPDB_het_dict_obsoleted_"+fetchtime+".txt") os.rename("reduce_wwPDB_het_dict_endlines.txt","reduce_wwPDB_het_dict.txt") os.remove("reduce_wwPDB_het_dict_download.txt") os.remove("reduce_wwPDB_het_dict_adjusted.txt") os.remove("reduce_wwPDB_het_dict_no_poh.txt") #----- End rename and remove -----