Bug fixes, docs, simulation tools

Author: Rob-Welch

README.md

@@ -4,14 +4,32 @@ A utility to automatically prepare structures from the PDB for molecular dynamic
 ## Features
 * [X] Automatically download structures, sequences and metadata from the PDB and UNIPROT
 * [X] Automatically fill missing loops with modeller
-* [X] Automatically add missing atoms and convert non-standard residues (using pdbfixer)
+* [X] Automatically add missing atoms and convert non-standard residues with pdbfixer
 * [X] Automatically propagate metadata through to finalised structure files
-* [X] Automatically resolve steric clashes
-* [ ] Automatically dock ligands
+* [X] Automatically resolve steric clashes and minimise structures
+* [X] A simple command-line interface for running molecular dynamics with OpenMM
 * [ ] AIIDA integration
 
 ## Installation
-* Download this repo and enter the root folder
+* Install [Conda](https://conda-forge.org/download/) (if you don't already have it)
+* Clone this repo and enter the folder: `git clone git@github.com:HECBioSim/prepmd.git && cd prepmd` 
 * Run `conda env create --name prepmd --file environment.yaml && conda activate prepmd && pip install .`
-* For the modeller part of the workflow to work, you need to get a [modeller license key](https://salilab.org/modeller/registration.html) and add it to modeller's config.py file. If you use conda, the key will be in envs/prepmd/lib/modeller-10.7/modlib/modeller/config.py relative to your conda install directory.
-* Run `prempmd` from the prepmd environment. A basic example: `prepmd 6xov 6xov_processed.pdb` will download the structure for PDB entry 6xov, process it and write it to `6xov_processed.pdb`.
+* For the modeller part of the workflow to work, you need to get a [modeller license key](https://salilab.org/modeller/registration.html) and add it to modeller's config.py file. If you use conda, the key will be in `envs/prepmd/lib/modeller-10.7/modlib/modeller/config.py` relative to the path where conda is installed.
+
+## Preparing structures from the PDB for simulation
+* A basic example: `prepmd 6xov 6xov_processed.pdb` will download the structure for PDB entry `6xov`, process it and write it to `6xov_processed.pdb`.
+* If you already have a pdb file, you can instead run: `prepmd --structure 6xov_input.pdb 6xov 6xov_processed.pdb`. You still need to supply a PDB code, as the various file formats used by prepmd require one to be present.
+* `prepmd` will attempt to guess the correct file formats from the filenames it's given. It won't perform implicit conversions, so make sure to start and end with the same file type.
+* By default, `prepmd` will leave intermediate files in a randomly-named temporary directory. You can set the name of this directory: `prepmd --wdir 6xov_temp 6xov 6xov.cif`.
+* By default, `prepmd` will read missing residues from the pdb/mmcif metadata, attempt to align the missing residues with the currently present residues, and then build missing loops. You can manually provide a FASTA file containing the alignment data with `--fasta`. You can also ask mdprep to get the sequence data from UNIPROT instead, with `--download`, though this is not recommended, as the raw sequence data can be very different from the PDB and cause the alignment to fail.
+* Note: while both pdb and mmCif are supported, using the mmCif format is strongly recommended, as the pdb format has been deprecated since 2024.
+* Use `prepmd --help` for a full list of parameters. 
+
+## Running MD simulations
+* `runmd` can run MD simulations using OpenMM.
+* A basic example: `runmd structure.cif --min_out structure_minimised.cif --traj_out traj.xtc --md_steps 5000 --step 100` will minimise and run a simulation of structure.cif, writing a trajectory to `traj_out.xtc`, for 5000 steps, saving one trajectory frame every 100 steps.
+* If you already have a minimised structure, you can skip minimisation: `runmd structure.cif --traj_out traj.xtc --md_steps 5000 --step 100 -nomin -notest`
+* Solvate the simulation box: `runmd structure.cif -o structure_minimised.cif --traj_out traj.xtc --md_steps 500 --step 10 -solv tip4pew`. tip3p, tip4pew and spce are supported. You can also add pressure coupling with `--pressure 1.0` (for 1 bar)
+* Run with different force fields: `runmd structure.cif -o structure_minimised.cif --traj_out traj.xtc --md_steps 500 --step 50 -ff amber14` runs with amber14. AMOEBA is also available, and amber19 is available if you have a recent version of OpenMM.
+* Finally, you may wish to fix the backbone in place and just equilibrate the side chains: `runmd structure.cif -o structure_minimised.cif --fix_backbone -solv tip4pew --notest`
+* Use `runmd --help` for a full list of parameters. 

prepmd/download.py

@@ -20,13 +20,23 @@ def get_structure(pdb_id, directory, file_format="mmCif"):
     returns:
         path to the downloaded file.
     """
-    if file_format == "mmCif":
+    
+    if file_format == "mmCif" or file_format == "cif":
         format_str = "cif"
     if file_format == "pdb":
         format_str = "pdb"
-    urllib.request.urlretrieve("https://files.rcsb.org/download/" +
-                               pdb_id+"."+format_str, directory+sep+pdb_id +
-                               "."+format_str)
+    try:
+        url = "https://files.rcsb.org/download/"+pdb_id+"."+format_str
+        destination = directory+sep+pdb_id+"."+format_str
+        urllib.request.urlretrieve(url, destination)
+    except urllib.error.HTTPError as e:
+        r = requests.get(url.replace(".pdb", ".cif"))
+        if r.status_code == 200:
+            msg = "No PDB for "+pdb_id+" exists (but an mmcif structure does). "
+            "Run with --fmt cif to use it."
+            raise IOError(msg)
+        else:
+            raise e
     return directory+sep+pdb_id+"."+format_str
 
 

prepmd/model.py

@@ -187,7 +187,7 @@ def fix_missing_residues(code, fastafile, alignmentout, inmodel, outmodel,
             print("Trying to use the sequence data in the input file...")
 
     if not fastafile:
-        print("No fasta file. Getting missing residues from input file...")
+        print("Getting missing residues from input file...")
         pdb_res = get_residues.get_residues_pdb(inmodel, code)
         pdb_fullseq = get_residues.get_fullseq_pdb(inmodel, code)
         with open(code+".seq", "w") as file:

prepmd/prep.py

@@ -11,6 +11,7 @@
 import random
 import string
 import sys
+import shutil
 
 from prepmd import download
 from prepmd import run
@@ -34,10 +35,10 @@
 parser.add_argument("-a", "--alignmentout",
     help="Alignment output file from sequences aligned for loop filling",
     default="alignment_out.fasta")
-parser.add_argument("-fmt", "--format",
+parser.add_argument("-fmt", "--dlformat",
     help="Structure format to download (only used when no structure file is "
     "provided)",
-    default="mmCif")
+    default=None)
 parser.add_argument(
     "-q", "--quiet", help="Do not print debug info", action="store_true")
 parser.add_argument("-e", "--fixstart",
@@ -62,17 +63,8 @@ def id_generator(size=6, chars=string.ascii_uppercase + string.digits):
     return ''.join(random.choice(chars) for _ in range(size))
 
 
-#def run_test(code, wdir):
-#    params = {"code": code,
-#              'outmodel': "/home/rob/temp/missingresiduestestrun/"+code+"_ready.pdb",
-#              "workingdir": "/home/rob/temp/missingresiduestestrun/"+code+"_"+id_generator(6)
-#              }
-#    prep(params["code"], params["outmodel"],
-#         params["workingdir"], structure_format="pdb")
-
-
 def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None,
-         alignmentout="alignment_out.fasta", structure_format="mmCif",
+         alignmentout="alignment_out.fasta", download_format="mmCif",
          quiet=False, fix_after=True, download_sequence=False,
          fix_missing_atoms=True):
     """
@@ -89,7 +81,7 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None,
         inmodel: path to input structure file, a string.
         alignmentout: output file from the alignment of the sequence from the
         structure file and the original sequence.
-        structure_format: format of the structure file, a string. Can be
+        download_format: format of the structure file, a string. Can be
         'mmCif' or 'pdb'.
         quiet: if true, don't print anything. boolean
         fix_after: if true, fix the pdb file after loop building. If false,
@@ -102,10 +94,36 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None,
     Returns:
         nothing, but writes out a file to outmodel.
     """
+    
+    # infer download format from output format
+    if not download_format:
+        if (".pdb") in outmodel:
+            download_format = "pdb"
+        elif (".cif") in outmodel or ".mmcif" in outmodel or (
+                ".pdbx") in outmodel:
+            download_format = "mmCif"
+    
+    # check that input/output are specified in the same format
+    # i'm not against converting the files but it shouldn't happen implicitly
+    in_string = lambda substr, text : text == None or substr in text.lower()
+    if in_string(".pdb", inmodel) and in_string(
+            ".pdb", outmodel) and in_string("pdb", download_format):
+        pass
+    elif in_string(".cif", inmodel) and in_string(
+            ".cif", outmodel) and in_string("cif", download_format):
+        pass
+    else:
+        raise IOError("Inputs and outputs are in different formats! Please "
+                      "use only one format (ideally cif)" )
 
     if not os.path.isdir(workingdir):
         pathlib.Path(workingdir).mkdir(parents=True, exist_ok=True)
 
+    if inmodel:
+        # todo: copy input file to working directory
+        shutil.copyfile(inmodel, workingdir+os.path.sep+inmodel)
+
+    run_dir = os.getcwd()
     os.chdir(workingdir)
 
     # check folder for strucutre/sequence files
@@ -129,7 +147,7 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None,
     if code and inmodel == None and folder == None:
         print("Downloading structure file")
         inmodel = download.get_structure(
-            code, workingdir, file_format=structure_format)
+            code, workingdir, file_format=download_format)
 
     # fix
     if not fix_after:
@@ -174,9 +192,12 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None,
     if ".cif" in inmodel or ".mmcif" in inmodel:
         print("Metadata restoration not implemented for mmCif (yet)")
 
-    print("Testing "+code)
+    print("Simulating "+code)
     run.test_sim(outmodel)
-    print("Wrote "+outmodel)
+    print("Done.")
+    
+    if not os.path.isabs(outmodel):
+        shutil.copyfile(outmodel, run_dir+os.path.sep+outmodel)
 
 def entry_point():
     "CLI entry point function. Uses sys.argv and argparse args object."
@@ -189,7 +210,7 @@ def entry_point():
         args.wdir = args.code+"_"+id_generator(6)
     prep(args.code, args.out, os.getcwd()+os.path.sep+args.wdir,
          folder=args.directory, fastafile=args.fasta, inmodel=args.structure,
-         alignmentout=args.alignmentout, structure_format=args.format,
+         alignmentout=args.alignmentout, download_format=args.dlformat,
          quiet=args.quiet, fix_after=fix_after,
          download_sequence=args.download, fix_missing_atoms=args.leavemissing)
 
@@ -199,7 +220,7 @@ def test_mmcif_support():
     prep("6xov",
          os.getcwd()+os.path.sep+"6xov"+"_"+genid+".cif",
          os.getcwd()+os.path.sep+"testout"+os.path.sep+"6xov"+"_"+genid,
-         structure_format="mmCif")
+         download_format="mmCif")
 
 def test_minimise():
     genid = id_generator(6)
@@ -221,29 +242,48 @@ class style():
         RESET = '\033[0m'
 
     #testinputs = ["6xov", "9CS5", "8RM8", ]#"8VV2", "9B8B", "8CAE", "8QZA", "8RTL", "8RTO", "9A9W"] # regular pdbs
-    testinputs = ["6TY4"] # clash
+    #testinputs = ["6TY4", "6XOV", "9CS5", "8CAE", "8QZA", "8RTO"]
+    
+    tests = [
+        {"id": "6TY4", "format":"pdb"},
+        {"id": "6XOV", "format":"pdb"},
+        {"id": "9CS5", "format":"pdb"},
+        {"id": "8CAE", "format":"pdb"},
+        {"id": "8QZA", "format":"pdb"},
+        {"id": "8RTO", "format":"mmCif"},
+        {"id": "7IB8", "format":"mmCif"},
+        {"id": "9A9G", "format":"mmCif"},
+        {"id": "9I3U", "format":"pdb"},
+        #test_minimise,
+        #test_mmcif_support
+    ]
+    #testinputs = ["8rto"]
 
     results = {}
     state = 0
     cwd = os.getcwd()
 
-    for test in range(len(testinputs)):
+    for test in range(len(tests)):
         try:
             os.chdir(cwd)
-            code = testinputs[test]
-            print(f"Testing {code} ({test}/{len(testinputs)})")
+            code = tests[test]["id"]
+            curr_format =  tests[test]["format"]
+            print(f"Testing {code} ({test}/{len(tests)})")
             genid = id_generator(6)
             pathlib.Path(os.getcwd()+os.path.sep+"testout").mkdir(
                 parents=True, exist_ok=True)
-            prep(code,
-                 os.getcwd()+os.path.sep+testinputs[test]+"_"+genid+".pdb",
-                 os.getcwd()+os.path.sep+"testout"+os.path.sep+code+"_"+genid,
-                 structure_format="pdb")
+            if type(tests[test]) == dict:
+                prep(code,
+                     os.getcwd()+os.path.sep+code+"_"+genid+"."+"cif",
+                     os.getcwd()+os.path.sep+"testout"+os.path.sep+code+"_"+genid,
+                     download_format=curr_format)
+            elif callable(tests[test]):
+                test()
             print(f"{style.GREEN}PASSED: {test} {style.RESET}")
-            results[testinputs[test]] = "PASS"
+            results[code] = "PASS"
         except Exception as e:
-            print(f"{style.RED}FAILED: {testinputs[test]}{style.RESET}")
-            results[testinputs[test]] = e
+            print(f"{style.RED}FAILED: {test}{style.RESET}")
+            results[code] = e
             state = 1
     print("")
     print("RESULTS:")