Merge pull request #10 from davidkastner/failure-modes

Updated the code to handle a variety of failure modes

qp/cli.py

@@ -98,9 +98,9 @@ def cli(i, o, modeller, protoss, coordination, skip):
                 show_default=False,
             )
         )
-        charge = click.confirm("> Compute charges (requires Protoss)")
-        count = click.confirm("> Count residues")
-        xyz = click.confirm("> Write XYZ files")
+        charge = click.confirm("> Compute charges (requires Protoss)", default=True)
+        count = click.confirm("> Count residues", default=True)
+        xyz = click.confirm("> Write XYZ files", default=True)
         if capping or charge:
             protoss = True
         click.echo("")
@@ -170,6 +170,7 @@ def cli(i, o, modeller, protoss, coordination, skip):
             try:
                 if protoss:
                     add_hydrogens.adjust_active_sites(path, metals)
+                    add_hydrogens.rename_nterminal(path)
                 clusters = coordination_spheres.extract_clusters(
                     path, f"{o}/{pdb}", metals,
                     limit, ligands, capping, charge, count, xyz

qp/cluster/coordination_spheres.py

@@ -343,7 +343,7 @@ def compute_charge(spheres):
         "GLU": ["HE2", "HOE1"],
         "CYS": ["HG"],
         "TYR": ["HH"],
-        "OCS": []
+        "OCS": [],
     }
 
     charge = []
@@ -358,6 +358,9 @@ def compute_charge(spheres):
             # Check for C-terminus
             if res.has_id("OXT"):
                 c -= 1
+            # Check for N-terminus
+            if res.has_id("NT"):
+                c += 1
         charge.append(c)
     return charge
 

qp/structure/add_hydrogens.py

@@ -154,6 +154,33 @@ def accept_atom(self, atom):
     io.save(path, AtomSelect())
 
 
+def rename_nterminal(path):
+    """Changes N-terminal nitrogen atom name to NT in the PDB file."""
+
+    with open(path, 'r') as f:
+        pdb_content = f.read()
+
+    lines = pdb_content.split('\n')
+    new_lines = []
+
+    seen_chains = set()
+
+    for line in lines:
+        if line.startswith("ATOM") or line.startswith("HETATM"):
+            current_chain = line[21]
+            atom_name = line[12:16].strip()  # Extracts the atom name
+
+            # Rename N-terminal nitrogen
+            if current_chain not in seen_chains and atom_name == "N":
+                line = line[:12] + " NT " + line[16:]
+                seen_chains.add(current_chain)
+
+        new_lines.append(line)
+
+    with open(path, 'w') as f:
+        f.write('\n'.join(new_lines))
+
+
 def compute_charge(path):
     """
     Computes the total charge of each ligand

qp/structure/missing_loops.py

@@ -71,7 +71,7 @@
     "GLU": "E",
     "TYR": "Y",
     "MET": "M",
-    "MSE": "M",  # other non-standard residues not recognized by MODELLER
+    "MSE": "M",
 }
 
 
@@ -164,7 +164,8 @@ def get_residues(path):
     for chain_res in residues:
         while chain_res and chain_res[0][2] == "R":
             chain_res.pop(0)
-        while chain_res and chain_res[-1][2] == "R":
+        # Now, check for unresolved residues
+        while chain_res and chain_res[-1][2] == "R" and chain_res[-1][1] in AA.values():
             chain_res.pop()
 
     return residues
@@ -199,6 +200,42 @@ def write_alignment(residues, pdb, path, out):
         f.write(f">P1;{pdb}_fill\nsequence:::::::::\n{seq_fill}*\n")
 
 
+def fix_pdb_numbering(pdb_content):
+    """Handles insertion codes from Modeller."""
+
+    lines = pdb_content.split('\n')
+    new_lines = []
+
+    prev_chain = None
+    prev_residue_number = None
+    seen_residues = set()
+
+    for line in lines:
+        if line.startswith("ATOM") or line.startswith("HETATM"):
+            current_chain = line[21]
+            residue_number_with_letter = line[22:27].strip()  # Residue number along with potential insertion code
+            has_char = line[26] != ' '
+
+            # If we encounter a new chain or a residue we haven't seen, reset/update values
+            if current_chain != prev_chain or residue_number_with_letter not in seen_residues:
+                prev_chain = current_chain
+                if has_char:
+                    prev_residue_number = prev_residue_number + 1 if prev_residue_number else int(residue_number_with_letter[:-1])
+                else:
+                    prev_residue_number = int(residue_number_with_letter)
+                seen_residues.add(residue_number_with_letter)
+
+            # If there's an insertion code, replace the residue number and remove the insertion code
+            if has_char:
+                line = line[:22] + f"{prev_residue_number:4} " + line[27:]
+
+            new_lines.append(line)
+        else:
+            new_lines.append(line)
+
+    return '\n'.join(new_lines)
+
+
 def build_model(residues, pdb, ali, out, optimize=1):
     """
     Runs MODELLER for the given alignment file
@@ -281,6 +318,15 @@ def special_patches(self, aln): # renumber residues to avoid hybrid-36 notation
 
         # Transfer the residue and chain ids and write out the modified MODEL:
         mdl_built.res_num_from(mdl_reference, aln)
-        mdl_built.write(file=os.path.basename(out))
+        file=os.path.basename(out)
+        mdl_built.write(file)
+
+        # Example usage:
+        with open(file, 'r') as f:
+            content = f.read()
+
+        corrected_content = fix_pdb_numbering(content)
+        with open(file, 'w') as f:
+            f.write(corrected_content)
 
     os.chdir(cwd)