Fix a few bugs to do with counting hetatm residues

Author: Rob-Welch

prepmd/get_residues.py

@@ -9,6 +9,7 @@
 except:
     NO_MODELLER = True
 from prepmd import util
+from prepmd import ligand
 
 
 def get_residues_pdb(pdb, code, get_hetatms=False):
@@ -76,13 +77,23 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False):
                         sequence = line.split()[2]
                         seqres[chain] .append(sequence)
                 if line.startswith("#"):
-                    reading_seq = False # TODO: option to add marker for hetatms
+                    reading_seq = False
                 if line.startswith("HET") and not hetatms_found and get_hetatms:
                     hetatms_found = True
 
     if hetatms_found:
+        # count hetatm residues and add an equivalent number of "." entries in
+        # the fasta sequence
+        universe = ligand.load_universe(pdb)
+        ligands = universe.select_atoms('not protein and not water')
+        num_hetatm_residues = len(ligands.split("residue"))
         last_key = sorted(seqres.keys())[-1]
-        seqres[last_key] += ["..."] # not '.h.'?
+        seqres[last_key] += [num_hetatm_residues*"."]
+
+        
+    if get_hetatms and not hetatms_found:
+        raise ValueError("Was told to retrieve hetatms from "+pdb+" but none "
+                         "were found.")
 
     # convert to fasta
     fastas = []

prepmd/ligand.py

@@ -369,3 +369,9 @@ def create_ligand_system(
     # intergrator - variable langevin, tolerance: 0.0005 (or timestep < 0.001ps)
     # minimisation - tolerance = 2.5 kilojoule/(nanometer*mole)
     return openmm_system, openmm_topology, openmm_positions
+
+
+def get_ligand_centroid(pdb, traj):
+    universe = load_universe(pdb)
+    
+    ligand = u.select_atoms('not protein and not water')

prepmd/prep.py

@@ -274,7 +274,6 @@ def in_string(substr, text): return text == None or substr in text.lower()
     no_sim_output = False
     if split_hetatms:
         ligands = ligand.split_pdb_ligand(inmodel)
-    #if ligands and split_hetatms:
         print("Wrote hetatms/ligands to "+", ".join(ligands))
         no_sim_output = True # don't minimise the ligandless structure without
         # the ligands!

prepmd/util.py

@@ -28,7 +28,6 @@
     "XAA": "X",
     "TYR": "Y",
     "GLX": "Z",
-    ".h.": ".h." # leave hetatms
 }
 
 
@@ -40,7 +39,7 @@ def is_residue(resid):
     Returns:
         a boolean that is true if the string is a valid residue
     """
-    return resid in codes
+    return resid in codes or list(set(resid)) == ["."]
 
 
 def is_residue_sequence(sequence):
@@ -68,6 +67,8 @@ def three_to_one(resid, ignore_non_standard=False):
     Returns:
         the fasta residue, a string
     """
+    if list(set(resid)) == ["."]: # hetatms
+        return resid
     if is_residue(resid):
         return codes[resid]
     if ignore_non_standard:
@@ -85,11 +86,13 @@ def three_to_one_sequence(resids):
     Returns:
         the residue sequence in FASTA format
     """
+    print(resids)
     pdb_sequence = ""
     non_standard = []
     for resid in resids:
         try:
-            pdb_sequence += codes[resid]
-        except KeyError:
+            pdb_sequence += three_to_one(resid)
+        except ValueError:
             non_standard.append(resid)
+    print(pdb_sequence)
     return pdb_sequence